Exact nonequilibrium hole dynamics, magnetic polarons and string excitations in antiferromagnetic Bethe lattices

To progress on solving for the coefficients of the nonequilbrium wave function, we introduce the retarded and advanced wave functions

allowing us to regularize the Fourier transformation to frequency space with the positive infinitesimal $\eta>0$, and the Heaviside step function $\theta(\tau)$. Consequently, we may express the Schrödinger equation, $i\partial_{\tau}\ket{\Psi(\tau)}=\hat{H}\ket{\Psi(\tau)}$, as

As these equations are complex conjugate of each other, it follows that $\ket{\Psi^{\rm A}(\omega)}=[\ket{\Psi^{\rm R}(\omega)}]^{*}$. We denote the coefficients of the retarded state $R^{(d)}(\omega)$. The dynamical coefficients in Eq. (10) are then retrieved from $R^{(d)}(\omega)$ by a Fourier transformation

Using Eq. (12), the equations of motion for the retarded state coefficients become

The term $i$ in the equation for $R^{(0)}(\omega)$ comes directly from the term $i\ket{\Psi(\tau=0)}$ in Eq. (12). Furthermore, the magnetic energy cost $E_{J}^{(d)}$ comes from applying the spin-spin interaction Hamiltonian in Eq. (8) to the $d$th term in the wave function. Finally, the hopping Hamiltonian simply relates the coefficient at depth $d$ to the coefficients at depth $d-1$ and $d+1$. At $d=0$, there are $q$ nearest neighbors at depth $d=1$. All subsequent depths $d\leq d_{\rm tot}-1$, however, only has $q-1$ nearest neighbors at depth $d+1$ and a single neighbor one depth up, $d-1$. This leads to the terms $(q-1)\cdot t\cdot R^{(d+1)}(\omega)$ and $1\cdot t\cdot R^{(d-1)}(\omega)$ for $1\leq d\leq d_{\rm tot}$. Finally, by assumption the bottom depth is $d_{\rm tot}$, and so $R^{(d)}=0$ for $d>d_{\rm tot}$.

The value of $E_{J}^{(d)}$ may be understood recursively, starting from $d=0$. Here, the hole is located at the central site, and there are $q$ broken antiferromagnetic spin bonds of strength $J_{z}/2$, as indicated in Fig. 2(a). Note that Eq. (8) implies that a hole at a site has the same magnetic energy cost as a spin excitation at that site. Therefore, $E_{J}^{(0)}=E_{0}(q,d)+qJ_{z}/2$, where $E_{0}(q,d)$ is the ground state energy of the Ising antiferromagnet [Eq. (4)]. When the hole is at depth $d=1$, as shown in Fig. 2(b), there are $q-1$ additional broken bonds. So $E_{J}^{(1)}=E_{J}^{(0)}+(q-1)J_{z}/2$. At every subsequent depth $d\geq 2$ [Figs. 2(c) and 2(d)], the spin bonds along the path of the hole are repaired. Therefore, each hop is associated with $q-1$ new broken bonds and $1$ repaired bond. The magnetic energy cost for $d\geq 2$ is, thus, $(q-2)J_{z}/2$ for every step in depth. In this way,

for $1\leq d\leq d_{\rm tot}-1$. Here, we let $E_{J}^{(d)}\to E_{J}^{(d)}-E_{J}^{(0)}$, measuring all energies with respect to $E_{J}^{(0)}$. Finally, at depth $d_{\rm tot}$ the absence of neighbors at $d_{\rm tot}+1$ in principle results in a surface effect, in which there is no additional magnetic energy cost for the hole to hop to the surface, $E_{J}^{(d_{\rm tot})}=E_{J}^{(d_{\rm tot}-1)}$. While this presumably has an effect on the dynamics for small lattices, we focus on the behavior for large total depths, $d_{\rm tot}\gg 1$, and ignore the surface effects. With the magnetic energy costs for the hole to hop to the $d$th depth in place [Eq. (15)], we can now solve the equations of motion in Eq. (14) iteratively. Specifically, the equation for $R^{(d_{\rm tot})}$ simply yields

Continuing this process for $d_{\rm tot}\geq d\geq 1$, we get the recursion relation

defining $G^{(d_{\rm tot}+1)}(\omega)=0$. We refer to $G^{(d)}(\omega)$ as the retarded Green’s function at depth $d$. Finally, inserting the solution from $d=1$ into the equation for $R^{(0)}$, we get $R^{(0)}(\omega)=iG^{(0)}(\omega)$, with

Using the recursion relation in Eq. (17) in reverse, we obtain the general result for the coefficients of the retarded wave function in frequency space

We end this subsection by connecting the coefficients of the wave function to a specific set of many-body correlation functions. Specifically, we have that

This more clearly demonstrates that the $d$th coefficient of the wave function is the probability amplitude of the hole to be at depth $d$ at site ${\bf j}_{d}$ with $d$ spin excitations after a time $\tau$.

In the thermodynamic limit, $d_{\rm tot}\to\infty$, the recursion relation for the retarded Green’s functions in Eq. (17) can be written as a self-consistency equation for a single Green’s function $G_{1}(\omega)$. Specifically, we can write $G^{(d)}(\omega)=G_{1}(\omega-E_{J}^{(d)})$ for $d\geq 1$, where

For brevity, in this subsection we let $\omega+i\eta\to\omega$. Similarly, we get the Green’s function for the central site, $d=0$, to be

Inspired by Ref. Chernyshev and Leung (1999), we anticipate that the $G_{1}$ Green’s function can be expressed in terms of Bessel functions. In fact, using the ansatz $G_{1}(\omega)=-1/(\sqrt{q-1}t)\cdot Y(\omega)/Y(\omega+(q-2)J/2)$, the resulting recursive relation for the $Y$ functions coincide with that for the Bessel functions of the first kind. This results in

for $q\geq 3$, with $\Omega(\omega)=-2\omega/(q-2)J_{z}$. This generalizes the result in Ref. Chernyshev and Leung (1999), where the $q=4$ Bethe lattice is effectively used as an approximate description of a two-dimensional square lattice, and coincides with what is found in Ref. Wrzosek and Wohlfeld (2021). We can now use the simple fraction of Bessel functions in Eq. (23) together with $G^{(d)}(\omega)=G_{1}(q,\omega-E_{J}^{(d)})$ for $d\geq 1$, to find the explicit expression for the coefficients of the nonequilibrium many-body wave function in the thermodynamic limit

where the magnetic energy cost for a hole at depth $d$, $E_{J}^{(d)}$, is given in Eq. (15). In the limit of infinitely strong interactions, $J_{z}/t\to 0^{+}$, Eq. (21) leads to a second order equation for $G_{1}$, yielding

This results in the hole Green’s function

Equations (25) and (26) together with Eq. (19) describes a free quantum walk of a single particle in Bethe lattices Eckstein et al. (2005), as one might expect in this extreme limit, where spin interactions play no role.

The static Schrödinger equation,

simply corresponds to removing the $i\ket{\Psi(\tau=0)}$ term in Eq. (12), and replacing the frequency $\omega+i\eta$ with the eigenstate energy $\varepsilon_{n}$. As a result, we simply have to remove the term $i$ in the equations of motion in Eq. (14) for the nonequilibrium case, and the same recursive relation in Eq. (17) for the nonequilibrium coefficients, therefore, applies to these eigenstates. Explicitly,

with the Green’s functions $G^{(d)}$ evaluated at the quasiparticle pole $\varepsilon_{n}$. Inserting the result at depth $d=1$ into the equation of motion at depth $d=0$, we get

Here, $\Sigma^{(0)}(\omega)=qt^{2}\cdot G^{(1)}(\omega)$ is the self-energy associated with the Green’s function $G^{(0)}(\omega)$ in Eq. (18). This equation simply states the fact that the poles of $G^{(0)}$ are eigenstate energies of the system, the lowest of which gives the true ground state of the system. To find the lowest order coefficient of the $n$th state, $C^{(0)}_{n}$, we compute the full norm of the wave function, $\ket{\Psi_{n}}$. Repeated use of the recursion relation in Eq. (29) yields

Investigating the quasiparticle residue, $Z_{n}=[1-\partial_{\omega}\Sigma^{(0)}(\omega)|_{\omega=\varepsilon_{n}}]^{-1}$, shows that the expression within the square brackets above is simply $Z^{-1}_{n}$. In turn, $C^{(0)}_{n}=\sqrt{Z_{n}}$, whereby we may write the eigenstate coefficients as

In the thermodynamic limit, these coefficients can, in complete analogy to the nonequilibrium case investigated in Sec. IV.2, also be written as fractions of Bessel functions. This has previously been established for the ground state in Ref. Bieniasz et al. (2019), where the Bethe lattice description is used as a variational ansatz for the square lattice case. In Fig. 3, we plot the spectral function $A(\omega)=-2{\rm Im}\,G_{0}(\omega)$, computed in the thermodynamic limit from Eq. (22). This reveals the poles of $G_{0}$, and thereby the eigenstate energies $\varepsilon_{n}$. In the special case of $q=2$ Batista and Ortiz (2000), the system is one-dimensional and the spectral function splits into a quasiparticle pole at

and a continuum of states residing between $J_{z}/2\pm 2t$; see also Fig. 3(a). This describes the effect known as spin-charge separation, which has been intensely studied both theoretically and experimentally in one-dimensional systems Tomonaga (1950); Kim et al. (1996); Segovia et al. (1999); Auslaender et al. (2005); Kim et al. (2006); Jompol et al. (2009); Vijayan et al. (2020). In the present case, as the charge degree of freedom moves (the hole), a single spin excitation is created and remains at the original site of the hole. Due to intense previous studies of these effects, we will not go any further with it in the present paper.

For more nearest neighbors, $q\geq 3$, a series of quasiparticles peaks arises [see Figs. 3(b)–3(d)]. In Sec. V.2, we show that at strong interactions, their energies all scale with $(J_{z}/t)^{2/3}$. The origin of this scaling is that the system in this limit can be rephrased as a continuum model, in which the hole moves in a linear potential with strength $J_{z}$. While such a rephrasing in a regular lattice is approximate Bulaevskii et al. (1968); Kane et al. (1989); Grusdt et al. (2018a), this becomes exact in Bethe lattices as we shall unfold in Sec. V.2. In the limit of weak interactions on the other hand, $J_{z}\gg t$, the hole becomes immobile and the $n$th eigenstate is completely localized at the $n$th depth with an energy simply given by the magnetic energy cost, $\varepsilon_{n}\to E_{J}^{(n)}=[(q-1)+(n-1)(q-2)]J_{z}/2$ [Eq. (15)] for $n\geq 1$. In fact, to order $\mathcal{O}(t^{2}/J_{z})$, only the magnetic polaron ground state and the first string excitation are affected and attain the energies

Here, $\varepsilon_{0}=\varepsilon^{(-)}$ and $\varepsilon_{1}=\varepsilon^{(+)}$ are the ground and lowest string excitation energies respectively.

The method used here to find the explicit coefficients of these eigenstates is similar to Refs. Ramšak and Horsch (1998); Reiter (1994), where the approximate magnetic polaron eigenstates within the SCBA Kane et al. (1989); Martinez and Horsch (1991); Liu and Manousakis (1991) are found. We finally note that these many-body eigenstates can in principle be used to completely characterize the full nonequilibrium dynamics. Specifically, $\braket{\Psi_{n}}{\Psi(\tau)}=\bra{\Psi_{n}}{\rm e}^{-i\hat{H}\tau}\hat{h}^{\dagger}_{0}\ket{{\rm AF}}={\rm e}^{-i\varepsilon_{n}\tau}\!\bra{\Psi_{n}}\hat{h}^{\dagger}_{0}\ket{{\rm AF}}=\sqrt{Z_{n}}\cdot{\rm e}^{-i\varepsilon_{n}\tau}$. Therefore, the dynamics is the quantum interference of the polaron ground state and the string excitations. While this lends insight into the conceptual nature of the dynamics, it does not actually simplify the characterization of the underlying hole motion.

In this subsection, we describe in detail the continuum limit arising for $J_{z}/t\to 0^{+}$ and a number of nearest neighbors $q\geq 3$. In addition to the dominant $(J_{z}/t)^{2/3}$ contribution to the energy found previously Kane et al. (1989); Grusdt et al. (2018a); Zhong and Sorella (1995), we also variationally determine a term linear in $J_{z}$, which allows us to accurately describe the limiting behavior of the quasiparticle residue for $J_{z}/t\to 0^{+}$. Additionally, we explicitly compare the exact eigenstates found in Sec. V.1 to the strong-coupling limit investigated here.

We begin by transforming the equations of motion to an effective one-dimensional setting valid for the depth symmetric eigenstates found in the previous Sec. V.1. To this end, we define $\psi(0)=C^{(0)}$, and $\psi(d)=(-1)^{d}\sqrt{q(q-1)^{d-1}}\cdot C^{(d)}$ for $d\geq 1$, similar to the approach in Ref. Grusdt et al. (2018a). In this manner, the wave function fulfills a one-dimensional normalization condition: $\sum_{d}|\psi(d)|^{2}=1$. With this in hand, the equations of motion become

where the lower equation is valid for $d\geq 2$. We may now formulate the continuum model by taking this lower equation and rephrase it as a continuous equation in the limit of $J_{z}\ll t$. This is possible, because the wave function spreads out over an ever increasing number of lattice sites in this limit, as we shall explicitly see in Sec. VI. Here, it is first beneficial to rewrite this equation according to

using Eq. (15), and defining $\sqrt{q-1}t\cdot\tilde{\varepsilon}=\left[\varepsilon+2\sqrt{q-1}t\right]$. Next, we let $d=x/\lambda$ and introduce a rescaled wave function $\phi(x)=\psi(x/\lambda)/\sqrt{\lambda}$. Inserting this in Eq. (35) and dividing out $\lambda^{2}$, we get

To eliminate $J_{z}/t$, we set $\lambda=[(q-2)J_{z}/(2\sqrt{q-1}t)]^{1/3}$. For $J_{z}/t\to 0^{+}$, the fraction in the second line of Eq. (36) approaches the second order derivative. Setting $a=\tilde{\varepsilon}/\lambda^{2}$, we, hereby, obtain

which is accurate up to order $\lambda^{2}\propto(J_{z}/t)^{2/3}$. The wave function $\phi(x)$ fulfills the continuous normalization condition: $\int_{0}^{\infty}\!{\rm d}x|\phi(x)|^{2}=1$. Letting $y=x-a$, we get the Airy equation

for $y\geq-a$, where $f(y)=\phi(y+a)$. The solutions to this equation can thus be written as a superposition of the Airy functions: $f(y)=A\cdot{\rm Ai}(y)+B\cdot{\rm Bi}(y)$. Since ${\rm Bi}(y)$ blows up for $y\to\infty$, $B=0$. Furthermore, the continuous one-dimensional description, which is exact for $J_{z}/t\to 0^{+}$, entails that $\phi(x)$ must vanish for $x<0$, and therefore also at $x=0$. In turn, $-a$ must be a zero of the Airy function ${\rm Ai}(y)$. In this way, we obtain the sought set of eigenstates in the strongly interacting limit

The $n$th eigenstate is, thus, defined by the $n$th order zero of the Airy function $-a_{n}$. Also, $A_{n}^{-2}=\int_{0}^{\infty}{\rm d}x\,|{\rm Ai}\left(x-a_{n}\right)|^{2}$ ensures normalized eigenstates. For any nonzero value of $J_{z}/t$, the eigenstates remain nonzero at the origin, $d=0$. This turns out to yield a correction linear in $J_{z}/t$. To accommodate for this, we let $d\to d-d_{0}$ in Eq. (39) and use $d_{0}$ as a variational parameter. A rather lengthy calculation (see Appendix A) yields the variational energy for the $n$th eigenstate