compat=1.0.0 compat=1.0.0

# Wave function and spatial structure of polarons in an antiferromagnetic
bilayer

J. H. Nyhegn  Center for Complex Quantum Systems, Department of Physics and
Astronomy, Aarhus University, Ny Munkegade, DK-8000 Aarhus C, Denmark  G. M.
Bruun  Center for Complex Quantum Systems, Department of Physics and
Astronomy, Aarhus University, Ny Munkegade, DK-8000 Aarhus C, Denmark  K.
Knakkergaard Nielsen  Max-Planck Institute for Quantum Optics, Hans-
Kopfermann-Str. 1, D-85748 Garching, Germany

###### Abstract

Adding a dopant to an antiferromagnetic spin background disturbs the magnetic
order and leads to the formation of a quasiparticle coined the magnetic
polaron, which plays a central role in understanding strongly correlated
materials. Recently, remarkably detailed insights into the spatial properties
of such polarons have been obtained using atoms in optical lattices. Motivated
by this we develop a nonperturbative scheme for calculating the wave function
of the magnetic polaron in a bilayer antiferromagnet using the self-consistent
Born approximation. The scheme includes an infinite number of spin waves,
which is crucial for an accurate description of the most interesting regime of
strong correlations. Utilizing the developed wave function, we explore the
spatial structure of the polaron dressing cloud consisting of magnetically
frustrated spins surrounding the hole. Mimicking the nonmonotonic behavior of
the antiferromagnetic order, we find that the dressing cloud first decreases
and then increases in size with increasing interlayer hopping. The increase
reflects the decrease in the magnetic order as a quantum phase transition to a
disordered state is approached for large interlayer hopping. We, furthermore,
find that the symmetry of the ground state dressing cloud changes as the
interlayer coupling increases. Our results should be experimentally accessible
using quantum simulation with optical lattices.

##  I Introduction

The properties of dopants in antiferromagnetic (AFM) environments have been
intensely studied both experimentally and theoretically for decades. The
reasons for this are at least two-fold. First, the dopants are believed to be
the charge carriers and fundamental building blocks for understanding
unconventional superconducting phases in cuprates [1]. Second, the
superconducting phases in cuprates as well as in pnictides [2], organic layers
[3], and twisted bilayer graphene [4] all reside close to the
antiferromagnetically ordered phase. The mobile charge carriers can,
therefore, provide important clues to the underlying physics behind
unconventional superconducting phases. The essential properties of these
systems are captured by the Fermi-Hubbard model and the associated $t$–$J$
model. Breakthroughs using quantum simulation of the Fermi-Hubbard model with
ultracold atoms [5, 6, 7, 8, 9, 10, 11, 12, 13, 12, 14, 15, 16] have
stimulated new interest in this fundamental topic [17, 18, 19, 20, 21, 22]. In
particular, the impressive single-site resolution [23, 24, 25, 26] of such
experiments makes it possible to measure spatial correlations, giving insights
into the spatial structure of the emergent quasiparticles [14] and their
dynamical formation [15], which are well beyond the reach of condensed matter
experiments. Recently, progress towards studying bilayer systems using optical
lattices has been reported [27]. This is not only interesting due to the
relation to cuprates with superconducting blocks of two CuO2 planes, but also
because bilayer systems exhibit a quantum phase transition to a disordered
state of interlayer spin singlets [28, 29, 30, 31, 32, 33].

Figure 1:  The magnetic dressing cloud of the polaron ground state for
$t_{\perp}/t=0.5$ (left) and $t_{\perp}/t=3.25$ (right) with a hole in the
center of layer 1. Here, $\tilde{M}=\pm 1$ corresponds to the average
magnetization predicted by linear spin wave theory (LSWT). For
$t_{\perp}/t=0.5$, the ground state momentum is ${\bf p}=(\pi/2,\pi/2)$
indicated by the arrow, and the cloud has $C_{2v}$ symmetry. Increasing
$t_{\perp}/t$ to $3.25$ changes the ground state momentum to ${\bf p}={\bf
0}$, where the magnetization cloud recovers the full $C_{4v}$ symmetry of the
antiferromagnet. In this limit, the spins in the dressing clouds in the two
layers are approximately antiparallel.

On the theory side, the self-consistent Born approximation (SCBA) has proven
to yield accurate descriptions of single dopants in these systems, whether it
be their equilibrium [34, 35, 36, 37, 38, 39, 40, 21] or nonequilibrium
properties [22]. The ability to describe the spatial structure of the emergent
magnetic polarons [37, 21] hinges in this framework on the analytical
expression for the many-body wave function [41, 22].

Motivated by the exciting progress regarding the experimental study of bilayer
systems, we generalize the SCBA approach for a single layer [41, 37, 21], and
derive an analytical expression for the magnetic polaron many-body wave
function for the bilayer geometry. The wave function includes an infinite
number of spin waves, which is important in the strongly correlated regime.
This, in turn, allows us to calculate the magnetic dressing cloud around the
hole by self-consistently solving a set of Dyson-like equations. We use this
to explore how the dressing cloud is modified by the interlayer interaction.
We find that the size of the ground state dressing cloud exhibits a
nonmonotonic behavior with increasing hopping between the two layers and that
it, additionally, changes its spatial symmetry (see Fig. 1).

The paper is organized as follows. In Sec. II, we introduce the $t$-$J$ model
and discuss how we retrieve an effective model in terms of the hole and spin
waves. The SCBA is then described and the quasiparticle wave function is
constructed in Sec. III. Using the quasiparticle wave function, we derive a
set of Dyson-like equations used to calculate the spin frustration surrounding
the hole in Sec. IV. In Sec. V, we explore the spatial structure of the
quasiparticle as the interlayer coupling is varied before we conclude in Sec.
VI.

##  II Bilayer $t$–$J$ model

To model the behavior of a single hole in an AFM bilayer, we use the $t$–$J$
model with the Hamiltonian

$\hat{H}=\hat{H}_{t}+\hat{H}_{J},$ (1)

where

$\hat{H}_{t}=-t\sum_{l,\braket{{\bf i},{\bf
j}},\sigma}\tilde{c}^{\dagger}_{l,{\bf i},\sigma}\tilde{c}_{l,{\bf
j},\sigma}-t_{\perp}\sum_{{\bf i},\sigma}\tilde{c}^{\dagger}_{1,{\bf
i},\sigma}\tilde{c}_{2,{\bf i},\sigma}+\text{H.c.}$ (2)

is the hopping Hamiltonian, and

$\displaystyle\hat{H}_{J}=$ $\displaystyle\ J\sum_{l,\braket{{\bf i},{\bf
j}}}\left[\hat{\mathbf{S}}_{l,{\bf i}}\cdot\hat{\mathbf{S}}_{l,{\bf
j}}-\frac{1}{4}\hat{n}_{l,{\bf i}}n_{l,{\bf j}}\right]$ $\displaystyle\
+J_{\perp}\sum_{{\bf i}}\left[\hat{\mathbf{S}}_{1,{\bf
i}}\cdot\hat{\mathbf{S}}_{2,{\bf i}}-\frac{1}{4}\hat{n}_{1,{\bf
i}}\hat{n}_{2,{\bf i}}\right]$ (3)

describes spin coupling within and between the layers. The subscript $l$ is an
index for the top ($l=1$) and bottom ($l=2$) layers, $J,J_{\perp}>0$ are the
antiferromagnetic couplings between the spins, and $t,t_{\perp}$ are the
intra- and interlayer hopping amplitudes. The restricted fermionic operators,
$\tilde{c}^{\dagger}_{i,\sigma}=\hat{c}^{\dagger}_{j,\sigma}(1-\hat{n}_{i,\bar{\sigma}})$,
are constructed so that no double occupancies are allowed. The isotropic spin
coupling in Eq. (3) naturally arises as an effective low-energy description of
superexchange processes in the Fermi-Hubbard model at large onsite repulsions
$U$, for which $J=4t^{2}/U$ and $J_{\perp}=4t_{\perp}^{2}/U$. As this is
directly relevant to current experiments using atoms in optical lattices as
well as strongly correlated two-dimensional (2D) materials, we will assume
these relations in the following such that $J_{\perp}/J=(t_{\perp}/t)^{2}$.

In the absence of holes, the bilayer system undergoes a quantum phase
transition at $J_{\perp}/J\simeq 2.5$ [29] corresponding to $t_{\perp}/t\simeq
1.6$, from an AFM ordered phase at low $J_{\perp}/J$ to a disordered phase at
large $J_{\perp}/J$, in which neighboring spins in the two layers form spin
singlets [30, 42, 43]. Deep in the AFM phase, the low energy states can be
quantitatively described using spin-wave theory [33, 44], and we will,
furthermore, use this framework to qualitatively describe magnetic polarons in
the region close to the quantum critical point (QCP).

###  II.1 Linear spin-wave theory

The AFM order defines two sublattices $A$ and $B$ where the spins
predominantly point up and down. Using a Holstein-Primakoff transformation
generalized to take the presence of holes into account [1, 34, 35, 45], the
mapping for sublattice $A$ is

$\displaystyle\hat{S}^{z}_{l,{\mathbf{i}}}$
$\displaystyle=(1-\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{h}_{l,{\mathbf{i}}})/2-\hat{s}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}},$
$\displaystyle\hat{S}_{l,{\mathbf{i}}}^{-}$
$\displaystyle=\hat{S}^{x}_{l,{\mathbf{i}}}-i\hat{S}^{y}_{l,{\mathbf{i}}}=\hat{s}^{\dagger}_{l,{\mathbf{i}}}(1-\hat{s}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}}-\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{h}_{l,{\mathbf{i}}})^{1/2},$
$\displaystyle\tilde{c}_{i,{\downarrow}}$
$\displaystyle=\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}},\
\tilde{c}_{i,{\uparrow}}=\hat{h}^{\dagger}_{l,{\mathbf{i}}}(1-\hat{s}^{\dagger}_{l,{\mathbf{i}}}\hat{s}_{l,{\mathbf{i}}}-\hat{h}^{\dagger}_{l,{\mathbf{i}}}\hat{h}_{l,{\mathbf{i}}})^{1/2},$
(4)

with $\hat{h}^{\dagger}_{l,{\mathbf{i}}}$ being a fermionic creation operator
of a hole and $\hat{s}^{\dagger}_{l,{\mathbf{i}}}$ a bosonic creation operator
of a spin excitation. Similar expressions apply to sublattice B. Applying the
generalized Holstein-Primakoff transformation, keeping only the linear terms,
Fourier transforming, and using the mapping

$\displaystyle\begin{bmatrix}\hat{s}_{1,{\bf k}}\\\ \hat{s}^{\dagger}_{1,-{\bf
k}}\\\ \hat{s}_{2,{\bf k}}\\\ \hat{s}^{\dagger}_{2,-{\bf
k}}\end{bmatrix}=\frac{1}{\sqrt{2}}\begin{bmatrix}\mathds{1}&-\mathds{1}\\\
\mathds{1}&\mathds{1}\end{bmatrix}\begin{bmatrix}U_{+,{\bf k}}&0\\\
0&U_{-,{\bf k}}\end{bmatrix}\begin{bmatrix}\hat{b}_{+,{\bf k}}\\\
\hat{b}^{\dagger}_{+,-{\bf k}}\\\ \hat{b}_{-,{\bf k}}\\\
\hat{b}^{\dagger}_{-,-{\bf k}}\end{bmatrix},$ (5)

it has previously been shown [45] that the Hamiltonian transforms into

$\hat{H}_{J}=\ E_{0}+\sum_{{\bf k},\mu=\pm}\omega_{\mu,{\bf
k}}\hat{b}^{\dagger}_{\mu,{\bf k}}\hat{b}_{\mu,{\bf k}}$ (6)

and

$\displaystyle\hat{H}_{t}=$ $\displaystyle\sum_{{\bf k},{\bf
p},l}\hat{h}_{l,{\bf p}+{\bf k}}^{\dagger}\hat{h}_{l,{\bf p}}[g_{+}({\bf
p},{\bf k})\hat{b}^{\dagger}_{+,-{\bf k}}+(-1)^{l}g_{-}({\bf p},{\bf
k})\hat{b}^{\dagger}_{-,-{\bf k}}]$ $\displaystyle+$ $\displaystyle\sum_{{\bf
k},{\bf p}}\hat{h}_{2,{\bf p}+{\bf k}}^{\dagger}\hat{h}_{1,{\bf p}}[f_{+}({\bf
k})(\hat{b}_{+,{\bf k}}+\hat{b}^{\dagger}_{+,-{\bf k}})$
$\displaystyle\phantom{\sum_{{\bf k},{\bf p}}\hat{h}_{2,{\bf p}+{\bf
k}}^{\dagger}\hat{h}_{1,{\bf p}}[\\!\\!\\!\\!\\!\\!}+f_{-}({\bf
k})(\hat{b}_{-,{\bf k}}-\hat{b}^{\dagger}_{-,-{\bf k}})]+\text{H.c.}$ (7)

Here, we only retain the lowest order terms in the operators
$\hat{b}^{\dagger}_{\mu,{\bf k}}$, describing bosonic spin waves in the branch
$\mu=\pm$ with dispersions

$\omega_{\pm,{\bf
k}}=\frac{1}{2}\sqrt{\left(Jz+J_{\perp}\right)^{2}-\left(Jz\gamma_{{\bf k}}\pm
J_{\perp}\right)^{2}}.$ (8)

The intra- and interlayer interaction vertices are

$\displaystyle g_{\pm}({\bf p},{\bf k})$
$\displaystyle=\frac{zt}{\sqrt{2N}}\left[u_{\pm,{\bf k}}\gamma_{{\bf p}+{\bf
k}}-v_{\pm,{\bf k}}\gamma_{\bf p}\right],$ $\displaystyle f_{\pm}({\bf k})$
$\displaystyle=\frac{t_{\perp}}{\sqrt{2N}}\left[u_{\pm,{\bf k}}\mp v_{\pm,{\bf
k}}\right],$ (9)

with $z=4$ being the lattice coordinate number, $\gamma_{{\bf
k}}=2(\cos{k_{x}}+\cos{k_{y}})/z$ the structure factor, and

$\displaystyle u_{\pm,{\bf k}}$
$\displaystyle=\sqrt{\frac{1}{2}\left(\frac{zJ+J_{\perp}}{2\omega_{\pm,{\bf
k}}}+1\right)}$ $\displaystyle v_{\pm,{\bf k}}$ $\displaystyle={\rm
sgn}\left[zJ\gamma_{\bf k}\pm
J_{\perp}\right]\sqrt{\frac{1}{2}\left(\frac{zJ+J_{\perp}}{2\omega_{\pm,{\bf
k}}}-1\right)}$ (10)

the coherence factors obtained from diagonalizing $\hat{H}_{J}$. The lattice
constant is taken to be unity. In the derivation of the quasiparticle wave
function, the exact expressions for the interaction vertices are not
important, but their symmetries will greatly simplify the expressions. These
symmetries are related to the AFM ordering vector $\bf{Q}=(\pi,\pi)$, and for
the coherence factors and the spin wave dispersions, we have

$\displaystyle u_{+,{\bf k}+{\bf Q}}$ $\displaystyle=u_{-,{\bf k}},\;v_{+,{\bf
k}+{\bf Q}}=-v_{-,{\bf k}},\;\omega_{+,{\bf k}+{\bf Q}}=\omega_{-,{\bf k}},$
(11)

whereas the interaction vertices fulfill

$\displaystyle g_{\pm}({\bf p},{\bf k}+{\bf Q})$ $\displaystyle=-g_{\mp}({\bf
p},{\bf k}),\;g_{\pm}({\bf p}+{\bf Q},{\bf k})=-g_{\pm}({\bf p},{\bf k})$
$\displaystyle f_{\pm}({\bf k}+{\bf Q})$ $\displaystyle=f_{\mp}({\bf k}).$
(12)

Equation (7) constitutes our effective Hamiltonian found within LSWT. For the
bilayer Heisenberg model this approximation scheme has been shown to neglect a
longitudinal mode [33]. To consider this mode, one could instead tend to the
bond-operator representation [44, 46, 47, 48]. That said, the result presented
in Sec. V, when approaching the QCP, reflects the decrease in AFM order which
is qualitatively well captured within LSWT, but overestimates the transition
point to lie around $t_{\perp}/t\simeq 3.67$ as can generally be expected from
a mean-field theory. With the effective Hamiltonian in place, we turn to the
quantum field theory approach applied previously [45] to describe the
quasiparticle properties of the emergent magnetic polarons.

###  II.2 Quantum field theory

With the effective Hamiltonian in Eq. (7), the Matsubara Green’s function of a
hole takes on a $2\times 2$ matrix structure [45]

$G_{lm}({\bf p},\tau)=-\braket{T_{\tau}[\hat{h}_{l,{\bf
p}}(\tau)\hat{h}^{\dagger}_{m,{\bf p}}]},$ (13)

where $T_{\tau}$ is the time ordering in imaginary time $\tau$ and $l,m=1,2$
denote the two layers. In frequency space, the Dyson equation reads
$G(p)=\mathds{1}G_{0}(p)+G_{0}(p)\Sigma(p)G(p)$, where $\Sigma=[\Sigma_{lm}]$
is the self-energy matrix and $p=({\bf p},i\omega_{p})$ with $i\omega_{p}$ a
fermionic Matsubara frequency. Finally, $G_{0}(p)=1/i\omega_{p}$ is the
noninteracting hole Green’s function. Due to the layer symmetry, there are
only two independent matrix elements for the Green’s function: $G_{{\rm
d}}=G_{11}=G_{22}$ and $G_{o}=G_{12}=G_{21}$ or equivalently
$\Sigma_{d}=\Sigma_{11}=\Sigma_{22}$ and $\Sigma_{o}=\Sigma_{12}=\Sigma_{21}$.
Using this, it follows from the Dyson equation that

$\displaystyle G_{\rm d}(p)=\frac{G_{0}^{-1}(p)-\Sigma_{\rm
d}(p)}{\left(G_{0}^{-1}(p)-\Sigma_{\rm d}(p)\right)^{2}-\left(\Sigma_{\rm
o}(p)\right)^{2}}$ $\displaystyle G_{\rm o}(p)=\frac{\Sigma_{\rm
o}(p)}{\left(G_{0}^{-1}(p)-\Sigma_{\rm d}(p)\right)^{2}-\left(\Sigma_{\rm
o}(p)\right)^{2}}.$ (14)

In order to calculate the self-energies, we employ the SCBA to sum the class
of so-called rainbow diagrams as illustrated in Fig. 2. This includes an
infinite number of spin waves that are absorbed in the same order as they are
emitted and turns out to provide an accurate description of hole dynamics in
the strong coupling limit $t\gg J$ that naturally arises when $U\gg t$ [21].

Figure 2: Feynman diagrams for $\Sigma_{\text{d}}$ within the SCBA. The green
color is associated with layer $l=1$ and yellow with $l=2$, such that the all
green and all yellow propagators represent $G_{{\rm d}}$ while the mixed green
and yellow propagator represent $G_{{\rm o}}$. The green (yellow) vertex
describes a hole in layer $l=1$ $(l=2)$ interacting with a spin wave and
continuing in the same layer. Instead, the red vertex describes interlayer
interactions where the hole jumps from one layer to the other. The off-
diagonal self-energy $\Sigma_{\rm o}$ is given by similar diagrams [45].

Using the vertex symmetries in Eq. (12), we obtain [45]

$\displaystyle\Sigma({\bf p},\omega)=$ $\displaystyle\ \sum_{{\bf
k},\mu}V_{\mu}({\bf p},{\bf k})G({\bf p}+{\bf k},\omega-\omega_{\mu,{\bf
k}})V^{\dagger}_{\mu}({\bf p},{\bf k})$ $\displaystyle=$ $\displaystyle\
2\sum_{{\bf k}}V_{+}({\bf p},{\bf k})G({\bf p}+{\bf k},\omega-\omega_{+,{\bf
k}})V^{\dagger}_{+}({\bf p},{\bf k}),$ (15)

with the coupling matrix

$\displaystyle V_{\mu}({\bf p},{\bf k})=$ $\displaystyle\begin{bmatrix}\mu
g_{+}({\bf p},{\bf k})&\mu f_{+}({\bf k})\\\ f_{+}({\bf k})&g_{+}({\bf p},{\bf
k})\end{bmatrix}.$ (16)

We iteratively solve Eqs. (14) and (15) self-consistently starting from
$\Sigma=0$. This yields the hole Green’s function with poles determining the
energy of the magnetic polarons, which due to the layer symmetry come in two
classes: either symmetric or antisymmetric under layer exchange [45].

In Fig. 3, we plot the crystal momentum of the lowest energy polaron as a
function of $t_{\perp}/t$. Here and in the rest of the paper we use $J/t=0.3$.
When the layers are decoupled, i.e. $t_{\perp}/t=0$, the magnetic polaron
ground state of the system has crystal momentum ${\bf p}=(\pi/2,\pi/2)$. For
increasing interlayer hopping, however, Fig. 3 shows that the ground state
momentum moves from ${\bf p}=(\pi/2,\pi/2)$ along the diagonal and settles at
${\bf p}={\bf 0},(\pi,\pi)$ for $t_{\perp}/t\gtrsim 2.3$ [45, 44]. As we will
show in Sec. V, this is essential for understanding the spatial structure of
the ground state as the coupling between the layers changes. The insets in
Fig. 3 show the energy dispersion of the magnetic polaron in the Brillouin
zone for different values of $t_{\perp}/t$ with the momentum of the ground
state indicated.

Figure 3: The ground state crystal momentum versus $t_{\perp}/t$ for
$J/t=0.3$. The yellow and green markers correspond to the symmetric and
antisymmetric solution respectively and they both represent the degenerate
ground state. The insets show the dispersion of the polaron in the first
quadrant of the Brillouin zone at written values of $t_{\perp}$ and with the
color bars indicating $\omega/t$. The energy minima are indicated by the
yellow and green markers which represent the symmetric and antisymmetric
solution respectively. Around $t_{\perp}/t\simeq 2.3$, the crystal momentum of
the ground state settles at ${\bf p}=\mathbf{0}$ for the antisymmetric
solution and ${\bf p}=(\pi,\pi)$ for the symmetric.

##  III Quasiparticle wave function

While quantities such as the quasiparticle energy and residue are naturally
obtained from the field theory described above, it is much less suited to
calculate observables related to the real-space structure of the magnetic
polaron. This is a serious limitation given the detailed spatial information
provided by the optical lattice experiments [14, 15]. Motivated by this, we
now generalize a nonperturbative scheme to calculate the many-body wave
function of the magnetic polaron within the SCBA from a single layer [41, 21]
to the bilayer case.

The wave function of the magnetic polaron can in general be expanded in the
number of spin waves it has excited in the antiferromagnet. This gives

$\displaystyle\ket{\Psi_{\bf p}}=$ $\displaystyle\sum_{l}\Big{[}a^{(l)}({\bf
p})\hat{h}^{\dagger}_{l,{\bf p}}+\sum_{\mu_{1},{\bf
k}_{1}}a^{(l)}_{\mu_{1}}({\bf p},{\bf k}_{1})\hat{h}^{\dagger}_{l,{\bf p}+{\bf
k}_{1}}\hat{b}^{\dagger}_{\mu_{1},-{\bf k}_{1}}$
$\displaystyle+\dots\Big{]}\ket{{\rm AF}},$ (17)

where $\ket{{\rm AF}}$ is the AFM ground state and $\ldots$ represents states
with multiple spin waves excited. The coefficient
$a^{(l)}_{\\{\mu_{n}\\}}({\bf p},\\{{\bf k}_{n}\\})$ is the amplitude for
having $n$ spin waves excited in the branches
$\\{\mu_{n}\\}=\mu_{1},\dots,\mu_{n}$ with crystal momenta $\\{{\bf
k}_{n}\\}=-{\bf k}_{1},\dots,-{\bf k}_{n}$. The $i$’th spin wave is, hereby,
excited into the branch $\mu_{i}$ with crystal momentum $-{\bf k}_{i}$. The
superscript, $l=1,2$ defines the layer in which the hole resides.

###  III.1 Recursive determination of amplitudes

To find explicit expressions for the amplitudes, we use the Schrödinger
equation $(\omega-\hat{H})\ket{\Psi_{{\bf p}}}=0$ for the polaron with energy
$\omega$. The Hamiltonian given by Eq. (7) can only excite and annihilate a
single spin wave at a time, any coefficient is, therefore, only coupled to the
coefficients for having one more or one less spin wave excited. In addition,
when taking matrix elements of the Schrödinger equation the SCBA dictates that
the order of the spin waves is kept fixed so that only a single Wick
contraction for a particular coefficient is maintained. This reduces the
complexity of the many-body wave function while retaining the possibility of
describing a highly correlated state with many spin waves excited.

With this approach, we derive a set of equations of motion, which at order
$n\geq 1$ has the matrix form

$\displaystyle\left[\omega-\sum_{i=1}^{n}\omega_{\mu_{i},{\bf
k}_{i}}\right]{\bf a}_{\\{\mu_{n}\\}}({\bf p},\\{{\bf k}_{n}\\})$
$\displaystyle=V^{\dagger}_{\mu_{n}}({\bf K}_{n-1},{\bf k}_{n}){\bf
a}_{\\{\mu_{n-1}\\}}({\bf p},\\{{\bf k}_{n-1}\\})$
$\displaystyle+\\!\sum_{\begin{subarray}{c}{\bf k}_{n+1}\\\
\mu_{n+1}\end{subarray}}\\!V_{\mu_{n+1}}({\bf K}_{n},{\bf k}_{n+1}){\bf
a}_{\\{\mu_{n+1}\\}}({\bf p},\\{{\bf k}_{n+1}\\}).$ (18)

Here, ${\bf a}=[a^{(1)},a^{(2)}]$ is a vector containing the wave-function
coefficients for the hole in each layer and ${\bf K}_{n}\equiv{\bf
p}+\sum_{i=1}^{n}{\bf k}_{i}$. In close analogy to the single-layer case [41,
37, 21], we find that the recursive relation

$\displaystyle{\bf a}_{\\{\mu_{n+1}\\}}({\bf p},\\{{\bf
k}_{n+1}\\})=G\big{(}{\bf K}_{n+1},\omega-\sum_{i=1}^{n+1}\omega_{\mu_{i},{\bf
k}_{i}}\big{)}$ $\displaystyle\times V^{\dagger}_{\mu_{n+1}}({\bf K}_{n},{\bf
k}_{n+1}){\bf a}_{\\{\mu_{n}\\}}({\bf p},\\{{\bf k}_{n}\\})$ (19)

self-consistently solves the equations of motion as can straightforwardly be
checked by inserting Eq. (19) into Eq. (18) (see Appendix A for details).
Equation (19) gives for $n=0$ the matrix equation

$\displaystyle\left[\omega\cdot\mathds{1}-\Sigma({\bf p},\omega)\right]{\bf
a}({\bf p})=0,$ (20)

which shows that $\omega$ must coincide with one of the two poles
$\epsilon^{\pm}_{\bf p}$ of the Green’s functions in Eq. (14) in order to have
a nonzero solution of the coefficients. Thus, the polaron energies obtained
from the wave function are the same as from the Green’s function. This
explicitly demonstrates the consistency of the two approaches. The zero-order
coefficients for the symmetric or antisymmetric polaron with wave function
$\ket{\Psi^{\pm}_{{\bf p}}}$ are $a^{(2)}({\bf p})=\pm a^{(1)}({\bf p})$
solving $a^{(1)}({\bf p})[\omega-\Sigma_{\rm d}({\bf p},\omega)\mp\Sigma_{\rm
o}({\bf p},\omega)]=0$. All higher-order coefficients can now be calculated
iteratively in terms of the zero-order coefficients.

To each order, there will be $2^{n+1}$ terms and an equal number of
coefficients. With the vertex symmetries in Eq. (12), we can simplify the
number of coefficients to a single unique coefficient for each order, which
reduces Eq. (19) to a scalar relation (see Appendix C). After normalizing the
wave function (see Appendix B) the polaron wave function to first order is

$\displaystyle\ket{\Psi^{\pm}_{{\bf p}}}=\sqrt{\frac{{Z_{{\bf
p}}^{\pm}}}{2}}\left(\hat{h}^{\dagger}_{1,{\bf p}}\pm\hat{h}^{\dagger}_{2,{\bf
p}}+\sum_{{\bf k}_{1}}R^{\pm}({\bf p},{\bf
k}_{1})\left[\left(\hat{h}^{\dagger}_{1,{\bf p}+{\bf
k}_{1}}\pm\hat{h}^{\dagger}_{2,{\bf p}+{\bf
k}_{1}}\right)\hat{b}^{\dagger}_{+,-{\bf
k}_{1}}+\left(\hat{h}^{\dagger}_{1,{\bf p}+{\bf k}_{1}+{\bf
Q}}\mp\hat{h}^{\dagger}_{2,{\bf p}+{\bf k}_{1}+{\bf
Q}}\right)\hat{b}^{\dagger}_{-,-{\bf k}_{1}-{\bf
Q}}\right]+...\right)\ket{{\rm AF}},$ (21)

where

$\displaystyle R^{\pm}({\bf p},\\{{\bf k}_{n}\\})=$ $\displaystyle[G_{\rm
d}({\bf K}_{n},\Omega^{\pm}_{\\{{\bf k}_{n}\\}})\pm G_{\rm o}({\bf
K}_{n},\Omega^{\pm}_{\\{{\bf k}_{n}\\}})]$ $\displaystyle\times[g_{+}({\bf
K}_{n-1},{\bf k}_{n})\pm f_{+}({\bf k}_{n})]$ (22)

and

$\Omega^{\pm}_{\\{{\bf k}_{n}\\}}({\bf p})=\epsilon^{\pm}_{\bf
p}-\sum_{i=1}^{n}\omega_{+,{\bf k}_{n}}.$ (23)

Here, $Z^{\pm}_{{\bf p}}=[1-\partial_{\omega}\left(\Sigma_{{\rm d}}({\bf
p},\omega)\pm\Sigma_{{\rm o}}({\bf
p},\omega)\right)|_{\omega=\epsilon^{\pm}_{{\bf p}}}]^{-1}$ is the
quasiparticle residue of the polaron.

##  IV Magnetic dressing cloud

Having derived the SCBA expression for the polaron wave function in a bilayer,
we now turn to its spatial structure. Specifically, we will compute the
magnetic dressing cloud, i.e. the magnetic frustration around the hole.
Following Refs. [37, 21], we define the magnetization at position ${\bf
r}+{\bf d}$ in layer $l$ given the hole is at position ${\bf r}$ in layer $1$
as

$\displaystyle M_{l}^{\pm}({\bf d},{\bf p})$
$\displaystyle=\frac{\braket{\hat{h}^{\dagger}_{1,{\bf r}}\hat{h}_{1,{\bf
r}}\hat{S}^{(z)}_{l,{\bf r}+{\bf d}}}_{{\bf
p}}^{\pm}}{\braket{\hat{h}^{\dagger}_{1,{\bf r}}\hat{h}_{1,{\bf
r}}}^{\pm}_{{\bf p}}\\!\braket{\hat{S}^{(z)}_{l,{\bf r}+{\bf d}}}_{{\bf
p}}^{\pm}}=\frac{1-2M_{l}^{(2),\pm}({\bf d};{\bf p},\epsilon^{\pm}_{{\bf
p}})}{M_{{\rm AF}}},$ (24)

where $\braket{\dots}_{{\bf p}}^{\pm}=\bra{\Psi^{\pm}_{{\bf
p}}}\dots\ket{\Psi^{\pm}_{{\bf p}}}$ is the average with respect to the
symmetric or antisymmetric polaron wave function given by Eq. (21). The hole-
spin correlator is

$\displaystyle M_{l}^{(2),\pm}({\bf d};{\bf p},\epsilon^{\pm}_{{\bf p}})$
$\displaystyle=2N\braket{\hat{h}^{\dagger}_{1,{\bf r}}\hat{h}_{1,{\bf
r}}\hat{s}^{\dagger}_{l,{\bf r}+{\bf d}}\hat{s}_{l,{\bf r}+{\bf d}}}_{{\bf
p}}^{\pm}$ $\displaystyle=2N\braket{\hat{h}^{\dagger}_{l,{\bf
r}}\hat{h}_{l,{\bf r}}\hat{s}^{\dagger}_{1,{\bf r}+{\bf d}}\hat{s}_{1,{\bf
r}+{\bf d}}}_{{\bf p}}^{\pm}$ (25)

where we have used that the hole is evenly spread out over the entire lattice
for crystal momentum states, $\braket{\hat{h}^{\dagger}_{1,{\bf
r}}\hat{h}_{1,{\bf r}}}_{{\bf p},\pm}=1/(2N)$. This also means that the hole
has negligible impact on the average local magnetization,
$\braket{\hat{S}^{z}_{m,{\bf r}+{\bf d}}}_{{\bf p}}^{\pm}=\bra{{\rm
AF}}\hat{S}^{z}\ket{{\rm AF}}+\mathcal{O}(1/N)$. Spin wave theory gives
$2|\bra{{\rm AF}}\hat{S}^{z}\ket{{\rm AF}}|\equiv M_{{\rm AF}}=1-2M_{\rm fl}$,
where $M_{\rm
fl}=\sum_{\mathbf{k}}\left(v_{\mathbf{k},+}^{2}+v_{\mathbf{k},-}^{2}\right)/(2N)$
describes the reduction of magnetic order due to quantum fluctuations in the
antiferromagnetic ground state. Additionally, we have used the layer symmetry
of the system to swap the indices of the holes and spin excitations. Using Eq.
(5), we can write

$\displaystyle M_{l}^{(2),\pm}({\bf d};{\bf p},\epsilon^{\pm}_{{\bf
p}})=M_{\text{fl}}+B^{\pm}_{l,{\bf p}}({\bf d})+C^{\pm}_{l,{\bf p}}({\bf d}).$
(26)

where the $B$ and $C$ functions describe the frustration induced by the hole.
Their explicit expressions read

$\displaystyle B^{\pm}_{l,{\bf p}}({\bf d})$
$\displaystyle=\frac{1}{N}\sum_{{\bf q}_{1},{\bf q}_{2}}e^{i({\bf q}_{1}-{\bf
q}_{2})\cdot{\bf d}}B^{\pm}_{l}({\bf q}_{1},{\bf q}_{2};{\bf
p},\epsilon^{\pm}_{{\bf p}}),$ $\displaystyle B^{\pm}_{l}({\bf q}_{1},{\bf
q}_{2};{\bf p},\omega)$ $\displaystyle=\sum_{{\bf
k},\mu_{1},\mu_{2}}b_{\mu_{1},\mu_{2}}({\bf q}_{1},{\bf
q}_{2})\braket{\hat{h}^{\dagger}_{l,{\bf k}+{\bf q}_{1}-{\bf
q}_{2}}\hat{h}_{l,{\bf k}}\hat{b}^{\dagger}_{\mu_{1},-{\bf
q}_{1}}\hat{b}_{\mu_{2},-{\bf q}_{2}}}_{{\bf p},\omega}^{\pm},$ (27)

and

$\displaystyle C^{\pm}_{l,{\bf p}}({\bf d})$
$\displaystyle=-\frac{1}{2N}\sum_{{\bf q}_{1},{\bf q}_{2}}e^{i({\bf
q}_{1}+{\bf q}_{2})\cdot{\bf d}}\cdot C^{\pm,l}({\bf q}_{1},{\bf q}_{2};{\bf
p},\epsilon^{\pm}_{{\bf p}})+\text{H.c.},$ $\displaystyle C^{\pm}_{l}({\bf
q}_{1},{\bf q}_{2};{\bf p},\omega)$ $\displaystyle=\sum_{{\bf
k},\mu_{1},\mu_{2}}c_{\mu_{1},\mu_{2}}({\bf q}_{1},{\bf
q}_{2})\braket{\hat{h}^{\dagger}_{l,{\bf k}+{\bf q}_{1}-{\bf
q}_{2}}\hat{h}_{l,{\bf k}}\hat{b}_{\mu_{1},-{\bf q}_{1}}\hat{b}_{\mu_{2},-{\bf
q}_{2}}}_{{\bf p},\omega}^{\pm}.$ (28)

Here $\langle\ldots\rangle_{{\bf p},\omega}^{\pm}$ means the average with
respect to the polaron wave function $|\Psi_{\mathbf{p}}^{\pm}\rangle$ given
by Eq. (21) but with the energy $\epsilon^{\pm}_{{\bf p}}$ replaced by
$\omega$ when calculating the coefficients. The vertex functions
$b_{\mu_{1},\mu_{2}}({\bf q}_{1},{\bf q}_{2})=\mu_{1}\mu_{2}(v_{\mu_{1},{\bf
q}_{1}}v_{\mu_{2},{\bf q}_{2}}+u_{\mu_{1},{\bf q}_{1}}u_{\mu_{2},{\bf
q}_{2}})$ and $c_{\mu_{1},\mu_{2}}({\bf q}_{1},{\bf
q}_{2})=\mu_{1}\mu_{2}(u_{\mu_{1},{\bf q}_{1}}v_{\mu_{2},{\bf
q}_{2}}+v_{\mu_{1},{\bf q}_{1}}u_{\mu_{2},{\bf q}_{2}})$ arise due to the
transformation to the spin wave operators $\hat{b}$. We have $M_{l}^{\pm}=1$
for $B^{\pm}_{l,{\bf p}}=C^{\pm}_{l,{\bf p}}=0$, meaning that the
magnetization is that of the underlying antiferromagnet in the absence of
holes. It follows that $M_{l}^{\pm}<1$ reflects a hole-induced suppression of
AFM order while $M_{l}^{\pm}>1$ gives an increased order.

Calculating the magnetization around the hole is now reduced to computing the
$B^{\pm}_{l,{\bf p}}$ and $C^{\pm}_{l,{\bf p}}$ functions in Eqs. (27) and
(28). For strong coupling $t/J>1$, it is important to do this including terms
in the polaron wave function $|\Psi_{\mathbf{p}}^{\pm}\rangle$ with an
arbitrary number of spin waves. This was first achieved recently for a single
layer in Ref. [21], and here we generalize this to the bilayer case at hand.
As detailed in Appendixes D and E, using a diagrammatic approach we derive
self-consistency equations for $B^{\pm}_{l,{\bf p}}$ and $C^{\pm}_{l,{\bf
p}}$, which we then solve numerically. This allows us to include terms in the
wave function with an arbitrary number of spin waves, which is crucial in the
strong coupling regime.

As described in Sec. III, the bilayer supports magnetic polaron eigenstates
that are either symmetric or antisymmetric with respect to layer exchange. It
follows from the AFM symmetries in Eqs. (11) and (12) that the associated
magnetization clouds are related by the mapping (see Appendix G)

$\displaystyle M_{l}^{-}({\bf d},{\bf p}+{\bf Q})$
$\displaystyle=M_{l}^{+}({\bf d},{\bf p}).$ (29)

In the following, we, therefore, focus only on antisymmetric eigenstates.
Also, we can in analogy with the single-layer case use time-reversal and
inversion symmetry to show that the dressing cloud around the polaron is
$C_{2}$ symmetric, even when the momentum of the polaron is nonzero. Indeed,
we have $M_{l}^{\pm}({\bf d},{\bf p})=M_{l}^{\pm}({\bf d},-{\bf
p})=M_{l}^{\pm}(-{\bf d},{\bf p})$ where the first equality follows from time-
reversal symmetry and the second from inversion symmetry [21]. It follows from
these symmetries that all degenerate polaron states have the same dressing
cloud.

##  V Results

We now discuss our numerical results for the dressing cloud of magnetic
frustration around the hole. Compared to the single layer case discussed in
Ref. [21], the bilayer system has a qualitatively new feature in the sense
that it undergoes a quantum phase transition to a disordered state with
increasing interlayer coupling [30, 42, 43]. We, therefore, focus on how this
interlayer coupling affects the magnetic dressing cloud of the polaron. The
calculations are performed on a $16\times 16$ lattice with periodic boundary
conditions by first solving the self-consistent equations for the Green’s
functions, after which the $B^{\pm}_{l,{\bf p}}$ and $C^{\pm}_{l{\bf p}}$
series are calculated.

###  V.1 Dressing cloud for fixed momentum

Figure 4: The dressing cloud of magnetic frustration in layer 1 (left) and
layer 2 (right) around the hole in layer 1 for the polaron with crystal
momentum ${\bf p}=(\pi/2,\pi/2)$ indicated by the arrow. The cloud has a
$C_{2v}$ symmetry in both layers with the spins primarily frustrated parallel
to the crystal momentum. For $t_{\perp}/t=0.5,1.5$ the frustration cloud is
clearly more prominent in the same layer as the hole. Instead for
$t_{\perp}/t=3$, we see clouds of similar sizes in the two layers.

Figure 4 shows the magnetization $M_{l}^{-}({\bf d},{\bf p})$ at sites
separated by the in-layer vector ${\mathbf{d}}=(x,y)$ from the hole for the
antisymmetric polaron. We take momentum ${\mathbf{p}}=(\pi/2,\pi/2)$, which is
the ground state for $t_{\perp}/t\lesssim 0.6$ (see Fig. 3). The left column
shows the magnetization in layer $1$ where the hole is located, and the right
column shows the magnetization in layer $2$. Different rows correspond to
different values of $t_{\perp}/t$. The value $M_{l}^{-}({\bf d},{\bf p})=1$
corresponds to a magnetization equal to that of the AFM state without the
hole. Figure 4 clearly shows an elongated shape of frustration in both layers
and for all values of $t_{\perp}/t$ with the dressing cloud primarily spread
out parallel to the diagonal $x=y$. This is analogous to what was previously
shown to be the case for crystal momenta along the magnetic Brillouin zone
(MBZ) boundary $|p_{x}|+|p_{y}|=\pi$ for a single layer. We also see that, in
addition to the $C_{2}$ symmetry, the dressing cloud has mirror symmetries
along the diagonals. One can in fact show that for momenta along the MBZ
boundary, the magnetization cloud has the symmetries of the group $C_{2v}$
[21]. We, furthermore, see that with increasing $t_{\perp}$ the dressing cloud
steadily increases in spatial size and magnitude in layer $2$. Contrary to
this, the dressing cloud in layer $1$ where the hole resides is initially
slightly decreasing with increasing $t_{\perp}$. It eventually starts
increasing for larger values of $t_{\perp}$ and the dressing cloud has a
similar size in the two layers for $t_{\perp}/t=3$, where the hole has even
reversed the sign of the magnetic order at neighboring lattice sites.

Figure 5: The magnetization in layer (a) 1 and (b) layer 2 in the vicinity of
a hole in layer 1 as a function of interlayer hopping for ${\bf
p}=(\pi/2,\pi/2)$. Deviations from $M=1$ signals frustration. (a) In layer 1,
we see a nonmonotonic behavior as $t_{\perp}/t$ is increased. This can be
understood through the strength of the AFM order, which is shown as a black
line. As the order increases the cloud diminishes in size before expanding as
the order decreases. (b) In layer 2 on the other hand, the magnetic
frustration monotonically increases.

To analyze this behavior further, we plot in Fig. 5 the magnetization in the
vicinity of the hole as a function of $t_{\perp}/t$. We also plot the
magnetization $2|\bra{{\rm AF}}\hat{S}^{z}\ket{{\rm AF}}|$ in the absence of
holes, which exhibits a nonmonotonic behavior: It initially increases with
$t_{\perp}$ since the Fermi-Hubbard relation $J_{\perp}\propto t_{\perp}^{2}$
makes the AFM ordering across the planes stronger, reaching a maximum at
$t_{\perp}/t\simeq 0.82$ after which it decreases as the quantum phase
transition to the disordered state is approached [45]. This nonmonotonic
behavior of the magnetic order is reflected in the size and magnitude of the
dressing cloud in layer 1, which initially decreases with $t_{\perp}$, i.e.
$M$ increases towards $1$ as can be seen in Fig. 5, because the increased size
of the order parameter reduces the mobility of the hole. On the other hand,
when $t_{\perp}>0.82t$ and the local AF order starts decreasing, it becomes
easier for the hole to delocalize which boosts the magnetic frustration around
the hole (see Fig. 5). For small values of $t_{\perp}/t$, the second layer is
essentially not affected by the presence of a hole in layer 1, and the
magnetization is hardly changed from its background value, $M^{-}_{2}=1$.
However, as the coupling between the layers is increased, frustration appears
here as well. This happens both because the hole can now make virtual hops to
the second layer via $t_{\perp}$, and because the interlayer spin coupling
$J_{\perp}\propto t_{\perp}^{2}$ favors the spins in the second layer to
antialign with the spins in layer 1. Again we see from Fig. 5 that for large
values of $t_{\perp}$ the presence of the hole even inverts the sign of the
magnetic order. For values other than $J/t=0.3$, the situation is
qualitatively the same. Quantitatively, lower values of $J/t<0.3$ mean a
lowered frustration cost and result in an increase in frustration, and vice
versa for $J/t>0.3$. This is detailed in Appendix I.

###  V.2 Dressing cloud of ground state

As we discussed in Sec. II.2, the momentum of the ground state polaron
smoothly changes from ${\bf p}=(\pi/2,\pi/2)$ at $t_{\perp}=0$ to ${\bf
p}={\bf 0}$ [or ${\bf p}=(\pi,\pi)$] as the QCP is approached with increasing
$t_{\perp}$ (see Fig. 3 [45, 44]). Figure 6 plots the magnetization in the
vicinity of the hole as a function of $t/t_{\perp}$ for the ground state, i.e.
the crystal momentum varies according to Fig. 3. In addition to the increase
of the dressing cloud with increasing $t_{\perp}$ also seen for the case of
the fixed momentum discussed in Sec. V.1, Fig. 6 shows that the _symmetry_ of
the dressing cloud changes. The reason for this change is that the crystal
momentum of the ground state varies with $t_{\perp}$. For small $t_{\perp}/t$,
the momentum ${\bf p}=(\pi/2,\pi/2)$ is at the magnetic Brillouin zone
boundary so that the dressing cloud has $C_{2v}$ symmetry [21]. As seen in
Fig. 3, the momentum of the ground state decreases with $t_{\perp}$, becoming
zero for $t_{\perp}/t\gtrsim 2.3$. At this point, the dressing cloud recovers
the full $C_{4v}$ symmetry of the AFM order. It follows that the symmetry of
the magnetic polaron undergoes the transition $C_{2v}\rightarrow C_{4v}$ with
increasing $t_{\perp}/t$.

Figure 6: The magnetization cloud in layer 1 (left) and layer 2 (right) in the
vicinity of a hole in layer 1 for the ground state of the system. The ground
state momentum changes from (a) ${\bf p}=(\pi/2,\pi/2)$, over (b) ${\bf
p}=(\pi/4,\pi/4)$, to (c) ${\bf p}={\bf 0}$ as illustrated by the black
arrows.

Figure 6 also shows that the dressing cloud in layer 2 closely mimics that in
layer 1 when $t_{\perp}$ is large. To analyze this quantitatively, we plot in
Fig. 7 the difference $M_{1}^{-}({\bf r})-M_{2}^{-}({\bf r})$ between the
magnetization in the two layers for ${\bf p}=(\pi/2,\pi/2)$ and ${\bf p}={\bf
0}$ for a few selected sites as a function of $t_{\perp}/t$. This clearly
shows that the difference vanishes rapidly with increasing $t_{\perp}$ for the
polaron with momentum ${\bf p}={\bf 0}$ [Fig. 7(a)]. The reason is that the
spins in the dressing cloud in the two layers antialign for strong interlayer
coupling making them approximately mirror images of each other. Notably, this
antialignment of the spins inside the dressing clouds in the two layers is
much less pronounced for the polaron with momentum ${\bf p}=(\pi/2,\pi/2)$,
and this yields spatial insight into why the ground state momentum changes
from ${\bf p}=(\pi/2,\pi/2)$ to ${\bf p}={\bf 0}$ for increasing $t_{\perp}$.
In fact, for $t_{\perp}\gg t$, it comes with a high magnetic energy cost
scaling as $J_{\perp}\propto t_{\perp}^{2}$ if spins do not antialign between
the layers, explaining why the ${\bf p}=(\pi/2,\pi/2)$ polaron has higher
energy than the ${\bf p}={\bf 0}$ polaron in this regime. In Appendix H, we
elaborate on how the interlayer hopping and interlayer spin-spin individually
promote this mirroring effect of the ground state by taking $J_{\perp}\neq
4t^{2}_{\perp}/U$. Here, we see that the hopping primarily causes the
mirroring effect, while the spin-spin interaction penalizes states not having
AFM order between the layers. As a result, both couplings are needed to ensure
that the ground state features mirrored layers for increasing $t_{\perp}/t$.

Figure 7: Difference $M^{-}_{1}-M^{-}_{2}$ in magnetization between the layers
versus the interlayer hopping $t_{\perp}/t$. (a) For ${\bf p}={\bf 0}$, the
magnetization rapidly becomes identical in the two layers. (b) For ${\bf
p}=(\pi/2,\pi/2)$, they remain different for the shown values of $t_{\perp}/t$
but will eventually become identical (Appendix D).

Finally, we plot in Fig. 8 the size of the magnetic dressing cloud of the
ground state in the two layers defined as

$r_{c}^{l}({\bf p})=\sqrt{\sum_{{\bf d}}|{\bf d}|^{2}\cdot|M^{-}_{l}({\bf
d},{\bf p})-1|}$ (30)

as a function of $t_{\perp}/t$. We subtract $1$ from the frustration since
$M^{-}_{l}({\bf d},{\bf p})-1=0$ corresponds to no frustration. This
quantifies how the dressing cloud reflects the size of the order parameter of
the environment and, in particular, how it grows as the phase transition to
the disordered state, occurring at $t_{\perp}/t\simeq 3.67$ within LSWT, is
approached. We believe this increase is a reliable result even though spin-
wave theory cannot describe the region close to the QCP in a quantitatively
accurate way. The reason is that the increase is caused by the increase in the
hole mobility, which in turn is due to the decrease in the AF order as the QCP
is approached. This is a robust physical effect and not an artifact of the
theoretical approximations. In Refs. [47, 48] they argue for the same
expansion of the dressing cloud when the QCP is approached from the disordered
phase, supporting our results. The same effect is shown in Fig. 1, which
clearly shows how the magnetic dressing cloud expands as $t_{\perp}$ increases
and the phase transition to the disordered state is approached. One can also
see how the dressing clouds in the two layers become mirror images of each
other with their spins antiparallel for large $t_{\perp}/t$.

Figure 8: The size of the ground state dressing cloud in the two layers with
the hole residing in layer $1$ as a function of $t_{\perp}/t$. The size of the
magnetic order parameter in the absence of the hole is shown as a black line.

##  VI Conclusions

Inspired by the impressive ability to experimentally probe the spatial
properties of strongly correlated fermions with single-site resolution, not
only in single-layer but also bilayer optical lattices [27], we developed a
nonperturbative scheme for calculating the polaron wave function for a bilayer
antiferromagnet. The scheme is the wave function version of the SCBA [41, 37,
21] and allows one to include an infinite number of spin waves, which is
crucial for describing the strongly correlated regime. With the developed wave
function formalism, it is possible to calculate the spatial structure of the
bilayer magnetic polarons via spin-hole correlators. In particular, we mapped
out the associated dressing cloud in the vicinity of a hole and investigated
how this was influenced by the interlayer coupling strength. Here, we observed
the cloud to mimic the nonmonotonic behavior of the antiferromagnetic order,
where it first contracts due to an increase in the AFM order, whereafter it
starts to expand as the quantum phase transition to the disordered phase is
approached. We, furthermore, found that the symmetry of the ground state
dressing cloud undergoes the transition $C_{2v}\rightarrow C_{4v}$ with
increasing $t_{\perp}/t$. This happens because the ground state momentum
smoothly changes along the diagonal from ${\bf p}=(\pi/2,\pi/2)$ for small
$t_{\perp}/t$ to ${\bf p}=\bf{0}$ for large $t_{\perp}/t$. In the latter
limit, we find that the strong interlayer coupling makes the spins in the two
layers antiparallel in the dressing cloud so that they become mirror images of
each other. These effects should be amenable to direct observation in quantum
simulation experiments. Here, however, it is important to note that one must
project onto a single crystal momentum eigenstate to observe something else
than the full $C_{4v}$ symmetry of the underlying system.

Our scheme for nonperturbatively calculating the polaron wave function in a
bilayer geometry within SCBA offers a tool that can help to elucidate the
behavior of dopants in antiferromagnetic environments. For a single layer, the
wave function has been used to investigate the nonequilibrium physics
associated with the creation of the magnetic polaron [22]. The similarities
between the expressions for the single-layer and bilayer system suggest that
it is possible to probe the same dynamics for a bilayer system with the here-
derived wave function. It also offers the possibility of a detailed comparison
with novel ultracold-atom experiments [27].

In the future, we hope to explore extensions to finite temperatures, as well
as comparisons with sophisticated numerical methods, using, e.g., density
matrix renormalization techniques [49, 50]. This will allow a more
quantitative description as the quantum phase transition from the AFM ordered
state at low interlayer coupling to the disordered phase at large interlayer
coupling is crossed.

###### Acknowledgements.

This work has been supported by the Danish National Research Foundation
through the Center of Excellence “CCQ” (Grant Agreement No.: DNRF156), as well
as by the Carlsberg Foundation through a Carlsberg Internationalisation
Fellowship.

##  Appendix A Recursion relation

In this appendix, we will explicitly show that the recursion relation
presented in Eq. (19) self-consistently solves Eq. (18). Inserting the
recursion relation into Eq. (18), we find

$\displaystyle\left[\left(\omega-\sum_{i=1}^{n}\omega_{\mu_{i},{\bf
k}_{i}}\right)1\kern-2.5pt\text{l}-\Sigma\big{(}{\bf
K}_{n},\omega-\sum_{i=1}^{n}\omega_{\mu_{i},{\bf k}_{i}}\big{)}\right]{\bf
a}_{\\{\mu_{n}\\}}({\bf p},\\{{\bf k}_{n}\\})$
$\displaystyle=V^{\dagger}_{\mu_{n}}({\bf K}_{n-1},{\bf k}_{n}){\bf
a}_{\\{\mu_{n-1}\\}}({\bf p},\\{{\bf k}_{n-1}\\}),$ (31)

where Eq. (15) has been used. Looking at the Green’s functions in Eq. (14), we
find

$\displaystyle G({\bf p},\omega)^{-1}=$
$\displaystyle\begin{bmatrix}\omega-\Sigma_{\rm d}({\bf
p},\omega)&-\Sigma_{\rm o}({\bf p},\omega)\\\ -\Sigma_{\rm o}({\bf
p},\omega)&\omega-\Sigma_{\rm d}({\bf p},\omega)\end{bmatrix},$ (32)

which corresponds to the parentheses on the left-hand side in Eq. (31).
Inserting this, we retrieve

$\displaystyle{\bf a}_{\\{\mu_{n}\\}}({\bf p},\\{{\bf k}_{n}\\})=G\big{(}{\bf
K}_{n},\omega-\sum_{i=1}^{n}\omega_{\mu_{i},{\bf k}_{i}}\big{)}$
$\displaystyle\times V^{\dagger}_{\mu_{n}}({\bf K}_{n-1},{\bf k}_{n}){\bf
a}_{\\{\mu_{n-1}\\}}({\bf p},\\{{\bf k}_{n-1}\\}),$ (33)

showing the recursion relation to be self-consistent.

##  Appendix B Normalization

In this Appendix, we show how the normalization of the wave function links the
lowest-order amplitude to the quasiparticle residue. We have,

$\displaystyle 1=\ $ $\displaystyle\langle\Psi^{\pm}_{{\bf
p}}|\Psi^{\pm}_{{\bf p}}\rangle=2\left(a^{\pm}_{0}({\bf
p})\right)^{2}+2^{2}\sum_{{\bf k}_{1}}\left(a^{\pm}_{1}({\bf p},{\bf
k}_{1})\right)^{2}+$ $\displaystyle 2^{3}\sum_{{\bf k}_{1},{\bf
k}_{2}}\left(a^{\pm}_{2}({\bf p},\\{{\bf k}_{2}\\})\right)^{2}+...,$ (34)

where the factor $2^{n+1}$ in front comes from the number of different states
possible to each order. Each amplitude is related to the previous order
through Eq. (47). Inserting this yields

$\displaystyle\langle\Psi^{\pm}_{{\bf p}}|\Psi^{\pm}_{{\bf
p}}\rangle=2\left(a^{\pm}_{0}({\bf p})\right)^{2}\left(1+2\sum_{{\bf
k}_{1}}\left[R^{\pm}({\bf p},\\{{\bf k}_{1}\\})\right]^{2}+\right.$
$\displaystyle\left.2^{2}\sum_{{\bf k}_{1},{\bf k}_{2}}\left[R^{\pm}({\bf
p},\\{{\bf k}_{2}\\})\right]^{2}\left[R^{\pm}({\bf p},\\{{\bf
k}_{1}\\})\right]^{2}+...\right).$ (35)

Using the expressions for the Green’s functions in Eq. (14) and the self-
energies in Eq. (15), we obtain

$\displaystyle 1-$ $\displaystyle\partial_{\omega}\left[\Sigma_{{\rm d}}({\bf
p},\omega)\pm\Sigma_{{\rm o}}({\bf p},\omega)\right]=1+2\sum_{{\bf
k}_{1}}\left[R^{\pm}({\bf p},\\{{\bf k}_{1}\\})\right]^{2}$
$\displaystyle\left(1-\partial_{\omega}\left[\Sigma_{{\rm d}}({\bf
k}_{1},\omega-\omega_{+,{\bf k}})\pm\Sigma_{{\rm o}}({\bf
k}_{1},\omega-\omega_{+,{\bf k}})\right]\right).$ (36)

Evaluating this at $\omega=\varepsilon^{\pm}_{\bf p}$ and comparing it to Eq.
(35), we see that they are identical up to the factor of $2(a^{\pm}_{0}({\bf
p}))^{2}$. Therefore,

$\displaystyle a^{\pm}_{0}({\bf p})$
$\displaystyle=\frac{1}{\sqrt{2}}\frac{1}{\sqrt{1-\partial_{\omega}\left(\Sigma_{{\rm
d}}({\bf p},\omega)\pm\Sigma_{{\rm o}}({\bf
p},\omega)\right)|_{\omega=\epsilon^{\pm}_{{\bf p}}}}},$ (37)

Finally, we relate this to the quasiparticle residue. To make the residue
coincide with that for a single layer in the limit of $t_{\perp}/t=0$, we
define this as half the residue of the poles in the Green’s function whereby

$\displaystyle Z^{\pm}_{\bf p}/2={\rm Res}\left(G_{{\rm
d}},\epsilon^{\pm}_{\bf p}\right)$ $\displaystyle=\left(a^{\pm}_{0}({\bf
p})\right)^{2}$ $\displaystyle Z^{\pm}_{\bf p}/2={\rm Res}\left(G_{{\rm
o}},\epsilon^{\pm}_{\bf p}\right)$ $\displaystyle=\pm\left(a^{\pm}_{0}({\bf
p})\right)^{2}.$ (38)

Consequently, the quasiparticle residue, hereby, describes the total
probability of finding the magnetic polaron as a bare hole

$\displaystyle|\bra{\Psi^{\pm}_{{\bf p}}}\hat{h}^{\dagger}_{1,{\bf
p}}\ket{{\rm AF}}|^{2}+|\bra{\Psi^{\pm}_{{\bf p}}}\hat{h}^{\dagger}_{2,{\bf
p}}\ket{{\rm AF}}|^{2}=Z^{\pm}_{{\bf p}}.$ (39)

Ultimately, this leaves the normalization condition to be

$\displaystyle a^{\pm}_{0}({\bf p})=\sqrt{\frac{Z^{\pm}_{{\bf p}}}{2}},$ (40)

##  Appendix C Reducing coefficients

In this appendix, we will see that it is possible to define a single unique
coefficient for each order using the recursion relation in Eq. (19), the
symmetries of the interaction vertices, Eq. (12), and those of the Green’s
functions

$\displaystyle G_{\rm{d}}(\mathbf{p}+\mathbf{Q},\omega)$
$\displaystyle=G_{\rm{d}}(\mathbf{p},\omega)$ $\displaystyle
G_{\rm{o}}(\mathbf{p}+\mathbf{Q},\omega)$
$\displaystyle=-G_{\rm{o}}(\mathbf{p},\omega).$ (41)

The recursion relations can be restated as

$\displaystyle\begin{bmatrix}a^{(1)}_{\\{\mu_{n}\\}}({\bf p},\\{{\bf
k}_{n+1}\\})\\\ a^{(2)}_{\\{\mu_{n}\\}}({\bf p},\\{{\bf
k}_{n+1}\\})\end{bmatrix}=\
\prod^{n-1}_{i=0}\mathcal{V}_{\\{\mu_{n-i}\\}}(\\{{\bf k}_{n-i}\\};{\bf
p},\epsilon^{\pm}_{{\bf p}})\cdot\begin{bmatrix}a^{\pm}_{0}({\bf p})\\\ \pm
a^{\pm}_{0}({\bf p})\end{bmatrix}$ (42)

with

$\displaystyle\mathcal{V}_{\\{\mu_{n}\\}}(\\{{\bf k}_{n}\\};{\bf
p},\omega)=G({\bf K}_{n},\Omega^{\\{{\bf
k}_{n}\\}}_{\\{\mu_{n}\\}})V^{\dagger}_{\mu_{n}}({\bf K}_{n-1},{\bf k}_{n}),$
(43)

and $a^{(1)}({\bf p})\equiv a^{\pm}_{0}({\bf p})$ being the zeroth order
coefficient. To relate the different coefficients we effectively need to
change all subscripts of $\mathcal{V}$ to $\mu=+$. If translating the last
excited spin waves crystal momentum by ${\bf Q}$, we find

$\mathcal{V}_{\\{\mu_{n-1}\\},-}(\\{{\bf k}_{n-1}\\},{\bf k}_{n}+{\bf Q};{\bf
p},\omega)=\\\ \begin{bmatrix}1&0\\\
0&-1\end{bmatrix}\cdot\mathcal{V}_{\\{\mu_{n-1}\\},+}(\\{{\bf k}_{n-1}\\},{\bf
k}_{n};{\bf p},\omega).$ (44)

Hence, we see that by translating the momentum by ${\bf Q}$ it is possible to
change the subscripts at the expense of a possible sign. This relation also
explains the signs in the expression for the wave function in Eq. (21). If
instead, we wish to change not the last subscript but any other, then we find

$\mathcal{V}_{\\{\mu_{n}\\},-,\\{\mu_{m}\\}}(\\{{\bf k}_{n}\\},{\bf
k}_{n+1}+{\bf Q},\\{{\bf k}_{m}\\};{\bf p},\omega)=\\\ \begin{bmatrix}1&0\\\
0&-1\end{bmatrix}\cdot\mathcal{V}_{\\{\mu_{n}\\},+,\\{\mu_{m}\\}}(\\{{\bf
k}_{n}\\},{\bf k}_{n+1},\\{{\bf k}_{m}\\};{\bf
p},\omega)\cdot\begin{bmatrix}-1&0\\\ 0&1\end{bmatrix}.$ (45)

If defining the coefficient with the hole present in layer 1 with excited spin
waves only in the $+$ branch as

$\displaystyle a^{(1)}_{+,+,...,+}({\bf p},\\{{\bf k}_{n}\\})$
$\displaystyle\equiv a^{\pm}_{n}({\bf p},\\{{\bf k}_{n}\\}),$ (46)

then it becomes possible with Eq. (49) and (44) to relate any $n$th order
coefficient to $a^{\pm}_{n}$ by translating all momenta relating to a spin
wave of the type $\mu=-$ by ${\bf k}\rightarrow{\bf k}+{\bf Q}$. They will be
related up to a sign which depends on the layer the hole is in, the types of
excited spin waves, and the order in which they are excited. As an example we
find $a^{(1)}_{+,-}({\bf k}_{1},{\bf k}_{2}+{\bf Q};{\bf
p},\epsilon^{\pm}_{{\bf p}})=-a^{(1)}_{-,+}({\bf k}_{1}+{\bf Q},{\bf
k}_{2};{\bf p},\epsilon^{\pm}_{{\bf p}})=-a^{\pm}_{2}({\bf k}_{1},{\bf
k}_{2};{\bf p},\epsilon^{\pm}_{{\bf p}})$. This also leaves it possible to
write the recursion relation as the scalar relation

$\displaystyle a^{\pm}_{n+1}({\bf p},\\{{\bf k}_{n+1}\\})=$ $\displaystyle\
R^{\pm}({\bf p},\\{{\bf k}_{n+1}\\})a^{\pm}_{n}({\bf p},\\{{\bf k}_{n}\\}),$
(47)

with the recursive coupling

$\displaystyle R^{\pm}({\bf p},\\{{\bf k}_{n}\\})=$ $\displaystyle\left(G_{\rm
d}({\bf K}_{n},\Omega^{\pm}_{\\{{\bf k}_{n}\\}})\pm G_{\rm o}({\bf
K}_{n},\Omega^{\pm}_{\\{{\bf k}_{n}\\}})\right)$
$\displaystyle\left(g_{+}\left({\bf K}_{n-1},{\bf k}_{n}\right)\pm
f_{+}\left({\bf k}_{n}\right)\right)$ $\displaystyle\Omega^{\pm}_{\\{{\bf
k}_{n}\\}}({\bf p})=$ $\displaystyle\ \epsilon^{\pm}_{\bf
p}-\sum_{i=1}^{n}\omega_{+,{\bf k}_{n}}.$ (48)

From the above relations, we will now derive relations between amplitudes in
the wave function that will be necessary for the derivations of the B- and
C-series. If changing the order of two different types of spin waves we get

$a^{(l)}_{\\{\mu_{n}\\},+,-,\\{\mu_{m}\\}}(\\{{\bf k}_{n}\\},{\bf k},{\bf
q}+{\bf Q},\\{{\bf k}_{m}\\};{\bf p},\omega)=\\\
-a^{(l)}_{\\{\mu_{n}\\},-,+,\\{\mu_{m}\\}}(\\{{\bf k}_{n}\\},{\bf k}+{\bf
Q},{\bf q},\\{{\bf k}_{m}\\};{\bf p},\omega).$ (49)

If we instead wish to change only one subscript from $-$ to $+$, we find

$a^{(l)}_{\\{\mu_{n}\\},-,\\{\mu_{m}\\}}(\\{{\bf k}_{n}\\},{\bf k}+{\bf
Q},\\{{\bf k}_{m}\\};{\bf p},\omega)=\\\
(-1)^{m+l}a^{(l)}_{\\{\mu_{n}\\},+,\\{\mu_{m}\\}}(\\{{\bf k}_{n}\\},{\bf
k},\\{{\bf k}_{m}\\};{\bf p},\omega).$ (50)

At last, we investigate what happens if changing the first or two first types
of spin waves. If changing just the first

$a^{l}_{-,\\{\mu_{n-1}\\}}({\bf k}_{1},\\{ k_{n-1}\\};{\bf p},\omega)=\\\
(-1)^{n+l}a^{l}_{+,\\{\mu_{n-1}\\}}({\bf k}_{1}+{\bf Q},\\{ k_{n-1}\\};{\bf
p},\omega).$ (51)

For the first two

$a^{(l)}_{-,-,\\{\mu_{n-2}\\}}({\bf k}_{1},{\bf k}_{2},\\{ k_{n-2}\\};{\bf
p},\omega)=\\\ -a^{(l)}_{+,+,\\{\mu_{n-2}\\}}({\bf k}_{1}+{\bf Q},{\bf
k}_{2}+{\bf Q},\\{ k_{n-2}\\};{\bf p},\omega)$ (52)
$a^{(l)}_{+,-,\\{\mu_{n-2}\\}}({\bf k}_{1},{\bf k}_{2},\\{ k_{n-2}\\};{\bf
p},\omega)=\\\ (-1)^{n+l}a^{(l)}_{+,+,\\{\mu_{n-2}\\}}({\bf k}_{1},{\bf
k}_{2}+{\bf Q},\\{ k_{n-2}\\};{\bf p},\omega)$ (53)
$a^{(l)}_{-,+,\\{\mu_{n-2}\\}}({\bf k}_{1},{\bf k}_{2},\\{ k_{n-2}\\};{\bf
p},\omega)=\\\ -(-1)^{n+l}a^{(l)}_{+,+,\\{\mu_{n-2}\\}}({\bf k}_{1}+{\bf
Q},{\bf k}_{2},\\{ k_{n-2}\\};{\bf p},\omega).$ (54)

##  Appendix D B-series

To understand the structure of this series, we will start by writing the wave
function in a diagrammatic form. This form will not explicitly state the full
complexity of the wave function and extra rules are needed to translate it
into the different terms. This is though to work as a guide in understanding
the structure The double line combined with the vertex translates to
$R^{\pm}({\bf p};{\bf k})$ while the straight line to $\sqrt{Z^{\pm}_{{\bf
p}}}$. Looking at Eq. (21) we see that this is not the whole story since
exciting $\mu=-$ spin wave leaves $R$ to be evaluated at $R({\bf p},{\bf
k}+{\bf Q})$ and then we have a minus sign for specific vertices. We will come
back to these features later, but for now, let us focus on the structure of
the $B$ series. Looking at Eq. (27) we see that the diagrammatic structure of
$B^{\pm,l}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)$ becomes where
$B^{\pm,l}=Z^{\pm}_{{\bf p}}\sum_{i=1}^{\infty}B^{\pm,l}_{i}$, and the orange
vertex describes the $b_{\mu_{1},\mu_{2}}$ factor. The factor $Z^{\pm}_{{\bf
p}}$ comes from omitting the lines at the end in the diagrams. There will also
be diagrams that are not symmetric around the center, $B^{asym}$, it is though
seen in Appendix F that $B^{asym}=0$. To handle these diagrams we collect all
the diagrams of the last type in $B^{\pm,l}_{0}$ such that

(57)

With these definitions, we see that it is possible to rewrite the $B$ series
as

(58)

With these relations, we will start translating the diagrammatic
representation. We will begin with $F^{\pm}_{B}$ where we have to remember,
that the missing vertex in the center means that the momentum in the two
propagators need not be the same. We also need to worry about the relation in
Eq. (51). This says, that it actually matters whether or not the diagrams
attached to $F_{B}$ contain two spin-wave of the same kind or opposite. If
they are opposite then the diagrams need to come with alternating signs
according to Eq. (51). For now, we will disregard the factor of $(-1)^{l}$ to
find the following for the same kinds of spin waves

$\displaystyle F^{\pm,s}_{B}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=\sum_{{\bf k}}\big{(}R^{\pm}({\bf q}_{1},{\bf k};{\bf
p},\omega)R^{\pm}({\bf q}_{2},{\bf k};{\bf p},\omega)\times$ $\displaystyle
\left[1+F^{\pm,s}_{B}({\bf q}_{1},{\bf q}_{2};{\bf p}+{\bf
k},\omega-\omega^{+}_{{\bf k}})\right]+R^{\pm}({\bf q}_{1},{\bf k}+{\bf
Q};{\bf p},\omega)\times$ $\displaystyle R^{\pm}({\bf q}_{2},{\bf k}+{\bf
Q};{\bf p},\omega)\left[1+F^{\pm,s}_{B}({\bf q}_{1},{\bf q}_{2};{\bf p}+{\bf
k},\omega-\omega^{-}_{{\bf k}})\right]\big{)},$ (59)

and for two different

$\displaystyle F^{\pm,a}_{B}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=-\sum_{{\bf k}}\big{(}R^{\pm}({\bf q}_{1},{\bf k};{\bf
p},\omega)\times$ $\displaystyle R^{\pm}({\bf q}_{2},{\bf k};{\bf
p},\omega)\left[1+F^{\pm,a}_{B}({\bf q}_{1},{\bf q}_{2};{\bf p}+{\bf
k},\omega-\omega^{+}_{{\bf k}})\right]+$ $\displaystyle R^{\pm}({\bf
q}_{1},{\bf k}+{\bf Q};{\bf p},\omega)R^{\pm}({\bf q}_{2},{\bf k}+{\bf Q};{\bf
p},\omega)\times$ $\displaystyle\left[1+F^{\pm,a}_{B}({\bf q}_{1},{\bf
q}_{2};{\bf p}+{\bf k},\omega-\omega^{-}_{{\bf k}})\right]\big{)},$ (60)

with the recursive coupling now being defined as

$\displaystyle R^{\pm}(\\{{\bf k}_{n}\\};{\bf p},\omega)=$
$\displaystyle\left(G_{\rm d}({\bf K}_{n},\Omega^{\pm}_{\\{{\bf k}_{n}\\}})\pm
G_{\rm o}({\bf K}_{n},\Omega^{\pm}_{\\{{\bf k}_{n}\\}})\right)$
$\displaystyle\left(g_{+}\left({\bf K}_{n-1},{\bf k}_{n}\right)\pm
f_{+}\left({\bf k}_{n}\right)\right)$ $\displaystyle\Omega^{\pm}_{\\{{\bf
k}_{n}\\}}(\omega)=$ $\displaystyle\ \omega-\sum_{i=1}^{n}\omega_{+,{\bf
k}_{n}},$ (61)

in order to encompass being evaluated away from the quasi-particle energies.
As a sanity check, we find that in the limit of $t_{\perp}=0$ the results are
equivalent to those obtained for a single layer [21].

With this, we can now approach $B^{\pm,l}_{0}$, where we need to remember that
the two boson lines in the first term in Eq. (57) need not be the same and the
vertex $b_{\mu_{1},\mu_{2}}({\bf q}_{1},{\bf q}_{2})$ depends on the types of
spin waves. The fact that this diagram also allows for two different kinds of
spin waves makes it important to think about which layer the hole is present
in. Looking at Eq. (17), we see that if the hole is present in layer 1, then
the two types of spin waves come with the same sign while for layer 2 this is
not the case. This means that $B^{\pm,l}_{1}$ comes with a $-(-1)^{l}$ in
front when the two spin waves are of the opposite kind. Having this in mind we
find

$\displaystyle B^{\pm,l}_{0}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=\frac{Z^{\pm}_{{\bf p}}}{2}\big{[}R^{\pm}({\bf q}_{1};{\bf
p},\omega)R^{\pm}({\bf q}_{2};{\bf p},\omega)\times$ $\displaystyle
b_{+,+}({\bf q}_{1},{\bf q}_{2})\cdot\big{[}1+F^{\pm,s}_{B}({\bf q}_{1},{\bf
q}_{2};{\bf p},\omega)\big{]}+$ $\displaystyle R^{\pm}({\bf q}_{1}+{\bf
Q};{\bf p},\omega)R^{\pm}({\bf q}_{2}+{\bf Q};{\bf p},\omega)b_{-,-}({\bf
q}_{1},{\bf q}_{2})\times$ $\displaystyle \big{[}1+F^{\pm,s}_{B}({\bf
q}_{1}+{\bf Q},{\bf q}_{2}+{\bf Q};{\bf p},\omega)\big{]}-(-1)^{l}R^{\pm}({\bf
q}_{1};{\bf p},\omega)\times$ $\displaystyle R^{\pm}({\bf q}_{2}+{\bf Q};{\bf
p},\omega)b_{+,-}({\bf q}_{1},{\bf q}_{2})\cdot\big{[}1+F^{\pm,a}_{B}({\bf
q}_{1},{\bf q}_{2}+{\bf Q};{\bf p},\omega)\big{]}$
$\displaystyle-(-1)^{l}R^{\pm}({\bf q}_{1}+{\bf Q};{\bf p},\omega)R^{\pm}({\bf
q}_{2};{\bf p},\omega)b_{-,+}({\bf q}_{1},{\bf q}_{2})\big{]}\times$
$\displaystyle\big{[}1+F^{\pm,a}_{B}({\bf q}_{1}+{\bf Q},{\bf q}_{2};{\bf
p},\omega)\big{]}$ (62)

Here, we find $B^{\pm,2}=0$ in the limit of $t_{\perp}/t=0$ and $B^{\pm,1}$ to
coincide with the result for a single layer. This is expected since in Eq.
(25) we choose the spin to be measured in layer 1, it will therefore make
sense that if choosing $l=1$ and taking the limit of $t_{\perp}=0$ we expect
to find $M^{1,2}=0$, $B^{\pm,1}=0$ and $C^{\pm,1}=0$. Another feature found in
Eq. (62) is that the only difference between the two layers comes from the
terms mixing the spin wave branches. If we take the extreme limit of
$t_{\perp}\gg t$ we have $b_{+,-}=0=b_{-,+}$. The same holds for the
$C$-series, such that all states, as expected, will show perfect antialignment
between the layers in this limit. It is not trivial what happens for
intermediate values [see Fig. 7].

With these relations, we can now present the translation of Eq. (58)

$\displaystyle B^{\pm,l}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=B^{\pm,l}_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)+$
$\displaystyle\sum_{{\bf k}}\big{[}R^{\pm}({\bf k};{\bf p},\omega)R^{\pm}({\bf
k};{\bf p},\omega)\cdot B^{\pm,l}({\bf q}_{1},{\bf q}_{2};{\bf p}+{\bf
k},\omega-\omega^{+}_{{\bf k}})+$ $\displaystyle R^{\pm}({\bf k}+{\bf Q};{\bf
p},\omega)R^{\pm}({\bf k}+{\bf Q};{\bf p},\omega)\cdot B^{\pm,l}({\bf
q}_{1},{\bf q}_{2};{\bf p}+{\bf k},\omega-\omega^{-}_{{\bf k}})\big{]}$ (63)

We can now find $B^{\pm,l}$ by first solving Eq. (59) and (60) self-
consiestently starting from $F^{\pm}_{B}=0$. The results are then inserted
into Eq. (62) which is then used to solve Eq. (63) self-consistently starting
from $B^{\pm,l}=B^{\pm,l}_{0}$.

##  Appendix E C-series

Finding an expression for the $C$ series is much like finding the one for the
$B$ series. The $C$ series differs in that it will connect different orders of
the wave function while the $B$ series connects the same orders. Again, we
define $C=\sum_{n=0}^{\infty}C_{n,n+2}$, where $C_{i,j}$ represents the
overlap between the $i^{th}$ order in $\bra{\Psi^{\pm}_{{\bf p}}}$ and the
$j^{th}$ order in $\ket{\Psi^{\pm}_{{\bf p}}}$. $C_{n,n+2}$ can
diagrammatically be represented as We also have other diagrams, but we find
them summing to zero (See Appendix F). If collecting the last type of diagrams
in the quantity $C_{0}$, we see

(65)

The derivation of $F_{C}$ follows that of $F_{B}$ with the only difference
being, that in the expression for $C^{\pm,l}_{0}$ we have two spin waves
excited to the right of $F_{C}$ and none to the left, for $F_{B}$ we had one
on each side. We, therefore, need to use the relations presented in Eq. (54).
Again we need to define two types of $F_{C}$ where one has alternating signs

$\displaystyle F^{\pm,s}_{C}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=\sum_{{\bf k}}\big{[}R^{\pm}({\bf k};{\bf p},\omega)R^{\pm}({\bf
q}_{1},{\bf q}_{2},{\bf k};{\bf p},\omega)\times$
$\displaystyle\left[1+F^{\pm,s}_{C}({\bf q}_{1},{\bf q}_{2};{\bf p}+{\bf
k},\omega-\omega^{+}_{{\bf k}})\right]+$ $\displaystyle R^{\pm}({\bf k}+{\bf
Q};{\bf p},\omega)R^{\pm}({\bf q}_{1},{\bf q}_{2},{\bf k}+{\bf Q};{\bf
p},\omega)\big{[}1+$ $\displaystyle F^{\pm,s}_{C}({\bf q}_{1},{\bf q}_{2};{\bf
p}+{\bf k},\omega-\omega^{-}_{{\bf k}})\big{]}$ $\displaystyle
F^{\pm,a}_{C}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=-\sum_{{\bf
k}}\big{[}R^{\pm}({\bf k};{\bf p},\omega)R^{\pm}({\bf q}_{1},{\bf q}_{2},{\bf
k};{\bf p},\omega)\times$ $\displaystyle\left[1+F^{\pm,a}_{C}({\bf q}_{1},{\bf
q}_{2};{\bf p}+{\bf k},\omega-\omega^{+}_{{\bf k}})\right]+$ $\displaystyle
R^{\pm}({\bf k}+{\bf Q};{\bf p},\omega)R^{\pm}({\bf q}_{1},{\bf q}_{2},{\bf
k}+{\bf Q};{\bf p},\omega)\big{[}1+$ $\displaystyle F^{\pm,a}_{C}({\bf
q}_{1},{\bf q}_{2};{\bf p}+{\bf k},\omega-\omega^{-}_{{\bf k}})\big{]}.$ (66)

Here, we again use the definition for the recursive coupling given in Eq.
(61). To find $C^{\pm,l}_{0}$ we have to be a bit careful since the purple
vertex is seen as only acting on $\ket{\Psi^{\pm}_{{\bf p}}}$. In the
correlator, Eq. (65), we have $\hat{b}_{\mu_{1},-{\bf
q}_{1}}\hat{b}_{\mu_{2},-{\bf q}_{2}}$ such that $\hat{b}_{\mu_{1},-{\bf
q}_{1}}$ can either annihilate the first or the second spin wave. Using that
$c_{\mu_{1},\mu_{2}}({\bf q}_{1},{\bf q}_{2})=c_{\mu_{2},\mu_{1}}({\bf
q}_{2},{\bf q}_{1})$ this will yield a factor of two and we thereby find

$\displaystyle C^{\pm,l}_{0}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=\
Z^{\pm}_{{\bf p}}\big{(}R^{\pm}({\bf q}_{1};{\bf p},\omega)R^{\pm}({\bf
q}_{1},{\bf q}_{2};{\bf p},\omega)\times$ $\displaystyle c_{+,+}({\bf
q}_{1},{\bf q}_{2})\cdot\big{[}1+F^{\pm,s}_{C}({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)\big{]}$ $\displaystyle-R^{\pm}({\bf q}_{1}+{\bf Q};{\bf
p},\omega)R^{\pm}({\bf q}_{1}+{\bf Q},{\bf q}_{2}+{\bf Q};{\bf
p},\omega)c_{-,-}({\bf q}_{1},{\bf q}_{2})\times$
$\displaystyle\big{[}1+F^{\pm,s}_{C}({\bf q}_{1}+{\bf Q},{\bf q}_{2}+{\bf
Q};{\bf p},\omega)\big{]}$ $\displaystyle-(-1)^{l}R^{\pm}({\bf q}_{1};{\bf
p},\omega)R^{\pm}({\bf q}_{1},{\bf q}_{2}+{\bf Q};{\bf p},\omega)c_{+,-}({\bf
q}_{1},{\bf q}_{2})\times$ $\displaystyle\big{[}1+F^{\pm,a}_{C}({\bf
q}_{1},{\bf q}_{2}+{\bf Q};{\bf p},\omega)\big{]}$
$\displaystyle+(-1)^{l}R^{\pm}({\bf q}_{1}+{\bf Q};{\bf p},\omega)R^{\pm}({\bf
q}_{1}+{\bf Q},{\bf q}_{2};{\bf p},\omega)c_{-,+}({\bf q}_{1},{\bf
q}_{2})\big{]}\times$ $\displaystyle\big{[}1+F^{\pm,a}_{C}({\bf q}_{1}+{\bf
Q},{\bf q}_{2};{\bf p},\omega)\big{]}\big{)}.$ (67)

Here, the factors in front are a bit ambiguous, but they can trivially be
derived from in Eq. (54). As for the $B$-series we find $c_{+,-}=0=c_{-,+}$ in
the extreme limit of $t_{\perp}\gg t$, thereby, leaving the spins perfectly
antialigned between the layers.

With the definition of $C^{\pm,l}_{0}$ we can write the $C$ series as

(68)

This translates to

$C^{\pm,l}({\bf q}_{1},{\bf q}_{2};{\bf p},\omega)=\ C^{\pm,l}_{0}({\bf
q}_{1},{\bf q}_{2};{\bf p},\omega)+\\\ \sum_{{\bf k}}\big{[}R^{\pm}({\bf
k};{\bf p},\omega)R^{\pm}({\bf k};{\bf p},\omega)\cdot C^{\pm,l}({\bf
q}_{1},{\bf q}_{2};{\bf p}+{\bf k},\omega-\omega^{+}_{{\bf k}})+\\\
R^{\pm}({\bf k}+{\bf Q};{\bf p},\omega)R^{\pm}({\bf k}+{\bf Q};{\bf
p},\omega)\cdot C^{\pm,l}({\bf q}_{1},{\bf q}_{2};{\bf p}+{\bf
k},\omega-\omega^{-}_{{\bf k}})\big{]}$ (69)

As a sanity check, we find that for $t_{\perp}\rightarrow 0$ this yields the
result presented in [21] for a single layer.

##  Appendix F Vanishing diagrams

In this section, we will show that the nonsymmetric diagrams for the $B$
series and the rest of the diagrams for the $C$ series vanish. The structure
of the nonvanishing diagrams is such that for the $B$ series they will not be
symmetric around the center. For the $C$ series, the two spin waves
annihilated by the purple vertex will not be excited one after the other. An
example from both series is stated here

(70)

The feature determining whether the diagrams vanish or not is which terms in
the diagrams come with the same vertex factor, so $b_{\mu_{1},\mu_{2}}$ or
$c_{\mu_{1},\mu_{2}}$. This is not obvious, we, therefore, examine the
diagrams from Eq. (70).

For $B^{\pm,l,van.}_{2}$ the diagram translates to

$\displaystyle B^{\pm,l,van.}_{2}$ $\displaystyle({\bf q}_{1},{\bf q}_{2};{\bf
p},\omega)=\ \frac{Z^{\pm}_{{\bf p}}}{2}\big{[}$ $\displaystyle+b_{+,+}({\bf
q}_{1},{\bf q}_{2})R^{\pm}({\bf q}_{1};{\bf p},\omega)\sum_{{\bf
k}}R^{\pm}({\bf q}_{1},{\bf k};{\bf p},\omega)R^{\pm}({\bf k},{\bf q}_{2};{\bf
p},\omega)R^{\pm}({\bf k};{\bf p},\omega)$ $\displaystyle $ $\displaystyle-
b_{+,+}({\bf q}_{1},{\bf q}_{2})R^{\pm}({\bf q}_{1};{\bf p},\omega)\sum_{{\bf
k}}R^{\pm}({\bf q}_{1},{\bf k}+{\bf Q};{\bf p},\omega)R^{\pm}({\bf k}+{\bf
Q},{\bf q}_{2};{\bf p},\omega)R^{\pm}({\bf k}+{\bf Q};{\bf p},\omega)$
$\displaystyle $ $\displaystyle-(-1)^{l}b_{+,-}({\bf q}_{1},{\bf
q}_{2})R^{\pm}({\bf q}_{1};{\bf p},\omega)\sum_{{\bf k}}R^{\pm}({\bf
q}_{1},{\bf k};{\bf p},\omega)R^{\pm}({\bf k},{\bf q}_{2}+{\bf Q};{\bf
p},\omega)R^{\pm}({\bf k};{\bf p},\omega)$ $\displaystyle $
$\displaystyle+(-1)^{l}b_{+,-}({\bf q}_{1},{\bf q}_{2})R^{\pm}({\bf
q}_{1};{\bf p},\omega)\sum_{{\bf k}}R^{\pm}({\bf q}_{1},{\bf k}+{\bf Q};{\bf
p},\omega)R^{\pm}({\bf k}+{\bf Q},{\bf q}_{2}+{\bf Q};{\bf
p},\omega)R^{\pm}({\bf k}+{\bf Q};{\bf p},\omega)$ $\displaystyle $
$\displaystyle+(-1)^{l}b_{-,+}({\bf q}_{1},{\bf q}_{2})R^{\pm}({\bf
q}_{1}+{\bf Q};{\bf p},\omega)\sum_{{\bf k}}R^{\pm}({\bf q}_{1}+{\bf Q},{\bf
k};{\bf p},\omega)R^{\pm}({\bf k},{\bf q}_{2};{\bf p},\omega)R^{\pm}({\bf
k};{\bf p},\omega)$ $\displaystyle $ $\displaystyle-(-1)^{l}b_{-,+}({\bf
q}_{1},{\bf q}_{2})R^{\pm}({\bf q}_{1}+{\bf Q};{\bf p},\omega)\sum_{{\bf
k}}R^{\pm}({\bf q}_{1}+{\bf Q},{\bf k}+{\bf Q};{\bf p},\omega)R^{\pm}({\bf
k}+{\bf Q},{\bf q}_{2};{\bf p},\omega)R^{\pm}({\bf k}+{\bf Q};{\bf p},\omega)$
$\displaystyle $ $\displaystyle-b_{-,-}({\bf q}_{1},{\bf q}_{2})R^{\pm}({\bf
q}_{1}+{\bf Q};{\bf p},\omega)\sum_{{\bf k}}R^{\pm}({\bf q}_{1}+{\bf Q},{\bf
k};{\bf p},\omega)R^{\pm}({\bf k},{\bf q}_{2}+{\bf Q};{\bf
p},\omega)R^{\pm}({\bf k};{\bf p},\omega)$ $\displaystyle $
$\displaystyle+b_{-,-}({\bf q}_{1},{\bf q}_{2})R^{\pm}({\bf q}_{1}+{\bf
Q};{\bf p},\omega)\sum_{{\bf k}}R^{\pm}({\bf q}_{1}+{\bf Q},{\bf k}+{\bf
Q};{\bf p},\omega)R^{\pm}({\bf k}+{\bf Q},{\bf q}_{2}+{\bf Q};{\bf
p},\omega)R^{\pm}({\bf k}+{\bf Q};{\bf p},\omega)\big{]}$ (71)

Here we see that by taking ${\bf k}\rightarrow{\bf k}+{\bf Q}$ for every
second term this diagram yields zero. The recurring feature leading to these
terms being zero, in comparison to the symmetric version, is that the terms
coming with the same vertex amplitude, $b_{\mu_{1},\mu_{2}}$, can be put into
pairs of the kind $a_{\\{\mu_{n}\\},\mu_{1},\mu,\\{\mu_{m}\\}}\cdot
a_{\\{\mu_{n}\\},\mu,\mu_{2},\\{\mu_{m}\\}}$ and
$a_{\\{\mu_{n}\\},\mu_{1},\bar{\mu},\\{\mu_{m}\\}}\cdot
a_{\\{\mu_{n}\\},\bar{\mu},\mu_{2},\\{\mu_{m}\\}}$, where $\bar{\mu}$ is the
opposite type of spin wave to $\mu$, and the associated momentum is summed
over. If looking at Eq. (49), we see that such terms come with an opposite
sign. Now, since the moment associated with $\mu_{2}$ is summed over, we can
translate them by ${\bf Q}$ such that the two terms will sum up to zero. The
equivalent symmetric diagrams will instead come with the same sign in front.

For the C series, the vanishing diagrams yield zero for similar reasons. Here,
terms coming with the same vertex factor can be grouped in pairs of
$a_{\\{\mu_{n}\\},\mu,\\{\mu_{m}\\}}\cdot
a_{\\{\mu_{n},\mu_{1}\\},\mu,\mu_{2},\\{\mu_{m}\\}}$ and
$a_{\\{\mu_{n}\\},\bar{\mu},\\{\mu_{m}\\}}\cdot
a_{\\{\mu_{n},\mu_{1}\\},\bar{\mu},\mu_{2},\\{\mu_{m}\\}}$, where the momentum
associated with $\mu$ is summed over. One should understand
$\\{\mu_{n},\mu_{2}\\}$ as series of excited spin waves where $\mu_{2}$
excited at some unspecified order. Looking at Eq. (50), we see that these two
terms come with a sign difference. We can therefore translate the momentum by
${\bf Q}$ to have them cancel one another.

##  Appendix G Symmetries of the polaron cloud

In this appendix we use the expression for the $B$\- and $C$-series, Eq. (27)
and (28), to show the symmetry

$\displaystyle M_{l}^{-}({\bf d},{\bf p}+{\bf Q})$
$\displaystyle=M_{l}^{+}({\bf d},{\bf p}).$ (72)

Looking at the recursive coupling, Eq. (48), we find the following symmetry
related to the AF ordering vector

$\displaystyle R^{+}({\bf p}+{\bf Q},\\{{\bf k}_{n}\\})=-R^{-}({\bf p},\\{{\bf
k}_{n}\\}).$ (73)

In the expression for the $B$\- and $C$-series, Appendix D and E, we find that
all $R^{\pm}$ are multiplied by even numbers of $R^{\pm}$s. The same applies
to the residue $Z^{\pm}_{{\bf p}}$, such that this minus will change nothing
and we retrieve $M_{l}^{+}({\bf d},{\bf p})=M_{l}^{-}({\bf d},{\bf p}+{\bf
Q})$. Realizing that $\Psi^{+}_{{\bf p}+{\bf Q}}$ and $\Psi^{-}_{{\bf p}}$ are
degenerate states, see Ref. [45], we find that this manifests itself in the
same frustration of the magnetic environment.

##  Appendix H Polaron cloud for $J_{\perp}\neq 4t^{2}_{\perp}/U$

In this Appendix, we investigate the influence of $J_{\perp}$ and $t_{\perp}$
separately in relation to how the polaron clouds in the two layers mimic one
another for large $J_{\perp}/t$. To do so, we deviate from
$J_{\perp}=4t^{2}_{\perp}/U$ used in the main text, such that the model no
longer describes the low energy physics of the Fermi-Hubbard model in the
large $U$ limit. In analog to Fig. 7(a) in the main text, Fig. 9 show
$M^{-}_{1}-M^{-}_{2}$ for $t_{\perp}/t=0$ (a) and $J_{\perp}/t=0$ (b) with
${\bf p}=(0,0)$. The grey curves show the nearest neighbor results plotted in
the main text for $J_{\perp}=4t^{2}_{\perp}/U$. Comparing the orange and grey
curves in Fig. 7(a), we see that a vanishing interlayer hopping,
$t_{\perp}/t=0$, yields a less prominent mirroring effect between the layers.
In Fig. 7(b), for $J_{\perp}/t=0$, the orange and grey curves follow each
other more closely, though $M^{-}_{1}-M^{-}_{2}$ stays larger with no
interlayer spin-spin interaction.

Figure 9: Difference in magnetization between the layers for $J/t=0.3$ and
${\bf p}=(0,0)$. (a) Here shown as a function of the interlayer spin-spin
correlation with $t_{\perp}/t=0$. $\sqrt{J_{\perp}/J}$ is chosen since
$\sqrt{J_{\perp}/J}=t_{\perp}/t$ for $J_{\perp}=4t^{2}_{\perp}/U$, and,
therefore, enables easy comparison. In (b), we instead show it as a function
of the interlayer hopping $t_{\perp}/t$ with $J_{\perp}/J=0$. The grey curves
in both panels correspond to the nearest neighbor difference when
$J_{\perp}=4t^{2}_{\perp}/U$, and should be compared to the orange curve.
Doing so, we see that in (a) the mirroring effect is less prominent while in
(b) it follows more closely. The inset shows the energy of the ${\bf
p}=(\pi/2,\pi/2)$ state in black and ${\bf p}=(0,0)$ in red for
$J_{\perp}/J=0$. The ${\bf p}=(\pi/2,\pi/2)$ state, thus, remains the ground
state with increasing interlayer hopping for $t_{\perp}=0$.

The mirroring effect where the cloud upholds the interlayer AF order is,
therefore, primarily created by interlayer hopping. That said, the energy
decrease associated with keeping the AF order between the layers is absent in
the case of $J_{\perp}/t\neq 0$. As seen in the inset in Fig. 7(b), the system
will, consequently, not prefer ${\bf p}=(0,0)$ as its ground state. As a
result, in order for the system to prefer the ${\bf p}=(0,0)$ polaron state as
its ground state, we need the interlayer spin-spin interaction.

##  Appendix I Deviating from $J/t=0.3$

This appendix discusses the influence of the value $J/t$ on the magnetic
dressing cloud. Similar to Fig. 5 in the main text, we show in Fig. 10 the
development of the nearest neighbor and next-nearest neighbor frustrations in
layer 1 and 2, respectively, as a function of interlayer hopping for
$J/t=0.1,0.3,1.0$. Here, we see that they all follow the same behavior but
with increasing frustration for smaller $J/t$. This is because for smaller
$J/t$ the energy penalty of being localized increases and the hole will be
more delocalized, thereby, increasing the frustration. The features described
in the main text are independent of this battle between the two terms since it
reflects the behavior of the order parameter dictated by $J_{\perp}/J$. By
changing $J/t$ the behavior described in the main text remains qualitative the
same, but will quantitatively take place with more ($J/t<0.3$) or less
($J/t>0.3$) overall frustration.

Figure 10: This figure shows the frustration of the nearest neighboring site
in layer 1 where the hole is present and the opposite sites in layer 2. This
is plotted for $J/t=0.1,0.3,1.0$ with varying interlayer hopping,
$t_{\perp}/t$. Here, we see the same behavior as described in the main text
with the overall frustration being larger for smaller $J/t$.

## References

  * Kane  _et al._ [1989] C. L. Kane, P. A. Lee, and N. Read, Motion of a single hole in a quantum antiferromagnet, Physical Review B 39, 6880 (1989).
  * Wen and Li [2011] H.-H. Wen and S. Li, Materials and Novel Superconductivity in Iron Pnictide Superconductors, Annual Review of Condensed Matter Physics 2, 121 (2011).
  * Wosnitza [2012] J. Wosnitza, Superconductivity in Layered Organic Metals, Crystals 2, 248 (2012).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * Esslinger [2010] T. Esslinger, Fermi-hubbard physics with atoms in an optical lattice, Annual Review of Condensed Matter Physics 1, 129 (2010), https://doi.org/10.1146/annurev-conmatphys-070909-104059 .
  * Boll  _et al._ [2016] M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, Spin- and density-resolved microscopy of antiferromagnetic correlations in fermi-hubbard chains, 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, Observation of spatial charge and spin correlations in the 2d fermi-hubbard model, 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, A cold-atom fermi–hubbard antiferromagnet, 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, Revealing hidden antiferromagnetic correlations in doped hubbard chains via string correlators, 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, Spin-imbalance in a 2d fermi-hubbard system, Science 357, 1385 (2017).
  * Chiu  _et al._ [2018] C. S. Chiu, G. Ji, A. Mazurenko, D. Greif, and M. Greiner, Quantum state engineering of a hubbard system with ultracold fermions, 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, Bad metallic transport in a cold atom Fermi-Hubbard system, Science 363, 379 (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, String patterns in the doped Hubbard model, Science 365, 251 (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, Imaging magnetic polarons in the doped Fermi–Hubbard model, Nature 572, 358 (2019).
  * 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, Coupling a Mobile Hole to an Antiferromagnetic Spin Background: Transient Dynamics of a Magnetic Polaron, Physical Review X 11, 021022 (2021).
  * Koepsell  _et al._ [2021] 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, Science 374, 82 (2021).
  * Carlström  _et al._ [2016] J. Carlström, N. Prokof’ev, and B. Svistunov, Quantum Walk in Degenerate Spin Environments, Physical Review Letters 116, 247202 (2016).
  * Kanász-Nagy  _et al._ [2017] M. Kanász-Nagy, I. Lovas, F. Grusdt, D. Greif, M. Greiner, and E. A. Demler, Quantum correlations at infinite temperature: The dynamical Nagaoka effect, Physical Review B 96, 014303 (2017).
  * Grusdt  _et al._ [2018a] F. Grusdt, Z. Zhu, T. Shi, and E. Demler, Meson formation in mixed-dimensional t-J models, SciPost Physics 5, 057 (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, Parton Theory of Magnetic Polarons: Mesonic Resonances and Signatures in Dynamics, Physical Review X 8, 011046 (2018b).
  * Nielsen  _et al._ [2021] K. K. Nielsen, M. A. Bastarrachea-Magnani, T. Pohl, and G. M. Bruun, Spatial structure of magnetic polarons in strongly interacting antiferromagnets, Physical Review B 104, 155136 (2021).
  * Nielsen  _et al._ [2022] K. K. Nielsen, T. Pohl, and G. M. Bruun, Non-equilibrium hole dynamics in antiferromagnets: Damped strings and polarons, arXiv:2203.04789v1 [cond-mat.quant-gas] 10.48550/arXiv.2203.04789 (2022), arxiv:2203.04789v1 [cond-mat.quant-gas] .
  * Bakr  _et al._ [2009] W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, and M. Greiner, A quantum gas microscope for detecting single atoms in a hubbard-regime optical lattice, Nature 462, 74 (2009).
  * Sherson  _et al._ [2010] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Single-atom-resolved fluorescence imaging of an atomic mott insulator, 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, Single-atom imaging of fermions in a quantum-gas microscope, Nature Physics 11, 738 (2015).
  * Yang  _et al._ [2021] J. Yang, L. Liu, J. Mongkolkiattichai, and P. Schauss, Site-resolved imaging of ultracold fermions in a triangular-lattice quantum gas microscope, PRX Quantum 2, 020344 (2021).
  * Gall  _et al._ [2021] M. Gall, N. Wurz, J. Samland, C. F. Chan, and M. Köhl, Competing magnetic orders in a bilayer Hubbard model with ultracold atoms, Nature 589, 40 (2021).
  * Hida [1992] K. Hida, Quantum Disordered State without Frustration in the Double Layer Heisenberg Antiferromagnet –Dimer Expansion and Projector Monte Carlo Study–, Journal of the Physical Society of Japan 61, 1013 (1992).
  * Sandvik and Scalapino [1994] A. W. Sandvik and D. J. Scalapino, Order-disorder transition in a two-layer quantum antiferromagnet, Physical Review Letters 72, 2777 (1994).
  * Scalettar  _et al._ [1994] R. T. Scalettar, J. W. Cannon, D. J. Scalapino, and R. L. Sugar, Magnetic and pairing correlations in coupled Hubbard planes, Physical Review B 50, 13419 (1994).
  * Millis and Monien [1994] A. J. Millis and H. Monien, Spin gaps and bilayer coupling in YBa 2 Cu 3 O 7 - $\delta$ and YBa 2 Cu 4 O 8, Physical Review B 50, 16606 (1994).
  * Sandvik  _et al._ [1995] A. W. Sandvik, A. V. Chubukov, and S. Sachdev, Quantum critical behavior in a two-layer antiferromagnet, Physical Review B 51, 16483 (1995).
  * Chubukov and Morr [1995] A. V. Chubukov and D. K. Morr, Phase transition, longitudinal spin fluctuations, and scaling in a two-layer antiferromagnet, Physical Review B 52, 3521 (1995).
  * Martinez and Horsch [1991] G. Martinez and P. Horsch, Spin polarons in the _t_ \- _J_ model, Physical Review B 44, 317 (1991).
  * Liu and Manousakis [1991] Z. Liu and E. Manousakis, Spectral function of a hole in the _t_ \- _J_ model, Physical Review B 44, 2414 (1991).
  * Marsiglio  _et al._ [1991] F. Marsiglio, A. E. Ruckenstein, S. Schmitt-Rink, and C. M. Varma, Spectral function of a single hole in a two-dimensional quantum antiferromagnet, Physical Review B 43, 10882 (1991).
  * Ramšak and Horsch [1998] A. Ramšak and P. Horsch, Spatial structure of spin polarons in the t - J model, Physical Review B 57, 4308 (1998).
  * Chernyshev and Leung [1999] A. L. Chernyshev and P. W. Leung, Holes in the t - J z model: A diagrammatic study, Physical Review B 60, 1592 (1999).
  * Wohlfeld  _et al._ [2008] K. Wohlfeld, M. Daghofer, A. M. Oleś, and P. Horsch, Spectral properties of orbital polarons in Mott insulators, Physical Review B 78, 214423 (2008).
  * Diamantis and Manousakis [2021] N. G. Diamantis and E. Manousakis, Dynamics of string-like states of a hole in a quantum antiferromagnet: A diagrammatic Monte Carlo simulation, New Journal of Physics 23, 123005 (2021).
  * Reiter [1994] G. F. Reiter, Self-consistent wave function for magnetic polarons in the _t_ \- _J_ model, Physical Review B 49, 1536 (1994).
  * Kancharla and Okamoto [2007] S. S. Kancharla and S. Okamoto, Band insulator to Mott insulator transition in a bilayer Hubbard model, Physical Review B 75, 193103 (2007), arxiv:cond-mat/0703728 .
  * Bouadim  _et al._ [2008] K. Bouadim, G. G. Batrouni, F. Hébert, and R. T. Scalettar, Magnetic and transport properties of a coupled Hubbard bilayer with electron and hole doping, Physical Review B 77, 144527 (2008).
  * Vojta and Becker [1999] M. Vojta and K. W. Becker, Doped bilayer antiferromagnets: Hole dynamics on both sides of a magnetic ordering transition, Physical Review B 60, 15201 (1999).
  * Nyhegn  _et al._ [2022] J. H. Nyhegn, K. K. Nielsen, and G. M. Bruun, Equilibrium and non-equilibrium dynamics of a hole in a bilayer antiferromagnet 10.48550/ARXIV.2208.04653 (2022).
  * Rademaker  _et al._ [2012] L. Rademaker, K. Wu, and J. Zaanen, Dynamics of a single exciton in strongly correlated bilayers, New Journal of Physics 14, 083040 (2012).
  * Holt  _et al._ [2012] M. Holt, J. Oitmaa, W. Chen, and O. P. Sushkov, Fermi Surface Reconstruction by Dynamic Magnetic Fluctuations, Physical Review Letters 109, 037001 (2012).
  * Holt  _et al._ [2013] M. Holt, J. Oitmaa, W. Chen, and O. P. Sushkov, Fermi surface reconstruction by dynamic magnetic fluctuations and spin-charge separation near an O ( 3 ) quantum critical point, Physical Review B 87, 075109 (2013).
  * White and Scalapino [1997] S. R. White and D. J. Scalapino, Hole and pair structures in the t-J model, Physical Review B 55, 6504 (1997).
  * Maier and Scalapino [2011] T. A. Maier and D. J. Scalapino, Pair structure and the pairing interaction in a bilayer Hubbard model for unconventional superconductivity, Physical Review B 84, 180513 (2011).