Localized dopant motion across the 2D Ising phase transition

Kristian Knakkergaard Nielsen1$\star$

1 Max-Planck-Institut für Quantentoptik, D-85748 Garching, Germany

$\star$<EMAIL_ADDRESS>
## Abstract

I investigate the motion of a single hole in 2D spin lattices with square and
triangular geometries. While the spins have nearest neighbor Ising spin
couplings $J$, the hole is allowed to move only in 1D along a single line in
the 2D lattice with nearest neighbor hopping amplitude $t$. The non-
equilibrium hole dynamics is initialized by suddenly removing a single spin
from the thermal Ising spin lattice. I find that for any nonzero spin coupling
and temperature, the hole is _localized_. This is an extension of the
thermally induced localization phenomenon [1] to the case, where there is a
phase transition to a long-range ordered ferromagnetic phase. The dynamics
depends only on the ratio of the temperature to the spin coupling,
$k_{B}T/|J|$, and on the ratio of the spin coupling to the hopping $J/t$. I
characterize these dependencies in great detail. In particular, I find
universal behavior at high temperatures, common features for the square and
triangular lattices across the Curie temperatures for ferromagnetic
interactions, and highly distinct behaviors for the two geometries in the
presence of antiferromagnetic interactions due geometric frustration in the
triangular lattice.

Copyright attribution to authors.
This work is a submission to SciPost Physics.
License information to appear upon publication.
Publication information to appear upon publication. Received Date
Accepted Date
Published Date

###### Contents

  1. 1 Introduction
  2. 2 System and setup
  3. 3 Monte Carlo sampling of exact trajectories
  4. 4 Results
    1. 4.1 Infinite temperatures
    2. 4.2 Ferromagnetic couplings
    3. 4.3 Antiferromagnetic couplings
    4. 4.4 Spin coupling scaling dependency
  5. 5 Discussion
    1. 5.1 Beyond 1D dopant motion
    2. 5.2 Beyond the Ising model
    3. 5.3 Increasing the doping level
    4. 5.4 Connecting to an external bath
  6. 6 Conclusions and outlook
  7. A Short-range correlators: square lattice
  8. B Short-range correlators: triangular lattice
  9. C Localization length in the ferromagnetic phase for large $|J|/t$
  10. D Entropy for antiferromagnetic couplings

##  1 Introduction

Is an impurity or a dopant necessarily delocalized in a system with polarized
long-range order? As the system is almost perfectly homogeneous, the immediate
response would presumably be ’yes’! However, if the dopant’s motion is
correlated with the background in which it moves, then it modifies the
background as it moves. Therefore, this question is more subtle than one would
immediately expect. Such a scenario arises in its most simplistic form in an
Ising magnet with a doped hole. Here, the ability of the spins to hop only
onto vacant sites, means that the motion of holes is completely contingent on
a counterpropagating spin. While the hole in a ferromagnetic ground state will
certainly delocalize, the same may not remain true as soon as one heats up the
system by any infinitesimal amount. The general scenario turns out to be
tremendously complex to analyze, however. In particular, the exponential
growth in the number of possible configurations of the system as the hole
moves away from its origin means that analytical treatments are generically
out of the window. However, numerical investigations of full two-dimensional
(2D) motion of holes in thermal spin ensembles have been carried out [2, 3,
4]. These intriguing results are unfortunately limited by the underlying
exponential complexity of the dynamics in a generic spin state to short times
and/or small systems, and no robust conclusions have been found for long
timescales in the thermodynamic limit. A mean-field approach [5] has, however,
shown the possibility of rich phase diagram, supporting stripe formation [6]
at low temperatures.

In a broader context, dopants in magnetic lattices have been extensively
studied for decades [7, 8, 9, 10, 11] due to their intimate links to high-
temperature superconductivity [12, 13, 14]. In recent years, this line of
research has seen a fruitful revival thanks to advances in quantum simulation
experiments with ultracold atoms in optical lattices [15, 16, 17, 18, 19, 20,
21, 22]. Such experiments can also manipulate the models otherwise fixed in
the solid state. In particular, clear signatures of magnetically mediated
pairing of dopants [22] has been observed and boosted by only allowing holes
to move in one dimension (1D)_along_ the investigated ladder geometry.
Building on these successes, I recently investigated such a model idealised
further to support only Ising type spin couplings [1]. Contrary to earlier
setups [2, 3, 4], 1D motion in the an Ising magnet facilitates the
investigation of very large system sizes and essentially arbitrarily long
evolution times can be achieved. This allowed me to systematically show that
the hole is _localized_ for any nonzero temperature and any nonzero spin
coupling, and only asymptotically delocalize in the limits $\beta
J=J/k_{B}T\to-\infty$ and/or $|J|/t\to 0$. In other words, even though there
are perfect quasiparticle excitations in these limits, in which they behave
exactly as free particles, any nonzero temperature and spin coupling
immediately localizes the hole.

Figure 1: The mixed-dimensional $t$-$J_{z}$ model with 1D dopant motion (green
region) in a square (a) and triangular (b) spin lattice. Here, both
spin-$\uparrow$ (red spheres) and $-\downarrow$ particles (blue spheres) may
hop to vacant sites (green circles) – holes – with hopping amplitude $t$
_along_ the system. The spins are coupled via isotropic Ising interactions
with magnitude $J$ (black lines).

From a statistical mechanics perspective, the ladder geometry is, however, a
one-dimensional system. Consequently, there is no phase transition for the
underlying Ising spins until zero temperature is reached. One might think,
therefore, that the hole localization is linked directly to the disordered
nature of the finite temperature (Gibbs) state. This is supported in my
previous paper by the fact that the localization length becomes proportional
to the spin-spin correlation length at low temperatures. Since this length
scale diverges as the para- to ferromagnetic phase transition is crossed from
above, this would seem to indicate that for a fully 2D spin lattice, the hole
should delocalize as the transition to a ferromagnetic phase is crossed at the
Curie temperature. As a result, the present Article is concerned with
addressing exactly this question: how is the 1D motion of the hole affected in
the presence of the Ising phase transition for the 2D Ising model? To study
this carefully, I will consider two lattice geometries: the square and
triangular lattices with nearest neighbor Ising interactions, see Fig. 1.
These support such a para- to ferromagnetic phase transition at their
respective Curie temperatures. The comparison of the two lattice structures
allow me to extract the common features of the systems and highlight the
importance of the change in the number of nearest neighbor spin couplings
(from $4$ for the square lattice to $6$ in the triangular case). To my
surprise, the realization of a long-range ordered ferromagnetic phase is _not_
enough to delocalize the hole. Instead, I find that the exponentially small
probability of meeting spin flips in the ferromagnetic phase is more than
adequate to keep the hole localized.

The localization effect in the two-leg ladder was tied directly to thermal
spin fluctuations of the disordered spin lattice. This realizes a novel
variant of Anderson localization in the presence of _strong_ disorder, in
which the dopant back-scatter off the spin fluctuations, as the energy cost of
propagating further away from its origin will inevitably fluctuate to values
larger than the initial kinetic energy of the hole. This framework also
explains the localization quantitatively well in the disordered phase above
the Curie temperature in this present setup. However, in the ordered phase
this picture breaks down for intermediate to large values of $|J|/t$. In
particular, the thermal spin fluctuations in the long-range ordered phase on
short length scales tend to happen as singular spin flips. This results in a
crossover between large and small values of $|J|/t$. For large $|J|/t$, the
hole backscatters off these single spin flips. For intermediate values, it
tunnels through many such flips, but destructively interfere for different
pathways to a specific point because there is a large statistical variation in
how long the segments are between such spin flips. This is the physics of
Anderson localization in the presence of _weak_ disorder [23]. Finally, for
sufficiently small $|J|/t$ this effect is too weak. Instead the hole once
again backreflects once the build up of potential energy overcomes its initial
kinetic energy.

Moreover, for a square lattice its bipartite structure results in a symmetry
between the ferro- and antiferromagnetic case, such that the thermodynamics of
the underlying spin systems are equivalent. Therefore, the Neél critical
temperature is simply the same as the Curie temperature for ferromagnetic
couplings. However, the motion of the hole is markedly different in the two
scenarios. In the antiferromagnet, the buildup of the staggered Neél order
leads to a linearly increasing potential which becomes stronger and stronger
as zero temperature is approached. In this case, there is, therefore, a
crossover between localization driven by spin fluctuations at high
temperatures to localization exerted by an effective confining potential at
low temperatures. The latter eventually leads to heavy coherent oscillations
in the position of the hole due to interference of the low-energy states as
zero temperature is approached, analogous to what was found for the two-leg
ladder at zero temperature [24].

In stark constrast, in the triangular case, antiferromagnetic couplings hinder
any phase transition all the way down to zero temperature and results in a
residual entropy $S_{0}/N\simeq 0.323$ of the ground state manifold, due to an
exponentially large ground state degeneracy [25, 26]. The present setup,
hereby, also allows me to study the influence of this extensive ground state
degeneracy. We will see that while the behavior for the square antiferromagnet
is characterized by more and more coherent oscillations of the hole’s motion,
the same is _not_ true for the triangular antiferromagnetic case, and the
dynamics retains a thermal character with smooth and incoherent behavior, even
as zero temperature is approached.

From a traditional condensed matter point of view, it is hard to imagine how
one would actually investigate the dynamical phenomena detailed above, as it
requires one to track the motion of dopants in real time. However, quantum
simulation experiments with ultracold atoms have made substantial
breakthroughs in this regard. Not only do these systems allow for single site
detection [27, 28], but has enough sensitivity to actually track the motion of
holes in real time [16]. Moreover, the Ising type of interactions investigated
are e.g. facilitated with Rydberg-dressed atoms in optical lattices [29],
which crucially still allow for the motion of dopants. Finally, the 1D
restriction of the hole motion have been achieved for dopants in such spin
lattices [22]. In this manner, one can combine these well-established
experimental capabilities to observe the predicted effects, as I have also
previously pointed out in a suggestion for an explicit experimental protocol
[1].

The Article is organized as follows. In Sec. 2, the system and setup is
explained. In Sec. 3, I describe the employed Monte-Carlo simulation of exact
trajectories. In Sec. 4, a detailed analysis of the numerical results is given
at infinite temperature [Sec. 4.1], for ferromagnetic [Sec. 4.2] and
antiferromagnetic couplings [Sec. 4.2], and finally the general spin-coupling
dependency [Sec. 4.4]. Moreover, a detailed discussion on the extension of the
model is given in Sec. 5, before I conclude in Sec. 6.

##  2 System and setup

In this paper, I will consider the mixed-dimensional $t$-$J_{z}$ model

$\displaystyle{\hat{H}}$
$\displaystyle={\hat{H}}_{t}+{\hat{H}}_{J}=-t\sum_{\braket{{\bf i},{\bf
j}}_{\parallel},\sigma}\left[\tilde{c}^{\dagger}_{{\bf
i},\sigma}\tilde{c}^{\phantom{\dagger}}_{{\bf j},\sigma}+{\rm
h.c.}\right]+J\sum_{\braket{{\bf i},{\bf j}}}S^{(z)}_{\bf i}S^{(z)}_{\bf j},$
(1)

in the presence of a single dopant – a hole. Here, the correlated motion of
the dopant is allowed through the nearest neighbor hopping Hamiltonian
${\hat{H}}_{t}$, with constrained operators $\tilde{c}^{\dagger}_{{\bf
i},\sigma}=\hat{c}^{\dagger}_{{\bf i},\sigma}(1-\hat{n}_{{\bf
i}\bar{\sigma}})$ to ensure at most a single particle per site. Note that
$\bar{\sigma}$ designates the opposite spin of $\sigma$, i.e.
$\bar{\uparrow}=\downarrow$ and vice versa. Also, while the Ising spin
couplings ${\hat{H}}_{J}$ couple all nearest neighbors isotropically, the
hopping is _only_ allowed along a 1D line, indicated by $\braket{{\bf i},{\bf
j}}_{\parallel}$. In a recent paper [1], I studied this model in a two-leg
square ladder system. Here, I will extend the studies to a 2D spin lattice and
encompass both square and triangular lattices. The main idea is to study the
importance of the appearing Ising phase transition on the motion of the single
hole, as well as understanding the importance of magnetic frustration in the
triangular case for antiferromagnetic couplings.

To easily describe hole and spin degrees of freedom, I employ an exact
Holstein-Primakoff transformation with the ferromagnetic state $\ket{{\rm
FM}}=\ket{\cdots\uparrow\uparrow\cdots}$ with all spins pointin up as the
reference state. This leads to the alternate expressions for the hopping

$\displaystyle{\hat{H}}_{t}=t\sum_{\braket{{\bf i},{\bf j}}_{\parallel}}$
$\displaystyle\Big{[}\hat{h}^{\dagger}_{{\bf j}}F(\hat{h}_{{\bf
i}},\hat{s}_{{\bf i}})F(\hat{h}_{{\bf j}},\hat{s}_{{\bf
j}})\hat{h}^{\phantom{\dagger}}_{{\bf i}}+\hat{h}^{\dagger}_{{\bf
j}}\hat{s}^{\dagger}_{\bf i}F(\hat{h}_{{\bf i}},\hat{s}_{{\bf
i}})F(\hat{h}_{{\bf j}},\hat{s}_{{\bf j}})\hat{s}^{\phantom{\dagger}}_{\bf
j}\hat{h}^{\phantom{\dagger}}_{\bf i}\Big{]}+{\rm H.c.},$ (2)

and the spin coupling Hamiltonians

$\displaystyle\hat{H}_{J}=J\sum_{\braket{{\bf i},{\bf
j}}}\Big{[}\frac{1}{2}-\hat{s}^{\dagger}_{{\bf
i}}\hat{s}^{\phantom{\dagger}}_{{\bf
i}}\Big{]}\Big{[}\frac{1}{2}-\hat{s}^{\dagger}_{{\bf
j}}\hat{s}^{\phantom{\dagger}}_{{\bf
j}}\Big{]}\Big{[}1-\hat{h}^{\dagger}_{{\bf
i}}\hat{h}^{\phantom{\dagger}}_{{\bf
i}}\Big{]}\Big{[}1-\hat{h}^{\dagger}_{{\bf
j}}\hat{h}^{\phantom{\dagger}}_{{\bf j}}\Big{]}.$ (3)

Here, the spin excitation operator $\hat{s}^{\dagger}_{\bf i}$ is bosonic, and
creates a spin-$\downarrow$ on site ${\bf i}$. Also, the hole is created by
the operator $\hat{h}^{\dagger}_{\bf i}$, and maintains the statistics of the
underlying spins, be it fermionic _or_ bosonic [24]. In the hopping
Hamiltonian ${\hat{H}}_{t}$, the operator
$F(\hat{h},\hat{s})=\sqrt{1-\hat{s}^{\dagger}\hat{s}-\hat{h}^{\dagger}\hat{h}}$
keeps the single-occupancy constraint in check. This construction enables me
to succintly describe the motion of holes.

##  3 Monte Carlo sampling of exact trajectories

The non-equilibrium dynamics of the holes is initialized in the following
manner. The system starts out in the absence of holes in its thermal Gibs
state, $\hat{\rho}_{J}=e^{-\beta\hat{H}_{J}}/Z_{J}$. I assume nothing about
how this equilibrium is initially established, but from the time of the
quenched insertion of the hole at $\tau=0$, I assume the system to be closed.
The state of the system immediately after the removal of a spin from the
origin ${\bf i}={\bf 0}$ is

$\displaystyle\hat{\rho}(\tau=0)=\sum_{\sigma_{\bf
0}}\hat{c}^{\phantom{\dagger}}_{{\bf 0},\sigma_{\bf
0}}\hat{\rho}_{J}\hat{c}^{\dagger}_{{\bf 0},\sigma_{\bf
0}}=\hat{h}^{\dagger}_{\bf 0}\hat{\rho}_{J}\hat{h}^{\phantom{\dagger}}_{\bf
0}+\hat{h}^{\dagger}_{\bf 0}\hat{s}^{\phantom{\dagger}}_{\bf
0}\hat{\rho}_{J}\hat{s}^{\dagger}_{\bf 0}\hat{h}^{\phantom{\dagger}}_{\bf 0}$
(4)

After that, since the system is assumed to be closed, it evolves unitarily
under the full Hamiltonian ${\hat{H}}$ in Eq. (1),
$\hat{\rho}(\tau)=e^{-i{\hat{H}}\tau}\hat{\rho}(0)e^{+i{\hat{H}}\tau}$. Next,
I express the density operator in the Ising basis with spin configurations
${\boldsymbol{\sigma}}$. Using this, I can write the time-evolved density
matrix as the Boltzmann-weighted sum of pure-state evolutions

$\hat{\rho}(\tau)=\sum_{\sigma_{0},{\boldsymbol{\sigma}}}\frac{e^{-\beta
E_{J}(\sigma_{0},{\boldsymbol{\sigma}})}}{Z_{0}}\ket{\Psi_{{\boldsymbol{\sigma}}}(\tau)}\bra{\Psi_{{\boldsymbol{\sigma}}}(\tau)},$
(5)

where $E_{J}(\sigma_{0},{\boldsymbol{\sigma}})$ is the magnetic energy of the
spin realization $\sigma_{0},{\boldsymbol{\sigma}}$ before the hole is
introduced. With the hole and spin excitation operators at hand, we may
express the non-equilibrium pure states
$\ket{\Psi_{{\boldsymbol{\sigma}}}(\tau)}$ quite concisely. In particular, I
consider a system of size $(2N_{x}+1)\times(2N_{y}+1)$ with open boundary
conditions, such that the coordinates is written as ${\bf i}=x,y$ with
$x\in\\{-N_{x},-N_{x}+1,\dots,N_{x}\\}$ and
$y\in\\{-N_{y},-N_{y}+1,\dots,N_{y}\\}$, and with the hole moving along the
$y=0$ leg of the 2D system. Here, the triangular lattice is implemented as a
square lattice but with one additional spin coupling along one of the
diagonals of the lattice [30].

In terms of the spin excitation operators, a certain subset
$S^{y}_{\boldsymbol{\sigma}}$ in each leg $y$ will have spin flips. This means
that the initial wave function for such a realization can be expressed as

$\ket{\Psi_{{\boldsymbol{\sigma}}}(\tau=0)}=\hat{h}^{\dagger}_{0,0}\prod_{y=-N_{y}}^{N_{y}}\prod_{j\in
S_{{\boldsymbol{\sigma}}}^{y}}\hat{s}^{\dagger}_{j,y}\ket{{\rm FM}}.$ (6)

When the hole moves along the system, the spins in leg $y=0$ countermove by a
single lattice site, while all other spins remain static. Therefore, the state
at any later time $\tau$ is

$\displaystyle\ket{\Psi_{{\boldsymbol{\sigma}}}(\tau)}\\!=\\!\Bigg{[}$
$\displaystyle\sum_{x\geq
0}C_{\boldsymbol{\sigma}}(x,\tau)\hat{h}^{\dagger}_{x,0}\prod_{\begin{subarray}{c}j\in
S_{{\boldsymbol{\sigma}}}^{0}\\\ 0\leq j\leq
x\end{subarray}}\hat{s}^{\dagger}_{j-1,0}\prod_{\begin{subarray}{c}j\in
S_{{\boldsymbol{\sigma}}}^{0}\\\ j>x\end{subarray}}\hat{s}^{\dagger}_{j,0}$
$\displaystyle+$
$\displaystyle\sum_{x<0}C_{\boldsymbol{\sigma}}(x,\tau)\hat{h}^{\dagger}_{x,0}\prod_{\begin{subarray}{c}j\in
S_{{\boldsymbol{\sigma}}}^{0}\\\ x\leq
j<0\end{subarray}}\hat{s}^{\dagger}_{j+1,0}\prod_{\begin{subarray}{c}j\in
S_{{\boldsymbol{\sigma}}}^{0}\\\
j<x\end{subarray}}\hat{s}^{\dagger}_{j,0}\Bigg{]}\\!\prod_{y\neq 0}\prod_{j\in
S_{{\boldsymbol{\sigma}}}^{y}}\hat{s}^{\dagger}_{j,y}\ket{{\rm FM}}.$ (7)

The upper (lower) line describes the scenario in which the hole has moved
$|x|$ sites to the right (left), and how the spin excitations countermove by
one site to the left (right). Crucially, the probability amplitude to find the
hole at site $x$ and time $\tau$ for a given spin realization
${\boldsymbol{\sigma}}$ only depends on these three variables, since the spin
background is static apart from the countermotion due to the hole motion. This
also means that the probability to observe the hole at position $x$ after time
$\tau$ is the thermal average of $|C_{\boldsymbol{\sigma}}(x,\tau)|^{2}$,

$\displaystyle P(x,\tau)$ $\displaystyle={\rm
tr}\left[\hat{h}^{\dagger}_{x,0}\hat{h}_{x,0}\hat{\rho}(\tau)\right]=\sum_{\sigma_{0},{\boldsymbol{\sigma}}}\frac{e^{-\beta
E_{J}(\sigma_{0},{\boldsymbol{\sigma}})}}{Z_{J}}|C_{\boldsymbol{\sigma}}(x,\tau)|^{2}.$
(8)

In this way, we have to determine the probability amplitudes
$C_{\boldsymbol{\sigma}}(x,\tau)$ for a given spin realization
${\boldsymbol{\sigma}}$ and then perform sampling of the thermal average in
Eq. (8). To determine $C_{\boldsymbol{\sigma}}(x,\tau)$, we may realize, in
complete analogy to my previous paper on the two-leg ladder [1], that the
$C_{\boldsymbol{\sigma}}(x,\tau)$ amplitudes obey free-particle equations of
motion

$\displaystyle
i\partial_{\tau}C_{\boldsymbol{\sigma}}(x,\tau)=V_{{\boldsymbol{\sigma}}}(x)C_{\boldsymbol{\sigma}}(x,\tau)+t\left[C_{\boldsymbol{\sigma}}(x-1,\tau)+C_{\boldsymbol{\sigma}}(x+1,\tau)\right].$
(9)

with an emergent potential $V_{{\boldsymbol{\sigma}}}(x)$. In particular, the
potential may be divided into two parts:
$V_{{\boldsymbol{\sigma}}}=V_{{\boldsymbol{\sigma}},\parallel}+V_{{\boldsymbol{\sigma}},\perp}$.
The first term gives the contributions within the leg $y=0$ the hole is
propagating in,

$\displaystyle V_{{\boldsymbol{\sigma}},\parallel}(x)$
$\displaystyle=J[\sigma_{1,0}\sigma_{-1,0}-\sigma_{x,0}\sigma_{x+1,0}],\;x>0,$
$\displaystyle V_{{\boldsymbol{\sigma}},\parallel}(x)$
$\displaystyle=J[\sigma_{1,0}\sigma_{-1,0}-\sigma_{x,0}\sigma_{x-1,0}],\;x<0.$
(10)

Here $\sigma_{x,y}=\pm 1/2$ designates spin-$\uparrow$ $(+)$ and -$\downarrow$
$(-)$ at site ${\bf i}=x,y$ of the original sample, i.e. _before_ the hole has
started to move. The second term $V_{{\boldsymbol{\sigma}},\perp}$ describes
the trans-leg potential giving the contributions from the neighboring legs
$y=\pm 1$. This depends on the geometry of the couplings. For the square
lattice,

$\displaystyle V_{{\boldsymbol{\sigma}},\perp}(x)$ $\displaystyle=J\sum_{y=\pm
1}\sum_{j=+1}^{x}\sigma_{j,0}[\sigma_{j-1,y}-\sigma_{j,y}],\;x>0,$
$\displaystyle V_{{\boldsymbol{\sigma}},\perp}(x)$ $\displaystyle=J\sum_{y=\pm
1}\sum_{j=-1}^{x}\sigma_{j,0}[\sigma_{j+1,y}-\sigma_{j,y}],\;x<0.$ (11)

For the triangular lattice,

$\displaystyle V_{{\boldsymbol{\sigma}},\perp}(x)$
$\displaystyle=J\sum_{j=+1}^{x}\sigma_{j,0}\left\\{[\sigma_{j-1,+1}-\sigma_{j+1,+1}]+[\sigma_{j-2,-1}-\sigma_{j,-1}]\right\\},\;x>0,$
$\displaystyle V_{{\boldsymbol{\sigma}},\perp}(x)$
$\displaystyle=J\sum_{j=-1}^{x}\sigma_{j,0}\left\\{[\sigma_{j+2,+1}-\sigma_{j,+1}]+[\sigma_{j+1,-1}-\sigma_{j-1,-1}]\right\\},\;x<0.$
(12)

Here, the coordinates of the triangular lattices is defined by using the
standard embedding on a square lattice with an additional diagonal coupling
[25]. With these explicit constructions, the effective hole potential
$V_{\boldsymbol{\sigma}}(x)$ can easily be computed from each realized spin
sample ${\boldsymbol{\sigma}}$. Furthermore, for each ${\boldsymbol{\sigma}}$
Eq. (9) can be solved highly efficiently by defining the effective Hamiltonian
$\mathcal{H}_{\boldsymbol{\sigma}}$ with diagonal entries given by the
effective potential,
$\mathcal{H}_{\boldsymbol{\sigma}}(x,x)=V_{\boldsymbol{\sigma}}(x)$, and off-
diagonal entries given by the hopping
$\mathcal{H}_{\boldsymbol{\sigma}}(x,x\pm 1)=t$. By concatenating
$C_{\boldsymbol{\sigma}}(x,\tau)$ as a vector ${\bf
C}_{\boldsymbol{\sigma}}(\tau)$, the effective Schrödinger equation

$\displaystyle i\partial_{\tau}{\bf
C}_{\boldsymbol{\sigma}}(\tau)=\mathcal{H}_{\boldsymbol{\sigma}}{\bf
C}_{\boldsymbol{\sigma}}(\tau)$ (13)

can be solved using standard linear algebra packages. As
$\mathcal{H}_{\boldsymbol{\sigma}}$ is a sparse matrix, I use the
"expm_multiply" function, part of the "scipy.sparse.linalg" package in Python.
This allows me to go to systems sizes of at least 10.000 sites long.

To obtain the samples ${\boldsymbol{\sigma}}$ in the first place, I perform
Monte-Carlo sampling. In the presence of a phase transition, i.e. for FM
couplings in the triangular lattice as well as both FM and AFM couplings in
the square lattice, I use a combined Wolff [31, 32] and Metropolis-Hastings
[33, 34] algorithm to increase the accuracy around the critical temperature.
For AFM couplings in the triangular lattice, there is no phase transition, and
I instead simply use a standard Metropolis-Hasting algorithm with single-spin
flip updates. These algorithms are used to generate $2000$ samples for a range
of inverse temperatures $\beta J$. In this manner, the dynamics of the hole is
computed highly efficiently, whereby I easily go to very large system sizes,
long evolution times and effortlessly monitor the behavior across the phase
transition. The main observable in this regard will be the root-mean-square
(rms) distance, calculated as

$x_{\rm rms}(\tau)=\left[\sum_{x}x^{2}P(x,\tau)\right]^{1/2},$ (14)

evaluated from the hole probability distribution function $P(x,\tau)$ in Eq.
(8). This methodology is completely equivalent to the one I developed in Ref.
[1]. Only the explicit expressions for the trans-leg potentials in Eqs. (11)
and (12) differ, as well as the implementation of the Wolff algorithm. From
the rms dynamics, I extract a localization length as the long-time average

$l_{\rm loc}=\lim_{\tau\to\infty}\frac{1}{\tau}\int_{0}^{\tau}dsx_{\rm
rms}(s).$ (15)

Figure 2: Length scales associated with the effective hole potential for the
square (a) and triangular (b) lattices for ferromagnetic couplings, $J<0$, as
a function of the inverse temperature $\beta|J|$. The length scales are
defined through the linearity of the mean and variance of the potential:
$\braket{V_{\boldsymbol{\sigma}}(x)}=|J|\cdot|x|/{\color[rgb]{0.55078125,0.01171875,0.98828125}x_{\rm
ave}}$ and ${\rm
Var}[V_{\boldsymbol{\sigma}}(x)]=J^{2}\cdot|x|/{\color[rgb]{1,0.44921875,0}x_{\rm
fl}}$, respectively. The solid violet lines show the analytic solutions in Eq.
(17), while the black lines give the asymptotic behaviors, scaling as
$e^{3\beta|J|}$ and $e^{5\beta|J|}$ for the square and triangular cases. The
thin vertical black lines designate the position of the para- to ferromagnetic
phase transition at $\beta_{c}|J|=2\ln(1+\sqrt{2})$ and $\beta_{c}|J|=\ln(3)$,
respectively.

As a good check of the sampling, I compute the average and variance of the
effective hole potential. As for the two-leg ladder [1], these are both found
to be linear in the distance $|x|$,

$\displaystyle\braket{V_{\boldsymbol{\sigma}}(x)}=|J|\frac{|x|}{x_{\rm
ave}}+{\rm const.},\hskip 14.22636pt{\rm
Var}[V_{\boldsymbol{\sigma}}(x)]=J^{2}\frac{|x|}{x_{\rm fl}}+{\rm const.},$
(16)

with temperature-dependent length scales $x_{\rm ave}$ and $x_{\rm fl}$,
respectively. These are explicitly shown in Fig. 2. The former may also quite
easily be expressed in terms of short-range correlators as

$\displaystyle x_{\rm ave}^{\square}$
$\displaystyle=\frac{2}{C(1)-C(\sqrt{2})},\hskip 14.22636ptx_{\rm
ave}^{\triangle}=\frac{2}{C(1)-C(\sqrt{13/4})},$ (17)

for the square (left) and triangular (right) lattices, respectively. The
nearest-neighbor correlators
$C(1)=4\braket{\hat{S}^{(z)}_{0,0}\hat{S}^{(z)}_{1,0}}$, and next-nearest
neighbor correlators
$C(\sqrt{2})=C(\sqrt{2})=4\braket{\hat{S}_{0,0}\hat{S}_{1,1}}$ and
$C(\sqrt{13/4})=4\braket{\hat{S}_{1,0}\hat{S}_{0,1}}$ are computed explicitly
in Appendices A and B and are seen to agree perfectly with the numerical
results in Fig. 2.

##  4 Results

In this section, I present the results for the hole dynamics. The section is
split into three subsections, describing the universal infinite temperature
limit in Sec. 4.1, ferromagnetic couplings in Sec. 4.2, and finally
antiferromagnetic couplings in Sec. 4.3.

###  4.1 Infinite temperatures

For infinite temperatures, $\beta J=0$, the results are quite similar to the
two-leg ladder case [1]. In particular, since each spin is now an independent
random variable $\sigma=\pm 1/2$, it follows that the potential
$V_{\boldsymbol{\sigma}}(x)$ performs a random walk as a function of the
distance $|x|$. In particular, the mean value of the potential over all the
spin realization vanishes identically, while the standard deviation scales as
$\sqrt{|x|}$,

$\displaystyle\sigma[V_{\boldsymbol{\sigma}}(x)]=\frac{J}{2}\sqrt{|x|+1}.$
(18)

This gives a length scale $x_{\rm fl}=4$ at infinite temperatures in excellent
agreement with the Monte-Carlo sampling result shown in Fig. 2. As was also
realized for the two-leg ladder [1], the dynamics in this limit becomes
independent of the sign of $J$. This can be seen directly from the explicit
expressions in Eqs. (10)-(12). Here, a sign change in $J$ can be absorbed in
$\sigma_{j,0}=\pm 1/2$, as this random variable is independent from the rest
of the spins at infinite temperatures. More interestingly, since the
expressions for the square and triangular lattices, respectively Eqs. (11) and
(12), have the same number of terms with the same quadratic structure,
$\sigma\sigma^{\prime}$, they give identical potentials _at infinite
temperature_. Hence, the dynamics is not only independent of the sign of the
spin coupling in this limit, but also whether the geometry is square or
triangular. I note, however, that this does not hold for general 2D
structures. Indeed, for kagome and honeycomb lattices, one can find 1D lines,
in which the structure is the same as the two-leg ladder, and consequently
gives the same hole motion along these lines as in the two-leg ladder at
infinite temperatures.

Figure 3: Universal infinite temperature dynamics. (a) Root-mean-square (rms)
distance as a function of time $\tau$ in units of the hopping $t$ for
indicated spin couplings at infinite temperature. The dashed lines are the
asymptotic averages defining the localization length $l_{\rm loc}$. All lines
collapse at short times to an initial ballistic motion with speed $\sqrt{2}t$
(black line). (b) Localization length $l_{\rm loc}$ at infinite temperature as
a function of the spin coupling (red dots). For small $|J|/t$, the
localization length diverges as $(t/J)^{2}$ (red solid line). Inset: effective
hole potential $V(x)$ in units of the spin coupling for $50$ spin samples
(grey lines). As a hole with initial kinetic energy $\sim t$ travels in a
specific spin sample (colored lines), it will eventually back-scatter off the
potential (colored lines with arrows, shown for $t\simeq 5|J|$), because the
standard deviation of the potential grows as $|J|\sqrt{|x|+1}/2$ (black
lines), see also Eq. (18).

In Fig. 3(a), I plot the root-mean-square distance, calculated from Eq. (14),
versus time for indicated values of the spin coupling. As was found for the
two-leg ladder [1], the dynamics crosses over from an initial universal
ballistic behavior [35] of a free particle [with speed $\sqrt{2}t$], to
localized dynamics on long timescales. As the spin coupling is lowered, the
associated localization length, shown in Fig. 3(b), increases and eventually
diverges as $(t/J)^{2}$ in the limit of $|J|/t\ll 1$. As was realized for the
two-leg ladder, the localization can be understood by equating the initial
kinetic energy of the hole $\sim t$ to the fluctuations of the potential at a
length scale $l_{\rm fl}$, i.e. the standard deviation
$\sigma(V_{\boldsymbol{\sigma}}(l_{\rm fl}))\simeq|J|\sqrt{l_{\rm fl}}/2$.
This results in the fluctuation-induced localization length,

$l_{\rm fl}\simeq 4\left[\frac{t}{J}\right]^{2},$ (19)

which, apart from an overall factor of $2$, is in quantitative agreement with
the observed results in Fig. 3(b). The physical picture is, hereby, that the
dopant will eventually back-scatter off the emergent effective potential
$V(x)$ [inset of Fig. 3(b)]. Indeed, taking the enhancement factor of $2$ into
account, I check when $|V_{\boldsymbol{\sigma}}(x)|$ first exceeds $\sqrt{2}t$
for each spin realization ${\boldsymbol{\sigma}}$ and average over the
achieved mean-square distance $x^{2}$ over all the spin realizations. This
quantitatively recovers the full numerical solution both in terms of the
localization length $l_{\rm loc}$, and in terms of the standard deviation
around this value.

###  4.2 Ferromagnetic couplings

In this Section, I take a detailed look at the temperature dependency for
ferromagnetic couplings. In Figs. 4(a) and 4(b), I show the rms dynamics
across the phase transitions at the inverse Curie temperatures
$\beta_{c}^{\square}|J|=2\ln(1+\sqrt{2})\simeq 1.76$ and
$\beta_{c}^{\triangle}|J|=\ln(3)\simeq 1.10$, for the square and triangular
cases respectively. At short times, they again all collapse to a ballistic
expansion with speed $v=\sqrt{2}t$, as they should [35]. For lower
temperatures – higher $\beta|J|$ – the dynamics generally follows this
ballistic behavior for longer. Essentially, this is because the system becomes
more and more ferromagnetically ordered, allowing the hole to move more
freely. Note, however, that there are exceptions to this general rule. For
example, for $|J|=0.5t$, we see that the dynamics at $\beta|J|=1.2$ in the
square lattice and $\beta|J|=0.8$ in the triangular lattice is _more_
localized than at infinite temperature. We will return to this subtlety later
on.

Moreover, it is strikingly apparent that the asymptotes, i.e. the localization
lengths of the hole, remain _finite_ even in the phase with long-range
ferromagnetic order. This is surprising with the analysis of the two-leg
ladder in mind [1]. Here, it was shown that the localization length scales
with the spin-spin correlation length at low temperatures. As this length
scale diverges across the para- to ferromagnetic phase transition in the
present 2D system, the expectation from there would be that the hole should
also delocalize across the transition.

Figure 4: Ferromagnetic couplings. Root-mean-square distance (rms) dynamics
for $J<0$ and indicated inverse temperatures $\beta|J|$ in the square (left)
and triangular (right) lattices for $|J|=5t$ in (a) and $|J|=0.5t$ in (b). At
long times and any finite $\beta|J|$, the rms distance saturates at long
times. The associated localization length is plotted in (c) versus inverse
temperature for indicated values of the spin coupling for the square (left)
and triangular (right) lattices. The vertical black lines indicate the phase
transitions at $\beta_{c}^{\square}|J|=2\ln(1+\sqrt{2})$ (left) and
$\beta_{c}^{\triangle}|J|=\ln(3)$ (right). At low temperatures, the
localization length scales as $c_{1}\exp(c_{2}\beta|J|)$ (solid lines). The
growth rate $c_{2}$ is extracted and plotted in the insets of (c) as a
function of $|J|/t$. The horizontal black and grey lines show the expected
limiting behaviors for small and large values of $|J|/t$, corresponding to the
black and grey lines in the main panels.

To analyze this puzzling situation further, I next calculate the localization
length across the phase transitions for the square and triangular lattices in
Fig. 4(c). This manifestly shows that even though the localization has a sharp
increase around the phase transition, no divergence appears. Instead, I find
that the localization length scales as

$\displaystyle l_{\rm
loc}=c_{1}\left(\frac{|J|}{t}\right)\exp\left[c_{2}\left(\frac{|J|}{t}\right)\beta|J|\right],$
(20)

for low temperatures. Here, the coefficients $c_{1}$ and $c_{2}$ are functions
of $|J|/t$. The exponential growth rate, $c_{2}$, is plotted in the insets of
Fig. 4(c), and is seen to increase for decreasing $|J|/t$, between what seems
to be two limiting behaviors. This is again in stark contrast to the two-leg
ladder case [1]. Here, the localization length scales with the spin-spin
correlation length at low temperatures, which in turn increases exponentially
as $\exp[\beta|J|]$, i.e. with an exponential coefficient $c_{2}=1$
independent of $|J|/t$.

Let us, therefore, analyze these limits in detail. First, for a large mobility
of the hole, $|J|\ll t$, we can use a semi-classical argument of energetic
turning points. The basic idea is that the average and standard deviation of
the effective hole potential defines two length scales which compete to
localize the hole even as $|J|/t$ becomes small. In particular, I found in
Fig. 2 that the mean and standard deviation at all temperatures and
ferromagnetic couplings behave as
$\braket{V_{\boldsymbol{\sigma}}(x)}=|J||x|/x_{\rm ave}$ and
$\sigma(V_{\boldsymbol{\sigma}}(x))=|J|\sqrt{|x|/x_{\rm fl}}$. By equating
these to the initial kinetic energy of the hole $\sim t$, we obtain the semi-
classical turning points

$\displaystyle l_{\rm ave}$ $\displaystyle=\frac{t}{|J|}x_{\rm ave},\;l_{\rm
fl}=\left(\frac{t}{J}\right)^{2}x_{\rm fl},$ (21)

for the average and standard deviation of the potential, respectively.
Whichever of these two length scales is the shortest is expected to describe
the localization as $|J|/t$ becomes small – where it is well-defined to talk
about an initial kinetic energy of the hole. This actually also explains the
non-monotonic behavior of the localization length at high temperatures and low
$|J|/t$ mentioned previously. To understand why, note that $l_{\rm ave}$
diverges as the infinite temperature limit is approached, $\beta|J|\to 0$.
However, the fluctuation length scale $l_{\rm fl}$ remains finite as discussed
in Sec. 4.1. Now, at sufficiently low $|J|/t$, $l_{\rm ave}$ will eventually
drop _below_ the infinite temperature localization length going as $l_{\rm
fl}\propto(t/J)^{2}$, because it has a weaker $t/|J|$ scaling. Physically, the
potential achieves a confining bias, $\braket{V_{\boldsymbol{\sigma}}}(x)>0$,
which can localize the hole more strongly than the fluctutations at infinite
temperatures. Finally, since the hole delocalizes asymptotically as zero
temperature is reached in the ferromagnetic phase, the localization length is
non-monotonic for low enough $|J|/t$.

In the ferromagnetic phase, I found in Fig. 2 that the fluctuation length
scale $x_{\rm fl}$ has the same scaling behavior with decreasing temperature
as the average length scale $x_{\rm ave}$. This means that $l_{\rm ave}$ and
$l_{\rm fl}$ scale with temperature in the same manner, with exponential
growth rates of $c_{2}=3$ and $c_{2}=5$ for the square and triangular
lattices, respectively. This is seen to match the numerical findings in Fig.
4(c) in the regime of $|J|\ll t$.

Finally, we should understand why the scaling behavior is different for
intermediate to large $|J|/t$. Here, it is important to again stress the
difference with the two-leg ladder case [1]. There, as low temperatures are
approached, one also approaches the phase transition from the para- to
ferromagnetic phase. As a result, one can expect to only see a single length
scale appear in this limit: the spin-spin correlation length. However, here as
low temperatures are reached we do not approach a phase transition, because
the system is already in the ferromagnetic phase. Therefore, there can easily
be more than one length scale available. Indeed by analyzing the length scale
over which a single spin flip occurs – see Appendix C – I find that this on
average scales as

$l_{\rm flip}\propto\exp\left[\frac{z}{2}\beta|J|\right],$ (22)

for $z$ nearest neighbors. This gives _another_ length scale with growth rates
of $c_{2}=2$ and $c_{2}=3$ for the square and triangular lattices,
respectively. At these length scales, the effective hole potential will,
hereby, jump by $-|J|/2$ before jumping up again to $0$. When $|J|\gg t$, this
length scale will thus define the localization length, because the hole will
reflect back, as soon as it meets this energetic barrier. This explains why
$c_{2}\to 2$ for the square lattice and $c_{2}\to 3$ for the triangular
lattice, when $|J|/t$ becomes large. As $|J|/t$ diminishes, however, the hole
can start to tunnel through these barriers. Eventually, it can tunnel through
enough barriers so that the lower length scale is given by either $l_{\rm
ave}$ or $l_{\rm fl}$ in Eq. (21). This, consequently, explains the crossover
between the two behaviors.

Finally, it is worth pointing out that at intermediate values $|J|\sim 3t$,
the barriers are low enough that the hole can tunnel through many of them, but
still the localization length is exponentially small compared to what one
expects from Eq. (21). In this regime, it seems most accurate to think of the
localization in terms of Anderson localization in the presence of weak
disorder, i.e. when the disorder strength is smaller than or comparable to the
kinetic energy. In such a scenario, instead of simple back-reflection, a
particle localizes because it accumulates randomly varying phases for arriving
to a particular point [23], i.e.

$C(x)=c_{1}e^{i\varphi_{1}}+c_{2}e^{i\varphi_{2}}+\dots,$ (23)

in which the phases $\varphi_{1}$ are basically chosen at random. In the
present setup, these varying phases arise, because the distribution of the
spin flips, happening on average on the length scale of $l_{\rm flip}$, is
random. Indeed, the standard deviation on $l_{\rm flip}$ is on the same order
as $l_{\rm flip}$ itself, as shown in Appendix C. As a result, the hole will
travel wildly different length scales between each barrier. As the hole can
arrive between two barriers in many different ways, this gives rise to the
destructive interference in Eq. (23).

The asymptotic delocalization of the hole in the low temperature limit
describes a reversed metal-insulator _crossover_ , in which the system is
highly insulating at high temperatures, and becomes more and more metallic
once the phase transition to the ferromagnetic phase is crossed and zero
temperatures are approached.

###  4.3 Antiferromagnetic couplings

In this Section, we delve into the regime of antiferromagnetic couplings,
$J>0$. While ferromagnetic couplings led to qualitatively the same behavior
for the square and triangular lattices, we shall see that antiferromagnetic
couplings define highly distinct behaviors both for the underlying spin
lattice and for the dopant dynamics.

The square lattice is bipartite, and may therefore be divided into sublattices
A and B. This means that it is possible to rotate the local reference frame on
every second site, such that $\hat{S}^{(z)}_{\bf j}\to-\hat{S}^{(z)}_{\bf j}$
on sublattice B. As the spin couplings are between nearest neighbors only,
this also corresponds to flipping the sign of the spin coupling. This simple
analysis shows that the ferromagnetic and antiferromagnetic cases are
completely equivalent for a bipartite lattice. However, in the presence of
holes and, as here, nearest neighbor hopping of the spins onto such vacant
sites, the AFM and FM scenarios are no longer equivalent. Indeed, at low
temperatures the staggered magnetization appearing in the Neél ordered ground
state for AFM couplings gives rise to a _confined_ hole as it starts to move.
This naturally realizes an exact version of the retraceable path approximation
due to Brinkman and Rice [7]. The confimenent comes about, because the 1D
motion of the hole realigns spins that were otherwise antialigned, giving rise
to a linear potential increasing as $J/2\cdot|x|$ [24]. The crossover from the
thermally induced localization at high temperatures to the confined hole
motion at low temperatures is illustrated in Fig. 5(left). Where the dynamics
at high temperatures is mostly featureless, the motion at low temperatures is
characterized by strong coherent oscillations. These have previously been
shown to be due to interferences between the so-called string states that
define the low-energy eigenstates at zero temperature [24]. Moreover, in Fig.
5(c), we see that the approach to the zero-temperature limit happens as the
entropy of the spin lattice [inset in Fig. 5(c)] approaches $0$, around $\beta
J\gtrsim 3$. The sharp decrease around the critical temperature behavior
$\beta_{c}^{\square}J=2\ln(1+\sqrt{2})$ originates in this sense directly from
the sharp decreasing behavior in the entropy at the phase transition.

Figure 5: Antiferromagnetic couplings. Root-mean-square (rms) distance of hole
to its original site as a function of time $\tau$ in units of hopping $t$ for
$|J|=5t$ (a) and $|J|=0.5t$ (b) for indicated values of the inverse
temperature $\beta$. For the square lattice (left), lower temperatures –
higher $\beta|J|$ – results in more and more pronounced _coherent_
oscillations eventually approaching the zero-temperature behavior in grey
lines. For the triangular lattice (right), on the other hand, the dynamics
depend only mildly on temperature, especially at larger spin couplings as in
(a), and retains a thermal character even as zero temperature is approached.
(c) Localization length versus inverse temperature for indicated values of the
spin coupling for the square (left) and triangular (right) lattices. In the
square lattice (right), the localization length decreases rapidly across the
phase transition at $\beta_{c}^{\square}J=2\ln(1+\sqrt{2})\simeq 1.76$
(vertical black line). In the triangular lattice (right), there is no phase
transition and the localization length consequently has a much slower
dependency on temperature. The insets in (c) show the entropy per particle
$S/N$, both starting out at $S/N=k_{B}\ln(2)$ at high temperatures. The
saturation at low temperatures happens as the entropy of the system approaches
its zero-temperature limit ($0$ for the square lattice, $\simeq 0.323$ for the
triangular lattice).

The triangular lattice is, however, markedly different. In this case, the
system is frustrated and there is no mapping between the ferro- and
antiferromagnetic cases. Even more dramatically, the frustration of the
lattice leads to a non-vanishing entropy at zero temperature [25], shown in
the inset of the right figure in Fig. 5(c). This strongly affects the dynamics
as temperature is lowered. As mentioned in Sec. 4.1, the high-temperature
limit gives the same dynamics both for ferro- and antiferromagnetic
interactions _and_ for the square and triangular geometries. However, as
temperature is lowered there is, in stark contrast to the square lattice, no
appearance of strong oscillations, as can be seen to the right in Figs. 5(a)
and 5(b). The reason is that the ground state degeneracy is exponentially
large in system size, meaning that the dynamics is averaged over many
different spin realizations, even at the lowest temperatures, and this washes
out the coherent oscillations that would otherwise appear. Moreover, the lack
of a phase transition means that the localization length changes much slower
with temperature as seen from the figure to the right in Fig. 5(c). In point
of fact, we need to wait for the entropy per particle to be close to its zero-
temperature limit, $S_{0}/N\simeq 0.323k_{B}$, for the localization length of
the hole to saturate to its zero-temperature limit. This only happens for
inverse temperatures $\beta J\gtrsim 6$ for the triangular case.

What is also initially confounding about the triangular lattice is that there
are ground state configurations, in which the effective hole potential is
completely flat. Such configurations are all variations of the Neél states
shown in Fig. 6(a). So if such configurations are there at low temperatures,
why is the hole even localized? Surely, their presence must make it possible
for the hole to escape its origin – even ballistically fast. Well, the answer
to this conundrum turns out again to lie in the entropy of such states. As
originally pointed out by Wannier [25], there are on the order of
$2^{\sqrt{N}}$ states with the structure in Fig. 6(a). This scaling comes from
realizing that the only alteration one can make to such a state is to shift
the rows of the lattice by 1. And since there are $\sqrt{N}$ rows this gives
$2^{\sqrt{N}}$ states. However, there is a much much larger family of states
shown in Fig. 6(b). Here all the purple spins, which is every third, can be
either spin-$\uparrow$ or -$\downarrow$ without a change in the energy. There
are, therefore, at least a staggering $2^{N/3}$ of these111Wannier [25] came
up with an even stronger lower bound of $2^{5N/12}$ based on simple
geometrical arguments. This gives a ground state entropy
$>5\ln(2)/12k_{B}N\simeq 0.289k_{B}N$, pretty close to the exact value of
$S_{0}\simeq 0.323k_{B}N$.. This also explains why the entropy of the ground
state manifold is nonzero. For just $N=100$ spins, the relative abundance of
the latter type of states to the former type is $2^{N/3}/2^{\sqrt{N}}\sim
10^{7}$, for $400$ spins the ratio is at $\sim 10^{34}$! As a result, even
though there are _in principle_ states available in the ground state manifold
in which the hole could delocalize, they have _zero_ statistical weight.

Figure 6: Origin of localization for triangular lattice at low temperatures.
(a) A perfect Neél structure is a part of the ground state manifold. As the
hole moves through such a state, it experiences a completely flat potential,
$V(x)=1/2$ for $x\neq 0$, since the perpendicular part vanishes
$V_{\perp}(x)=0$. This suggests that the hole should be able to delocalize for
antiferromagnetic couplings. However, the number of perfect Neél ordered
states scales only as $2^{\sqrt{N}}$, and there are configurations (b) that
have a _much_ larger weight. Here, all purple spins can be chosen freely
between $\ket{\uparrow},\ket{\downarrow}$, leading to at least $2^{N/3}$
states. For these, the hole always experiences an overall growing potential,
as it is guaranteed to increase for every third hop. As the latter type in (b)
completely outnumbers the first type in (a), the hole remains localized.

###  4.4 Spin coupling scaling dependency

Before concluding, this Section is concerned with describing in detail the
dependency of the localization length on the ratio between the spin coupling
and the hopping amplitude, $J/t$. Examples of this dependency are shown in
Figs. 7(a) and 7(b) for ferro- and antiferromagnetic couplings, respectively.
This is furthermore compared to the universal behavior found at infinite
temperatures, identical for the square and triangular lattices and for both
ferro- and antiferromagnetic spin couplings, scaling as $(t/J)^{2}$ for
$|J|/t\ll 1$. Analogous to the two-leg ladder [1], this scaling behavior turns
out to be highly specific to the infinite-temperature limit. In point of fact,
for any finite temperature we observe that the asymptotic behavior is rather
$t/|J|$ with a temperature dependent prefactor. This is because as $t/|J|$ is
lowered, the bias of the effective hole potential will eventually dominate
over the fluctuations, i.e. $l_{\rm ave}\sim t/|J|<l_{\rm fl}\sim(t/J)^{2}$
for sufficiently small $|J|/t$, no matter the temperature.

Moreover, we see that for ferromagnetic couplings in Fig. 7(a), equal values
of $\beta/\beta_{c}^{\square}$ and $\beta/\beta_{c}^{\triangle}$ – i.e.
temperatures in the same proportion to their respective critical temperatures
– are quantitatively very similar, especially for low values of $|J|/t$. On
the antiferromagnetic side in Fig. 7(b), we additionally observe that for the
square lattice, the localization length already at the critical temperature is
quantitatively close to the zero-temperature limit for low $J/t$, with a
factor of $1.5$ between them. For larger $J/t$, the localization length is
short, and it becomes important that the system locally has occasional spin
flips with respect to the perfect Neél-ordered state at zero temperature.
Finally, for hole motion in the triangular lattice the effective hole
potential in the ground state manifold only increases in every third hop, as
described in Fig. 6. As a result, the localization length is strictly larger
for the hole in the triangular lattice compared to the square lattice at
similar temperatures, here shown for $\beta J=1.8\simeq\beta_{c}^{\square}J$.
Also, at these temperatures the localization length of the hole in the
triangular lattice even follows the infinite temperature behavior for
$J\gtrsim 2t$.

Figure 7: Localization length dependency for ferromagnetic (a) and
antiferromagnetic (b) couplings as a function of $J/t$ on a log-log plot for
indicated inverse temperatures. Square and triangular symbols are for the
square and triangular lattices, respectively. Also,
$\beta_{c}^{\square}=2\ln(1+\sqrt{2})/|J|$ and
$\beta_{c}^{\triangle}=\ln(3)/|J|$ indicate the inverse critical temperatures
in the square and triangular lattice, respectively. Only at infinite
temperature is the asymptotic scaling $(t/J)^{2}$ (solid red line). For any
finite temperature, the asymptotic behavior (dashed lines) is $t/|J|$ with a
temperature dependent prefactor.

##  5 Discussion

In this Section, I will discuss generalizations to the mixed-dimensional
$t$-$J_{z}$ model considered in this Article.

###  5.1 Beyond 1D dopant motion

A central assumption in the present Article has been the one-dimensional
nature of the dopant motion. This not only enables the numerical computation
of the results to large system sizes and long times, but also clarifies
certain physical situations that are far more complex in higher dimensions.
For example, when the hole moves solely in 1D it has to scatter on any domain
walls it may find on its way. If the hole were to move in 2D, it could very
well be able to circumvent these, enabling a slow, but continuing,
propagation. In the ferromagnetic phase, where the system macroscopically
occupies one spin state, say spin-$\ket{\uparrow}$, the domains of
spin-$\ket{\downarrow}$ become smaller and rarer, and eventually consist of
singly flipped spins as discussed in Sec. 4.2. In this low-temperature
ferromagnetic limit, the system hereby resembles a perfect lattice with
occasional "impurities" of spin-$\ket{\downarrow}$ that the hole may scatter
on. This seems essentially equivalent to the propagation of electrons in
lattices with a low density of defects – or local impurities. As a result, in
the ferromagnetic phase we should at the very least expect the hole to perform
diffusion with _weak localization_ corrections [36], if not full-scale
ballistic motion.

Contrary to the usual case in the solid state, however, the effective defect
density – the domain walls – in the present setup increases rapidly with
temperature. As a result, the mean free path decreases to be on on the order
of the lattice spacing at high temperatures. In this limit, therefore,
scattering happens all the time and the dopant may very well localize
completely. These considerations illustrate that the inverted metal-insulator
_crossover_ found in the present analysis for 1D dopant motion may be replaced
by an actual _transition_ for 2D motion, and should probably happen at the
underlying Curie temperature of the spin lattice.

Detailed numerical analyses in the infinite temperature limit of vanishing
spin couplings for 2D motion of dopants [2, 3] illustrate that the situation
even in this restricted setup is quite complex. At the short to intermediate
timescales investigated, the motion is a lot slower than ballistic motion, but
remain faster than the predictions of the retraceable path approximation [2,
7], and faster than expected from the mapping to a Bethe lattice [3].
Conclusions on long timescales, however, remain to be drawn and must await
further analyses.

###  5.2 Beyond the Ising model

Another crucial assumption in the presented analysis is the Ising type spin
couplings. One could go beyond this realm by introducing spin flip-flop terms
into the Hamiltonian. The computational complexity in this case, however,
increases in a daunting manner, as the description even just of the underlying
spin lattice now becomes highly complex and – at best – approximate. At low
temperatures and antiferromagnetic couplings, there is good evidence that
quasiparticles – magnetic polarons – form [9]. Indeed, linear spin-wave theory
combined with the selfconsistent Born approximation has been shown to compare
well to exact diagonalization studies [10, 37] and Monte-Carlo simulations
[38]. Even more importantly, it has successfully explained [35] the
experimentally observed propagation of holes [16] through the formation and
propagation of such magnetic polarons. As temperature is increased, however,
there are currently only limited approaches available [39], and the
quasiparticle picture at low temperatures is even still debated [40, 41, 42,
43, 44, 45].

One possible path forward could be to go in the opposite regime of the XY
model, where the _only_ present spin couplings are flip-flop terms. In one
dimension, this supports a simple analytical solution via the Jordan-Wigner
transformation [46]. If one can also formulate the motion of dopants in an
efficient manner in this model, this would certainly be an interesting pathway
to pursue.

###  5.3 Increasing the doping level

In the present analysis, I have focused on the propagation of a single dopant.
One may wonder, what happens as the doping level rises. If their initial
mutual distances are much larger than the single dopant localization length,
one can expect them to retain their single-dopant characteristics uncovered in
the present analysis. However, if two dopants start out close to each another,
they may significantly alter each others motion. This could lead to novel
phenomena, as their motion may become strongly correlated even though they are
submerged in, e.g., an infinite temperature spin environment.

###  5.4 Connecting to an external bath

A final obvious extension of the ideas pursued in the present Article is to
investigate the same type of dynamics, but in the presence of a coupling to an
external heat bath. This should drive the system back to thermal equilibrium
by allowing spin flip dynamics to occur, and may be a way to investigate their
influence in a simpler manner than through spin flip-flop terms in the
Hamiltonian. It is clear, however, that the methodology must be changed
substantially in such an open system case, as now the propagation dynamics
cannot be calculated as the thermal – Boltzmann-weighted – average of pure
state evolutions.

##  6 Conclusions and outlook

In this Article, I have investigated the one-dimensional motion of a dopant in
two-dimensional square and triangular lattices of Ising coupled spins. The
thermally induced localization effect originally found in a two-leg ladder
geometry [1] has, hereby, been extended to the case in which there exists a
finite temperature transition to a long-range ordered ferromagnetic phase.
While the high-temperature limit features universal localized hole dynamics
across ferro- and antiferromagnetic couplings as well as the two investigated
lattice geometries, finite temperatures break this correspondence. On the
ferromagnetic side, the hole remains localized across the Curie temperature
and feature very similar behaviors for the two investigated geometries. While
the localization in the two-leg ladder was found to scale identically with
temperature across any (negative) value of $J/t$, this is no longer true for
the two-dimensional system. At small $|J|/t$, the localization can be
understood as the back-scattering off an effective hole potential that
fluctuates to large values as the hole moves away from its origin. At large
$|J|/t$, the hole instead back-reflects on singular spin flips happening on an
exponentially shorter length scale. Inbetween, increased tunneling through
these singular spin flips describes a crossover from one to the other behavior
as $|J|/t$ is lowered.

On the antiferromagnetic side, only the square lattice features a phase
transition to a long-range ordered AFM phase. As the antiferromagnetic
correlations grow, the hole experiences a stronger and stronger linear
potential, because its motion now starts to realign spins that were otherwise
antialigned. This leads to a crossover between (thermal) disorder induced
localization to confinement, with a particularly sharp decrease of the
localization length around the Curie temperature in the square lattice. The
characteristics of the dynamics in these two regimes is also markedly
different. At high temperatures, the motion is incoherent with a smooth,
featureless behavior of the root-mean-square distance of the hole to its
origin. As zero temperature is approached, and the entropy of the system
vanishes, stronger coherent oscillation occur due to quantum interference of
the low-lying energy states, the so-called string states [47, 24].

In the triangular lattice, no phase transition happens for antiferromagnetic
couplings due to frustration. This makes the approach to the zero-temperature
limit much slower, explained well by when the entropy drops to its nonzero
zero-temperature limit [25]. The associated exponentially large ground state
manifold means that the hole dynamics retains its featureless, thermal
behavior.

These detailed investigations show that the motion of dopants, even in these
highly simplistic models, host rich and diverse behaviors, in which an indepth
knowledge of the underlying spin lattice is crucial for understanding the
dopant dynamics. For ferromagnetc couplings, it describes an intriguing
_reversed_ metal-insulator crossover, from a highly insulating regime at high
temperatures to an increasingly metallic regime at low temperatures. Moreover,
there are several exciting research pathways that may be undertaken to expand
these considerations. First and foremost, it would be interesting to study the
stability of the localization effect. Here, one could study the influence of
coupling the spins to an external heat bath driving them towards
thermalization. In such a scenario, the thermal fluctuations might play a
similar role to flip-flop spin interactions, and such links could be pursued
further. One could also introduce such flip-flop terms in the Hamiltonian
directly, and finally one could pursue the understanding of less restrained
hole motion, where it is allowed to move not only along a one-dimensional line
in the lattice.

## Acknowledgements

The author thanks J. Ignacio Cirac, Pavel Kos, Dominik S. Wild, and Marton
Kanasz-Nagy for valuable discussions.

#### Funding information

This article was supported by the Carlsberg Foundation through a Carlsberg
Internationalisation Fellowship, grant number CF21_0410.

##  Appendix A Short-range correlators: square lattice

In this Appendix, I compute the nearest and next-nearest spin correlators for
the square lattice. The calculation is based Ref. [48]. Starting from the
Hamiltonian

${\hat{H}}_{J}=-|J|\sum_{\braket{{\bf i},{\bf j}}}\hat{S}^{(z)}_{\bf
i}\hat{S}^{(z)}_{\bf j}=-\frac{|J|}{4}\sum_{\braket{{\bf i},{\bf
j}}}\hat{s}_{\bf i}\hat{s}_{\bf j},$ (A.1)

we can express the desired correlators as

$\displaystyle
C(1)=4\braket{\hat{S}^{(z)}_{0,0}\hat{S}^{(z)}_{1,0}}=\braket{\hat{s}_{0,0}\hat{s}_{1,0}},\;C(\sqrt{2})=4\braket{\hat{S}^{(z)}_{0,0}\hat{S}^{(z)}_{1,1}}=\braket{\hat{s}_{0,0}\hat{s}_{1,1}}.$
(A.2)

Here, I define $\hat{s}=2\hat{S}^{(z)}$, such that it can take on the values
$\pm 1$. The nearest and next-nearest neighbor correlators are computed from
the expression

$a_{0}(\alpha_{1},\alpha_{2})=\int_{0}^{2\pi}\frac{d\theta}{2\pi}\left[\frac{(1-\alpha_{1}e^{i\theta})(1-\alpha_{2}e^{-i\theta})}{(1-\alpha_{1}e^{-i\theta})(1-\alpha_{2}e^{+i\theta})}\right]^{1/2}.$
(A.3)

Here, the only difference between the two correlators are in the choice of the
$\alpha_{i}$. Explicitly,

$\displaystyle C(1)$
$\displaystyle:\alpha_{1}=e^{-\beta|J|/2}\tanh\left(\frac{\beta|J|}{4}\right),\;\alpha_{2}=e^{-\beta|J|/2}{\rm
coth}\left(\frac{\beta|J|}{4}\right),$ $\displaystyle C(\sqrt{2})$
$\displaystyle:\alpha_{1}=0,\;\alpha_{2}=\frac{1}{\sinh^{2}\left(\frac{\beta|J|}{2}\right)}.$
(A.4)

This allows me to numerically compute these correlators. Also, a rather
tedious low-temperature expansion shows that

$\displaystyle C(1)$ $\displaystyle\to
1-4e^{-2\beta|J|}-\frac{47}{4}e^{-3\beta|J|},$ $\displaystyle C(\sqrt{2})$
$\displaystyle\to 1-4e^{-2\beta|J|}-16e^{-3\beta|J|}.$ (A.5)

The length scale arising from the mean value of the hole potential then
asymptotically scales as

$x_{\rm
ave}=\frac{2}{C(1)-C(\sqrt{2})}\to\frac{8}{64-47}e^{+3\beta|J|}=\frac{8}{17}e^{+3\beta|J|},$
(A.6)

having a fast $e^{+3\beta|J|}$ scaling behavior.

##  Appendix B Short-range correlators: triangular lattice

In this Appendix, I compute the nearest and next-nearest neighbor spin
correlators for ferromagnetic couplings in the triangular lattice. The
calculation is based on Refs. [49, 50]. The setup is identical to the one in
the previous Appendix, albeit with the diagonal coupling appropriate for the
triangular lattice. I also define $v=\tanh(\beta|J|/4)$. I need the nearest
and next-nearest neighbor correlators (at distance $1$ and $\sqrt{13/4}$)

$\displaystyle
C(1)=4\braket{\hat{S}^{(z)}_{0,0}\hat{S}^{(z)}_{1,0}}=\braket{\hat{s}_{0,0}\hat{s}_{1,0}},\;C(\sqrt{13/4})=4\braket{\hat{S}^{(z)}_{1,0}\hat{S}^{(z)}_{0,1}}=\braket{\hat{s}_{1,0}\hat{s}_{0,1}}.$
(B.1)

The nearest-neighbor correlator is explicitly computed in Ref. [49] to be

$C(1)=\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\left[\frac{a-be^{+i\omega}-ce^{-i\omega}}{a-be^{-i\omega}-ce^{+i\omega}}\right]^{1/2},$
(B.2)

with

$\displaystyle a=2v(1+v^{2}),\;b=v^{2}c=v^{2}(1-v)^{2}.$ (B.3)

The next-nearest correlator $C(\sqrt{13/4})$ does not seem to be explicitly
computed in Stephenson’s papers. However, we may use the results for 4-point
correlators to get $C(\sqrt{13/4})$. In particular, slightly rewritting Eq.
(2.15) in Ref. [50], I get

$\displaystyle\braket{\hat{s}_{0,0}\hat{s}_{1,0}\hat{s}_{p,q}\hat{s}_{p,q+1}}=$
$\displaystyle C(1)^{2}+(1+v)^{2}(1-v)^{2}$
$\displaystyle\times\left\\{[p-1,q]_{4,3}[p,q+1]_{1,6}-[p,q]_{1,3}[p-1,q+1]_{4,6}\right\\}.$
(B.4)

Here, the notation $[p,q]$ is short-hand for the $6\times 6$ matrix

$\displaystyle[p,q]=A^{-1}(p,q)=\int_{-\pi}^{\pi}\frac{d\varphi_{1}}{2\pi}\int_{-\pi}^{\pi}\frac{d\varphi_{2}}{2\pi}e^{-i(p\varphi_{1}+q\varphi_{2})}A^{-1}(\phi_{1},\phi_{2}).$
(B.5)

Here, $A(\phi_{1},\phi_{2})$ is a specific $6\times 6$ matrix depending on
$v,\varphi_{1},\varphi_{2}$, which I will return to in a moment. Inserting
$p=q=0$ in Eq. (B.4), the 4-point correlator collapses to
$\braket{\hat{s}_{0,0}\hat{s}_{1,0}\hat{s}_{0,0}\hat{s}_{0,1}}$
$=\braket{\hat{s}_{1,0}\hat{s}_{0,1}}=C(\sqrt{2})$, as the $\hat{s}$ operators
commute and $\hat{s}_{p,q}^{2}=1$. In this manner,

$\displaystyle
C(\sqrt{13/4})=C(1)^{2}+(1+v)^{2}(1-v)^{2}\left\\{[-1,0]_{4,3}[0,+1]_{1,6}-[0,0]_{1,3}[-1,+1]_{4,6}\right\\}.$
(B.6)

The matrix that we need to invert is given in Eq. (2.4) in Ref. [50]

$\displaystyle
A(\varphi_{1},\varphi_{2})=\begin{bmatrix}0&1&1&1-v\,e^{i\varphi_{1}}&1&1\\\
-1&0&1&1&1-v\,e^{i(\varphi_{1}+\varphi_{2})}&1\\\
-1&-1&0&1&1&1-v\,e^{i\varphi_{2}}\\\ -1+v\,e^{-i\varphi_{1}}&-1&-1&0&1&1\\\
-1&-1+v\,e^{-i(\varphi_{1}+\varphi_{2})}&-1&-1&0&1\\\
-1&-1&-1+v\,e^{-i\varphi_{2}}&-1&-1&0\end{bmatrix}.$ (B.7)

Here, I perform the inversion of the matrix in Mathematica. I express the
result in terms of the cofactor matrix $C$: $A^{-1}=C^{T}/\Delta$, where
$C^{T}$ is the transpose of $C$ and $\Delta=\det(A)$ is the determinant.
Explicitly,

$\displaystyle\Delta=(1+v)^{2}\left[\frac{(1+v^{2})^{3}+8v^{3}}{(1+v)^{2}}-2v(1-v)^{2}\left\\{\cos(\varphi_{1})+\cos(\varphi_{2})+\cos(\varphi_{1}+\varphi_{2})\right\\}\right]$
(B.8)

Transforming the variables as
$\theta=-\varphi_{2},\omega=\varphi_{1}+\varphi_{2}$, one can express the
determinant in the form

$\displaystyle\Delta(\theta,\omega)=(1+v)^{2}\left[A+B\cos(\theta)+C\sin(\theta)\right].$
(B.9)

Here,

$\displaystyle A$
$\displaystyle=\frac{(1+v^{2})^{3}+8v^{3}}{(1+v)^{2}}-2v(1-v)^{2}\cos(\omega),$
(B.10) $\displaystyle B$
$\displaystyle=-2v(1-v)^{2}\left[1+\cos(\omega)\right],$ (B.11) $\displaystyle
C$ $\displaystyle=2v(1-v)^{2}\sin(\omega).$ (B.12)

Moreover, we need combonents $C_{3,4},C_{6,1},C_{3,1}$ and $C_{6,4}$ to
compute the correlator. First,

$\displaystyle C_{3,4}(\varphi_{1},\varphi_{2})$
$\displaystyle=(1+v)\left\\{\left[1-v(1-2v)\right]-v\left[2-v(1-v)\right]e^{-i(\varphi_{1}+\varphi_{2})}-v(1-v)\left[e^{-i\varphi_{1}}+e^{-i\varphi_{2}}\right]\right\\}$
$\displaystyle=(1+v)\left\\{\left[1-v(1-2v)\right]-v\left[2-v(1-v)\right]e^{-i\omega}-v(1-v)\left[e^{-i(\omega+\theta)}+e^{i\theta}\right]\right\\}$
$\displaystyle=(1+v)\left\\{D+Ee^{-i\omega}+F\left[e^{-i(\omega+\theta)}+e^{i\theta}\right]\right\\},$
(B.13)

with

$\displaystyle D$
$\displaystyle=\left[1-v(1-2v)\right],\;E=-v\left[2-v(1-v)\right],\;F=-v(1-v).$
(B.14)

Second,

$\displaystyle C_{6,1}(\varphi_{1},\varphi_{2})$
$\displaystyle=-(1+v)\left\\{D+Ee^{+i(\varphi_{1}+\varphi_{2})}+F\left[e^{+i\varphi_{1}}+e^{+i\varphi_{2}}\right]\right\\}=-C_{3,4}^{*}(\varphi_{1},\varphi_{2})$
$\displaystyle=-(1+v)\left\\{D+Ee^{+i\omega}+F\left[e^{+i(\omega+\theta)}+e^{-i\theta}\right]\right\\}.$
(B.15)

Third,

$\displaystyle C_{3,1}(\varphi_{1},\varphi_{2})$
$\displaystyle=-(1+v)\left\\{-(1-v)+v^{2}(1-v)e^{+i(\varphi_{1}-\varphi_{2})}+v(1+v)\left[e^{+i\varphi_{1}}+e^{-i\varphi_{2}}\right]\right\\}$
$\displaystyle=-(1+v)\left\\{-(1-v)+v^{2}(1-v)e^{+i(\omega+2\theta)}+v(1+v)e^{i\theta}\left[e^{+i\omega}+1\right]\right\\}.$
(B.16)

And finally,

$\displaystyle C_{6,4}(\varphi_{1},\varphi_{2})$
$\displaystyle=-(1+v)\left\\{-(1-v)+v^{2}(1-v)e^{-i(\varphi_{1}-\varphi_{2})}+v(1+v)\left[e^{-i\varphi_{1}}+e^{+i\varphi_{2}}\right]\right\\}=C_{3,1}^{*}(\varphi_{1},\varphi_{2})$
$\displaystyle=-(1+v)\left\\{-(1-v)+v^{2}(1-v)e^{-i(\omega+2\theta)}+v(1+v)e^{-i\theta}\left[e^{-i\omega}+1\right]\right\\}.$
(B.17)

Now, I will the terms in Eq. (B.6) explicitly. First,

$\displaystyle[-1,0]_{4,3}=\int_{-\pi}^{\pi}\frac{d\varphi_{1}}{2\pi}\int_{-\pi}^{\pi}\frac{d\varphi_{2}}{2\pi}e^{i\varphi_{1}}\frac{C_{3,4}(\varphi_{1},\varphi_{2})}{\Delta(\varphi_{1},\varphi_{2})}=\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}e^{i(\omega+\theta)}\frac{C_{3,4}(\omega,\theta)}{\Delta(\omega,\theta)}$
$\displaystyle=\frac{1}{1+v}\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\frac{F+(De^{i\omega}+E)e^{i\theta}+Fe^{i\omega}e^{2i\theta}}{A+B\cos(\theta)+C\sin(\theta)}$
$\displaystyle=\frac{1}{1+v}\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\left[FI_{0}(\omega)+(De^{i\omega}+E)I_{1}(\omega)+Fe^{i\omega}I_{2}(\omega)\right].$
(B.18)

Here, I follow Stephenson [49] and define

$I_{n}(\omega)=\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\frac{e^{in\theta}}{A+B\cos(\theta)+C\sin(\theta)}=\frac{\alpha^{n}}{(A^{2}-B^{2}-C^{2})^{1/2}}.$
(B.19)

Importantly, this integral is solved exactly. Here,
$\alpha=[(A^{2}-B^{2}-C^{2})^{1/2}-A]/[B-iC]$. I, furthermore, define
$I_{n,m}=\int_{-\pi}^{\pi}d\omega{\rm e}^{im\omega}I_{n}(\omega)/(2\pi)$. Then

$\displaystyle[-1,0]_{4,3}$
$\displaystyle=\frac{1}{1+v}\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\left[FI_{0}(\omega)+(De^{i\omega}+E)I_{1}(\omega)+Fe^{i\omega}I_{2}(\omega)\right]$
$\displaystyle=\frac{1}{1+v}\left[DI_{1,1}+EI_{1,0}+F(I_{0,0}+I_{2,1})\right].$
(B.20)

Likewise,

$\displaystyle[0,1]_{1,6}=\int_{-\pi}^{\pi}\frac{d\varphi_{1}}{2\pi}\int_{-\pi}^{\pi}\frac{d\varphi_{2}}{2\pi}e^{-i\varphi_{2}}\frac{C_{6,1}(\varphi_{1},\varphi_{2})}{\Delta(\varphi_{1},\varphi_{2})}=\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}e^{i\theta}\frac{C_{6,1}(\omega,\theta)}{\Delta(\omega,\theta)}$
$\displaystyle=-\frac{1}{1+v}\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\frac{F+(D+Ee^{i\omega})e^{i\theta}+Fe^{i\omega}e^{2i\theta}}{A+B\cos(\theta)+C\sin(\theta)}$
$\displaystyle=-\frac{1}{1+v}\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\left[FI_{0}(\omega)+(D+Ee^{i\omega})I_{1}(\omega)+Fe^{i\omega}I_{2}(\omega)\right]$
$\displaystyle=-\frac{1}{1+v}\left[DI_{1,0}+EI_{1,1}+F(I_{0,0}+I_{2,1})\right]$
(B.21)

which is very similar to $[-1,0]_{4,3}$. Moreover,

$\displaystyle[0,0]_{1,3}=\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\frac{C_{3,1}(\omega,\theta)}{\Delta(\omega,\theta)}$
$\displaystyle=-\frac{1}{1+v}\left[-(1-v)I_{0,0}+v^{2}(1-v)I_{2,1}+v(1+v)(I_{1,1}+I_{1,0})\right].$
(B.22)

And finally,

$\displaystyle[-1,+1]_{4,6}=\int_{-\pi}^{\pi}\frac{d\omega}{2\pi}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}e^{i(\omega+2\theta)}\frac{C_{6,4}(\omega,\theta)}{\Delta(\omega,\theta)}$
$\displaystyle=-\frac{1}{1+v}\left[-(1-v)I_{2,1}+v^{2}(1-v)I_{0,0}+v(1+v)(I_{1,1}+I_{1,0})\right].$
(B.23)

So, to compute $C(\sqrt{13/4})$, I need to compute the four terms
$I_{0,0},I_{1,0},I_{1,1},I_{2,1}$.

##  Appendix C Localization length in the ferromagnetic phase for large
$|J|/t$

In this Appendix, I derive the average distance between spin flips in the
ferromagnetic phase. I then use this to calculate what the asymptotic
localization length is at low temperatures and large $|J|/t$.

First, in a lattice with $z$ nearest neighbor interactions, the probability to
have a single spin flip is proportional to the Boltzmann factor $p_{\rm
flip}=e^{-z\beta|J|/2}$. Here $z=4,6$ correspond to the square and triangular
lattices, respectively. This means that at low temperatures, the probability
to find a strip of length $l\geq 1$ with exactly one spin flip at the end is
proportional to $p_{\rm flip}[1-p_{\rm flip}]^{l-1}$. Since the normalization
constant is simply unity, $A=\sum_{l=1}^{\infty}p_{\rm flip}[1-p_{\rm
flip}]^{l-1}$ $=p_{\rm flip}\sum_{l=0}^{\infty}[1-p_{\rm flip}]^{l}=1$, the
probability to find such a segment of length $l$ is

$P(l)=p_{\rm flip}[1-p_{\rm flip}]^{l-1}.$ (C.1)

The average length of such a segment gives the mean distance between spin
flips (along a line) at low temperatures

$\displaystyle l_{\rm flip}=\braket{l}$
$\displaystyle=\sum_{l=1}^{\infty}lP(l)=p_{\rm
flip}\sum_{l=1}^{\infty}l[1-p_{\rm flip}]^{l-1}=p_{\rm
flip}\frac{d}{dr}\sum_{l=0}^{\infty}r^{l}\Big{|}_{r=1-p_{\rm flip}}=p_{\rm
flip}\frac{d}{dr}\frac{1}{1-r}\Big{|}_{r=1-p_{\rm flip}}$
$\displaystyle=p_{\rm flip}\frac{1}{p_{\rm flip}^{2}}=p_{\rm
flip}^{-1}=e^{z\beta|J|/2}.$ (C.2)

Now, spin flips along any of the three lines $y=-1,0,+1$ will give a change of
$-|J|/2$ in the potential. However, at low temperatures we may treat the legs
as independent. As a result, taking the other legs into account will not alter
this asymptotic result.

Second, I now use this average distance to calculate the localization length
for $|J|/t\gg 1$. In this limit, the hole is effectively a single particle in
a one-dimensional infinite square well potential of length $l_{\rm flip}$.
Now, the initial state is centered at $x=0$. Since it is on a single lattice
site, in this continuum limit I will take it to be a constant with width $1$
(in units of the lattice spacing)

$\displaystyle\Psi(x,t=0)=1,\;-1/2\leq x\leq+1/2.$ (C.3)

Since this is even, there is only an overlap with the even eigenfunctions,
$\psi_{2n+1}(x)=\sqrt{2/l_{\rm flip}}$ $\cos[(2n+1)\pi x/l_{\rm flip}]$, where
$n=0,1,2,\dots$. The overlaps are then

$c_{2n+1}=\braket{\psi_{2n+1}}{\Psi(t=0)}=\int_{-1/2}^{+1/2}dx\psi_{2n+1}(x)=\frac{2\sqrt{2l_{\rm
flip}}}{(2n+1)\pi}\sin\left[\frac{(2n+1)\pi}{2l_{\rm flip}}\right].$ (C.4)

To compute the asymptotic mean-square distance,
$\braket{x^{2}}=\sum_{n}|c_{2n+1}|^{2}\bra{\psi_{2n+1}}x^{2}\ket{\psi_{2n+1}}$,
we need the mean-square distance for each of the contributing eigenfunctions.
These are

$\bra{\psi_{2n+1}}x^{2}\ket{\psi_{2n+1}}=\frac{2}{l_{\rm flip}}\int_{-l_{\rm
flip}/2}^{+l_{\rm flip}/2}dxx^{2}|\psi_{2n+1}(x)|^{2}=\frac{l_{\rm
flip}^{2}}{4}\left[\frac{1}{3}-\frac{2}{(2n+1)^{2}\pi^{2}}\right].$ (C.5)

The resulting asymptotic mean-square distance is

$\braket{x^{2}}=\sum_{n}|c_{2n+1}|^{2}\bra{\psi_{2n+1}}x^{2}\ket{\psi_{2n+1}}=\sum_{n}|c_{2n+1}|^{2}\frac{l_{\rm
flip}^{2}}{4}\left[\frac{1}{3}-\frac{2}{(2n+1)^{2}\pi^{2}}\right]\to\frac{l_{\rm
flip}^{2}}{12}.$ (C.6)

Here, I use that the $1/(2n+1)^{2}$ term will contribute with a term linear in
$l_{\rm flip}$, which is exponentially small compared to $l_{\rm flip}$ at low
temperatures. Note that this asymptotic behavior is actually independent of
the particular choice of the initial state in Eq. (C.3), because only the
common $l_{\rm flip}^{2}/12$ part for the mean-square distance is retained.
Hence, the asymptotic rms distance is

$x_{\rm rms}\to\sqrt{\braket{x^{2}}}=\frac{l_{\rm
flip}}{2\sqrt{3}}=\frac{e^{z\beta|J|/2}}{2\sqrt{3}}.$ (C.7)

This is valid in the limit of strong spin couplings, $|J|/t\gg 1$, and at low
temperatures, $\beta|J|\gg 1$.

##  Appendix D Entropy for antiferromagnetic couplings

In this Appendix, I calculate the entropy as a function of temperature for
antiferromagnetic couplings, both for the square and triangular lattice.

The calculation is carried through by using the thermodynamic relation
$F=U-TS$, between the free energy $F$, the average energy
$U=\braket{{\hat{H}}_{J}}$, and the entropy, $S$. In particular, both for the
square and triangular lattice there are explicit expressions for $F$ and $U$,
whereby the entropy per particle may readily be computed as

$\displaystyle\frac{S}{k_{B}N}=\frac{\beta}{N}[U-F],$ (D.1)

with $\beta=1/(k_{B}T)$ the inverse temperature. For the square lattice in
particular [51, 52],

$\displaystyle-\beta\frac{F}{N}=\ln\left[2\cosh\left(\frac{\beta|J|}{2}\right)\right]+\frac{1}{\pi}\int_{0}^{\pi/2}d\varphi\ln\left[f(k,\varphi)\right],$
(D.2)

with
$f(k,\varphi)=\frac{1}{2}\left\\{1+(1-k^{2}(\beta|J|)\sin^{2}\varphi)^{1/2}\right\\}$
and $k(\beta|J|)=2\sinh(\beta|J|/2)/\cosh^{2}(\beta|J|/2)$. From the relation
$U=\braket{{\hat{H}}_{J}}=-\partial_{\beta}\ln(Z)$, where $Z={\rm
tr}[e^{-\beta{\hat{H}}_{J}}]=e^{-\beta F}$ is the partition function, it
follows that $F=-\ln(Z)/\beta$, and hereby

$\displaystyle\frac{U}{N|J|}=-\partial_{\beta|J|}\left[-\frac{\beta
F}{N}\right]=$ $\displaystyle\frac{1}{2}{\rm
tanh}[\frac{\beta|J|}{2}]+\frac{1-2{\rm
tanh}(\beta|J|/2)}{\cosh(\beta|J|/2)}k(\beta|J|)\frac{1}{\pi}\int_{0}^{\pi/2}d\varphi
f^{-1}(k,\varphi)(1-k^{2}\sin^{2}\varphi)^{-1/2}.$ (D.3)

The integrals appearing in Eqs. (D.2) and (D.3) are easily solved numerically.
By insertion in Eq. (D.1), we hereby have the entropy per particle for the
square Ising lattice.

For the triangular lattice, there are also exact expressions. The free energy
can be calculated from [25]

$\displaystyle-\beta\frac{F}{N}=$ $\displaystyle\ln\left[2e^{-\beta
J/4}\cosh\left(\frac{\beta J}{2}\right)\right]$
$\displaystyle+\frac{1}{2\pi^{2}}\int_{0}^{\pi}\\!\\!d\omega_{1}\\!\\!\int_{0}^{\pi}\\!\\!d\omega_{2}\ln\left[1+4\kappa\cos(\omega_{1})\cos(\omega_{2})-4\kappa\cos^{2}(\omega_{2})\right],$
(D.4)

where $\kappa(\beta J)=[e^{-\beta J}-1]/[e^{-\beta J}+1]^{2}$. Moreover, there
is a closed form expression for the average energy [26],

$\displaystyle\frac{U}{NJ}=$
$\displaystyle\frac{1}{2(1-\mu)}\left[1-\frac{4\mu(3-\mu)}{4\sqrt{|\mu|}+\sqrt{(|\mu|+1)^{3}(3-|\mu|)}}\frac{2}{\pi}K(x)\right],$
(D.5)

where $\mu=1-2{\rm tanh}(-\beta J/2)$, $K(x)$ is the complete elliptic
integral of the first kind, and $x(\beta
J)=[4\sqrt{|\mu|}-\sqrt{(|\mu|+1)^{3}(3-|\mu|)}\,]/[4\sqrt{|\mu|}+\sqrt{(|\mu|+1)^{3}(3-|\mu|)}]$.
By inserting Eqs. (D.4) and (D.5) into Eq. (D.1), we thus obtain the entropy
per particle for the triangular Ising lattice.

## References

  * [1] K. K. Nielsen,  _Thermally induced localization of dopants in a magnetic spin ladder_ (2023), 2310.11193.
  * [2] J. Carlström, N. Prokof’ev and B. Svistunov,  _Quantum walk in degenerate spin environments_ ,  Phys. Rev. Lett. 116, 247202 (2016),  10.1103/PhysRevLett.116.247202.
  * [3] 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_ ,  Phys. Rev. B 96, 014303 (2017),  10.1103/PhysRevB.96.014303.
  * [4] L. Hahn, A. Bohrdt and F. Grusdt,  _Dynamical signatures of thermal spin-charge deconfinement in the doped ising model_ ,  Phys. Rev. B 105, L241113 (2022),  10.1103/PhysRevB.105.L241113.
  * [5] F. Grusdt and L. Pollet,  _${\mathbb{z}}_{2}$ parton phases in the mixed-dimensional $t-{J}_{z}$ model_,  Phys. Rev. Lett. 125, 256401 (2020),  10.1103/PhysRevLett.125.256401.
  * [6] J. Zaanen, O. Y. Osman, H. V. Kruis, Z. Nussinov and J. Tworzydlo,  _The geometric order of stripes and luttinger liquids_ ,  Philosophical Magazine B 81(10), 1485 (2001),  10.1080/13642810108208566,  https://doi.org/10.1080/13642810108208566.
  * [7] W. F. Brinkman and T. M. Rice,  _Single-particle excitations in magnetic insulators_ ,  Phys. Rev. B 2, 1324 (1970),  10.1103/PhysRevB.2.1324.
  * [8] S. Schmitt-Rink, C. M. Varma and A. E. Ruckenstein,  _Spectral function of holes in a quantum antiferromagnet_ ,  Phys. Rev. Lett. 60, 2793 (1988),  10.1103/PhysRevLett.60.2793.
  * [9] C. L. Kane, P. A. Lee and N. Read,  _Motion of a single hole in a quantum antiferromagnet_ ,  Phys. Rev. B 39, 6880 (1989),  10.1103/PhysRevB.39.6880.
  * [10] G. Martinez and P. Horsch,  _Spin polarons in the t-j model_ ,  Phys. Rev. B 44, 317 (1991),  10.1103/PhysRevB.44.317.
  * [11] G. F. Reiter,  _Self-consistent wave function for magnetic polarons in the t-j model_ ,  Phys. Rev. B 49, 1536 (1994),  10.1103/PhysRevB.49.1536.
  * [12] J. R. Schrieffer and J. S. Brooks,  _Handbook of high-temperature superconductivity_ ,  Springer (2007).
  * [13] E. Dagotto, R. Joynt, A. Moreo, S. Bacci and E. Gagliano,  _Strongly correlated electronic systems with one hole: Dynamical properties_ ,  Phys. Rev. B 41, 9049 (1990),  10.1103/PhysRevB.41.9049.
  * [14] E. Dagotto,  _Correlated electrons in high-temperature superconductors_ ,  Rev. Mod. Phys. 66, 763 (1994),  10.1103/RevModPhys.66.763.
  * [15] C. Gross and I. Bloch,  _Quantum simulations with ultracold atoms in optical lattices_ ,  Science 357(6355), 995 (2017),  10.1126/science.aal3837.
  * [16] 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_ ,  Phys. Rev. X 11, 021022 (2021),  10.1103/PhysRevX.11.021022.
  * [17] Y. Wang, A. Bohrdt, S. Ding, J. Koepsell, E. Demler and F. Grusdt,  _Higher-order spin-hole correlations around a localized charge impurity_ ,  Phys. Rev. Research 3, 033204 (2021),  10.1103/PhysRevResearch.3.033204.
  * [18] M. L. Prichard, B. M. Spar, I. Morera, E. Demler, Z. Z. Yan and W. S. Bakr,  _Directly imaging spin polarons in a kinetically frustrated hubbard system_ (2023), 2308.12951.
  * [19] 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(7769), 358 (2019),  10.1038/s41586-019-1463-1.
  * [20] 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(6563), 82 (2021),  10.1126/science.abe7165.
  * [21] M. Lebrat, M. Xu, L. H. Kendrick, A. Kale, Y. Gang, P. Seetharaman, I. Morera, E. Khatami, E. Demler and M. Greiner,  _Observation of nagaoka polarons in a fermi–hubbard quantum simulator_ ,  Nature 629(8011), 317 (2024),  10.1038/s41586-024-07272-9.
  * [22] S. Hirthe, T. Chalopin, D. Bourgund, P. Bojović, A. Bohrdt, E. Demler, F. Grusdt, I. Bloch and T. A. Hilker,  _Magnetically mediated hole pairing in fermionic ladders of ultracold atoms_ ,  Nature 613(7944), 463 (2023),  10.1038/s41586-022-05437-y.
  * [23] A. Lagendijk, B. v. Tiggelen and D. S. Wiersma,  _Fifty years of Anderson localization_ ,  Physics Today 62(8), 24 (2009),  10.1063/1.3206091.
  * [24] K. K. Nielsen,  _Exact dynamics of two holes in two-leg antiferromagnetic ladders_ ,  Phys. Rev. B 108, 085125 (2023),  10.1103/PhysRevB.108.085125.
  * [25] G. H. Wannier,  _Antiferromagnetism. the triangular ising net_ ,  Phys. Rev. 79, 357 (1950),  10.1103/PhysRev.79.357.
  * [26] G. H. Wannier,  _Antiferromagnetism. the triangular ising net_ ,  Phys. Rev. B 7, 5017 (1973),  10.1103/PhysRevB.7.5017.
  * [27] 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(7269), 74 (2009),  10.1038/nature08482.
  * [28] 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(7311), 68 (2010),  10.1038/nature09378.
  * [29] J. Zeiher, R. van Bijnen, P. Schauß, S. Hild, J.-y. Choi, T. Pohl, I. Bloch and C. Gross,  _Many-body interferometry of a rydberg-dressed spin lattice_ ,  Nature Physics 12(12), 1095 (2016),  10.1038/nphys3835.
  * [30] G. F. Newell,  _Crystal statistics of a two-dimensional triangular ising lattice_ ,  Phys. Rev. 79, 876 (1950),  10.1103/PhysRev.79.876.
  * [31] U. Wolff,  _Collective monte carlo updating for spin systems_ ,  Phys. Rev. Lett. 62, 361 (1989),  10.1103/PhysRevLett.62.361.
  * [32] E. Luijten,  _Introduction to Cluster Monte Carlo Algorithms_ , pp. 13–38,  Springer Berlin Heidelberg, Berlin, Heidelberg,  ISBN 978-3-540-35273-0,  10.1007/3-540-35273-2_1 (2006).
  * [33] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller,  _Equation of state calculations by fast computing machines_ ,  The Journal of Chemical Physics 21(6), 1087 (1953),  10.1063/1.1699114.
  * [34] W. K. Hastings,  _Monte carlo sampling methods using markov chains and their applications_ ,  Biometrika 57(1), 97 (1970),  10.1093/biomet/57.1.97.
  * [35] K. K. Nielsen, T. Pohl and G. M. Bruun,  _Nonequilibrium hole dynamics in antiferromagnets: Damped strings and polarons_ ,  Phys. Rev. Lett. 129, 246601 (2022),  10.1103/PhysRevLett.129.246601.
  * [36] H. Bruus and K. Flensberg,  _Many-body quantum theory in condensed matter physics – an introduction_ ,  Oxford University Press (2016).
  * [37] Z. Liu and E. Manousakis,  _Dynamical properties of a hole in a heisenberg antiferromagnet_ ,  Phys. Rev. B 45, 2425 (1992),  10.1103/PhysRevB.45.2425.
  * [38] N. 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 (2021).
  * [39] K. Shen, K. Sun, M. F. Gelin and Y. Zhao,  _Finite-temperature hole–magnon dynamics in an antiferromagnet_ ,  The Journal of Physical Chemistry Letters 15(2), 447 (2024),  10.1021/acs.jpclett.3c03298.
  * [40] D. N. Sheng, Y. C. Chen and Z. Y. Weng,  _Phase string effect in a doped antiferromagnet_ ,  Phys. Rev. Lett. 77, 5102 (1996),  10.1103/PhysRevLett.77.5102.
  * [41] K. Wu, Z. Y. Weng and J. Zaanen,  _Sign structure of the $t\text{$-$}j$ model_,  Phys. Rev. B 77, 155102 (2008),  10.1103/PhysRevB.77.155102.
  * [42] Z. Zhu and Z.-Y. Weng,  _Quasiparticle collapsing in an anisotropic $t-j$ ladder_,  Phys. Rev. B 92, 235156 (2015),  10.1103/PhysRevB.92.235156.
  * [43] S. R. White, D. J. Scalapino and S. A. Kivelson,  _One hole in the two-leg $t\text{$-$}j$ ladder and adiabatic continuity to the noninteracting limit_,  Phys. Rev. Lett. 115, 056401 (2015),  10.1103/PhysRevLett.115.056401.
  * [44] R.-Y. Sun, Z. Zhu and Z.-Y. Weng,  _Localization in a $t\text{$-$}j$-type model with translational symmetry_,  Phys. Rev. Lett. 123, 016601 (2019),  10.1103/PhysRevLett.123.016601.
  * [45] J.-Y. Zhao, S. A. Chen, H.-K. Zhang and Z.-Y. Weng,  _Two-hole ground state: Dichotomy in pairing symmetry_ ,  Phys. Rev. X 12, 011062 (2022),  10.1103/PhysRevX.12.011062.
  * [46] P. Jordan and E. Wigner,  _Über das paulische äquivalenzverbot_ ,  Zeitschrift für Physik 47(9), 631 (1928),  10.1007/BF01331938.
  * [47] K. K. Nielsen,  _Omnipresent bound state of two holes in antiferromagnetic bethe lattices_ ,  Phys. Rev. B 108, 085122 (2023),  10.1103/PhysRevB.108.085122.
  * [48] B. M. McCoy and J.-M. Maillard,  _The Importance of the Ising Model_ ,  Progress of Theoretical Physics 127(5), 791 (2012),  10.1143/PTP.127.791.
  * [49] J. Stephenson,  _Ising-Model Spin Correlations on the Triangular Lattice_ ,  Journal of Mathematical Physics 5(8), 1009 (1964),  10.1063/1.1704202.
  * [50] J. Stephenson,  _Ising Model Spin Correlations on the Triangular Lattice. II. Fourth-Order Correlations_ ,  Journal of Mathematical Physics 7(6), 1123 (1966),  10.1063/1.1705003.
  * [51] C. K. Majumdar,  _Analytic properties of the onsager solution of the ising model_ ,  Phys. Rev. 145, 158 (1966),  10.1103/PhysRev.145.158.
  * [52] L. Onsager,  _Crystal statistics. i. a two-dimensional model with an order-disorder transition_ ,  Phys. Rev. 65, 117 (1944),  10.1103/PhysRev.65.117.