# Spatial structure of magnetic polarons in strongly interacting
antiferromagnets

K. K. Nielsen  Center for Complex Quantum Systems, Department of Physics and
Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark  M.
A. Bastarrachea-Magnani  Center for Complex Quantum Systems, Department of
Physics and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C,
Denmark  Departamento de Física, Universidad Autónoma Metropolitana-
Iztapalapa, Av. San Rafael Atlixco 186, CP 09340 CDMX, Mexico  T. Pohl  Center
for Complex Quantum Systems, Department of Physics and Astronomy, Aarhus
University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark  G. M. Bruun  Center
for Complex Quantum Systems, Department of Physics and Astronomy, Aarhus
University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark  Shenzhen Institute
for Quantum Science and Engineering and Department of Physics, Southern
University of Science and Technology, Shenzhen 518055, China

###### Abstract

The properties of mobile impurities in quantum magnets are fundamental for our
understanding of strongly correlated materials and may play a key role in the
physics of high-temperature superconductivity. Hereby, the motion of hole-like
defects through an antiferromagnet has been of particular importance. It
creates magnetic frustrations that lead to the formation of a quasiparticle,
whose complex structure continues to pose substantial challenges to theory and
numerical simulations. In this article, we develop a non-perturbative
theoretical approach to describe the microscopic properties of such magnetic
polarons. Based on the self-consistent Born approximation, which is provenly
accurate in the strong-coupling regime, we obtain a complete description of
the polaron wave function by solving a set of Dyson-like equations that permit
to compute relevant spin-hole correlation functions. We apply this new method
to analyze the spatial structure of magnetic polarons in the strongly
interacting regime and find qualitative differences from predictions of
previously applied truncation schemes. Our calculations reveal a remarkably
high spatial symmetry of the polaronic magnetization cloud and a surprising
misalignment between its orientation and the polaron crystal momentum. The
developed framework opens up an approach to the microscopic properties of
doped quantum magnets and will enable detailed analyses of ongoing experiments
based on cold-atom quantum simulations of the Fermi-Hubbard model.

##  I Introduction

The Fermi-Hubbard Hamiltonian is a paradigmatic model in condensed matter
physics, introduced to describe the behavior of electrons in a solid Hubbard
(1963). It supports a remarkably broad spectrum of quantum phases of matter,
and it is believed to capture the essential phenomenology of strongly
correlated materials including the cuprates Lee  _et al._ (2006). Yet, the
Fermi-Hubbard model has proven extremely difficult to analyze and continues to
challenge theoretical and numerical efforts for more than four decades
Schmitt-Rink  _et al._ (1988); Shraiman and Siggia (1988); Kane  _et al._
(1989); Martinez and Horsch (1991); Liu and Manousakis (1991); Emery (1987);
Schrieffer  _et al._ (1988); Dagotto (1994); Anderson (1987); Bonča  _et al._
(1989); Hasegawa and Poilblanc (1989); Dagotto  _et al._ (1990); Sachdev
(1989); Trugman (1990); Boninsegni and Manousakis (1992); Brunner  _et al._
(2000); Mishchenko  _et al._ (2001); Blomquist and Carlström (2020); White and
Affleck (2001); Zhu  _et al._ (2014); Wang  _et al._ (2021); Chen  _et al._
(2021); Bulaevskii  _et al._ (1968); Brinkman and Rice (1970); Trugman (1988);
Manousakis (2007); Grusdt  _et al._ (2018a, b, 2019); Bohrdt  _et al._ (2019,
2021); Soriano and Katsnelson (2020). An important case emerges close to half
filling where each lattice site is occupied by one fermion. Then, strong on-
site particle repulsion leads to the build-up of antiferromagnetic order of
the spins of the fermions, which competes with the delocalization of holes
that can be present in the lattice Chao  _et al._ (1977); Hirsch (1985);
Izyumov (1997). This results in a buildup of magnetic frustrations around such
holes and the formation of quasiparticles, termed magnetic polarons Schmitt-
Rink  _et al._ (1988); Shraiman and Siggia (1988); Kane  _et al._ (1989);
Martinez and Horsch (1991); Liu and Manousakis (1991). The emerging magnetic
dressing cloud induces effective interactions between two such holes that have
been conjectured to provide a mechanism for high-temperature superconductivity
Emery (1987); Schrieffer  _et al._ (1988); Dagotto (1994). Understanding and
characterizing magnetic polarons has therefore been of key interest for many
decades.

Owing to the shear complexity of the problem, only a few theoretical
approaches have been applied under different conditions and with varying
success. This includes exact diagonalization for small system sizes Bonča  _et
al._ (1989); Hasegawa and Poilblanc (1989); Dagotto  _et al._ (1990); Wang
_et al._ (2021), mean field approaches Anderson (1987), and variational
calculations Sachdev (1989); Trugman (1990). One has also analysed the string-
excitations caused by defect motion through the magnet Bulaevskii  _et al._
(1968); Brinkman and Rice (1970); Trugman (1988); Manousakis (2007); Grusdt
_et al._ (2018a, b, 2019); Bohrdt  _et al._ (2019), and employed numerical
techniques such as Monte-Carlo simulations Boninsegni and Manousakis (1992);
Brunner  _et al._ (2000); Mishchenko  _et al._ (2001); Blomquist and Carlström
(2020), machine learning methods Bohrdt  _et al._ (2019) as well as
renormalization group techniques White and Affleck (2001); Zhu  _et al._
(2014); Wang  _et al._ (2021); Chen  _et al._ (2021).

Figure 1: Magnetization around a hole in a Heisenberg antiferromagnet in a 16
by 16 square lattice for different interaction strengths and crystal momenta
along the magnetic Brillouin zone. As the hopping amplitude $t$ becomes large
compared to the spin-spin coupling $J$, moving from panels (a) to (b), the
magnetic frustrations around the hole increase in size and magnitude and the
magnetic order even flips sign at the nearest neighbor sites. The dashed lines
indicate the reflection symmetries of the magnetization, which along with a
$C_{2}$ rotation symmetry means that the dressing cloud has the remarkably
high symmetry group $C_{2v}$. As a consequence, the orientation of the
dressing cloud is in general misaligned with the crystal momentum leading to a
non-trivial behavior.

Recent experimental breakthroughs in manipulating ultracold atoms in optical
lattices have opened up the possibility to perform quantum simulations of the
Fermi-Hubbard model Esslinger (2010); Boll  _et al._ (2016); Cheuk  _et al._
(2016); Mazurenko  _et al._ (2017); Hilker  _et al._ (2017); Brown  _et al._
(2017); Chiu  _et al._ (2018); Brown  _et al._ (2019); Koepsell  _et al._
(2019); Chiu  _et al._ (2019); Koepsell  _et al._ (2020); Brown  _et al._
(2020); Vijayan  _et al._ (2020); Hartke  _et al._ (2020); Guardado-Sanchez
_et al._ (2020); Ji  _et al._ (2021); Gall  _et al._ (2021). In particular,
the ability to image individual atoms with single-site spatial resolution Bakr
_et al._ (2009); Sherson  _et al._ (2010); Haller  _et al._ (2015); Yang  _et
al._ (2021) makes it possible to probe the microscopic structure of the
magnetic polaron Koepsell  _et al._ (2019); Chiu  _et al._ (2019). Such
detailed insights offer stringent tests of the understanding of the Fermi-
Hubbard model, and enable a systematic improvement of theoretical approaches
to these quasiparticles.

A particularly successful approach has been the self-consistent Born
approximation (SCBA) Schmitt-Rink  _et al._ (1988); Kane  _et al._ (1989). The
SCBA permits a non-perturbative calculation of the Green’s function of the
hole and was shown to yield quantitatively accurate results Martinez and
Horsch (1991); Liu and Manousakis (1991, 1992); Chernyshev and Leung (1999) in
the strongly interacting regime. Despite this success, its applicability has
thus far been limited to single-particle observables such as the energy
dispersion and quasiparticle residue of the magnetic polaron, while the
extraction of finer structural information, such as that contained in spin-
hole correlations or the polaron’s magnetic dressing cloud, have proved
difficult. The construction of the polaron wave function in terms of spin-wave
excitations Reiter (1994) in principle offers a solution to this problem.
Determining correlation functions, however, entails an infinite series of
terms with increasing number of spin excitations Ramšak and Horsch (1993,
1998), whose truncation has restricted calculations to the weak coupling
domain Ramšak and Horsch (1993, 1998); Bała and Oleś (2000); Bała  _et al._
(2002). Moreover, these conditions violate the underlying assumptions of the
SCBA treatment Schmitt-Rink  _et al._ (1988); Kane  _et al._ (1989) in the
context of the Fermi-Hubbard model.

Here, we develop a theoretical framework that makes it possible to overcome
this obstacle. The approach is based on a set of self-consistency equations
that are reminiscent of the Dyson equation, and enable the inclusion of _all_
terms up to infinite numbers of spin excitations in the SCBA wave function for
the magnetic polaron. It thereby extends its application into the regime of
strong interactions. We use this new approach to explore the microscopic
spatial structure of the magnetic polaron in the regime of strong coupling. In
general, the obtained magnetic dressing cloud has an elongated shape that
increases in size and magnitude with the strength of interactions, and, in the
strong-coupling regime, differs qualitatively from previous calculations based
on a truncated quasiparticle wave function. Our analysis reveals that the
symmetries of the antiferromagnetic spin lattice decisively determine the form
of the dressing cloud, and lead to a remarkably high symmetry of the magnetic
dressing cloud for polaron momenta along the edge of the magnetic Brillouin
zone (MBZ). Surprisingly and in contrast to previous expectation, this can
lead to a misalignment between the dressing cloud and the crystal momentum,
such that the spatial structure of the polaron is generally not oriented along
its direction of motion.

Our theoretical framework moreover permits to explore the transition of the
underlying quantum magnet from the isotropic Heisenberg spin-lattice to the
Ising model. We find that the gap opening in the spin-wave excitation spectrum
results in a shrinking of the polaronic magnetization cloud with increasing
anisotropy of the effective spin-spin interaction. In the Ising limit, the
full symmetry of the antiferromagnetic order is restored, which makes it
possible to determine the dressing cloud of the magnetic polaron analytically.
More generally, the developed framework may open up a new approach for
microscopic explorations of doped quantum magnets in the strong-coupling
regime, including finite temperature effects and non-equilibrium dynamics to
induced interactions between multiple defects.

This article is organised as follows. Section II provides the Fermi-Hubbard
Hamiltonian and the $t$-$J$ model that derives from it. Based on the $t$-$J$
model, we describe the transformation into a magnetic polaron Hamiltonian
within linear spin wave theory. In Sec. III, we summarize the quasiparticle
properties of the magnetic polaron, including the calculation of the polaron
Green’s and wave functions within the SCBA. The magnetization in the vicinity
of a hole is explored in Sec. IV for a two-dimensional square lattice as a
function of interaction strength and anisotropy. In Sec. V, we demonstrate the
non-perturbative effects predicted by our developed formalism, and Sec. VI
gives a detailed derivation of the self-consistency equations for the local
magnetization. Finally, in Sec. VII, we describe the prospects of testing our
theory experimentally.

##  II The anisotropic $t$-$J$ model

The Fermi-Hubbard model

$\hat{H}_{\rm FH}=-t\sum_{\braket{{\bf i},{\bf
j}},\sigma}\left[\hat{c}^{\dagger}_{{\bf i},\sigma}\hat{c}_{{\bf
j},\sigma}+{\rm H.c.}\right]+U\sum_{\bf i}\hat{n}_{{\bf
i},\uparrow}\hat{n}_{{\bf i},\downarrow}$ (1)

describes spin-1/2 fermions moving in a lattice with hopping amplitude $t$ and
onsite repulsive interactions $U>0$. Here, $\braket{{\bf i},{\bf j}}$ denotes
nearest neighbor lattice sites, $\hat{c}^{\dagger}_{{\bf i},\sigma}$ creates a
fermion at site $\bf i$ and spin $\sigma$, while $\hat{n}_{{\bf
i},\sigma}=\hat{c}^{\dagger}_{{\bf i},\sigma}\hat{c}_{{\bf i},\sigma}$ is the
corresponding counting operator. Despite its apparent simplicity, many open
questions remain concerning its properties. Near half filling, one can expand
the Hubbard model in the particle hopping for large onsite repulsion $U\gg t$
to derive an effective low-energy description given by the so-called $t$-$J$
model Chao  _et al._ (1977); Dagotto (1994); Izyumov (1997). The Hamiltonian
is $\hat{H}=\hat{H}_{t}+\hat{H}_{J}$, where

$\hat{H}_{t}=-t\sum_{\braket{{\bf i},{\bf
j}},\sigma}\left[\tilde{c}^{\dagger}_{{\bf i},\sigma}\tilde{c}_{{\bf
i},\sigma}+{\rm H.c.}\right],$ (2)

describes the restrained nearest neighbor particle hopping, where
$\tilde{c}^{\dagger}_{{\bf j},\sigma}=\hat{c}_{{\bf
j},\sigma}^{\dagger}(1-\hat{n}_{{\bf j},\bar{\sigma}})$, and the factor
$(1-\hat{n}_{{\bf j},\bar{\sigma}})$ with $\bar{\sigma}$ the opposite spin
restrains the Hilbert space of the model to states with maximally one particle
per site. Furthermore,

$\hat{H}_{J}=J\sum_{\braket{{\bf i},{\bf j}}}\\!\left[\hat{S}^{z}_{\bf
i}\hat{S}^{z}_{\bf j}\\!+\\!\frac{\alpha}{2}\\!\left(\hat{S}^{+}_{\bf
i}\hat{S}^{-}_{\bf j}\\!+\\!\hat{S}^{-}_{\bf i}\hat{S}^{+}_{\bf
j}\right)\\!-\\!\frac{\hat{n}_{\bf i}\hat{n}_{\bf j}}{4}\right]$ (3)

gives the antiferromagnetic ($J>0$) spin-spin interactions. The Schwinger-
fermion representation of spin $1/2$ reads as usual

$\mathbf{S}_{\bf j}=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}}\hat{c}_{{\bf
j},\sigma}^{\dagger}\boldsymbol{\sigma}_{\sigma\sigma^{\prime}}\hat{c}_{{\bf
j},\sigma^{\prime}}$ (4)

with $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ a vector of the
Pauli matrices. The $t$-$J$ model with $\alpha=1$ thus yields an accurate
description of the low-energy physics of the underlying Fermi-Hubbard model
for $t\gg J=4t^{2}/U$ close to half filling. More generally, other
experimental platforms Porras and Cirac (2004); Gorshkov  _et al._ (2011);
Britton  _et al._ (2012); Zeiher  _et al._ (2016, 2017) make it possible to
tune from the isotropic Heisenberg limit with $\alpha=1$ to an Ising magnet
with $\alpha=0$. In Eq. (3), we have therefore generalized the model to
include the case of anisotropic spin interactions by introducing the parameter
$0\leq\alpha\leq 1$.

At half filling with exactly one fermion per lattice site, the first term in
Eq. (2) is ineffective and a positive superexchange coupling $J>0$ between the
spins enforces antiferromagnetic ordering for any value of $\alpha$. Lattice
defects, or holes, in such an antiferromagnet tend to delocalize and thereby
lower their kinetic energy, as given by $\hat{H}_{t}$. The associated motion
of holes, on the other hand, leads to the buildup of magnetic frustration
which increases the energy of the system according to $\hat{H}_{J}$. The
competition between these two processes eventually gives rise to the magnetic
polaron, i.e. a mobile hole that is surrounded by a finite magnetization
cloud. Small ratios of $J/t$, thus, correspond to the strong coupling regime
in which a high hole mobility leads to a significant disturbance of its
magnetic environment and thereby generates strong spin-hole correlations. We
can accurately describe this process using spin-wave theory as outlined in the
next section.

###  II.1 Slave fermion representation

We begin by performing a Holstein-Primakoff transformation generalized to take
into account the presence of holes Schmitt-Rink  _et al._ (1988); Kane  _et
al._ (1989); Martinez and Horsch (1991); Liu and Manousakis (1991). The
antiferromagnetic state defines a bipartite lattice, whereby one sublattice
carries fermions in the spin-up state, while the other sublattice is formed by
particles in the spin-down state. In the former, we rewrite $\hat{S}_{\bf
i}^{-}=\hat{s}_{\bf i}^{\dagger}(1-\hat{s}_{\bf i}^{\dagger}\hat{s}_{\bf
i}-\hat{h}^{\dagger}_{\bf i}\hat{h}_{\bf i})^{1/2}$, $\tilde{c}_{{\bf
i},\downarrow}=\hat{h}^{\dagger}_{\bf i}\hat{s}_{\bf i}$, and $\tilde{c}_{{\bf
i},\uparrow}=\hat{h}^{\dagger}_{\bf i}(1-\hat{s}^{\dagger}_{\bf i}\hat{s}_{\bf
i}-\hat{h}^{\dagger}_{\bf i}\hat{h}_{\bf i})^{1/2}$ in terms of fermionic
operators $\hat{h}^{\dagger}_{\bf i}$ and bosonic operators
$\hat{s}^{\dagger}_{\bf i}$ that create a hole and a spin excitation at site
${\bf i}$ respectively. The factor $(1-\hat{s}^{\dagger}_{\bf i}\hat{s}_{\bf
i}-\hat{h}^{\dagger}_{\bf i}\hat{h}_{\bf i})^{1/2}$ ensures that there is at
most one hole or one spin excitation at each site. Finally, the spin-$z$
operator can be rewritten as $\hat{S}_{\bf i}^{z}=(1-\hat{h}_{\bf
i}^{\dagger}\hat{h}_{\bf i})/2-\hat{s}^{\dagger}_{\bf i}\hat{s}_{\bf i}$. The
representation of the spin and holes on the other sublattice of spin-down
fermions, proceeds analogously by swapping spin $\uparrow$ and $\downarrow$ in
the transformations given above. Using this so-called slave-fermion
representation in Eq. (3), keeping only the linear terms, and diagonalizing
the transformed Hamiltonian yields the spin wave Hamiltonian Schmitt-Rink  _et
al._ (1988); Kane  _et al._ (1989); Martinez and Horsch (1991); Liu and
Manousakis (1991)

$\displaystyle H_{J}=E_{0}+\sum_{\bf k}\omega_{\bf k}\hat{b}^{\dagger}_{\bf
k}\hat{b}_{\bf k}.$ (5)

This describes spin waves with the energy $\omega_{\bf
k}=zJ\sqrt{1-\alpha^{2}\gamma_{\bf k}^{2}}/2$, where the structure factor

$\displaystyle\gamma_{\bf k}=\frac{1}{z}\sum_{{\boldsymbol{\delta}}}{\rm
e}^{i{\bf k}\cdot{\boldsymbol{\delta}}},$ (6)

is the sum of the $z$ nearest neighbor phases. The associated bosonic spin
wave operators $\hat{b}_{\bf k}$ are related to the physical spin excitations
$\hat{s}_{\bf k}$ via a Bogoliubov transformation $\hat{b}_{\bf k}=u_{\bf
k}\hat{s}_{\bf k}+v_{\bf k}\hat{s}^{\dagger}_{-{\bf k}}$ with the
antiferromagnetic coherence factors given by $u_{\bf
k}=[(1/\sqrt{1-\alpha^{2}\gamma_{\bf k}^{2}}+1)/2]^{1/2}$ and $v_{\bf k}={\rm
sgn}(\gamma_{\bf k})[(1/\sqrt{1-\alpha^{2}\gamma_{\bf k}^{2}}-1)/2]^{1/2}$.
Using the slave fermion representation in Eq. (2) yields Schmitt-Rink  _et
al._ (1988); Kane  _et al._ (1989); Martinez and Horsch (1991); Liu and
Manousakis (1991)

$\displaystyle\hat{H}_{t}=\sum_{{\bf q},{\bf k}}\hat{h}_{{\bf q}+{\bf
k}}^{\dagger}\hat{h}_{{\bf q}}\left[g({\bf q},{\bf k})\hat{b}^{\dagger}_{-{\bf
k}}+g({\bf q}+{\bf k},-{\bf k})\hat{b}_{\bf k}\right],$ (7)

where the interaction vertex is $g({\bf q},{\bf k})=zt\cdot(u_{\bf
k}\gamma_{{\bf q}+{\bf k}}-v_{\bf k}\gamma_{\bf q})/\sqrt{N}$ with $N$ the
number of lattice sites. Here, we retain terms linear in the spin-wave
operators. Equation (7) explicitly shows how the hopping of a hole gives rise
to the emission/absorption of spin waves and directly represents the
competition between hole delocalization and magnetic order.

##  III Magnetic polaron

The above formulation in terms of interactions with magnetic spin-wave
excitations enables a non-perturbative description of the microscopic
structure of the magnetic polaron, as we shall outline in this section. The
approach is based on the self-consistent Born approximation (SCBA) for the
hole Green’s function Schmitt-Rink  _et al._ (1988); Kane  _et al._ (1989),
which has been shown to yield quantitatively accurate results for the polaron
energy in the Heisenberg limit across all interaction strengths Martinez and
Horsch (1991).

###  III.1 The Green’s function

The SCBA includes the so-called rainbow diagrams in the computation of the
hole Green’s function $G({\bf p},\omega)=1/(\omega-\Sigma({\bf
p},\omega)+i\eta)$ Schmitt-Rink  _et al._ (1988); Kane  _et al._ (1989);
Martinez and Horsch (1991); Liu and Manousakis (1991), where $\eta=0^{+}$ is a
positive infinitesimal. Using the spin wave Green’s function $G_{b}({\bf
k},\omega)=1/(\omega-\omega_{\bf k}+i\eta)$, the diagrammatic structure shown
in Fig. 2(a) leads to the self-consistent equation for the self-energy

$\Sigma({\bf p},\omega)=\sum_{\bf k}\frac{g^{2}({\bf p},{\bf
k})}{\omega-\omega_{\bf k}-\Sigma({\bf p}+{\bf k},\omega-\omega_{\bf
k})+i\eta},$ (8)

which can be solved iteratively starting from $\Sigma=0$. Knowing the self-
energy, one can determine several important quantities, such as the
quasiparticle residue

$\displaystyle Z_{\bf p}=\frac{1}{1-\partial_{\omega}\Sigma({\bf
p},\omega)|_{\omega=\epsilon_{\bf p}}},$ (9)

which is the overlap $Z_{\bf p}=|\bra{{\rm AF}}\hat{h}_{{\bf p}}\ket{\Psi_{\bf
p}}|^{2}$ of the polaron many-body wave function $\ket{\Psi_{\bf p}}$ with the
state of a bare hole in an otherwise unperturbed antiferromagnetic state,
$\hat{h}_{{\bf p}}^{\dagger}\ket{{\rm AF}}$, for a given crystal momentum
${\bf p}$ of the hole. Hereby, the antiferromagnetic quantum Néel state is
defined as $\hat{b}_{{\bf k}}\ket{{\rm AF}}=0$ for any spin-wave momentum
${\bf k}$. The numerical solution of Eq. (8) under strong-coupling conditions,
$J/t=0.3$, in a 36 by 36 square lattice is shown in Fig. 2. The hole spectral
function, depicted in Fig. 2(b), exhibits a clear quasiparticle peak at
$\omega\simeq-2.4$t giving the energy of the magnetic polaron. The
corresponding quasiparticle residue is $Z\simeq 0.3$ and there is a continuum
of many-body states at higher energies with large spectral weight, reflecting
the strongly interacting nature of the problem. The dispersion of the magnetic
polaron in the first Brillouin zone is shown in Fig. 2(d). It features four
degenerate ground states at the crystal momenta ${\bf p}=(\pm\pi/2,\pm\pi/2)$,
given in units of the inverse lattice constant. The surprising predictive
power of the SCBA result for the hole spectral function in the Heisenberg
limit Martinez and Horsch (1991) compared to exact diagonalization studies
Dagotto  _et al._ (1990) is attributed to small vertex corrections to the SCBA
even for strong coupling Liu and Manousakis (1992). The dependence of
quasiparticle residue on $J/t$ [Fig. 2(c)] clearly illustrates the necessity
of a non-perturbative theory for $J/t\ll 1$, since small values of $Z_{{\bf
p}}$ indicate a large number of spin-wave excitations and strong correlations
between the generated spin fluctuations and the motion of the hole. The
Green’s function alone is, however, not well suited to study such correlations
and requires additional analysis as we will now discuss.

Figure 2: Quasiparticle properties. (a) Dyson equation for the hole Green’s
function, $G({\bf p},\omega)$, within SCBA (double black line), in terms of
the non-interacting Green’s function $G_{0}(\omega)=1/\omega$ (single black
line), and the spin wave Green’s function $G_{b}=1/(\omega-\omega_{\bf k})$
(blue wiggly line). (b) Hole spectral function $A({\bf p},\omega)=-2{\rm
Im}G({\bf p},\omega)$ in the Heisenberg limit in a 36 by 36 square lattice,
featuring a quasiparticle peak at $\varepsilon_{\bf p}\simeq-2.4t$ for
$J/t=0.3$ (vertical red line). For $J/t\ll 1$, the quasiparticle residue
[panel (c)] becomes very small. Consequently, an increasing number of terms
must be retained in the quasiparticle wave function [panel (e)]. The
diagrammatic rules for the construction of the wave function are: (1) Single
straight line: $Z_{\bf p}^{1/2}$. (2) Blue wiggly line with momentum $-{\bf
k}_{i}$: spin wave operator, $\hat{b}^{\dagger}_{-{\bf k}_{i}}$. (3) nth red
dot from the right in a diagram: interaction vertex, $g({\bf K}_{n},{\bf
k}_{n+1})$. Here, ${\bf K}_{0}={\bf p}$, ${\bf K}_{1}={\bf p}+{\bf k}_{1}$,
${\bf K}_{2}={\bf p}+{\bf k}_{1}+{\bf k}_{2}$, and so forth. (4) nth double
line from the right: $G\left({\bf K}_{n},\varepsilon_{\bf
p}-\sum_{i=1}^{n}\omega_{{\bf k}_{i}}\right)$. (5) nth open-ended hole line
from the right: hole operator, $\hat{h}^{\dagger}_{{\bf K}_{n}}$. (6) Sum over
all spin wave momenta, ${\bf k}_{i}$. (d) The quasiparticle dispersion
throughout the first Brillouin zone with the four degenerate ground states at
${\bf p}=(\pm\pi/2,\pm\pi/2)$ indicated by circles. The magnetic Brillouin
zone (MBZ) is indicated by black diagonal lines.

###  III.2 The polaron wave function

Our non-perturbative approach to determine spin-hole correlations exploits the
fact that it is formally possible to write the wave function $\ket{\Psi_{\bf
p}}$ of the magnetic polaron within the SCBA. Explicitly, the wave function
Reiter (1994); Ramšak and Horsch (1998)

$\displaystyle\ket{\Psi_{\bf p}}=\sum_{n=0}^{\infty}\sum_{\\{{\bf
k}_{i}\\}}a^{(n)}({\bf p},\\{{\bf
k}_{i}\\}_{i=1}^{n})\cdot\hat{h}^{\dagger}_{{\bf
K}_{n}}\prod_{i=1}^{n}\hat{b}_{-{\bf k}_{i}}^{\dagger}\ket{{\rm AF}}$
$\displaystyle=\sqrt{Z_{\bf p}}\Big{[}h^{\dagger}_{\bf p}\\!+\\!\sum_{{\bf
k}_{1}}g({\bf p},{\bf k}_{1})G\left({\bf p}\\!+\\!{\bf k}_{1},\varepsilon_{\bf
p}\\!-\\!\omega_{{\bf k}_{1}}\right)\hat{h}^{\dagger}_{{\bf p}+{\bf
k}_{1}}\hat{b}^{\dagger}_{-{\bf k}_{1}}$ $\displaystyle\phantom{=\sqrt{Z_{\bf
p}}\Big{[}}+\dots\Big{]}\ket{{\rm AF}},$ (10)

can be expressed as an expansion in the number of spin wave excitations on the
antiferromagnetic quantum Néel state $\ket{{\rm AF}}$. The lowest order
coefficient, $a^{(0)}({\bf p})=\sqrt{Z_{\bf p}}$, is given by the square root
of the quasiparticle residue ensuring the overall normalization of the wave
function, while the higher order coefficients can be computed from the
recurrence relation

$\displaystyle a^{(n+1)}({\bf p},\\{{\bf k}_{i}\\}_{i=1}^{n+1})=$
$\displaystyle g({\bf K}_{n},{\bf k}_{n+1})a^{(n)}({\bf p},\\{{\bf
k}_{i}\\}_{i=1}^{n})$ $\displaystyle G\left({\bf K}_{n+1},\varepsilon_{\bf
p}-\sum_{i=1}^{n}\omega_{{\bf k}_{i}}\right),$ (11)

with ${\bf K}_{n}={\bf p}+\sum_{i=1}^{n}{\bf k}_{i}$ for $n\geq 1$, and ${\bf
K}_{0}={\bf p}$. The structure of the first order term $a^{(1)}$ can be
understood by using the recursion relation in Eq. (11) in the coupling between
the two lowest order coefficients, $a^{(0)}$ and $a^{(1)}$, of the
quasiparticle wave function (10). This yields $\varepsilon_{\bf p}a^{(0)}({\bf
p})=\sum_{{\bf k}}g({\bf p},{\bf k})a^{(1)}({\bf p},{\bf k})=\Sigma({\bf
p},\varepsilon_{\bf p})a^{(0)}({\bf p})$, by using Eq. (8) for the self-energy
$\Sigma$. Hence, the construction of the wave function (10) relies on the
presence of a well-defined quasiparticle peak at $\varepsilon_{\bf p}$
determined by $\varepsilon_{\bf p}=\Sigma({\bf p},\varepsilon_{\bf p})$,
corresponding to the energy of the magnetic polaron.

The wave function is visualized diagrammatically Ramšak and Horsch (1998) in
Fig. 2(e), along with the diagrammatic rules for its construction. The
iterative structure of Fig. 2(e) together with precise diagrammatic rules is
reminiscent of the Dyson equation in quantum field theory [Fig. 2(a)], and is
at the heart of our non-perturbative framework developed below. A major
advantage of the polaron wave function is that it allows for the computation
of spin-hole correlation functions in a much more direct way than the hole
Green’s function. The problem is nevertheless still far from straightforward,
since the number of important terms in Eq. (10) increases with the
interaction, $t/J$. As a consequence, there is no controlled way to truncate
the series for the wave function, while still obtaining reliable results in
the strong coupling regime $J/t\ll 1$. This has so far limited the use of this
wave function to the weak coupling regime, where it is sufficient to include
only a small of number spin wave excitations Ramšak and Horsch (1993, 1998).
However, since the SCBA is least accurate precisely in the weak coupling
regime of $J/t\gtrsim 1$, the validity of such an approach is not clear. The
non-perturbative framework, developed in the present work, thus represents a
major step as it now makes it possible to utilize the SCBA for reliable
calculations of correlation functions in the strongly coupled regime.

##  IV Magnetization around a hole

Before presenting the derivation of our non-perturbative approach in Sec. VI,
in this section, we first illustrate its application by calculating the local
magnetization in the neighborhood of a hole. This is a fundamental property of
the magnetic polaron that determines its microscopic structure, and eventually
the form of the induced interaction between multiple holes.

Consider first the magnetization in the absence of holes. The local
magnetization at a given lattice site ${\bf d}$ is given by

$\\!\\!2\bra{{\rm AF}}\\!\hat{S}^{z}_{\bf d}\\!\ket{{\rm AF}}=(-1)^{l_{\bf
d}}M_{\text{AF}}=(-1)^{l_{\bf
d}}\\!\left(1\\!-\\!2M_{\text{fl}}\right),\\!\\!$ (12)

where $0\leq M_{\text{fl}}\leq 1$ quantifies the effect of quantum spin
fluctuations to suppress the magnetic order from its maximum value
$M_{\text{AF}}=1$ in the Ising limit $\alpha=0$. Without loss of generality,
we take $\bra{{\rm AF}}\hat{S}^{z}_{{\bf d}={\bf 0}}\ket{{\rm AF}}>0$ in the
spontaneously broken symmetry state $\ket{{\rm AF}}$. The site-dependent
parameter $l_{\bf d}$ is defined as the minimal number of lattice links
between the two sites ${\bf r}=\mathbf{0}$ and ${\bf r}={\bf d}$. From linear
spin-wave theory, we have $\hat{S}^{z}_{{\bf d}}=(-1)^{l_{\bf
d}}(1/2-\hat{s}^{\dagger}_{{\bf d}}\hat{s}_{{\bf d}})$ and
$\hat{s}^{\dagger}_{{\bf d}}=\sum_{\bf k}{\rm e}^{-i{\bf k}\cdot{\bf
d}}\hat{s}_{{\bf k}}^{\dagger}/\sqrt{N}$, which gives $M_{\text{fl}}=\sum_{\bf
k}\bra{{\rm AF}}\hat{s}_{{\bf k}}^{\dagger}\hat{s}_{{\bf k}}\ket{{\rm
AF}}/N=\sum_{{\bf k}}v_{\bf k}^{2}/N$.

With these definitions, we can now formulate the magnetization of a given
lattice site at a distance ${\bf d}$ from the hole

$\displaystyle M_{\bf p}({\bf d})=\frac{\braket{\hat{h}^{\dagger}_{{\bf
r}}\hat{h}_{{\bf r}}\hat{S}^{z}_{{\bf r}+{\bf d}}}_{\bf
p}}{\braket{\hat{h}^{\dagger}_{{\bf r}}\hat{h}_{{\bf r}}}_{\bf
p}\braket{\hat{S}^{z}_{{\bf r}+{\bf d}}}_{\bf p}}=\frac{1-2M_{\bf
p}^{(2)}({\bf d})}{M_{\text{AF}}},$ (13)

where $\braket{\ldots}_{\bf p}=\bra{\Psi_{\bf p}}\ldots\ket{\Psi_{\bf p}}$ is
the expectation value for the ground state Eq. (10) of a magnetic polaron with
crystal momentum ${\bf p}$, and

$\displaystyle M_{\bf p}^{(2)}({\bf d})=N\braket{\hat{h}^{\dagger}_{{\bf
r}}\hat{h}_{{\bf r}}\hat{s}^{\dagger}_{{\bf r}+{\bf d}}\hat{s}_{{\bf r}+{\bf
d}}}_{\bf p}.$ (14)

The translational symmetry of the system ensures that these correlation
functions only depend on the distance vector ${\bf d}$. In Eqs. (13) and (14),
we have used $\braket{\hat{h}^{\dagger}_{{\bf r}}\hat{h}_{{\bf r}}}_{\bf
p}=1/N$, reflecting the fact that the hole is equally distributed across the
lattice for a given momentum state. We can also omit corrections to the
average magnetization from the presence of a single hole, since they scale as
$\braket{\hat{S}^{z}_{{\bf d}}}_{\bf p}-\bra{{\rm AF}}\hat{S}^{z}_{{\bf
d}}\ket{{\rm AF}}\sim{\mathcal{O}}(1/N)$.

Fourier transforming and rotating to the bosonic spin wave operators
$\hat{b}_{\bf k}$, the remaining two-point correlator $M_{\bf p}^{(2)}$ given
by Eq. (14) can be decomposed as

$\displaystyle M^{(2)}_{\bf p}({\bf d})=M_{\text{fl}}+B_{\bf p}({\bf
d})+C_{\bf p}({\bf d}).$ (15)

Notice that the zero-point fluctuations $M_{\text{fl}}$ of the quantum
antiferromagnet appears explicitly. The corrections to the magnetization due
to the presence of the hole are thus described by the functions

$\displaystyle\\!\\!\\!\\!B_{\bf p}({\bf d})=\frac{1}{N^{2}}\\!\\!\sum_{{\bf
q}_{1},{\bf q}_{2}}{\rm e}^{-i({\bf q}_{2}-{\bf q}_{1})\cdot{\bf
d}}\left[u_{{\bf q}_{1}}u_{{\bf q}_{2}}+v_{{\bf q}_{1}}v_{{\bf q}_{2}}\right]$
$\displaystyle\phantom{\\!\\!\\!\\!B_{\bf p}({\bf
d})=\frac{1}{N^{2}}\\!\\!\sum_{{\bf q}_{1},{\bf q}_{2}}}\cdot B({\bf
q}_{1},{\bf q}_{2};{\bf p},\varepsilon_{\bf p}),$
$\displaystyle\\!\\!\\!\\!B({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=N\sum_{\bf
k}\braket{\hat{h}^{\dagger}_{{\bf k}+{\bf q}_{1}}\hat{h}_{{\bf k}+{\bf
q}_{2}}\hat{b}^{\dagger}_{-{\bf q}_{1}}\hat{b}_{-{\bf q}_{2}}}_{{\bf
p},\omega},$ (16)

and

$\displaystyle\\!\\!\\!\\!C_{\bf p}({\bf d})=-\frac{1}{2N^{2}}\\!\\!\sum_{{\bf
q}_{1},{\bf q}_{2}}{\rm e}^{-i({\bf q}_{2}+{\bf q}_{1})\cdot{\bf
d}}\left[u_{{\bf q}_{1}}v_{{\bf q}_{2}}+v_{{\bf q}_{1}}u_{{\bf q}_{2}}\right]$
$\displaystyle\phantom{\\!\\!\\!\\!C({\bf d})=-\frac{1}{N^{2}}\\!\\!\sum_{{\bf
q}_{1},{\bf q}_{2}}}\cdot C({\bf q}_{1},{\bf q}_{2};{\bf p},\epsilon_{\bf
p})+\textrm{c.c.}$ $\displaystyle\\!\\!\\!\\!C({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=N\sum_{{\bf k}}\braket{\hat{h}^{\dagger}_{{\bf k}-{\bf q}_{1}-{\bf
q}_{2}}\hat{h}_{\bf k}\hat{b}_{-{\bf q}_{1}}\hat{b}_{-{\bf q}_{2}}}_{{\bf
p},\omega},$ (17)

where c.c. stands for the complex conjugate. In the definition of $B_{\bf
p}(\cdot,\omega)$ and $C_{\bf p}(\cdot,\omega)$, we allow the energy $\omega$
of the magnetic polaron to vary, meaning that the expectation value
$\braket{\ldots}_{{\bf p},\omega}$ is taken with respect to the state
$\ket{\Psi_{\bf p}}$ given by Eq. (10), where the polaron energy
$\varepsilon_{\bf p}$ is replaced by $\omega$ in the appearing Green’s
functions. As we will describe in Sec. VI, this generalization combined with
the diagrammatic rules for the wave function given in Fig. 2 makes it possible
to derive self-consistency equations for the correlation functions in Eqs.
(16) and (17). These equations are similar to Eq. (8) for the self-energy, and
ultimately facilitate evaluation of correlation functions to _all_ orders in
the number of spin-wave excitations in the wave function Eq. (10).

###  IV.1 Heisenberg limit

In Figs. 1(a) and 1(b), we show the spatial structure of the local
magnetization $\tilde{M}_{\bf p}({\bf d})=2\braket{\hat{S}^{z}_{{\bf d}}}_{\bf
p}M_{\bf p}({\bf d})$ in the isotropic Heisenberg limit, $\alpha\to 1$, for
two different coupling strengths. The crystal momentum is ${\bf
p}=(\pi/2,\pi/2)$, which corresponds to one of the four degenerate ground
states of the polaron [Fig. 2(d)]. The results show that the hole
significantly reduces the magnetic order in its vicinity. This effect is
particularly prominent for strong coupling, $J/t=0.05$, where the size of the
magnetization cloud is increased, and where the nearest neighbor spins are
even flipped as a consequence of strong correlations between the motion of the
hole and the local magnetization of the lattice.

Figure 3: Magnetization in the neighborhood of a hole for the ground state
magnetic polaron with ${\bf p}=(\pi/2,\pi/2)$ in a Heisenberg antiferromagnet
as a function of interaction strength, $J/t$, in a 16 by 16 square lattice.
(a) Red squares are for the nearest neighbor sites, blue dots are the next
nearest neighbors, and green triangles are for the perpendicular direction as
illustrated in the inset. The lines are guides to the eye. Note that the
magnetization at the green points are larger than in the absence of holes for
intermediate interaction strengths, $J/t>0.3$. (b) Heat maps of the
magnetization for a strong and intermediate interaction strength, $J/t=0.05$
(left) and $J/t=0.5$ (right) respectively.

The emergence of this sign flip is illustrated in Fig. 3(a), where we show the
magnetization $M_{{\bf p}}({\bf d})$ as a function of the inverse interaction
strength $J/t$, indicating that this strong magnetic disturbance extends to
larger and larger distances as we enter the strong coupling regime and further
increase $t/J$. The results also show that the magnetization of the next-
nearest neighbors is anisotropic and reflects the direction of the crystal
momentum of the moving polaron. Perpendicular to this direction, the
magnetization can even be larger than in the absence of the hole. This
surprising effect results from the coherent addition of the generated spin
waves to produce a net increase in the magnetization for intermediate
interaction strengths.

The elongated shape of the magnetization cloud is shown more directly in Fig.
3(b). It appears that the magnetization cloud is oriented along the crystal
momentum of the polaron, as also reported previously based on truncated wave
function calculations. There, the observed alignment has been attributed to
the semiclassical idea that the hole will predominantly disturb the
magnetization in the direction of its motion Ramšak and Horsch (1993, 1998).
However, as we will discuss below, this is _generally not the case_ , as the
symmetry properties of the underlying antiferromagnetic competes with the
directed motion of the polaron and yields a nontrivial orientation of the
magnetization cloud with respect to the polaron crystal momentum.

While the SCBA can in principle lead to an unphysical spin state in the
vicinity of the hole, we find no evidence for this in the entire investigated
region. More precisely, the physical limits of the magnetization in Eq. (13)
is $\pm 1/(1-2M_{\rm fl})\simeq\pm 1.67$ corresponding to having exactly $0$
or $1$ spin excitations at a given site. The most extreme value is associated
with the nearest neighbor magnetization $M(d=1)\simeq-0.81$ at $J/t=0.01$.
This corresponds to a mean value of spin excitations of $M^{(2)}(d=1)=0.74$.

###  IV.2 Symmetries

The preceding discussion suggests an underlying symmetry of the magnetization
around the hole, which we will explore for general crystal momentum of the
polaron in this section. The antiferromagnetic spin lattice exhibits several
point symmetries, namely mirror symmetries with respect to the two principal
axes and the diagonals of the lattice, as well as $C_{4}$ rotations. This
combines to the symmetry group $C_{4v}$.

Figure 4: Magnetization around the hole for $J/t=0.3$ vs indicated crystal
momenta [panel (a)] in the Heisenberg limit. For crystal momenta on the MBZ
boundary [panels (b)–(f)], i.e., $|p_{x}|+|p_{y}|=\pi$, the system retains the
reflection symmetries around the two diagonals (green-dashed lines), even
though the total crystal momentum of the magnetic polaron might point off
these axes. For crystal momenta marked with green (purple) symbols, the
magnetization cloud is oriented along the $x=y$ ($x=-y$) diagonal. When ${\bf
p}=(\pm\pi,0),(0,\pm\pi)$ [black squares, panel (d)], the dressing cloud
recovers the full $C_{4v}$ symmetry of the background antiferromagnetic order.
This symmetry is again reduced when the crystal momentum does not lie on the
MBZ boundary, as can be seen by comparing panels (d) [in the cases of ${\bf
p}=(\pm\pi,0)$] and (g) [${\bf p}=(\pm\pi/2,0)$], which have different
symmetries even though the momenta are parallel. For a general momentum in
panel (h), only the $C_{2}$ symmetry is retained and four crystal momentum
states [black stars] show the same magnetization cloud pattern.

The symmetry group for the magnetic polaron must therefore be a descendant of
$C_{4v}$. It turns out that the magnetic dressing cloud of the polaron has a
remarkably high symmetry and that the full $C_{4v}$ is recovered for certain
crystal momenta ${\bf p}$ of the polaron. First, it follows from time-reversal
symmetry that the dressing cloud is inversion symmetric, i.e. $M_{\bf p}({\bf
d})=M_{\bf p}(-{\bf d})$, for all crystal momenta, corresponding to the
$C_{2}$ point symmetry group in two dimensions. The two-point hole-spin
operator $\hat{M}^{(2)}({\bf
d})\propto\hat{h}^{\dagger}_{\mathbf{0}}\hat{h}_{\mathbf{0}}\hat{s}^{\dagger}_{{\bf
d}}\hat{s}_{{\bf d}}$ [see Eq. (14)] gives the spatial structure of the
magnetization around the hole, and we can therefore argue for the spatial
symmetries from here. Since it is hermitian, it fulfills the identity
$\bra{\Psi_{\bf p}}\hat{M}^{(2)}({\bf d})\ket{\Psi_{\bf
p}}=\bra{\tilde{\Psi}_{\bf p}}\hat{\mathcal{T}}\hat{M}^{(2)}({\bf
d})\hat{\mathcal{T}}^{-1}\ket{\tilde{\Psi}_{\bf p}}$ Sakurai (1994). Here,
$\hat{\mathcal{T}}$ is the anti-unitary time-reversal operator and
$\ket{\tilde{\Psi}_{\bf p}}=\hat{\mathcal{T}}\ket{\Psi_{\bf p}}$ is the time-
reversed polaron state. Now, the magnetization operator is invariant under
$\hat{\mathcal{T}}$, i.e. $\hat{\mathcal{T}}\hat{M}^{(2)}({\bf
d})\hat{\mathcal{T}}^{-1}=\hat{M}^{(2)}({\bf d})$, since reversal of time does
not affect position operators. On the other hand, reversal of time flips the
crystal momentum of the polaron so that $\hat{\mathcal{T}}\ket{\Psi_{\bf
p}}=\ket{\Psi_{-{\bf p}}}$. Consequently, $M_{-{\bf p}}({\bf d})=M_{\bf
p}({\bf d})$ and using the total inversion symmetry of the system $M_{-{\bf
p}}({\bf d})=M_{{\bf p}}(-{\bf d})$, we finally arrive at the $C_{2}$
inversion symmetry of the magnetization for any crystal momentum

$\displaystyle M_{\bf p}(-{\bf d})=M_{{\bf p}}({\bf d}).$ (18)

Note that this symmetry holds in any state given by a real linear combination
of the crystal momentum eigenstates. This general $C_{2}$ point symmetry group
of the magnetic dressing cloud may come as a surprise, since ${\bf p}$ defines
a specific direction which stands at odds with inversion symmetry along the
direction of the crystal momentum.

Figure 5: Magnetization around the hole for $J/t=0.3$ and different crystal
momenta in the Heisenberg limit (top) compared to the scattering probability
$\bar{P}^{\rm scat}_{{\bf p}\to{\bf q}}=P^{\rm scat}_{{\bf p}\to+{\bf
q}}+P^{\rm scat}_{{\bf p}\to-{\bf q}}$ defined in Eq. (19) and normalized to
its maximal value (bottom). While the crystal momentum of the magnetic polaron
may be misaligned with the magnetization cloud, we find that the scattering
probability, characterizing in which direction the hole preferably moves under
the emission of spin waves, generally reflects the orientation of the magnetic
dressing cloud. For the cases illustrated in panels (d) and (e), the polaron
shows a strong tendency to scatter to very few momentum states, making the
scattering profile remarkably sharp.

Higher spatial symmetries emerge for special crystal momenta. When ${\bf p}$
is parallel to one of the lattice axes or one of the diagonals, the
magnetization is symmetric under reflection operations parallel or
perpendicular to the momentum. Combined with the general $C_{2}$ symmetry this
forms the $C_{2v}$ symmetry group. Remarkably, the reflection symmetries along
the lattice diagonals are retained for _all_ crystal momenta along the edge of
the magnetic Brillouin zone (MBZ) given by $|p_{x}|+|p_{y}|=\pi$, as indicated
by black lines in Fig. 4(a). This can be understood from the symmetry of the
magnetization under the translation ${\bf p}\to{\bf p}\pm{\bf Q}$, where ${\bf
Q}=(\pi,\pi)$ is the wave vector of the antiferromagnetic spin-density wave.
For crystal momenta along the line $(p_{x},p_{y})=(p_{x},\pi-p_{x})$, it
follows that states with momenta $(p_{x},\pi-p_{x})-{\bf
Q}=(p_{x}-\pi,-p_{x})=(-p_{y},-p_{x})$ must show the same magnetization
pattern. This gives the reflection symmetry around the $y=-x$ diagonal. Using
the $C_{2}$ symmetry described above then gives that $(p_{y},p_{x})$ and
$(p_{x},p_{y})$ show the same magnetization pattern, which corresponds to a
reflection symmetry around the $y=x$ diagonal. As a result, the dressing cloud
of the magnetic polaron with momenta along the edge of the MBZ is
characterized by the high point symmetry group $C_{2v}$. It even follows that
when the momentum is at the corners of the MBZ, i.e. ${\bf p}=(\pm\pi,0)$ or
${\bf p}=(0,\pm\pi)$, the dressing cloud of the polaron has the full symmetry
$C_{4v}$ of the AF state without the polaron. We note that this unusual
spatial symmetry of the moving polaron is a fundamental property of the
system, and holds generally irrespective of our SCBA treatment. Remarkably,
the highest spatial symmetry of the magnetic polaron does, hereby, not emerge
in its ground state at ${\bf p}=(\pm\pi/2,\pm\pi/2)$.

The described symmetries are summarized in Fig. 4 showing the magnetic
dressing cloud for different crystal momenta. Moving along the MBZ [Figs.
4(b)–4(f)], the magnetization cloud undergoes a discrete rotation where it is
always oriented along one of the diagonals, which leads to a misalignment of
the dressing cloud and the crystal momentum. The dressing cloud recovers the
full $C_{4v}$ symmetry of the background antiferromagnetic order when ${\bf
p}=(\pm\pi,0),(0,\pm\pi)$ [Fig. 4(d)], whereas the symmetry descends to
$C_{2v}$ when the momentum reduces to ${\bf p}=(\pm\pi/2,0)$ [Fig. 4(g)].
Finally, Fig. 4(h) shows how the magnetic dressing cloud exhibits the minimal
$C_{2}$ symmetry for a general momentum, while still being misaligned with the
crystal momentum of the hole.

We can develop a microscopic picture of the discussed symmetries by
considering the probability to find a hole at momentum ${\bf q}$ in a polaron
state with momentum ${\bf p}$. This probability is given by

$\displaystyle P^{\rm scat}_{{\bf p}\to{\bf q}}=\bra{\Psi_{\bf
p}}\hat{h}^{\dagger}_{{\bf q}}\hat{h}_{{\bf q}}\ket{\Psi_{\bf p}}-Z_{\bf
p}\cdot\delta_{{\bf p},{\bf q}},$ (19)

and describes the scattering between the hole and excitations of the AF in the
magnetic polaron state. Since the bare hole state $\hat{h}_{\bf
p}^{\dagger}\ket{{\rm AF}}$ carries no spin wave excitations, and therefore
does not contain any information about the magnetization cloud, we subtract
the probability, $Z_{\bf p}$, of remaining in that state. By itself, this
scattering probability is not symmetric under inversion of the scattered
momentum, ${\bf q}\to-{\bf q}$, but we can consider its symmetrized form
$\bar{P}^{\rm scat}_{{\bf p}\to{\bf q}}=P^{\rm scat}_{{\bf p}\to+{\bf
q}}+P^{\rm scat}_{{\bf p}\to-{\bf q}}$ which is shown in Fig. 5 for selected
crystal momenta. The depicted momentum distributions indicate the preferred
directionality of the motion of the hole for a given momentum ${\bf p}$ of the
magnetic polaron, and its maxima indeed reflect the orientation of the
magnetization cloud. In particular, for the ground state momentum, ${\bf
p}=(\pi/2,\pi/2)$, we see in Fig. 5(a) that the hole motion predominantly
remains along the $x=y$ diagonal, thereby reducing the magnetization the most
in this direction. The comparison between $M_{{\bf p}}({\bf d})$ and
$\bar{P}^{\rm scat}_{{\bf p}\to{\bf q}}$ for ${\bf p}$ clearly shows that it
is not the polaron crystal momentum but the preferred direction of the hole
momentum that determines the symmetry and orientation of of the magnetic
dressing cloud. Like the magnetization, $\bar{P}^{\rm scat}_{{\bf p}\to{\bf
q}}$ is obtained from a self-consistency equation, which we derive in Appendix
C.

###  IV.3 Anisotropic spin interactions and the Ising limit

The developed framework also permits to study the transition from the
Heisenberg to the Ising limit, by tuning the parameter $\alpha$ from $1$ to
$0$ in Eq. (3). While choosing $\alpha\neq 1$ looses the correspondence with
the original Fermi-Hubbard Hamiltonian, approaching the Ising limit makes the
$t$-$J$ model better accessible to approximate treatments, such as a
description in terms of defect-strings described in Bulaevskii  _et al._
(1968); Brinkman and Rice (1970); Trugman (1988); Manousakis (2007); Grusdt
_et al._ (2018a, b, 2019); Bohrdt  _et al._ (2019).

In Fig. 6(a), we plot the local magnetization around the hole as a function of
$\alpha$ for $J/t=0.3$. As one approaches the Ising limit, the magnetization
cloud deforms and becomes more symmetric. This is associated with the
appearance of a gap in the spin wave dispersion $\omega_{\bf
k}=zJ\sqrt{1-\alpha^{2}\gamma_{\bf k}^{2}}/2$, making it harder for the hole
to emit spin waves when $\alpha<1$. The direction perpendicular to ${\bf p}$
stands out since the magnetization decreases as the Ising limit is approached,
and the coherent increase of the magnetization in a Heisenberg magnet is lost.
The Ising limit, $\alpha=0$, restores the $C_{4v}$ symmetry of the
magnetization, whereby the magnetization of all next nearest neighbors becomes
identical.

Figure 6: Magnetization in the neighborhood of a hole for the polaron ground
state with ${\bf p}=(\pi/2,\pi/2)$ as a function of $\alpha$ determining the
anisotropy of the spin-spin interactions for $J/t=0.3$ in a 16 by 16 square.
(a) Near the Heisenberg limit, $\alpha\lesssim 1$, the magnetization of the
nearest neighbors (red squares) and of the next-nearest neighbors in direction
of the hole crystal momentum (blue dots) rapidly changes with $\alpha$,
whereas it only changes slightly for $0\leq\alpha\lesssim 0.8$. (b)
Magnetization heat maps for the Ising case ($\alpha=0$, left) and for
$\alpha=0.5$ (right). Note the difference in scale from Fig. 5(a), which
illustrates the shrinking dressing cloud as the Ising limit is approached.

This isotropy, shown in Fig. 6(b), makes it possible to derive an analytical
expression for the magnetization in the Ising limit, as we show explicitly in
Appendix D. Due to corrections to the hole Green’s function beyond the SCBA
Chernyshev and Leung (1999) and the so-called Trugman loops Trugman (1988),
this result is not exact. In fact, these small corrections have been shown
Trugman (1988); Chernyshev and Leung (1999) to slightly lift the massive
degeneracy of the crystal momentum states, favoring the crystal momentum ${\bf
p}=(0,0)$ as the ground state. Therefore, the reemergence of the $C_{4v}$
symmetry in the Ising limit Grusdt  _et al._ (2018b) can be expected to be an
exact result for the ground state of the $t$-$J_{z}$ model.

##  V Strong-coupling effects

Previous SCBA calculations of spin-hole correlations have required a
truncation of the SCBA polaron wave function restricting the number number of
spin wave excitations Ramšak and Horsch (1998); Bała and Oleś (2000) or
considered essentially flat excitation spectra Ramšak and Horsch (1998); Bała
_et al._ (2002) to simplify higher-order terms in the polaron wave function.

Since the hole-spin interaction vertex scales as $g({\bf q},{\bf p})\sim t$,
and the spin wave energy scales as $\sim J$, such a truncation of the wave
function can be understood as a perturbative series in $t/J$. This naturally
limits the accuracy of such calculations to the weak to intermediate coupling
regime $J/t\gtrsim 1$. This is illustrated in Fig. 7, where we compare the
results from this truncation approach to our non-perturbative theory. For the
interaction strength $J/t=1$, the results converge nicely, and the
magnetization calculated by including up to three spin waves in the wave
function is essentially identical to the non-perturbative result. However, we
see that the effects of the hole on the surrounding magnetization is
significantly underestimated by the truncated wave functions for $J/t=0.05$.
In particular, it completely misses the sign flip in the magnetization at the
nearest neighbor sites.

One should note that the $t$-$J$ model no longer describes the Fermi-Hubbard
model when $J\sim t$. The found discrepancies, therefore, render any such
truncation procedures practically inapplicable when comparing to the Fermi-
Hubbard model. In contrast, the non-perturbative framework developed here
allows for the inclusion of all terms in Eq. (10), whereby one can describe
the spatial correlation of holes and spins deep in the strongly correlated
regime of $J/t\ll 1$, in which the SCBA is expected to yield an efficient and
accurate description of the $t$-$J$ Hamiltonian as well as the Fermi-Hubbard
model around half filling.

Figure 7: Comparison of the full magnetization obtained from our non-
perturbative calculation (solid lines) with that obtained by truncating the
summation of the lowest terms (dashed lines) for $J=t$ (a) and $J=0.05t$ (b)
as a function of distance $d$ from the hole along the $y$-direction in a 16 by
16 square lattice. The dashed lines correspond to including 1, 2, and 3 spin
waves in the magnetic polaron wave function, Eq. (10). For $J=t$, keeping only
2 spin waves gives an accurate description, while the 3 spin wave result is
indistinguishable from the full calculation. For $J=0.05t$, however, the full
calculation yields quantitatively and qualitatively different results,
including a sign flip of the magnetization at the nearest neighbor, $d=1$.

##  VI Derivation of self-consistency equations

We now describe the theoretical framework to include all terms in the SCBA
wave function, Eq. (10), as used in the previous section. In particular, we
derive self-consistency equations for the $B$ and $C$ functions in Eqs. (16)
and (17), which makes it possible to sum the infinite series in terms of spin
wave excitation numbers.

Figure 8 shows the first few diagrams in the $B$ series. We construct these
diagrams in the following way. First, we take the polaron wave function
$\ket{\Psi_{\bf p}}$ and its adjoint $\bra{\Psi_{\bf p}}$, corresponding to
Fig. 2(b) and its mirror image respectively. Second, we place the operator
$\hat{B}=N\sum_{\bf k}\hat{h}^{\dagger}_{{\bf k}+{\bf q}_{1}}\hat{h}_{{\bf
k}+{\bf q}_{2}}\hat{b}^{\dagger}_{-{\bf q}_{1}}\hat{b}_{-{\bf q}_{2}}$ between
$\bra{\Psi_{\bf p}}$ and $\ket{\Psi_{\bf p}}$ to compute the expectation value
$B({\bf q}_{1},{\bf q}_{2};{\bf p},\varepsilon_{\bf p})$ in Eq. (16). A non-
zero contribution to this expectation value involves the annihilation of a
hole and a spin wave for both $\bra{\Psi_{\bf p}}$ and $\ket{\Psi_{\bf p}}$,
which yields the structure $B_{1}$ in Fig. 8(b). Finally, all spin waves from
$\bra{\Psi_{\bf p}}$ and $\ket{\Psi_{\bf p}}$ that are not annihilated by
$\hat{B}$ must be joined together. Visually, this means that the series can be
constructed from the norm series shown in Fig. 8(a). Specifically, we detach
all double lines from black and red dots, insert the likewise detached $B_{1}$
diagram and reassemble the double lines with the dots. As a result, an
infinite series of terms $B_{1},B_{2},\dots$ emerges Ramšak and Horsch (1993,
1998), as shown in Fig. 8(b). The $n$-th term,
$B_{n}=\sum_{i=1}^{n}B_{n}^{(i)}$, contains $n$ nonzero diagrams coming from
the $n$ spin wave term in the polaron wave function (10). In all diagrams, we
have suppressed the two overall single lines corresponding to the residue,
$Z_{\bf p}$, as depicted in Fig. 2(e). All diagrams, where the spin wave lines
of the wave function cross are not allowed within the SCBA, and have to be
omitted for consistency. This is at the heart of the SCBA, in which only
rainbow diagrams are included. Also, all diagrams that are left-right
asymmetric vanish, as we show explicitly in Appendix B.

(a)$\\!\\!\\!\\!\\!\mathcal{N}/Z_{\bf
p}=$$1\;\,+\;\;$$+$$\phantom{\;\dots}+\;\dots$(b)$\\!\\!\\!\\!\\!B_{1}=$$\\!\\!\\!\\!\\!B_{2}=$$+$$\\!\\!\\!\\!\\!B_{3}=$$+$$+$$\\!\\!\\!\\!\\!B_{4}=$$\dots$
Figure 8: Diagrammatic representation for the norm series (a) and the $B$
series (b). The $B$ series comes from inserting the basic diagram, $B_{1}$, in
the norm series (a). The shown diagrams are the ones that give a nonzero
contribution to the $B$ function. These come in an increasing order denoted
$B_{n}$, where $n$ denotes the number of interaction vertices (red dots) to
each side. At a given order $n$, there are $n$ nonzero diagrams.

The central idea to obtain the self-consistency equations is the following.
First, we take the last diagram from each order, $B_{n}^{(n)}$, and sum up
only these, to obtain

$B_{0}=\sum_{n=1}^{\infty}B_{n}^{(n)}.$ (20)

This leads to the diagrammatic structure shown in Figs. 9(a) and 9(b). Second,
we notice that all other terms in Fig. 8 are related to $B_{0}$ by putting a
number of spin wave lines around $B_{0}$. The result is the full $B$ function
shown in Fig. 9(c). Algebraically, the $B_{0}$ function is thus

$\displaystyle B_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=$
$\displaystyle\;Ng({\bf p},{\bf q}_{1})g({\bf p},{\bf q}_{2})$
$\displaystyle\cdot G({\bf p}\\!+\\!{\bf q}_{1},\omega\\!-\\!\omega_{{\bf
q}_{1}})G({\bf p}\\!+\\!{\bf q}_{2},\omega\\!-\\!\omega_{{\bf q}_{2}})$
$\displaystyle\cdot\left[1+F_{B}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)\right].$ (21)

This is written in terms of the function $F_{B}$, which fulfills the self-
consistency equation

$\displaystyle\\!\\!\\!\\!\\!F_{B}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=$
$\displaystyle\sum_{\bf k}g({\bf p}\\!+\\!{\bf q}_{1},{\bf k})g({\bf
p}\\!+\\!{\bf q}_{2},{\bf k})$ $\displaystyle\cdot G({\bf p}\\!+\\!{\bf
q}_{1}\\!+\\!{\bf k},\omega\\!-\\!\omega_{{\bf q}_{1}}\\!-\\!\omega_{\bf k})$
$\displaystyle\cdot G({\bf p}\\!+\\!{\bf q}_{2}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{{\bf q}_{2}}\\!-\\!\omega_{\bf k})$
$\displaystyle\cdot\left[1+F_{B}({\bf q}_{1},{\bf q}_{2};{\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})\right],$ (22)

that is depicted diagrammatically in Fig. 9(b). Both of these depend on the
general energy $\omega\leq\varepsilon_{\bf p}$. Finally, the self-consistency
equation shown in Fig. 9(c) can be written as

$\displaystyle B({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=$
$\displaystyle\;B_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)$
$\displaystyle+\sum_{\bf k}g^{2}({\bf p},{\bf k})G^{2}({\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})$ $\displaystyle\cdot B({\bf q}_{1},{\bf
q}_{2};{\bf p}\\!+\\!{\bf k},\omega\\!-\\!\omega_{\bf k}).$ (23)

Importantly, we can relate the structure of Eqs. (22) and (23) to that of the
self-energy, in Eq. (8). Specifically, taking the derivative of the self-
energy equation yields

$\displaystyle-\partial_{\omega}\Sigma({\bf p},\omega)=$
$\displaystyle\,\sum_{\bf k}g^{2}({\bf p},{\bf k})G^{2}({\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})$
$\displaystyle\,\cdot\left[1-\partial_{\omega}\Sigma({\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})\right],$ (24)

which corresponds to the norm series in Fig. 8(a), apart from the first term.
Comparing this to the self-consistency equation for the $B$ function, we see
that $-\partial_{\omega}\Sigma({\bf p},\omega)$ takes on the role of $B$,
while $\sum_{\bf k}g^{2}({\bf p},{\bf k})G^{2}({\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})$ corresponds to $B_{0}$. For each external
momenta ${\bf q}_{1}$ and ${\bf q}_{2}$, these more advanced self-consistency
equations, thus, show the exact same structure as that of the self-energy. In
fact, computing $F_{B}({\mathbf{0}},{\mathbf{0}};{\bf
p},\omega+\omega_{\mathbf{0}})$ in Eq. (22) and comparing it to Eq. (24) shows
that

$\displaystyle F_{B}({\mathbf{0}},{\mathbf{0}};{\bf
p},\omega+\omega_{\mathbf{0}})=-\partial_{\omega}\Sigma({\bf p},\omega),$ (25)

which in turn is related to the quasiparticle residue in Eq. (9). This link
provides a useful consistency check for the numerical calculations, and also
offers an alternative way of computing the residue. Furthermore, it shows that
the order of the $B$ and $C$ functions (including the overall factor of
$Z_{\bf p}$) is $Z_{\bf p}\cdot\mathcal{O}[-\partial_{\omega}\Sigma({\bf
p},\omega)]=Z_{\bf p}(1/Z_{\bf p}-1)=1-Z_{\bf p}$. This emphasizes that small
quasiparticle residues, $Z_{\bf p}\ll 1$, correspond to large changes in the
local magnetization cloud, $B,C\sim 1$, as was discussed in more general terms
in Sec. III.

(a)$\;B_{0}$$=$$+$(b)$=$$F_{B}$$+$$F_{B}$(c)$\;B\phantom{{}_{0}}$$=$$B_{0}\;+\;$$B_{0}$$\phantom{\dots}+\dots$$=$$B_{0}\;+\;$$B$
Figure 9: (a) Summation of all the last diagrams, $B_{0}=\sum_{n}B_{n}^{(n)}$.
(b) The $F_{B}$ function thus emerges. (c) Finally, the full $B$ function
comes about by putting $0,1,2\dots$ spin wave lines around $B_{0}$. This
results in a self-consistent equation for the $B$ function.

For the $C$ function, an equivalent derivation, see Appendix A, leads to the
feeding term

$\displaystyle C_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=N\big{[}$
$\displaystyle g({\bf p},{\bf q}_{1})g({\bf p}\\!+\\!{\bf q}_{1},{\bf
q}_{2})G({\bf p}\\!+\\!{\bf q}_{1},\omega\\!-\\!\omega_{{\bf q}_{1}})$
$\displaystyle+$ $\displaystyle g({\bf p},{\bf q}_{2})g({\bf p}\\!+\\!{\bf
q}_{2},{\bf q}_{1})G({\bf p}\\!+\\!{\bf q}_{2},\omega\\!-\\!\omega_{{\bf
q}_{2}})\big{]}$ $\displaystyle\cdot G({\bf p}\\!+\\!{\bf q}_{1}\\!+\\!{\bf
q}_{2},\omega\\!-\\!\omega_{{\bf q}_{1}}\\!-\\!\omega_{{\bf q}_{2}}),$
$\displaystyle\cdot\left[1+F_{C}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)\right].$ (26)

This is written in terms of the $F_{C}$ function

$\displaystyle\\!\\!\\!\\!\\!F_{C}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=$
$\displaystyle\sum_{\bf k}g({\bf p},{\bf k})g({\bf p}\\!+\\!{\bf
q}_{1}\\!+\\!{\bf q}_{2},{\bf k})$ $\displaystyle\cdot G({\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})$ $\displaystyle\cdot G({\bf p}\\!+\\!{\bf
q}_{1}\\!+\\!{\bf q}_{2}\\!+\\!{\bf k},\omega\\!-\\!\omega_{{\bf
q}_{1}}\\!-\\!\omega_{{\bf q}_{2}}\\!-\\!\omega_{\bf k})$
$\displaystyle\cdot\left[1+F_{C}({\bf q}_{1},{\bf q}_{2};{\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})\right],$ (27)

which has the same role as $F_{B}$ in the $B$ series, and which it also
closely resembles. The final self-consistency equation for the $C$ function is
identical to Eq. (23) with $C$ and $C_{0}$ playing the role of $B$ and $B_{0}$
respectively.

The full computation of the magnetization in the neighborhood of a hole can
then be performed as follows. For each pair of external spin wave momenta
${\bf q}_{1},{\bf q}_{2}$, we first solve the two self-consistency equations
for $F_{B}$ and $F_{C}$, Eqs. (22) and (27) respectively. From these, we
calculate $B_{0}$ and $C_{0}$, Eqs. (21) and (26). These are used in the self-
consistency equation (23). We then multiply $B$ and $C$ with the overall
residue $Z_{\bf p}$, and finally transform them to position space according to
Eqs. (16) and (17). The task of solving $N^{2}$ self-consistency equations may
seem daunting at first sight. A major simplification, however, comes from a
number of symmetry properties of the $B$ and $C$ functions in momentum space,
as outlined in Appendix E. Therefore, only a small fraction of the $N^{2}$
equations has to be solved explicitly, whereby all results presented in this
work could be obtained with modest computational resources.

##  VII Experiments

The spatial structure of magnetic polarons can play an important role for the
transport properties of electrons in a solid. While experiments show evidence
for their formation in the cuprates Wells  _et al._ (1995); Ronning  _et al._
(1998); Kim  _et al._ (1998); Yoshida  _et al._ (2003); Ronning  _et al._
(2005); Rösch  _et al._ (2005); Schrieffer and Brooks (2007), a detailed and
direct probing of the underlying correlations on a microscopic level has not
been possible in condensed matter measurements. This has changed with the
development of quantum simulation platforms based on cold atoms in optical
lattices Esslinger (2010); Gross and Bloch (2017); Schäfer  _et al._ (2020),
and has in recent years made it possible to implement near-perfect
realizations of the Fermi-Hubbard model Boll  _et al._ (2016); Cheuk  _et al._
(2016); Mazurenko  _et al._ (2017); Hilker  _et al._ (2017); Brown  _et al._
(2017); Chiu  _et al._ (2018); Brown  _et al._ (2019); Koepsell  _et al._
(2019); Chiu  _et al._ (2019); Koepsell  _et al._ (2020); Brown  _et al._
(2020); Vijayan  _et al._ (2020); Hartke  _et al._ (2020); Guardado-Sanchez
_et al._ (2020); Ji  _et al._ (2021).

The single-site resolution achievable in current experiments permits to image
any desired correlation function between particles, and in particular between
a dopant and its surrounding effective spins. This opens up the possibility to
probe the inner structure of the magnetic polaron and its motion, as explored
in the present work. So far, the lowest achieved temperatures are around
$k_{\rm B}T=0.5J$ Chiu  _et al._ (2019), at which there may be significant
thermal corrections to the pure polaron states investigated here. Reaching
lower temperatures will be an important step for direct comparisons and is
widely expected to yield key insights into the microscopic physics of magnetic
polarons and their role in high-$T_{c}$ superconductivity in strongly
correlated materials.

While such quantum simulators of the Fermi-Hubbard Hamiltonian naturally
realize the isotropic $t$-$J$ model, the versatile toolbox to control and
manipulate cold atoms also makes it possible to implement more general spin
Hamiltonians, such as considered in Eq. (3). This includes polar molecules
Gorshkov  _et al._ (2011), as well as Rydberg-dressed atoms in optical
lattices Glaetzle  _et al._ (2015); van Bijnen and Pohl (2015); Zeiher  _et
al._ (2016, 2017); Borish  _et al._ (2020); Guardado-Sanchez  _et al._ (2021),
which will make it possible to continuously tune between the $t$-$J$ and
$t$-$J_{z}$ model, and to realize low temperatures compared to the much larger
spin interactions achievable in these systems.

##  VIII Conclusions and outlook

Inspired by recent experimental breakthroughs, we explored the properties of
magnetic polarons that are formed by a hole and spin fluctuations in an
antiferromagnetic square lattice, as described by the $t$-$J$ model. By
combining the SCBA for the hole Green’s function with the many-body wave
function wave of the polaron, we developed a non-perturbative resummation
scheme that now makes it possible to determine spin-hole correlations in the
strongly interacting regime. This method thus enables broad explorations of
the microscopic structure of magnetic polarons, which, so far, has not been
possible within the SCBA. Given the proven accuracy of the SCBA for one-body
observables such as the energy and residue of the polaron, this constitutes a
significant step forward and will enable detailed analyses of ongoing
experiments based on cold-atom quantum simulators. To illustrate the power of
the approach, we have considered the magnetization in the vicinity of the
hole, which turns out to deviate considerable from previous perturbative
results under conditions where the $t$-$J$ model is valid. For a moving hole
the magnetization cloud has an elongated shape, which features a surprising
misalignment with the hole momentum that originates from the various
symmetries of the antiferromagnetic state of the underlying spin lattice.

The demonstrated possibility to explore correlations within the SCBA opens up
new perspectives for studying strongly correlated quantum matter, including
the physics of cuprates in the limit of small doping. It has been shown Kyung
and Ferrell (1996); Xiang and Wheatley (1996); Manousakis (2007) that these
systems can be modeled quantitatively by including next-nearest neighbor
hopping terms in the Hamiltonian, defining the so-called
$t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model. Within linear spin wave theory,
this amounts to the addition of kinetic energy terms for the hole and can,
thus, straightforwardly be including in our methodology.

The developed approach in general offers a promising starting point for
extensions of the method along several directions and will facilitate detailed
comparisons to recent cold-atom experiments on the Fermi-Hubbard model. While
we have focused here on two-dimensional square lattices, the developed
theoretical framework can be equally applied to any bipartite Bravais lattice
in one, two or three spatial dimensions. It will moreover be interesting to
assess corrections beyond linear spin wave theory Singh and Gelfand (1995);
Kim  _et al._ (1999); Sandvik and Singh (2001); Rønnow  _et al._ (2001);
Coldea  _et al._ (2001); Christensen  _et al._ (2004) and to explore higher-
order correlation functions Grusdt  _et al._ (2018b, 2019); Blomquist and
Carlström (2020), which have been observed in recent cold-atom experiments
Koepsell  _et al._ (2019); Chiu  _et al._ (2019); Koepsell  _et al._ (2020).
The SCBA approach as used here to describe spin-hole correlations may also be
employed to analyze correlations between two holes Dagotto (1994); Chernyshev
_et al._ (1994); Belinicher  _et al._ (1997), which will contribute to the
understanding of pairing and a potential mechanism for high temperature
superconductivity in the limit of small doping Schrieffer  _et al._ (1988);
Shraiman and Siggia (1989); Frenkel and Hanke (1990); Eder (1992); Izyumov
(1997); Riera and Dagotto (1998). While the present pure-state treatment
restricts our predictions to zero temperature, generalizing our framework to
finite temperatures will make it possible to characterize the impact of
temperatures that are currently achievable in optical-lattice experiments
Koepsell  _et al._ (2019); Chiu  _et al._ (2019); Koepsell  _et al._ (2020)
and enable tests of the SCBA framework based on direct comparisons to measured
correlation functions at finite temperature and strong interactions.
Circumventing current temperature limitations, recent experiments Ji  _et al._
(2021) have probed the transient dynamics following hole creation Bohrdt  _et
al._ (2020); Hubig  _et al._ (2020), which permits to trace the formation of
magnetic polarons. We anticipate that the SCBA approach developed in this work
can also provide an accurate framework to describe the non-equilibrium
dynamics of polarons in strongly interacting quantum magnets.

###### Acknowledgements.

We thank Annabelle Bohrdt and Fabian Grusdt for valuable feedback on our
manuscript. KKN would like to thank Simon Panyella Pedersen for setting up the
numerical calculation on the CSCAA Grendel cluster. This work has been
supported by the Danish National Research Foundation through the Center of
Excellence “CCQ” (Grant agreement no.: DNRF156).

##  Appendix A Derivation of $C$ function

The derivation of the $C$ function closely follows the derivation of the $B$
function in Sec. VI. The basic diagram for the $C$ series is denoted $C_{02}$
and shown in Fig. 10(a). Analogous to the $B$ series, the full $C$ series is
then obtained by placing this basic diagram in the norm series, Fig. 8(a).
This leads to the terms $C_{n,n+2}$, where $n$ is the number of spin waves
from the adjoint state $\bra{\Psi_{\bf p}}$ and $n+2$ is the number of spin
waves from $\ket{\Psi_{\bf p}}$. There are $n+1$ terms at order $n$:
$C_{n,n+2}=\sum_{i=0}^{n}C_{n,n+2}^{(i)}$. As for the $B$ series, we, then,
first sum up only the last diagrams $C_{n,n+2}^{(n)}$ at each order, defining

$C_{0}=\sum_{n}C_{n,n+2}^{(n)}.$ (28)

This leads to the diagrammatic structure shown in Figs. 10(b) and 10(c). The
result can be written in a similar form to $B_{0}$ and $F_{B}$, see Eqs. (21)
and (22). Explicitly,

$\displaystyle C_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=N\big{[}$
$\displaystyle g({\bf p},{\bf q}_{1})g({\bf p}\\!+\\!{\bf q}_{1},{\bf
q}_{2})G({\bf p}\\!+\\!{\bf q}_{1},\omega\\!-\\!\omega_{{\bf q}_{1}})$
$\displaystyle+$ $\displaystyle g({\bf p},{\bf q}_{2})g({\bf p}\\!+\\!{\bf
q}_{2},{\bf q}_{1})G({\bf p}\\!+\\!{\bf q}_{2},\omega\\!-\\!\omega_{{\bf
q}_{2}})\big{]}$ $\displaystyle\cdot G({\bf p}\\!+\\!{\bf q}_{1}\\!+\\!{\bf
q}_{2},\omega\\!-\\!\omega_{{\bf q}_{1}}\\!-\\!\omega_{{\bf q}_{2}}),$
$\displaystyle\cdot\left[1+F_{C}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)\right].$ (29)

Here, $F_{C}$ is given by Eq. (27), while the above expression is identical to
Eq. (26). By putting $0,1,2,\dots$ spin wave lines around $C_{0}$, as shown in
Fig. 10(c), the self-consistent equation for $C$ is achieved

$\displaystyle C({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=$
$\displaystyle\;C_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)$
$\displaystyle+\sum_{\bf k}g^{2}({\bf p},{\bf k})G^{2}({\bf p}\\!+\\!{\bf
k},\omega\\!-\\!\omega_{\bf k})$ $\displaystyle\cdot C({\bf q}_{1},{\bf
q}_{2};{\bf p}\\!+\\!{\bf k},\omega\\!-\\!\omega_{\bf k}),$ (30)

which has the exact same structure as Eq. (23) for the $B$ function.

(a)$\\!\\!\\!\\!\\!C_{02}=$$\\!\\!\\!\\!\\!C_{13}=$$+$$\\!\\!\\!\\!\\!C_{24}=$$+$$+$$\\!\\!\\!\\!\\!C_{35}=$$\dots$(b)$\\!\\!\\!\\!\\!C_{0}$$=$$+$$F_{C}$(c)$\\!\\!\\!\\!\\!C\phantom{{}_{0}}$$=$$C_{0}\;+\;$$C_{0}$$\phantom{\dots}+\dots$$=$$C_{0}\;+\;$$C$
Figure 10: Diagrammatic representation for the $C$ series. (a) The $C$ series
comes from inserting the basic diagram, $C_{02}$, in the norm series, Fig.
8(a). The shown diagrams are the ones that give a nonzero contribution to the
$C$ function. These come in an increasing order denoted $C_{n,n+2}$. At a
given order $n$, there are $n$ interaction vertices to the left, $n+2$ to the
right and a total of $n+1$ nonzero diagrams. (b) Summation of all the last
diagrams, $C_{0}=\sum_{n}C_{n,n+2}^{(n)}$. The appearing $F_{C}$ function is
shown in double spin wave lines. Diagrammatically, this looks identical to
$F_{B}$. However, because the exterior is different, the structure of the
self-consistent equation for $F_{C}$ [Eq. (27)] is different in this case. (c)
The self-consistent equation for the full $C$ function is achieved by putting
$0,1,2,\dots$ spin wave lines around $C_{0}$.

##  Appendix B Vanishing diagrams

In this section, we show that all left-right asymmetric diagrams in the $B$
series vanish, and that the corresponding diagrams in the $C$ series vanish as
well.

To understand how these diagrams vanish, we must first analyze a certain
symmetry of the interaction vertex $g({\bf p},{\bf k})$. Consider, therefore,
a Bravais lattice in which the lattice points can be written

${\bf r}=n_{1}{\bf a}_{1}+n_{2}{\bf a}_{2}+n_{3}{\bf a}_{3}.$ (31)

Here, the ${\bf a}_{i}$’s are the primitive vectors. In 2D, we simply set
$n_{3}=0$. In 1D, $n_{3}=n_{2}=0$. When the lattice is bi-partite, we can
choose the primitive vectors to be nearest neighbors to a given site, such
that ${\boldsymbol{\delta}}\in\\{\pm{\bf a}_{1},\pm{\bf a}_{2},\pm{\bf
a}_{3}\\}$. Because of the periodic boundary conditions ${\rm e}^{i{\bf
k}\cdot{\bf a}_{j}\cdot L_{j}}=1$ for a $L_{x}\times L_{y}\times L_{z}$
lattice. In turn,

${\bf k}\cdot{\bf a}_{j}=n_{j}\cdot\frac{2\pi}{L_{j}},$ (32)

where $n_{j}=-L_{j}/2+1,-L_{j}/2+2,\dots,L_{j}/2$ is an integer. Now, the wave
vector of the antiferromagnetic spin-density wave ${\bf Q}$ is defined by
letting $n_{j}=L_{j}/2$ for $j=x,y,z$. Then ${\rm e}^{i{\bf Q}\cdot{\bf
a}_{j}}={\rm e}^{i\pi}=-1$. In turn,

$\gamma_{{\bf k}+{\bf Q}}=\frac{1}{z}\sum_{{\boldsymbol{\delta}}}{\rm
e}^{i{\bf Q}\cdot{\boldsymbol{\delta}}}{\rm e}^{i{\bf
k}\cdot{\boldsymbol{\delta}}}=-\gamma_{\bf k},$ (33)

since ${\rm e}^{i{\bf Q}\cdot{\boldsymbol{\delta}}}=-1$, using that
${\boldsymbol{\delta}}\in\\{\pm{\bf a}_{1},\pm{\bf a}_{2},\pm{\bf a}_{3}\\}$.
As in Appendix B, the change in sign of $\gamma_{\bf k}$ also means that the
interaction changes sign

$g({\bf p},{\bf k}+{\bf Q})=g({\bf p}+{\bf Q},{\bf k})=-g({\bf p},{\bf k}).$
(34)

On the other hand, the self-energy, and thereby the Green’s function $G({\bf
p},\omega)$, is insensitive to this change in sign, because it scales with the
square of the interaction. Therefore,

$G({\bf p}+{\bf Q},\omega)=G({\bf p},\omega).$ (35)

As we shall now show this leads to the vanishing of all asymmetric $B$
diagrams, as well as the corresponding diagrams in the $C$ series.

Specifically, the asymmetric diagrams are all of the form shown in the top of
Fig. 11. The example shown evaluates to

$\displaystyle B_{2}^{\rm asym.}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)$
$\displaystyle=N\;g({\bf p},{\bf q}_{2})G({\bf p}\\!+\\!{\bf
q}_{2},\omega\\!-\\!\omega_{{\bf q}_{2}})$ $\displaystyle\cdot\sum_{{\bf
k}}g({\bf p},{\bf k})G({\bf p}\\!+\\!{\bf k},\omega\\!-\\!\omega_{\bf k})$
$\displaystyle\cdot g({\bf p}\\!+\\!{\bf k},{\bf q}_{1})g({\bf p}\\!+\\!{\bf
q}_{2},{\bf k})$ $\displaystyle\cdot G({\bf p}\\!+\\!{\bf k}\\!+\\!{\bf
q}_{1},\omega\\!-\\!\omega_{\bf k}\\!-\\!\omega_{{\bf q}_{1}})$
$\displaystyle\cdot G({\bf p}\\!+\\!{\bf k}\\!+\\!{\bf
q}_{2},\omega\\!-\\!\omega_{\bf k}\\!-\\!\omega_{{\bf q}_{2}}).\\!$ (36)

Here, the key point is that unlike the symmetric diagrams in Fig. 8, there is
an _odd_ number of terms with interaction vertices $g(\cdot,\cdot)$ that
depend on the summation index, ${\bf k}$. It is, therefore, sensitive to
changes in sign of $g$. As a result, the two terms ${\bf k}$ and ${\bf k}+{\bf
Q}$ in the sum in Eq. (36) have the same magnitude, but opposite signs.
Therefore, they cancel each other exactly. In this way, _all_ asymmetric
diagrams in the $B$ series vanish identically. Another way to understand this
vanishing is in terms of sublattice states. Every time the hole hops, it
changes sublattice. Therefore, the asymmetric diagrams like the one shown in
the top of Fig. 11 features overlaps of holes in opposite sublattices and thus
vanish. Indeed, the symmetry in Eq. (34) of the interaction is due to the
underlying sublattice symmetry of the system.

$B_{2}^{\rm asym.}=$$C_{13}^{\rm van.}=$ Figure 11: Examples of vanishing $B$
and $C$ diagrams from orders $B_{2}$ and $C_{13}$ respectively.

We finally comment on the vanishing of the remaining diagrams in the $C$
series. These diagrams have the structure shown in bottom part of Fig. 11,
where at least one spin wave line is between the external spin wave lines at
momenta $-{\bf q}_{1}$ and $-{\bf q}_{2}$ joining the black dot. The reason
that this vanishes is exactly the same as why the asymmetric diagrams in the
$B$ series all vanish. Specifically, a sum $\sim\sum_{\bf k}g({\bf p},{\bf
k})g({\bf p}+{\bf q}_{1},{\bf k})g({\bf p}+{\bf q}_{1}+{\bf k},{\bf q}_{2})$
appears, in which terms at ${\bf k}$ and ${\bf k}+{\bf Q}$ cancel.

##  Appendix C Self-consistency equation for the scattering probability

In this appendix, we derive a self-consistency equation for the scattering
probability

$P^{\rm scat}_{{\bf p}\to{\bf q}}=\bra{\Psi_{\bf p}}\hat{h}^{\dagger}_{{\bf
q}}\hat{h}_{{\bf q}}\ket{\Psi_{\bf p}}-Z_{\bf p}\cdot\delta_{{\bf p},{\bf q}}$
(37)

also defined in Eq. (19) of the main text. While this can be done
diagrammatically as for the $B$ and $C$ functions, it is just as simple to
write down the lowest order terms at a general energy $\omega$ and recognize
the pattern. As for the $B$ and $C$ series, we omit the overall factor of the
residue, $Z_{\bf p}$, in the following. This must be included in the end.
Since the final momentum must be ${\bf q}$, the total change in crystal
momentum is ${\bf q}-{\bf p}$. The lowest order term from the polaron wave
function is thus $g^{2}({\bf p},{\bf q}-{\bf p})G^{2}({\bf
q},\omega-\omega_{{\bf q}-{\bf p}})$. As required, the scattering probability
to a specific momentum state scales as $1/N$, since $g^{2}(\cdot,\cdot)\propto
1/N$. Inclusion of the next term yields

$\displaystyle P^{\rm scat}_{{\bf p}\to{\bf q}}(\omega)=\,$ $\displaystyle
g^{2}({\bf p},{\bf q}-{\bf p})G^{2}({\bf q},\omega-\omega_{{\bf q}-{\bf p}})$
$\displaystyle+\sum_{{\bf k}}g^{2}({\bf p},{\bf k})G^{2}({\bf p}+{\bf
k},\omega-\omega_{\bf k})$ $\displaystyle\cdot$ $\displaystyle g^{2}({\bf
p}+{\bf k},{\bf q}-{\bf p}-{\bf k})G^{2}({\bf q},\omega-\omega_{\bf
k}-\omega_{{\bf q}-{\bf p}-{\bf k}})$ $\displaystyle+\dots$

In the second term, there is a single free momentum, ${\bf k}$, but the final
momentum must again be ${\bf q}$. This gives the structure of the second term.
This series continues indefinitely, but we notice that it can be rewritten as
a self-consistency equation

$\displaystyle P^{\rm scat}_{{\bf p}\to{\bf q}}(\omega)=$
$\displaystyle\,g^{2}({\bf p},{\bf q}-{\bf p})G^{2}({\bf
q},\omega-\omega_{{\bf q}-{\bf p}})$ $\displaystyle\,+\sum_{{\bf k}}g^{2}({\bf
p},{\bf k})G^{2}({\bf p}+{\bf k},\omega-\omega_{\bf k})$ $\displaystyle\,\cdot
P^{\rm scat}_{{\bf p}+{\bf k}\to{\bf q}}(\omega-\omega_{\bf k}).$ (38)

Like the equations for the $B$\- and $C$-functions, this has the exact same
structure as the equation for the derivative of the self-energy [Eq. (24)]. By
solving this iteratively, evaluating at the quasiparticle peak,
$\varepsilon_{\bf p}$, and multiplying the result by the residue $Z_{\bf p}$,
we obtain the scattering probabilities plotted in Fig. 5. If we sum up all
contributions in Eq. (37), we obtain $\sum_{{\bf q}}P^{\rm scat}_{{\bf
p}\to{\bf q}}=1-Z_{\bf p}$, as one might expect. This summation rule gives a
good consistency check for the numerical calculations.

Figure 12: $C_{2}$ symmetrized scattering probability $\bar{P}^{\rm
scat}_{{\bf p}\to{\bf q}}=P^{\rm scat}_{{\bf p}\to+{\bf q}}+P^{\rm scat}_{{\bf
p}\to-{\bf q}}$ normalized to their maximal values for various indicated
crystal momenta and an interaction strength $J/t=0.3$ in a 16 by 16 square
lattice. Moving through panels (a) to (d), we increase $p_{y}$ from $4\pi/8$
to $7\pi/8$, keeping $p_{x}=7\pi/8$ constant. Panel (a) corresponds to Fig.
5(e).

In Figs. 5(d)–5(e), we noticed that the scattering profiles are remarkably
sharp. We look into this in Fig. 12, where we vary $p_{y}$ from $4\pi/8$ to
$7\pi/8$, keeping a constant $p_{x}=7\pi/8$ as in Fig. 5(e). This shows that
the sharp scattering profiles in Figs. 5(d)–5(e) are by no means exceptional.
In fact, at ${\bf p}=(7\pi/8,7\pi/8)$, the scattering probability is sharpened
further, and the hole dominantly scatters to only eight distinct momentum
states.

##  Appendix D Analytical result in the Ising limit

In this appendix, we derive an analytical formula for the magnetization in the
Ising limit, $\alpha=0$. This is possible due to a huge simplification in the
interaction, $g({\bf q},{\bf k})\to zt\cdot\gamma_{{\bf q}+{\bf k}}/\sqrt{N}$.
As a result of this simplification, the Green’s function within SCBA is
independent of crystal momentum, and fulfills the equation
$G^{-1}(\omega)=\omega-zJ/2-zt^{2}\cdot G(\omega\\!-\\!zJ/2)$. This allows for
an analytical solution in terms of a continued fraction Starykh and Reiter
(1996)

$G(\omega)=\frac{1}{\omega\\!-\\!zJ/2\\!-\\!zt^{2}\frac{1}{\omega-2\cdot
zJ/2-zt^{2}\frac{1}{\omega-3\cdot zJ/2\dots}}},$ (39)

facilitating a numerically simple implementation. We now show that a similar
description applies to the magnetization.

Since the antiferromagnetic coherence factors in the Ising limit are $u_{\bf
k}\to 1$ and $v_{\bf k}\to 0$, the $C$ series does not contribute to the
magnetization [Eq. (17)]. Focusing then on the $B$ series, we first calculate
$F_{B}$ from Eq. (22)

$\displaystyle F_{B}({\bf q}_{1}\\!-\\!{\bf q}_{2};\omega)=$
$\displaystyle\,zt^{2}\gamma_{{\bf q}_{1}-{\bf q}_{2}}G^{2}(\omega\\!-\\!zJ)$
$\displaystyle\cdot\left[1+F_{B}({\bf q}_{1}\\!-\\!{\bf
q}_{2},\omega\\!-\\!zJ/2)\right]$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}(zt^{2}\gamma_{{\bf q}_{1}-{\bf
q}_{2}})^{n}\prod_{k=1}^{n}G^{2}\\!\left(\omega\\!-\\!(k+1)\frac{zJ}{2}\right).$
(40)

To obtain the upper line, we assume that $F_{B}$ is independent of the polaron
crystal momentum ${\bf p}$, and use that $\sum_{{\bf k}}\gamma_{{\bf k}+{\bf
q}_{1}}\gamma_{{\bf k}+{\bf q}_{2}}/N=\gamma_{{\bf q}_{1}-{\bf q}_{2}}/z$.
Note that $F_{B}$ only depends on relative momentum, ${\bf q}_{1}-{\bf
q}_{2}$. Then, by repeatedly reinserting $F_{B}$ as described by the upper
line, we obtain the infinite series in the lower line. From Eq. (21), it then
follows that

$\displaystyle B_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=$
$\displaystyle\;(zt)^{2}\gamma_{{\bf p}+{\bf q}_{1}}\gamma_{{\bf p}+{\bf
q}_{2}}G^{2}(\omega\\!-\\!zJ/2)$ $\displaystyle\cdot\left[1+F_{B}({\bf
q}_{1}\\!-\\!{\bf q}_{2};\omega)\right].$ (41)

This, therefore, still depends on the momentum of the polaron, ${\bf p}$. To
get rid of this momentum dependency in the self-consistency equations, we
define

$\displaystyle\\!\\!\\!\\!\\!\Delta B({\bf q}_{1},{\bf q}_{2};\omega)=B({\bf
q}_{1},{\bf q}_{2};{\bf p},\omega)-B_{0}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega),\\!\\!$ (42)

and make the ansatz that this is independent of ${\bf p}$. Using Eq. (23), we,
in fact, get

$\displaystyle\Delta B({\bf q}_{1},{\bf q}_{2};\omega)=$
$\displaystyle\;\Delta B_{0}({\bf q}_{1},{\bf q}_{2};\omega)\\!+\\!zt^{2}\cdot
G^{2}(\omega\\!-\\!zJ/2)$ $\displaystyle\cdot\Delta B({\bf q}_{1},{\bf
q}_{2};\omega\\!-\\!zJ/2),$ (43)

where

$\displaystyle\Delta B_{0}({\bf q}_{1},{\bf q}_{2};\omega)=$
$\displaystyle(zt^{2})^{2}\Big{[}2\gamma_{{\bf q}_{1}}\gamma_{{\bf
q}_{2}}+\left(1-\frac{2}{z}\right)\gamma_{{\bf q}_{1}-{\bf q}_{2}}$
$\displaystyle-\frac{\gamma_{{\bf q}_{1}+{\bf
q}_{2}}}{z}\Big{]}G^{2}(\omega\\!-\\!zJ/2)G^{2}(\omega\\!-\\!zJ)$
$\displaystyle\cdot\left[1+F_{B}({\bf q}_{1}\\!-\\!{\bf
q}_{2};\omega\\!-\\!zJ/2)\right]$ (44)

Here, we use that $\sum_{\bf k}\gamma_{{\bf k}}^{2}\gamma_{{\bf k}+{\bf
q}_{1}}\gamma_{{\bf k}+{\bf q}_{2}}/N=[2\gamma_{{\bf q}_{1}}\gamma_{{\bf
q}_{2}}+(1-2/z)\gamma_{{\bf q}_{1}-{\bf q}_{2}}-\gamma_{{\bf q}_{1}+{\bf
q}_{2}}/z]/z^{2}$. In the final step, we repeatedly insert $\Delta B$ on the
right-hand side, starting from the initial value of $\Delta B=0$. Thus,

$\displaystyle\Delta B({\bf q}_{1},{\bf q}_{2};\omega)=\Delta B_{0}({\bf
q}_{1},{\bf q},\omega)\\!+\\!\sum_{n=1}^{\infty}(zt^{2})^{n}$
$\displaystyle\cdot\prod_{k=1}^{n}\\!G^{2}\\!\left(\omega-k\frac{zJ}{2}\right)\Delta
B_{0}\\!\left({\bf q}_{1},{\bf q}_{2};\omega\\!-\\!n\frac{zJ}{2}\right).$ (45)

To get the sought analytical result for the $B$ series, we then combine Eqs.
(45) for $\Delta B$ and (42) for $B$ with Eq. (44) for $\Delta B_{0}$, Eq.
(41) for $B_{0}$, and Eq. (40) for $F_{B}$. In all of these expressions, we
use the continued fraction form for the Green’s function in Eq. (39). Finally,
the magnetization in the neigborhood of the hole is achieved by multiplying
$B$ with the (momentum independent) residue $Z$ and then transforming to
position space using Eqs. (15) and (16).

Since $B_{0}$ depends explicitly on the polaron crystal momentum, ${\bf p}$,
we might expect that the magnetization in position space will as well. We show
now, however, that the magnetization becomes independent of the polaron
momentum just as the Green’s function [Eq. (39)]. The magnetization, $M({\bf
d})=\Delta M({\bf d})+M_{0}({\bf p},{\bf d})$, can be separated into two terms
using Eq. (42)

$\displaystyle\Delta M({\bf d})$ $\displaystyle=\frac{1}{N^{2}}\sum_{{\bf
q}_{1},{\bf q}_{2}}{\rm e}^{i({\bf q}_{1}-{\bf q}_{2})\cdot{\bf d}}\Delta
B({\bf q}_{1},{\bf q}_{2};\varepsilon_{0})$ $\displaystyle M_{0}({\bf p},{\bf
d})$ $\displaystyle=\frac{1}{N^{2}}\sum_{{\bf q}_{1},{\bf q}_{2}}{\rm
e}^{i({\bf q}_{1}-{\bf q}_{2})\cdot{\bf d}}B_{0}({\bf q}_{1},{\bf q}_{2};{\bf
p},\varepsilon_{0}).$ (46)

Here, we explicitly evaluate the functions at the quasiparticle ground state
energy, $\varepsilon_{0}=\Sigma(\varepsilon_{0})$. While the first term,
$\Delta M({\bf d})$, is explicitly independent of ${\bf p}$, the second term,
$M_{0}({\bf p},{\bf d})$, might still depend on ${\bf p}$. We now insert Eq.
(41), writing $\Delta{\bf q}={\bf q}_{1}-{\bf q}_{2}$ and using that
$B_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\varepsilon_{0})=Z\gamma_{{\bf p}+{\bf
q}_{1}}\gamma_{{\bf p}+{\bf q}_{2}}b_{0}(\Delta{\bf q};\varepsilon_{0})$.
Here, $b_{0}(\Delta{\bf
q};\varepsilon_{0})=(zt)^{2}G^{2}(\omega\\!-\\!zJ/2)\left[1+F_{B}(\Delta{\bf
q};\omega\\!-\\!zJ/2)\right]$. Consequently, we get

$\displaystyle M_{0}({\bf p},{\bf d})=\frac{Z}{N^{2}}\\!\\!\sum_{\Delta{\bf
q},{\bf q}_{1}}{\rm e}^{i\Delta{\bf q}\cdot{\bf d}}b_{0}(\Delta{\bf
q};\varepsilon_{0})\gamma_{{\bf p}+{\bf q}_{1}}\gamma_{{\bf p}+{\bf
q}_{1}-\Delta{\bf q}}$ $\displaystyle=\frac{Z}{(zN)^{2}}\\!\\!\sum_{\Delta{\bf
q},{\bf q}_{1}}\\!\\!\\!{\rm e}^{i\Delta{\bf q}\cdot{\bf d}}b_{0}(\Delta{\bf
q};\varepsilon_{0})\\!\sum_{{\boldsymbol{\delta}}_{1},{\boldsymbol{\delta}}_{2}}\\!\\!{\rm
e}^{i({\bf p}+{\bf
q}_{1})\cdot({\boldsymbol{\delta}}_{1}+{\boldsymbol{\delta}}_{2})-i\Delta{\bf
q}\cdot{\boldsymbol{\delta}}_{2}}.$

Note that $Z$ is the quasiparticle residue, while $z$ is the coordination
number. The sum over ${\bf q}_{1}$ now enforces
${\boldsymbol{\delta}}_{2}=-{\boldsymbol{\delta}}_{1}$,

$\displaystyle M_{0}({\bf d})=$
$\displaystyle\,\frac{Z}{z^{2}N}\sum_{\Delta{\bf q}}{\rm e}^{i\Delta{\bf
q}\cdot{\bf d}}b_{0}(\Delta{\bf
q};\varepsilon_{0})\sum_{{\boldsymbol{\delta}}_{1}}{\rm e}^{-i\Delta{\bf
q}\cdot{\boldsymbol{\delta}}_{1}}$ $\displaystyle=$
$\displaystyle\,\frac{Z}{zN}\sum_{\Delta{\bf q}}{\rm e}^{i\Delta{\bf
q}\cdot{\bf d}}b_{0}(\Delta{\bf q};\varepsilon_{0})\gamma_{\Delta{\bf q}}.$
(47)

This shows explicitly that the magnetization is independent of the polaron
momentum, ${\bf p}$, in the Ising case.

##  Appendix E Symmetries of the $B$ and $C$ functions

From the equations for the $B$ and $C$ function [Eqs. (16) and (17)]

$\displaystyle\\!\\!\\!\\!B({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)$
$\displaystyle=N\sum_{\bf k}\braket{\hat{h}^{\dagger}_{{\bf k}+{\bf
q}_{1}}\hat{h}_{{\bf k}+{\bf q}_{2}}\hat{b}^{\dagger}_{-{\bf
q}_{1}}\hat{b}_{-{\bf q}_{2}}}_{{\bf p},\omega},$
$\displaystyle\\!\\!\\!\\!C({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)$
$\displaystyle=N\sum_{{\bf k}}\braket{\hat{h}^{\dagger}_{{\bf k}-{\bf
q}_{1}-{\bf q}_{2}}\hat{h}_{\bf k}\hat{b}_{-{\bf q}_{1}}\hat{b}_{-{\bf
q}_{2}}}_{{\bf p},\omega}.\\!$ (48)

and the spatial symmetries of the interactions, it follows that the $B$ and
$C$ functions have 3 essential symmetry properties ($X=B,C$)

  1. 1.

Total exchange symmetry:
$X({\bf q}_{2},{\bf q}_{1};{\bf p},\omega)=X({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)$

  2. 2.

Any exchange of $(x,y,z)$ coordinates:
$X({\bf q}_{1}^{*},{\bf q}_{2}^{*};{\bf p}^{*},\omega)=X({\bf q}_{1},{\bf
q}_{2};{\bf p},\omega)$.
Example: ${\bf q}^{*}=(q_{x},q_{y},q_{z})^{*}=(q_{y},q_{x},q_{z})$.

  3. 3.

Sign flip of individual coordinates:
$X(\bar{{\bf q}}_{1},\bar{{\bf q}}_{2};\bar{{\bf p}},\omega)=X({\bf
q}_{1},{\bf q}_{2};{\bf p},\omega)$.
Example $\bar{{\bf q}}=(-q_{x},q_{y},q_{z})$.

For the $C$ function, the total exchange symmetry (1.) follows directly from
the symmetric form of $C$ in Eq. (48). For the $B$ function, we additionally
use that it is real so that $B({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=B^{*}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=N\sum_{\bf
k}\braket{(\hat{h}^{\dagger}_{{\bf k}+{\bf q}_{1}}\hat{h}_{{\bf k}+{\bf
q}_{2}}\hat{b}^{\dagger}_{-{\bf q}_{1}}\hat{b}_{-{\bf
q}_{2}})^{\dagger}}_{{\bf p},\omega}=B({\bf q}_{2},{\bf q}_{1};{\bf
p},\omega)$. The reality of $B({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)$ is a
result of the fact that all coefficients in the polaron wave function
expansion in Eq. (10) are real. This, in turn, is a consequence of the fact
that all appearing Green’s functions are always evaluated _below_ the
quasiparticle peak. The second and third symmetries reflect that there is no
preferred direction of the system. Therefore, we can swap the coordinates as
we wish (2.) and reverse 1, 2 or all 3 spatial directions (3.).

## References

  * Hubbard (1963) J. Hubbard, Proc. R. Soc. London A, Containing Papers of a Mathematical and Physical Character 276, 238 (1963).
  * Lee  _et al._ (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * Schmitt-Rink  _et al._ (1988) S. Schmitt-Rink, C. M. Varma, and A. E. Ruckenstein, Phys. Rev. Lett. 60, 2793 (1988).
  * Shraiman and Siggia (1988) B. I. Shraiman and E. D. Siggia, Phys. Rev. Lett. 61, 467 (1988).
  * Kane  _et al._ (1989) C. L. Kane, P. A. Lee, and N. Read, Phys. Rev. B 39, 6880 (1989).
  * Martinez and Horsch (1991) G. Martinez and P. Horsch, Phys. Rev. B 44, 317 (1991).
  * Liu and Manousakis (1991) Z. Liu and E. Manousakis, Phys. Rev. B 44, 2414 (1991).
  * Emery (1987) V. J. Emery, Phys. Rev. Lett. 58, 2794 (1987).
  * Schrieffer  _et al._ (1988) J. R. Schrieffer, X.-G. Wen, and S.-C. Zhang, Phys. Rev. Lett. 60, 944 (1988).
  * Dagotto (1994) E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
  * Anderson (1987) P. W. Anderson, Science 235, 1196 (1987).
  * Bonča  _et al._ (1989) J. Bonča, P. Prelovšek, and I. Sega, Phys. Rev. B 39, 7074 (1989).
  * Hasegawa and Poilblanc (1989) Y. Hasegawa and D. Poilblanc, Phys. Rev. B 40, 9035 (1989).
  * Dagotto  _et al._ (1990) E. Dagotto, R. Joynt, A. Moreo, S. Bacci, and E. Gagliano, Phys. Rev. B 41, 9049 (1990).
  * Sachdev (1989) S. Sachdev, Phys. Rev. B 39, 12232 (1989).
  * Trugman (1990) S. A. Trugman, Phys. Rev. B 41, 892 (1990).
  * Boninsegni and Manousakis (1992) M. Boninsegni and E. Manousakis, Phys. Rev. B 45, 4877 (1992).
  * Brunner  _et al._ (2000) M. Brunner, F. F. Assaad, and A. Muramatsu, Phys. Rev. B 62, 15480 (2000).
  * Mishchenko  _et al._ (2001) A. S. Mishchenko, N. V. Prokof’ev, and B. V. Svistunov, Phys. Rev. B 64, 033101 (2001).
  * Blomquist and Carlström (2020) E. Blomquist and J. Carlström, Communications Physics 3, 172 (2020).
  * White and Affleck (2001) S. R. White and I. Affleck, Phys. Rev. B 64, 024411 (2001).
  * Zhu  _et al._ (2014) Z. Zhu, H.-C. Jiang, D. N. Sheng, and Z.-Y. Weng, Scientific Reports 4, 5419 (2014).
  * Wang  _et al._ (2021) Y. Wang, A. Bohrdt, S. Ding, J. Koepsell, E. Demler, and F. Grusdt, Phys. Rev. Research 3, 033204 (2021).
  * Chen  _et al._ (2021) B.-B. Chen, C. Chen, Z. Chen, J. Cui, Y. Zhai, A. Weichselbaum, J. von Delft, Z. Y. Meng, and W. Li, Phys. Rev. B 103, L041107 (2021).
  * Bulaevskii  _et al._ (1968) L. Bulaevskii, E. Nagaev, and D. Khomskii, JETP 27, 836 (1968).
  * Brinkman and Rice (1970) W. F. Brinkman and T. M. Rice, Phys. Rev. B 2, 1324 (1970).
  * Trugman (1988) S. A. Trugman, Phys. Rev. B 37, 1597 (1988).
  * Manousakis (2007) E. Manousakis, Phys. Rev. B 75, 035106 (2007).
  * Grusdt  _et al._ (2018a) F. Grusdt, Z. Zhu, T. Shi, and E. Demler, SciPost Phys. 5, 57 (2018a).
  * Grusdt  _et al._ (2018b) F. Grusdt, M. Kánasz-Nagy, A. Bohrdt, C. S. Chiu, G. Ji, M. Greiner, D. Greif, and E. Demler, Phys. Rev. X 8, 011046 (2018b).
  * Grusdt  _et al._ (2019) F. Grusdt, A. Bohrdt, and E. Demler, Phys. Rev. B 99, 224422 (2019).
  * Bohrdt  _et al._ (2019) A. Bohrdt, C. S. Chiu, G. Ji, M. Xu, D. Greif, M. Greiner, E. Demler, F. Grusdt, and M. Knap, Nature Physics 15, 921 (2019).
  * Bohrdt  _et al._ (2021) A. Bohrdt, Y. Wang, J. Koepsell, M. Kánasz-Nagy, E. Demler, and F. Grusdt, Phys. Rev. Lett. 126, 026401 (2021).
  * Soriano and Katsnelson (2020) D. Soriano and M. I. Katsnelson, Phys. Rev. B 101, 041402(R) (2020).
  * Chao  _et al._ (1977) K. A. Chao, J. Spalek, and A. M. Oles, Journal of Physics C: Solid State Physics 10, L271 (1977).
  * Hirsch (1985) J. E. Hirsch, Phys. Rev. Lett. 54, 1317 (1985).
  * Izyumov (1997) Y. A. Izyumov, Physics-Uspekhi 40, 445 (1997).
  * Esslinger (2010) T. Esslinger, Annual Review of Condensed Matter Physics 1, 129 (2010).
  * Boll  _et al._ (2016) M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, Science 353, 1257 (2016).
  * Cheuk  _et al._ (2016) L. W. Cheuk, M. A. Nichols, K. R. Lawrence, M. Okan, H. Zhang, E. Khatami, N. Trivedi, T. Paiva, M. Rigol, and M. W. Zwierlein, Science 353, 1260 (2016).
  * Mazurenko  _et al._ (2017) A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature 545, 462 (2017).
  * Hilker  _et al._ (2017) T. A. Hilker, G. Salomon, F. Grusdt, A. Omran, M. Boll, E. Demler, I. Bloch, and C. Gross, Science 357, 484 (2017).
  * Brown  _et al._ (2017) P. T. Brown, D. Mitra, E. Guardado-Sanchez, P. Schauß, S. S. Kondov, E. Khatami, T. Paiva, N. Trivedi, D. A. Huse, and W. S. Bakr, Science 357, 1385 (2017).
  * Chiu  _et al._ (2018) C. S. Chiu, G. Ji, A. Mazurenko, D. Greif, and M. Greiner, Phys. Rev. Lett. 120, 243201 (2018).
  * Brown  _et al._ (2019) P. T. Brown, D. Mitra, E. Guardado-Sanchez, R. Nourafkan, A. Reymbaut, C.-D. Hébert, S. Bergeron, A.-M. S. Tremblay, J. Kokalj, D. A. Huse, P. Schauß, and W. S. Bakr, Science 363, 379 (2019).
  * Koepsell  _et al._ (2019) J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Nature 572, 358 (2019).
  * Chiu  _et al._ (2019) C. S. Chiu, G. Ji, A. Bohrdt, M. Xu, M. Knap, E. Demler, F. Grusdt, M. Greiner, and D. Greif, Science 365, 251 (2019).
  * Koepsell  _et al._ (2020) J. Koepsell, D. Bourgund, P. Sompet, S. Hirthe, A. Bohrdt, Y. Wang, F. Grusdt, E. Demler, G. Salomon, C. Gross, and I. Bloch, “Microscopic evolution of doped mott insulators from polaronic metal to fermi liquid,” (2020), arXiv:2009.04440 [cond-mat.quant-gas] .
  * Brown  _et al._ (2020) P. T. Brown, E. Guardado-Sanchez, B. M. Spar, E. W. Huang, T. P. Devereaux, and W. S. Bakr, Nature Physics 16, 26 (2020).
  * Vijayan  _et al._ (2020) J. Vijayan, P. Sompet, G. Salomon, J. Koepsell, S. Hirthe, A. Bohrdt, F. Grusdt, I. Bloch, and C. Gross, Science 367, 186 (2020).
  * Hartke  _et al._ (2020) T. Hartke, B. Oreg, N. Jia, and M. Zwierlein, Phys. Rev. Lett. 125, 113601 (2020).
  * Guardado-Sanchez  _et al._ (2020) E. Guardado-Sanchez, A. Morningstar, B. M. Spar, P. T. Brown, D. A. Huse, and W. S. Bakr, Phys. Rev. X 10, 011042 (2020).
  * Ji  _et al._ (2021) G. Ji, M. Xu, L. H. Kendrick, C. S. Chiu, J. C. Brüggenjürgen, D. Greif, A. Bohrdt, F. Grusdt, E. Demler, M. Lebrat, and M. Greiner, Phys. Rev. X 11, 021022 (2021).
  * Gall  _et al._ (2021) M. Gall, N. Wurz, J. Samland, C. F. Chan, and M. Köhl, Nature 589, 40 (2021).
  * Bakr  _et al._ (2009) W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M. Greiner, Nature 462, 74 (2009).
  * Sherson  _et al._ (2010) J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Nature 467, 68 (2010).
  * Haller  _et al._ (2015) E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Nature Physics 11, 738 (2015).
  * Yang  _et al._ (2021) J. Yang, L. Liu, J. Mongkolkiattichai, and P. Schauss, PRX Quantum 2, 020344 (2021).
  * Liu and Manousakis (1992) Z. Liu and E. Manousakis, Phys. Rev. B 45, 2425 (1992).
  * Chernyshev and Leung (1999) A. L. Chernyshev and P. W. Leung, Phys. Rev. B 60, 1592 (1999).
  * Reiter (1994) G. F. Reiter, Phys. Rev. B 49, 1536 (1994).
  * Ramšak and Horsch (1993) A. Ramšak and P. Horsch, Phys. Rev. B 48, 10559 (1993).
  * Ramšak and Horsch (1998) A. Ramšak and P. Horsch, Phys. Rev. B 57, 4308 (1998).
  * Bała and Oleś (2000) J. Bała and A. M. Oleś, Phys. Rev. B 61, 6907 (2000).
  * Bała  _et al._ (2002) J. Bała, A. M. Oleś, and P. Horsch, Phys. Rev. B 65, 134420 (2002).
  * Porras and Cirac (2004) D. Porras and J. I. Cirac, Phys. Rev. Lett. 92, 207901 (2004).
  * Gorshkov  _et al._ (2011) A. V. Gorshkov, S. R. Manmana, G. Chen, J. Ye, E. Demler, M. D. Lukin, and A. M. Rey, Phys. Rev. Lett. 107, 115301 (2011).
  * Britton  _et al._ (2012) J. W. Britton, B. C. Sawyer, A. C. Keith, C. C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Nature 484, 489 (2012).
  * Zeiher  _et al._ (2016) J. Zeiher, R. van Bijnen, P. Schauß, S. Hild, J.-y. Choi, T. Pohl, I. Bloch, and C. Gross, Nature Physics 12, 1095 (2016).
  * Zeiher  _et al._ (2017) J. Zeiher, J.-y. Choi, A. Rubio-Abadal, T. Pohl, R. van Bijnen, I. Bloch, and C. Gross, Phys. Rev. X 7, 041063 (2017).
  * Sakurai (1994) J. J. Sakurai,  _Modern Quantum Mechanics_ (Addison-Wesley, Reading, MA, 1994) p. 274.
  * Wells  _et al._ (1995) B. O. Wells, Z. X. Shen, A. Matsuura, D. M. King, M. A. Kastner, M. Greven, and R. J. Birgeneau, Phys. Rev. Lett. 74, 964 (1995).
  * Ronning  _et al._ (1998) F. Ronning, C. Kim, D. L. Feng, D. S. Marshall, A. G. Loeser, L. L. Miller, J. N. Eckstein, I. Bozovic, and Z.-X. Shen, Science 282, 2067 (1998).
  * Kim  _et al._ (1998) C. Kim, P. J. White, Z.-X. Shen, T. Tohyama, Y. Shibata, S. Maekawa, B. O. Wells, Y. J. Kim, R. J. Birgeneau, and M. A. Kastner, Phys. Rev. Lett. 80, 4245 (1998).
  * Yoshida  _et al._ (2003) T. Yoshida, X. J. Zhou, T. Sasagawa, W. L. Yang, P. V. Bogdanov, A. Lanzara, Z. Hussain, T. Mizokawa, A. Fujimori, H. Eisaki, Z.-X. Shen, T. Kakeshita, and S. Uchida, Phys. Rev. Lett. 91, 027001 (2003).
  * Ronning  _et al._ (2005) F. Ronning, K. M. Shen, N. P. Armitage, A. Damascelli, D. H. Lu, Z.-X. Shen, L. L. Miller, and C. Kim, Phys. Rev. B 71, 094518 (2005).
  * Rösch  _et al._ (2005) O. Rösch, O. Gunnarsson, X. J. Zhou, T. Yoshida, T. Sasagawa, A. Fujimori, Z. Hussain, Z.-X. Shen, and S. Uchida, Phys. Rev. Lett. 95, 227002 (2005).
  * Schrieffer and Brooks (2007) J. R. Schrieffer and J. S. Brooks,  _Handbook of High-Temperature Superconductivity_ (Springer, New York, 2007).
  * Gross and Bloch (2017) C. Gross and I. Bloch, Science 357, 995 (2017).
  * Schäfer  _et al._ (2020) F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu, and Y. Takahashi, Nature Reviews Physics 2, 411 (2020).
  * Glaetzle  _et al._ (2015) A. W. Glaetzle, M. Dalmonte, R. Nath, C. Gross, I. Bloch, and P. Zoller, Phys. Rev. Lett. 114, 173002 (2015).
  * van Bijnen and Pohl (2015) R. M. W. van Bijnen and T. Pohl, Phys. Rev. Lett. 114, 243002 (2015).
  * Borish  _et al._ (2020) V. Borish, O. Marković, J. A. Hines, S. V. Rajagopal, and M. Schleier-Smith, Phys. Rev. Lett. 124, 063601 (2020).
  * Guardado-Sanchez  _et al._ (2021) E. Guardado-Sanchez, B. M. Spar, P. Schauss, R. Belyansky, J. T. Young, P. Bienias, A. V. Gorshkov, T. Iadecola, and W. S. Bakr, Phys. Rev. X 11, 021036 (2021).
  * Kyung and Ferrell (1996) B. Kyung and R. A. Ferrell, Phys. Rev. B 54, 10125 (1996).
  * Xiang and Wheatley (1996) T. Xiang and J. M. Wheatley, Phys. Rev. B 54, R12653 (1996).
  * Singh and Gelfand (1995) R. R. P. Singh and M. P. Gelfand, Phys. Rev. B 52, R15695 (1995).
  * Kim  _et al._ (1999) Y. J. Kim, A. Aharony, R. J. Birgeneau, F. C. Chou, O. Entin-Wohlman, R. W. Erwin, M. Greven, A. B. Harris, M. A. Kastner, I. Y. Korenblit, Y. S. Lee, and G. Shirane, Phys. Rev. Lett. 83, 852 (1999).
  * Sandvik and Singh (2001) A. W. Sandvik and R. R. P. Singh, Phys. Rev. Lett. 86, 528 (2001).
  * Rønnow  _et al._ (2001) H. M. Rønnow, D. F. McMorrow, R. Coldea, A. Harrison, I. D. Youngson, T. G. Perring, G. Aeppli, O. Syljuåsen, K. Lefmann, and C. Rischel, Phys. Rev. Lett. 87, 037202 (2001).
  * Coldea  _et al._ (2001) R. Coldea, S. M. Hayden, G. Aeppli, T. G. Perring, C. D. Frost, T. E. Mason, S.-W. Cheong, and Z. Fisk, Phys. Rev. Lett. 86, 5377 (2001).
  * Christensen  _et al._ (2004) N. Christensen, D. McMorrow, H. Rønnow, A. Harrison, T. Perring, and R. Coldea, Journal of Magnetism and Magnetic Materials 272-276, 896 (2004).
  * Chernyshev  _et al._ (1994) A. L. Chernyshev, A. V. Dotsenko, and O. P. Sushkov, Phys. Rev. B 49, 6197 (1994).
  * Belinicher  _et al._ (1997) V. I. Belinicher, A. L. Chernyshev, and V. A. Shubin, Phys. Rev. B 56, 3381 (1997).
  * Shraiman and Siggia (1989) B. I. Shraiman and E. D. Siggia, Phys. Rev. B 40, 9162 (1989).
  * Frenkel and Hanke (1990) D. M. Frenkel and W. Hanke, Phys. Rev. B 42, 6711 (1990).
  * Eder (1992) R. Eder, Phys. Rev. B 45, 319 (1992).
  * Riera and Dagotto (1998) J. Riera and E. Dagotto, Phys. Rev. B 57, 8609 (1998).
  * Bohrdt  _et al._ (2020) A. Bohrdt, F. Grusdt, and M. Knap, New Journal of Physics 22, 123023 (2020).
  * Hubig  _et al._ (2020) C. Hubig, A. Bohrdt, M. Knap, F. Grusdt, and J. I. Cirac, SciPost Phys. 8, 21 (2020).
  * Starykh and Reiter (1996) O. A. Starykh and G. F. Reiter, Phys. Rev. B 53, 2517 (1996).