# Continuous similarity transformation for critical phenomena: easy-axis
antiferromagnetic XXZ model

Matthias R. Walther<EMAIL_ADDRESS>Department of Physics, Friedrich-
Alexander-Universität Erlangen-Nürnberg, Staudtstraße 7, 91058 Erlangen,
Germany  Dag-Björn Hering<EMAIL_ADDRESS>Condensed Matter Theory,
Technische Universität Dortmund, Otto-Hahn-Straße 4, 44221 Dortmund, Germany
Götz S. Uhrig<EMAIL_ADDRESS>Condensed Matter Theory, Technische
Universität Dortmund, Otto-Hahn-Straße 4, 44221 Dortmund, Germany  Kai P.
Schmidt<EMAIL_ADDRESS>Department of Physics, Friedrich-
Alexander-Universität Erlangen-Nürnberg, Staudtstraße 7, 91058 Erlangen,
Germany

###### Abstract

We apply continuous similarity transformations (CSTs) to the easy-axis
antiferromagnetic XXZ-model on the square lattice. The CST flow equations are
truncated in momentum space by the scaling dimension $d$ so that all
contributions with $d\leq 2$ are taken into account. The resulting quartic
magnon-conserving effective Hamiltonian is analyzed in the zero-, one-, and
two-magnon sector. In this way, a quantitative description of the ground-state
energy, the one-magnon dispersion and its gap as well as of two-magnon bound
states is gained for anisotropies ranging from the gapped Ising model to the
gapless Heisenberg model. We discuss the critical properties of the gap
closing as well as the evolution of the one-magnon roton mininum. The
excitation energies of two-magnon bound states are calculated and their decay
into the two-magnon continuum is determined via the inverse participation
ratio.

##  I Introduction

The collective behavior of interacting quantum matter is an important topic in
condensed matter physics as well as quantum optics over the last decades,
since more and more intriguing facets of quantum systems are discovered giving
rise to unexpected properties of quantum materials and fascinating
perspectives in quantum technological applications. Such quantum behavior
manifests itself most prominently at low temperatures where quantum
fluctuations dominate and therefore zero-temperature quantum phases and phase
transition with emergent properties are a relevant subject of current research
[1]. Here two-dimensional systems are most challenging, since, unlike in one
dimension, no analytical nor generic numerical solutions are available and,
unlike in three dimensions, quantum fluctuations are still important. At the
same time two-dimensional correlated quantum systems are known for exciting
physical properties such as high-$T_{\rm c}$ superconductivity, the fractional
quantum Hall effect as well as highly entangled topological order.

One paradigmatic microscopic model is the antiferromagnetic
spin-$\tfrac{1}{2}$ square lattice Heisenberg model, which has been
intensively studied, in particular since it represents the relevant low-energy
description of the quantum magnetism in the undoped cuprate superconductors
[2]. While properties of the long-range ordered Néel state and the gapless
low-energy magnon excitations are well known [3], a quantitative understanding
of the high-energy part of the one-magnon dispersion displaying a
characteristic roton minimum as well as of dynamical correlation functions
including important contributions of magnon continua have only been achieved
in recent years [4, 5, 6, 7, 8]. In parallel, the latter quantities become
more and more important experimentally due to the continuously increasing
resolution of various spectroscopic techniques such as inelastic neutron
scattering, terahertz spectroscopy, or resonant inelastic x-ray scattering.

A quantitative description of the magnon excitation spectrum in the Heisenberg
model on the square lattice requires a proper treatment of the magnon-magnon
interaction. Technically, this can be achieved by continuous similarity
transformations (CSTs) using the scaling dimension of operators as truncation
criterion of the associated flow equations [4, 5]. The CST yields an effective
Hamiltonian in momentum space which is block-diagonal in the dressed magnon
quasi-particles and which allows one to quantitatively access one- and multi-
magnon excitation energies. It is a natural next step to ask whether CSTs can
be applied successfully when modifying the model away from the gapless
Heisenberg point. The easy-axis antiferromagnetic XXZ model on the square
lattice represents one such modification which tunes continuously from the
isotropic Heisenberg model to the classical Ising model with local and gapped
magnon excitations. The magnon excitation spectrum of the XXZ model has been
investigated by several means, in particular with high-order series expansions
(SE) [9, 10, 11], diagrammatic spin wave theory [12, 13], density matrix
renormalization group (DMRG) [7], and quantum Monte Carlo simulations (QMC)
[14, 15]. Apart from the gap closing transition, a notable feature is the
appearance of two-magnon bound states in large parts of the phase diagram,
which are not present at the Heisenberg point where there is only a resonance
of the Higgs-like amplitude mode within the magnon continuum [16, 5].
Interestingly, in coupled XXZ spin ladder compounds the finite easy-axis spin
anisotropy has lead to the experimental discovery of the Higgs amplitude mode
near a quantum critical point [17, 18].

This paper is structured as follows. In Sec. II we describe the relevant
properties of the easy-axis XXZ model on the square lattice. A description of
the CST as well as of technical aspects are provided in Sec. III; this section
can be skipped by those who want to focus on the results. Results for the
ground-state energy, the one-magnon dispersion, and the two-magnon sector are
presented and discussed in Sec. IV.We conclude our study in Sec. V.

##  II Model

We study the antiferromagnetic easy-axis XXZ model on the square lattice as
prototypical extension of the isotropic Heisenberg model. Here, we set the
lattice constant $a$ to unity. The anisotropy induces a finite energy gap
which vanishes in the isotropic limit where the Goldstone theorem applies [3].
Hence, we can study how the gap closes upon approaching the isotropic point
and whether binding effects occur. To establish a deeper understanding of the
performance and the applicability of continuous basis transformations [4, 5]
for vanishing energy gaps we investigate the XXZ model with easy-axis
anisotropies from the Ising limit to the isotropic Heisenberg limit. The
Hamiltonian reads

$\displaystyle\mathcal{H}$
$\displaystyle=J\sum_{\Braket{i,j}}\left[S_{i}^{z}S_{j}^{z}+\lambda\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}\right)\right]$
(1)

where $J>0$ is the exchange coupling between spins on site ${i}$,
$\Braket{i,j}$ represents the sum over pairs of nearest neighbors on the
square lattice where each bond is counted once. The spin anisotropy is
introduced by $\lambda\in[0,1]$ with the $z$-axis being the easy axis. For
$\lambda=0$ the XXZ model reduces to the Ising model which is exactly solvable
[19]. It displays a long-range ordered ground-state, a gapped, completely flat
dispersion and a two-particle sector with four nearest-neighbor bound states
per site in one sublattice. The Ising model is a frequently used starting
point for perturbative treatments of Heisenberg and XXZ models, see for
instance Ref. [10, 20, 11]. This Ising phase continues adiabatically, i.e.,
without a quantum phase transition, up to $\lambda=1$. At this point the
symmetry of Eq. 1 changes from $\mathds{Z}_{2}\times\text{U}(1)$ to
$\text{SU}(2)$ of all spin rotations. Here, at $\lambda=1$, the XXZ model
becomes the isotropic Heisenberg model still displaying long-range magnetic
order which breaks its _continuous_ symmetry spontaneously so that the
Goldstone theorem implies the existence of gapless elementary excitations, the
magnons [3]. Due to the absence of an energy gap, multi-magnon continua exist
starting right at the energies of single magnons extending to higher energies.
No bound states are known at the isotropic point, but the attractive
interaction between magnons leads to considerable shift of spectral weight to
lower energies [4, 5]. We focus on the following physical quantities:

##### Ground-state energy

The ground-state in the Ising limit ($\lambda=0$) is precisely the Néel state
with a ground-state energy per site $e=-J/2$ that spontaneously breaks the
discrete $\mathds{Z}_{2}$ symmetry. The other limiting case is $\lambda=1$
where the ground-state is the one of the isotropic Heisenberg model. In the
thermodynamics limit, it spontaneously breaks the SU(2) symmetry. The ground-
state energy per site obtained from QMC takes the value
$e^{\text{QMC}}=-0.669437(5)J$ [14]. This value serves as a first testbed for
our approach.

##### Dispersion $\omega(\bm{k})$ and single-particle gap $\Delta$

For $\lambda=0$ the dispersion is completely flat $\omega(\bm{k})=2J$ because
the elementary excitations, spin flips, are completely immobile. For
$\lambda>0$ more and more features emerge. In the magnetic Brillouin zone
(MBZ), the dispersion displays a minimum at $\bm{k}=(0,0)$ which defines the
spin gap $\Delta\coloneqq\omega(\bm{k}=\bm{0})$. It closes at the isotropic
Heisenberg point $\lambda=1$ according to Goldstone’s theorem. The closure of
the spin gap follows a square root power law in spin wave theory [10] although
series expansions would be consistent for slightly higher exponents as well
[9]. Additionally, the dispersion along the line from $\bm{k}=(\pi,0)$ to
$\bm{k}=(\frac{\pi}{2},\frac{\pi}{2})$ in the Brillouin zone develops a
distinct local minimum at $\bm{k}=(\pi,0)$, called roton minimum, and the
global maximum at $\bm{k}=(\frac{\pi}{2},\frac{\pi}{2})$. The roton minimum
occurs only in third order spin wave theory [21]. It is understood to be due
to magnon-magnon interactions [4, 5] and a peculiar cancellations for motion
along the diagonals [7]. We will analyse the spin gap $\Delta$, the roton
minimum, and dispersion maximum in dependence on the anisotropy
$\lambda\in[0,1]$ with a special focus on the critical behavior of $\Delta$
for $\lambda\to 1$.

##### Two-particle sector

The eigenstates which are built from two elementary excitations, i.e., two
magnons, are particularly interesting because they reflect the degree of
magnon-magnon interaction. For strong interaction relative to the kinetic
energy of the magnons, bound states occur in the gapped phase. The two-
particle spectrum of the Ising model features degenerate bound states at $3J$
stemming from adjacent pairs of spin flips. All other states of two magnons
have energy $4J$ and are thus highly degenerate. For finite $\lambda$, these
states evolve into energy continua of scattering states. The four bound states
exist also for finite $\lambda$, but merge with the continuum for values
$\lambda$ close to but smaller than one [22, 13, 11] consistent with quantum
Monte Carlo simulations of coupled XXZ spin ladders [18]. We track these bound
states for all $\lambda$ and determine at which values of $\lambda$ the bound
states vanish in the two-magnon continuum by analyzing the locality of bound
states with the help of the inverse participation ratio (IPR) [23, 24].

##  III Methods and Techniques

Here we recall the fundamental idea of continuous basis transformations in
general and the continuous similarity transformation (CST) in particular which
we employ here. Further methodological aspects are also recapitulated such as
the self-consistent mean-field approach to obtain a good starting point for
the CST and the IPR as a measure of locality. Technical aspects such as the
discretization of the magnetic Brillouin zone (MBZ) and various extrapolations
are addressed as well.

###  III.1 Continuous Similarity Transformation (CST)

We use the CST method in momentum space as introduced in Powalski_2015,
Powalski_2018,[0]Ref. [4, 5]. As in the flow equation approach [25, 26, 27],
we introduce an auxiliary flow variable $\ell$ that parametrizes a continuous
similarity transformation from $\mathcal{H}(\ell=0)=\mathcal{H}_{0}$ to an
effective Hamiltonian $\mathcal{H}(\ell=\infty)=\mathcal{H}_{\text{eff}}$. The
flow is governed by the flow equation
$\partial_{\ell}\mathcal{H}(\ell)=\left[\eta(\ell),\mathcal{H}(l)\right]$,
where the generator $\eta(\ell)$ determines the course of the flow and its end
point, the fixed point. We use the quasi-particle conserving generator
$\eta_{ij}=(q_{ii}-q_{jj})h_{ij}$ [26, 28, 29], where $q_{ii}$ are the eigen
values of the operator $\hat{Q}$ counting the number of quasi-particles, here
magnons, present in the system. Then the final effective Hamiltonian
$\mathcal{H}_{\text{eff}}$ has a block diagonal form, where the subspace on
which each block acts has a given number of quasi-particles. Subspaces with
different numbers of magnons do not influence each other; they are
disentangled. The coefficient in the zero-magnon block is the ground-state
energy, the one-magnon block contains the dispersion, and the two-magnon
energies are obtained by diagonalizing the two-magnon block.

Computing $\left[\eta(\ell),\mathcal{H}(\ell)\right]$ for the flow equation
generally yields an infinite number of terms which cannot be handled
numerically. Therefore, one has to truncate the flow equations in a systematic
and controlled way. We follow here the idea to use the scaling dimension as
truncation criterion. The higher this dimension the less important are the
corresponding terms for the physics at low energies. For the sine-Gordon
model, the operator product expansion implements this idea and allows one to
systematically truncate the flow according to Wegner’s generator [30, 31].
Recently, the truncation according to scaling dimension justified to truncate
the flow to the bilinear and quadrilinear terms in a bosonic Dyson-Maleev
representation [32, 33] of the spins [4, 5]. At the isotropic point
$\lambda=1$, the dispersion encoded in the bilinear terms has dimension 1
while the quadrilinear terms have dimension 2 for essentially constant
prefactors. The neglected hexalinear terms would have dimension 4. Here, we
focus on the gapped spectrum of the anisotropic easy-axis so that the
dimension of the bilinear dispersion is 0. Hence the same truncation remains
justified and we do not change it because the limit $\lambda\to 1$ shall be
captured as well. Eventually, the flow equations are solved numerically for a
discretized MBZ. Details of the numerical methods are presented in Secs.
III.4, III.5 and III.7. We will see that the results for the dispersion and
the energy gap in Sec. IV.3 strongly corroborate our approach. In addition, we
discern and discuss deviations from exact series results [13, 11] for the
energies of the bound states, see Sec. IV.4.

###  III.2 Self-consistent mean-field approach

Starting truncated transformations from spin representations in real space has
turned out to be difficult and prone to large truncation errors [34]. Hence,
we use the self-consistent mean-field solution in momentum space as initial
Hamiltonian even though a variation of this solution can sometimes perform
better [35]. The crucial advantage of the self-consistent mean-field approach
is that it captures the Goldstone bosons of the isotropic Heisenberg model. To
be self-contained we sketch the necessary analytical steps applied to
Hamiltonian (1). A detailed derivation with intermediate results is given in
Appendix A.

The ground-state for $0\leq\lambda\leq 1$ displays long-range Néel order.
Thus, we choose the same classical Néel state with alternating spin-up and
spin-down on sublattices A and B, respectively, as the reference state for all
$\lambda$. We introduce bosonic degrees freedom that represent deviations from
this state by the non-Hermitian Dyson-Maleev representation [32, 33]. Then we
perform a standard mean-field decoupling [36] based on Wick’s theorem in real
space. A Fourier transformation and a subsequent Bogoliubov transformation to
the bosonic operators $\alpha^{(\dagger)}_{\bm{k}}$ and
$\beta^{(\dagger)}_{\bm{k}}$ diagonalize the Hamiltonian on the bilinear level
where the $\alpha$ bosons refer to the A sublattice and the $\beta$ bosons to
the B sublattice. The mean-field Hamiltonian reads

$\displaystyle\mathcal{H}=E_{0}+\sum_{\bm{k}}\omega_{0}(\bm{k})\big{(}:{\alpha^{\dagger}_{\bm{k}}\alpha^{\phantom{\dagger}}_{\bm{k}}}:+:{\beta^{\dagger}_{\bm{k}}\beta^{\phantom{\dagger}}_{\bm{k}}}:\big{)}+\Gamma_{0}+\mathcal{V}_{0}.$
(2)

The wave vectors $\bm{k}$ are taken from the MBZ depicted in Fig. 1. The
colons $:{...}:$ indicate normal ordered operators with respect to the mean-
field ground-state. The exact expressions of the ground-state energy $E_{0}$,
the dispersion $\omega_{0}(\bm{k})$, and the two-particle interactions
$\mathcal{V}_{0}$ are given in Appendix A. The off-diagonal bilinear part
$\Gamma_{0}$ is zero after the self-consistent solution of the mean-field
decoupling.

After the mean-field decoupling and the CST flow the effective Hamiltonian is

$\displaystyle\mathcal{H}_{\mathrm{eff}}=E+\sum_{\bm{k}}\omega(\bm{k})\big{(}:{\alpha^{\dagger}_{\bm{k}}\alpha^{\phantom{\dagger}}_{\bm{k}}}:+:{\beta^{\dagger}_{\bm{k}}\beta^{\phantom{\dagger}}_{\bm{k}}}:\big{)}+\mathcal{V}.$
(3)

with the dispersion $\omega(\bm{k})$ of the quasi-particles in the new basis
and the quasiparticle conserving interaction term $\mathcal{V}$.

Figure 1: Magnetic Brillouin zone (MBZ) with different kinds of discretizing
mesh. In panel (a) important points such as $\bm{k}=(0,0)$ are sampled
directly; this kind of mesh is denoted by $N_{0}$. For $N_{0}$, the
calculations for $\lambda<1$ or the two-magnon sector do not require
additional interpolations. But at $\bm{k}=(0,0)$ divergences in the flow occur
for $\lambda\lesssim 1$. They can be avoided by a mesh of the type shown in
panel (b) for odd $L$; this mesh is denoted by $N_{\emptyset}$. We use it
especially for the analysis of the critical behavior of the spin gap. Due to
the point group symmetries only lattice sites in the marked blue area need to
be taken into account.

###  III.3 Inverse Participation Ratio

In Sec. IV.4 we determine for each of the four bound states at which value of
the anisotropy $\lambda$ it enters the two-magnon continuum, i.e., where it
ceases to exist. This merging into a continuum generically occurs slowly in
the sense that the slope of the lower band edge of the continuum and of the
energy of the bound state are the same [37, 38, 39]. Hence, the parameters
where the bound state ceases to exist are extremely difficult to determine
numerically. Therefore, we use another criterion, namely the locality of the
bound state in its relative coordinate. A true bound state is local in the
sense that it has a finite extension. If the extension is given by the system
size no binding occurs anymore. A good measure of the locality of a state is
the inverse participation ratio (IPR) [23, 24]

$I=\sum_{\bm{r}}|\Psi(r)|^{4}=\sum_{\bm{r}}|\rho(r)|^{2}.$ (4)

with density $\rho(r)$ at site $r$ of the normalized wave function $\Psi(r)$.
For an extended state the IPR $I$ scales proportional to $1/L^{2}=1/N$ if $L$
is the number of sites in one direction in two dimensions, This is easily seen
in the limiting case of $\rho(r)=1/N$ where $I=1/N$ holds. The other limiting
case is a completely local state with $\rho(0)=1$ and zero elsewhere so that
$I=1$ holds. Hence one has to determine whether the IPR vanishes $\propto 1/N$
for $N\to\infty$ or stays finite in this limit.

###  III.4 Flow equation and residual off-diagonality

The flow equations are a set of ordinary differential equations (ODE) for the
coefficients of the contributing monomials of operators which depend on wave
vectors. We solve these ODEs numerically for a finite lattice with a standard
Dormand-Prince 5 method as implemented by the odeint project [40] up to linear
lattice size $L=22$. In order to avoid the accumulation of numerical errors
and to enhance the performance we use the point group symmetries in the MBZ as
depicted in Fig. 1. In addition, we use an unambiguous symmetrized
representation of the coefficients. We stop the flow when the residual off-
diagonality (ROD), i.e., the square root of the sum over all squared entries
of the generator $\eta$ [29], has dropped to values below $10^{-6}J$.

###  III.5 Discretization of the magnetic Brillouin zone

One has to pay special attention to the wave vector $\bm{k}=(0,0)$ at
$\lambda\lesssim 1$. The gapless mean-field solution (2) at the isotropic
point displays diverging coefficients whenever one of the wave vectors in the
quadrilinear monomials vanishes, see for instance Uhrig_2013,[0]Ref. [35].
Hence, this solution cannot be used for a discretized MBZ. Note that these
divergences are integrable in the thermodynamic limit where the wave vector
dependence is continuous. So they do not constitute a physical, but rather a
numerical issue. For $\lambda\to 1$ and especially for small $L$ we
encountered signs of problematic divergences in form of a non-converging flow.
One way to circumvent this problem consists in using a reciprocal lattice
where $\bm{k}=(0,0)$ does not occur, for instance by using odd values of $L$
with antiperiodic boundary conditions as illustrated in right panel of Fig. 1.
In real space, the antiperiodic boundary conditions are achieved by giving
spin flip terms across the boundary of the system’s cluster a minus sign. In
this case, the cluster is a rhombus with side length of $\sqrt{2}L$. In the
following, we denote lattice discretizations sampling $\bm{k}=(0,0)$ by
$N_{0}$ and discretizations not sampling $\bm{k}=(0,0)$ by $N_{\emptyset}$. We
also implemented the approach used by Powalski et al. [4, 5] where
coefficients are set to zero, for which at least one wave vector vanishes,
i.e., they are neglected. The results extrapolated to $L\rightarrow\infty$
agree well with the results presented below, but they show a significantly
slower convergence in $1/L$. For this reason we do not consider them any
further.

Since the approach $N_{\emptyset}$ excludes the wave vector at which the spin
gap $\Delta=\omega(\bm{k}=(0,0))$ can be read off by construction, the gap
value must be retrieved otherwise. To this end, we adapt the fit

$\displaystyle\omega(k)^{2}=\sum_{m=0}^{M}A_{m}\cos(mk),$ (5)

used in Ref. [41] for a one-dimensional model to our two-dimensional case. We
fit the discrete dispersion values obtained by CST along the line from $(0,0)$
to $(\pi,0)$ in order to determine the coefficients $A_{m}$. A very fast
convergence of these coefficients is found for increasing $M$ so that we can
interpolate the dispersion very reliably in the whole $N_{\emptyset}$.

###  III.6 Determination of correlation lengths

As mentioned in Sec. III.2, we perform the mean-field decoupling to determine
the initial conditions in the thermodynamic limit, while we solve the flow
equations for a discretized lattice as discussed in the previous section.
Hence, our results combine elements of calculations in the thermodynamic limit
with calculations on a finite lattice. In order to judge if the discretization
in the MBZ induces finite-size effects, we calculate the ground-state
correlation length $\xi$. Whenever it is significantly larger than the linear
length scale of the system, finite-size effects are likely to occur. We expect
the correlation length $\xi$ to be the same in all directions for $\xi\gg a$
where $a$ is the lattice constant. For $\lambda\rightarrow 1$, $\xi$ diverges.
Thus, it is sufficient to perform the same one-dimensional fit for a single
direction as in Sec. III.5, Eq. 5 and calculate $\xi$ from the one-magnon
dispersion $\omega$ via

$\displaystyle\omega(\mathrm{i}\kappa)=0\quad\text{with}\quad\mathrm{Re}=1/\xi$
(6)

adapted from Okunishi_2001,[0]Ref. [41].

###  III.7 Extrapolation in the linear system size $L$

For all results shown an extrapolation in $1/L\to 0$ is performed. For small
values of $\lambda$ up to $\lambda\approx 0.8$, the computed values do not
show any relevant finite-size effects. There is no discernible difference
between the two kinds of mesh $N_{\emptyset}$ and $N_{0}$. Hence, all
extrapolating schemes yield essentially the same result, see panel (a) in Fig.
2. For $\lambda>0.8$, we observe two trends for increasing $L$, see panels (b)
and (c) in Fig. 2. First, the values for $N_{0}$ decrease monotonically
decreasing while those for $N_{\emptyset}$ increase monotonically. Second, the
values from the mesh $N_{\emptyset}$ show an almost linear relation as
function of $1/L$ for $\lambda\lesssim 1$. The values we provide as the
estimate of the extrapolation are the results of monotonic quadratic fits

$\displaystyle x(1/L)=a+b(1/L+c)^{2}$ (7)

for the mesh $N_{\emptyset}$ with $c\geq 0$, evaluated at $1/L=0$. The error
estimate results from the difference between the values from the meshes
$N_{\emptyset}$ and $N_{0}$, see Fig. 2. For the correlation length $\xi$ a
linear extrapolation in $1/L$ was used; the error estimate remains the same.

The locality of bound states is measured by the IPR given in Eq. 4. This
quantity is extrapolated linearly in $1/N$ which works remarkably well
yielding positive values for bound states and values in the close vicinity of
zero otherwise.

Figure 2: Generic extrapolations in $1/L$ for $\lambda=0.7$, $\lambda=0.96$,
and $\lambda=0.995$. For small $\lambda$, calculations on the meshes
$N_{\emptyset}$ and $N_{0}$ converge, see panel (a). At higher values of
$\lambda$, see panel (b) and (c), values from $N_{\emptyset}$ are
monotonically increasing while values from $N_{0}$ are monotonically
decreasing. For the whole range, a monotonic quadratic fit for the
$N_{\emptyset}$ data is used to determine the values for $L\rightarrow\infty$.
The error estimate depicted in yellow results from the difference between the
values from $N_{\emptyset}$ and $N_{0}$ for the highest reached value of $L$.

##  IV Results

###  IV.1 Convergence

Generally, flow-equation approaches based on the quasi-particle conserving
generator transform the original Hamiltonian to a block-diagonal form, where
each block describes the action of the effective Hamiltonian on the subspace
of a given, constant number of quasi-particles, here magnons, determined by
$Q$. If, however, this quasi-particle picture breaks down, for example by a
second-order phase transition or by level crossings of modes with a different
number of quasi-particles, this separation into blocks is not possible
anymore. As a consequence, the flow does not converge and no description in
terms of the chosen quasi-particles is obtained. This feature makes sense
because it reflects that the physics can no longer be captured in terms of the
particular choice of quasi-particles. Conversely, if the flow converges this
indicates that the underlying quasi-particles picture does not break down and
that there are no inter-block level crossings in the chosen truncation scheme.
Hence, the description in terms of the chosen quasi-particles reflects the
underlying physics correctly.

The flow of the CST calculations with the truncation based on scaling
dimension converges for both types of meshes $N_{\emptyset}$ and $N_{0}$ for
all $\lambda$. Only for $\lambda\geq 0.99925$ and the mesh $N_{0}$ including
$\mathbf{k}=(0,0)$ we find a divergent flow, that can be traced back to
numerical artifacts, see Sect. III.5. In general, an analysis of the ground-
state energy, the magnon dispersion, and bound states in the two-magnon block
is feasible and well-founded.

###  IV.2 Ground-state energy

The first quantity of the XXZ-model we inspect is the ground-state energy per
site $e(\lambda)=E(\lambda)/(2L^{2})$. It is shown in Fig. 3 as computed by
CST, series expansion up to order 14 in $\lambda$ [9], and QMC calculations
[14].

Figure 3: Ground-state energy per site $e$ as a function of $\lambda$. Red
crosses are CST data which are compared to results of a series expansion(SE)
[9] drawn as blue line; the green star marks the QMC results for the isotropic
Heisenberg model [14]. As expected from second-order perturbation theory the
ground-state energy displays a monotonically decreasing behavior with negative
curvature.

Starting from the classical Néel state in the Ising limit, increasing
$\lambda$ induces more and more quantum fluctuations resulting in a ground-
state energy which continuously decreases. The results of the CST are
identical to the series expansion results within line width for the whole
$\lambda$ range. For $\lambda=1$, the CST result ($e=$-0.669\,38$\,J$) agrees
well with QMC calculations ($e=$-0.669\,436(7)$\,J$), and the extrapolated
series expansion results ($e=$0.6693(1)$\,J$). We recall, however, that the
self-consistent spin-wave theory at $\lambda=1$ already yields
$e=$-0.670\,421$\,J$, fairly close to the above values in spite of its mean-
field character. This shows that the ground-state energy is a very robust
quantity, accessible by many approaches. This is different for the dispersion
which we consider next.

###  IV.3 Dispersion

The results of the CST for the dispersion along a high-symmetry path through
the MBZ for various $\lambda$ and $L=16$ are shown in Fig. 4. At $\lambda=0$,
the XXZ model is reduced to the antiferromagnetic Ising model with a
completely flat dispersion $\omega(\bm{k})=2J$. For $\lambda=1$, the XXZ model
coincides with the antiferromagnetic Heisenberg model which was already
studied using the CST method in Refs. [4, 5]. At this isotropic point, the
dispersion features gapless Goldstone bosons and a distinct minimum at
$\bm{k}=(\pi,0)$, called the roton minimum. For $0<\lambda<1$, a continuous
evolution between these two limits takes place. Spin gap and roton minimum are
discussed below.

Figure 4: CST result of the dispersion $\omega(\bm{k})$ for $N_{0}$ and $L=16$
for various $\lambda$. For $\lambda=0$, the XXZ model is the Ising model with
a flat dispersion while it is the antiferromagnetic Heisenberg model for
$\lambda=1$ with a gapless spectrum and a distinct roton minimum at
$\bm{k}=(\pi,0)$ [4, 7]. For $0<\lambda<1$ the CST results interpolate
smoothly between these two limits.

####  IV.3.1 Spin gap

Fig. 5 shows the extrapolated CST results for the one-magnon gap $\Delta$
compared with previous results. For small $\lambda$, we compare to high-order
series expansion calculation [10, 11] and to data obtained by the Coupled
Cluster Method (CCM) [42]. For large $\lambda\lesssim 1$, the asymptotic power
law in the critical region derived from series expansion results [10] and
spin-wave theory [12] are presented. The CST results interpolate smoothly
between these two limits, capturing both limits quantitatively. The inset in
Fig. 5 shows a double logarithmic plot of the gap as a function of
$1-\lambda^{2}$ in the critical region, i.e., for $\lambda\rightarrow 1$. The
last data point for $N_{0}$ and $L=22$ is $\lambda=0.99925$; for higher values
no reliable convergence was achieved. We attribute this to the fact that for
values of $\lambda$ even closer to the isotropic point the cluster size $L$ is
not large enough to capture the relevant physics, see also Fig. 6 in Sec.
IV.3.2.

Figure 5: CST results for the one-magnon gap $\Delta=\omega(\bm{k}=0)$ of the
XXZ model for $0\leq\lambda\leq 1$ compared to third order spin wave theory
(SWT) [12, 21], data from Coupled Cluster Method (CCM) [42], results from
series expansion (SE) about the Ising limit [10, 11], and the critical power
law extracted from SE [10]. A quantitative comparison of the critical behavior
is shown in the inset underlining very good agreement.

Our error estimate (see also inset of Fig. 5) grows for larger values of
$\lambda$ implying larger correlation lengths in the vicinity of the isotropic
point. For a quantitative analysis, we fit

$\displaystyle\Delta(\lambda\approx 1)=c\left(1-\lambda^{2}\right)^{\nu z},$
(8)

which is shown in the inset. We obtain $c=$1.30(2)$$ and $\nu z=$0.498(2)$$
which agrees within 5% with the value $c=$1.265(2)$$ obtained from series
expansion [10] if the critical exponent $\nu z=1/2$ was fixed. For fit with
fixed exponent of $\nu z=1/2$ in our results we see a slight increase to
$c=$1.313(1)$$ which still agrees within 5% to the literature. The agreement
between all approaches is very good for almost all values of $\lambda$. For
smaller values of $\lambda$ the various approaches agree excellently. Only
close to the vanishing of the gap, the CCM data deviates from the other
results which is likely to be linked to the attainable cluster sizes. The
agreement between the extrapolated series expansion and the CST results is
remarkable since the latter are obtained on a finite lattice of moderate size
$L\approx 20$ and extrapolated to the thermodynamic limit without imposing a
certain power law. This finding corroborates the applicability of CST in
momentum space to describe phase transitions with closing gaps.

####  IV.3.2 Correlation length $\xi$

Fig. 6 shows inverse correlation lengths $1/\xi$ along the used fit-axes from
$(0,0)$ to $(\pi,0)$ for different values $\lambda$. Note that we use units
$1/a=1$, since the lattice constant $a$ is set to $1$. Additionally, we show a
horizontal dashed line for the relevant length scale of our largest system
$L=22$, namely the shortest occurring wrap-around $\sqrt{2}L$, in this case,
in diagonal direction. In the inset of Fig. 6 a zoom to the domain
$\lambda>$0.993$$ is shown.

For values of $\lambda$ up to $0.995$ we find finite correlation lengths, well
below $\sqrt{2}L$. For $\lambda>$0.995$$, $1/\xi$ starts to vanish, which is
expected since the system at $\lambda=1$ is gapless with infinite correlation
length. The error bar for values of $\lambda>0.9985$ crosses the threshold
$1/\sqrt{2}L$ arising from the finite system size. Additionally, for the
largest achieved value of $\lambda=0.99925$, $1/\xi$ itself approaches
$1/\sqrt{2}L$.

We conclude that for $\lambda\leq$0.995$$, finite-size effects are not
expected. For $\lambda\geq 0.9985$ the finite size of the system might play a
noticeable role. However, $\xi$ never exceeds these relevant length scales
dramatically. Additionally, the extrapolated results for closing of the gap
are in agreement with the literature and the CST results for the Heisenberg
model in Ref. [4, 5] performed for $L=16$ are consistent with results in the
literature. Hence, it appears that the influence of the finite size is not
severe in the CST approach.

The error estimate from $\lambda=0$ to $\lambda=0.95$ lies below
$0.1\text{\,}\%$, below $1\text{\,}\%$ for $\lambda\leq 0.97$, below
$5\text{\,}\%$ for $\lambda\leq 0.99$ and it is larger than $10\text{\,}\%$
for $\lambda\geq 0.995$.

Figure 6:  CST results for the ground-state correlation length $\xi$ for
$0.7\leq\lambda\leq 1$. A comparison of the correlation length with the
shortest wrap around $\sqrt{2}L$ in the largest cluster with $L=22$ (dashed
line) is shown.

####  IV.3.3 Roton mode

Fig. 7 depicts the CST results for both the roton minimum and the dispersion
maximum as a function of $\lambda$. Additionally, results from series
expansion [10] and DMRG [7] are shown. Also, QMC data [15] are available for
the isotropic point. First, we focus on the dispersion maximum at
$\bm{k}=(\pi/2,\pi/2)$ which displays a monotonic increase in $\lambda$ for
CST, the series expansion and the DMRG data. All data agrees very well. At the
isotropic point, similar values are obtained: CST ($2.369\,72(5)$ $J$), series
expansion ($2.385(1)$ $J$), QMC($2.39$ $J$), and DMRG($2.40$ $J$).

Figure 7: CST results for the roton mininum $\bm{k}=(\pi,0)$ and the
dispersion maximum $\bm{k}=(\pi/2,\pi/2)$ for $0\leq\lambda\leq 1$ compared to
data from series expansion (SE) [10], DMRG [7], and QMC [15]. All methods
agree well for dispersion maximum. For the roton mininum, all methods predict
an inflection point for $\lambda\gtrsim 0.8$. The values for the depth of the
roton mode differ slightly.

For the roton minimum, CST, series expansion, and DMRG show qualitatively
similar results: first the dispersion increases and for larger $\lambda\gtrsim
0.8$ an inflection point occurs. This inflection point is consistent with the
observation in Ref. [4, 5] that the hybridization between the one-magnon state
and the three-magnon scattering states causes the dip due to level repulsion.
This mechanism matters more and more for higher $\lambda$ when the dispersion
and the continuum approach each other. The numerical values for the roton dip
for $\lambda=1$ are $2.211(3)$ $J$ for the CST, $2.18(1)$ $J$ for series
expansion, $2.16$ $J$ for QMC, and $2.06(1)$ $J$ for DMRG. Hence, these values
are reasonably close, but further spread than the maximum values. We assume
that the SE and QMC results reflect the true value best while CST is still a
bit too high because of truncation beyond the quartic monomials. The DMRG data
results from a cylindrical geometry with fixed linear extension $L=10$. A
comprehensive finite-size scaling has therefore could not been performed. All
in all, we conclude that a consistent description for the high-energy
properties of the dispersion of the XXZ model is reached by CST except for a
slight deviation of about $2\text{\,}\mathrm{\char 37\relax}$ at the roton
minimum.

###  IV.4 Bound states and two-magnon continuum

A particular asset of the systematic basis changes by means of CST is that
bound states can be directly addressed and computed. The conservation of the
particle number, here the magnon number, is achieved by the basis change. It
allows one to compute the bound states in the two-particle sector. Recently,
even three-particle bound states could be addressed in antiferromagnetic spin
ladders which are induced by three-particle irreducible interactions, i.e., by
hexatic terms in second quantization [43]. Here, we focus on two-magnon bound
states as established perturbatively for not too large values of $\lambda$
[13, 11].

For bound states to be infinitely long-lived they may not decay into
scattering states. Hence, their energy may not overlap with the corresponding
continuum; it should stay below the lower boundary of the continuum.
Generically, the energetically lowest lying bound states occur at total
momentum zero $\mathbf{K}=(0,0)$ corresponding to the bound state at rest.
This is confirmed by series expansions [13, 11]. Hence, we compute four bound
states for zero center-of-mass momentum $\bm{K}=\bm{k}_{1}+\bm{k}_{2}=(0,0)$.
In addition, the strongest binding occurs for zero total $S^{z}$, i.e.,
between an $\alpha$\- and a $\beta$-magnon living on the A- and B-sublattice,
respectively. Again, this is confirmed perturbatively.

In order to distinguish the four states and to avoid numerical inaccuracies we
exploit the point group symmetry of rotations $R$ by $90^{\circ}$ about any
site of the lattice; the effective Hamiltonian is block-diagonal within each
eigen subspace of $R$. Four eigen values of $R$ are possible and do occur:
$\pm 1$ and $\pm i$. Since the XXZ model is time-reversal invariant the
imaginary eigen values are degenerate and their complex eigen vectors are
complex conjugates. The numerical diagonalization within each eigen subspace
of $R$ yields the energies of the bound states as depicted in Fig. 8.

Figure 8: The upper panel shows the energies of the four magnon-magnon bound
states $\tau_{i}$ calculated by CST as a function of the anisotropy $\lambda$
and the corresponding lower edge $2\Delta$ of the two-magnon continuum. The
inset compares data from truncated CST to data from perturbative CUT up to
terms $\propto\lambda^{4}$ [11]. One discerns a certain small deviation in
quadratic order. In the lower panel, the inverse participation ratio IPR of
the four bound states is plotted; only three curves are shown because the
eigen wave functions of the two degenerate bound states are complex conjugates
with the same IPR $I$. The vertical dashed lines mark the decay points for the
corresponding bound state determined by a crossing of $I$ with the threshold
$I=10^{-4}$ (horizontal dashed line).

Starting in the Ising limit, we find four bound states $\tau_{i}$ below the
lower edge of the two-magnon continuum. For increasing $\lambda$ they are
absorbed in the continuum one after the other before the isotropic Heisenberg
model is reached. We determine the decay point where each bound state merges
with the continuum by computing and extrapolating its IPR as described in Sec.
III.3 and Sec. III.7. We recall that a finite IPR indicates binding while its
vanishing indicates delocalization, i.e., the bound states ceases to exist.
The results are displayed in the lower panel of Fig. 8. It is advantageous to
determine the merging into the continuum from the vanishing of the IPR because
this vanishing occurs as very rapid drop by many orders of magnitude while the
crossing of the binding energies with the lower continuum edge is an
intersection with small, rigorously even vanishing, angle. We choose a
numerical threshold for the IPR to distinguish between a bound state and a
delocalized state in the continuum and set it to $1\text{\times}{10}^{-4}$,
i.e., only for IPR $>$1\text{\times}{10}^{-4}$$ we deduce binding. This value
is the typical value we found for the IPR of generic states well within the
continuum, extrapolated to the infinite thermodynamic limit $N\to\infty$. The
absorption points are marked as dashed vertical lines.

We find that $\tau(-1)$ disappears for $\lambda\approx$0.565(5)$$ which is a
deviation of $\approx$5\text{\,}\mathrm{\char 37\relax}$$ to
$\lambda\approx$0.5401(1)$$ from the pCUT [11]. The merging of $\tau(\pm i)$
occurs at $\lambda\approx$0.72(1)$$ in CST; no results were given for these
bound states in Ref. [11]. For the curves for $\tau(\pm i)$ in the inset we
analyzed the hopping amplitudes given in Dusuel_2010,[0]Ref. [11] to calculate
the series up to order $\lambda^{4}$. For $\tau(+1)$ we find
$\lambda\approx$0.975(5)$$ with a deviation of only
$\approx$0.6\text{\,}\mathrm{\char 37\relax}$$ to the extrapolated pCUT value
of $\lambda\approx$0.97$$. We conclude that for almost all values
$\lambda\in[0,1]$ at least one bound state exists. But even the lowest bound
state $\tau(+1)$ dives into the continuum very close to the isotropic
Heisenberg model. We note that this finding is very similar to the decay of
bound states in the XXZ model on coupled two-leg ladders. Here QMC simulations
reveal the decay of the lowest bound state at $\lambda\approx$0.96$$ for
different inter-ladder couplings [18].

Comparing the CST result with the results from the perturbative series near
the Ising limit, see inset of Fig. 8, we discern a sizable deviation in order
of $\lambda^{2}$ for all bound states, in contrast to what we found for the
ground-state energy and the dispersion. These deviations are a consequence of
the truncation scheme neglecting hexatic operators. Indeed, one finds that for
a perturbative CST in second order of $\lambda$ in position space,
coefficients for quartic interactions of the type
$a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}b^{\dagger}_{j}b^{\phantom{\dagger}}_{j}$
are coupled to hexatic terms during the flow. Hence, these contributions are
not considered rigorously in the CST as implemented in this paper. Therefore,
the quantitative analysis for properties at intermediate values of the
anisotropy has to be taken with a grain of salt. Nevertheless, we observe that
the series expansion and the CST results show very similar behavior. In
particular, the results close to the Heisenberg point agree very well.
Obviously, the benefits of the truncation based on the scaling dimension
outweighs the caveat of not capturing all perturbative terms exactly in the
Ising limit. This observation holds an attractive promise for future
application of the CST combined with a truncation based on scaling.

##  V Summary and Outlook

We extended the systematic basis change by a continuous similarity
transformation as applied to the spin isotropic Heisenberg model in two
dimensions in Refs. [4, 5] to the two-dimensional XXZ model. This model
displays the transition from a spontaneously broken discrete symmetry with
finite energy gap to a spontaneously broken continuous symmetry with vanishing
energy gap. Starting point of our calculation was the Hamiltonian in self-
consistent spin wave theory based on the Dyson-Maleev representation of the
spin operators. The calculation is performed by numerically integrating the
flow equations in the prefactors of terms in second quantization. To this end,
we discretized the magnetic Brillouin zone in various ways and checked that
the results converged within small error bars. In addition, the proliferation
of terms was limited by a truncation based on the scaling dimension. We kept
all terms of quadratic and quartic form, but neglected hexatic and higher
terms after normal-ordering them with respect to the magnon vacuum.

We verified that this approach can be applied to the paradigmatic transition
of the two-dimensional XXZ model including the vanishing of the spin gap. The
ground-state energy, the one-magnon dispersion as well as the two-magnon
interaction are derived and computed. The dispersion displays the expected
appearance of the roton minimum and agrees well with all available previous
results. We successfully computed the two-magnon bound states at zero momentum
of the center-of-mass in the longitudinal channel, i.e., for vanishing total
$S^{z}$ component. Four bound states are identified of which two with complex
eigen wave functions are degenerate. The latter states have not been addressed
so far. Upon approaching the isotropic Heisenberg model they merge
successively with the two-magnon continuum. The lowest bound state is
rotationally invariant in space, i.e., with respect to rotations by
$90^{\circ}$, and vanishes only very close to the isotropic point in nice
quantitative agreement with previous perturbative results [11] within less
than one percent. This agreement is even more remarkable because the energies
of the bound states show a larger deviation of the order of five percent to
the perturbative results for intermediate values of $\lambda$. We trace this
effect back to the truncations at beyond quartic level. The degree of binding
and the delocalization upon diving into the continuum are determined reliably
by the inverse participation ratio.

The above sketched achievements hold the promise that the systematic CST with
truncation based on the scaling argument can be applied to many more models
which up to date have not been studied in great detail. We emphasize that the
closing and opening of gaps can be captured as well as the occurrence of
binding phenomena. Models suggesting themselves comprise the bilayer model of
two stacked square lattices which displays a transition to from the long-range
ordered phase with magnons to triplons on the interlayer bonds upon increasing
interlayer coupling [44, 45]. The frustrated $J_{1}$-$J_{2}$ model on the
square lattice is another intriguing candidate to be investigated its two
limits displaying long-range order [46, 47, 48, 49, 50, 8]. Antiferromagnets
displaying non-collinear long-range order represent another wide field of
application.

###### Acknowledgements.

We thank Frank Pollmann and Ruben Verresen for providing the DMRG results for
the roton modes and Stefan Wessel and Andres Sandvik for useful discussions.
We gratefully acknowledge financial support by the German Research Foundation
(DFG) in projects UH 90-13/1 (GSU), UH 90-14/1 (DBH) and SCHM 2511/13-1
(MRW/KPS).

## References

  * Sachdev [1999] S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, Cambridge, UK, 1999).
  * Manousakis [1991] E. Manousakis, The Spin-$\frac{1}{2}$ Heisenberg Antiferromagnet on a Square Lattice and its Application to the Cuprous Oxides, Rev. Mod. Phys. 63, 1 (1991).
  * Auerbach [1994] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_, Graduate Texts in Contemporary Physics (Springer, New York, 1994).
  * Powalski  _et al._ [2015] M. Powalski, G. S. Uhrig, and K. P. Schmidt, Roton Minimum as a Fingerprint of Magnon-Higgs Scattering in Ordered Quantum Antiferromagnets, Phys. Rev. Lett. 115, 207202 (2015).
  * Powalski  _et al._ [2018] M. Powalski, K. P. Schmidt, and G. S. Uhrig, Mutually attracting spin waves in the square-lattice quantum antiferromagnet, SciPost Phys. 4, 001 (2018).
  * Shao  _et al._ [2017] H. Shao, Y. Qin, S. Capponi, S. Chesi, Z. Meng, and A. Sandvik, Nearly deconfined spinon excitations in the square-lattice spin-$1/2$ heisenberg antiferromagnet, Phys. Rev. X 7, 041072 (2017).
  * Verresen  _et al._ [2018] R. Verresen, F. Pollmann, and R. Moessner, Quantum dynamics of the square-lattice Heisenberg model, Phys. Rev. B 98, 155102 (2018).
  * Ferrari and Becca [2018] F. Ferrari and F. Becca, Spectral signatures of fractionalization in the frustrated Heisenberg model on the square lattice, Phys. Rev. B 98, 100405 (2018).
  * Weihong  _et al._ [1991] Z. Weihong, J. Oitmaa, and C. J. Hamer, Square-lattice Heisenberg antiferromagnet at $T=0$, Physical Review B 43, 8321 (1991).
  * Zheng  _et al._ [2005] W. Zheng, J. Oitmaa, and C. J. Hamer, Series studies of the spin-$\frac{1}{2}$ Heisenberg antiferromagnet at $T=0$: Magnon dispersion and structure factors, Phys. Rev. B 71, 184440 (2005).
  * Dusuel  _et al._ [2010] S. Dusuel, M. Kamfor, K. P. Schmidt, R. Thomale, and J. Vidal, Bound states in two-dimensional spin systems near the Ising limit: A quantum finite-lattice study, Phys. Rev. B 81, 064412 (2010).
  * Hamer  _et al._ [1992] C. J. Hamer, Z. Weihong, and P. Arndt, Third-order spin-wave theory for the Heisenberg antiferromagnet, Phys. Rev. B 46, 6276 (1992).
  * Hamer [2009] C. J. Hamer, Energy spectrum of the two-magnon bound states in the Heisenberg-Ising antiferromagnet on the square lattice, Phys. Rev. B 79, 212413 (2009).
  * Sandvik [1997] A. W. Sandvik, Finite-size scaling of the ground-state parameters of the two-dimensional Heisenberg model, Phys. Rev. B 56, 11678 (1997).
  * Sandvik and Singh [2001] A. W. Sandvik and R. R. P. Singh, High-Energy Magnon Dispersion and Multimagnon Continuum in the Two-Dimensional Heisenberg Antiferromagnet, Phys. Rev. Lett. 86, 528 (2001).
  * Weidinger and Zwerger [2015] S. A. Weidinger and W. Zwerger, Higgs mode and magnon interactions in 2D quantum antiferromagnets from Raman scattering, Eur. Phys. J B 88, 237 (2015).
  * Hong  _et al._ [2017] T. Hong, M. Matsumoto, Y. Qiu, W. Chen, T. Gentile, S. Watson, F. Awwadi, M. Turnbull, S. Dissanayake, H. Agrawal, R. Toft-Petersen, B. Klemke, K. Coester, K. Schmidt, and D. Tennant, Higgs amplitude mode in a two-dimensional quantum antiferromagnet near the quantum critical point, Nat. Phys. 13, 638 (2017).
  * Ying  _et al._ [2019] T. Ying, K. P. Schmidt, and S. Wessel, Higgs Mode of Planar Coupled Spin Ladders and its Observation in ${\mathrm{C}}_{9}{\mathrm{H}}_{18}{\mathrm{N}}_{2}{\mathrm{CuBr}}_{4}$, Phys. Rev. Lett. 122, 127201 (2019).
  * Onsager [1944] L. Onsager, Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition, Phys. Rev. 65, 117 (1944).
  * J. Oitmaa and Zheng [2006] C. H. J. Oitmaa and W. Zheng,  _Series Expansion Methods for Strongly Interacting Lattice Models_ (Cambridge University Press, Cambridge, UK, 2006).
  * Syromyatnikov [2010] A. V. Syromyatnikov, Spectrum of short-wavelength magnons in a two-dimensional quantum Heisenberg antiferromagnet on a square lattice: third-order expansion in $1/S$, Journal of Physics: Condensed Matter 22, 216003 (2010).
  * Oguchi and Ishikawa [1973] T. Oguchi and T. Ishikawa, Theory of two-magnon bound states in a two-dimensional antiferromagnet, Journal of the Physical Society of Japan 34, 1486 (1973).
  * Kramer and MacKinnon [1993] B. Kramer and A. MacKinnon, Localization: theory and experiment, Reports on Progress in Physics 56, 1469 (1993).
  * Malki and Uhrig [2019] M. Malki and G. S. Uhrig, Delocalization of edge states in topological phases, EPL (Europhysics Letters) 127, 27001 (2019).
  * Wegner [1994] F. Wegner, Flow-equations for Hamiltonians, Ann Phys 506, 77 (1994).
  * Knetter and Uhrig [2000] C. Knetter and S. Uhrig, Perturbation theory by flow equations: dimerized and frustrated $S=1/2$ chain, Eur. Phys. J. B 13, 209 (2000).
  * Kehrein [2006] S. Kehrein,  _The Flow Equation Approach to Many-Particle Systems_ , Springer Tracts in Modern Physics, Vol. 217 (Springer, Berlin, 2006).
  * Mielke [1998] A. Mielke, Flow equations for band-matrices, Eur. Phys. J. B 5, 605 (1998).
  * Fischer  _et al._ [2010] T. Fischer, S. Duffe, and G. S. Uhrig, Adapted continuous unitary transformation to treat systems with quasi-particles of finite lifetime, New J. Phys. 10, 033048 (2010).
  * Kehrein [1999] S. Kehrein, Flow Equation Solution for the Weak- to Strong-Coupling Crossover in the Sine-Gordon Model, Phys. Rev. Lett. 83, 4914 (1999).
  * Kehrein [2001] S. Kehrein, Flow equation approach to the sine-Gordon model, Nuclear Physics B 592, 512 (2001).
  * Dyson [1956] F. J. Dyson, General Theory of Spin-Wave Interactions, Phys. Rev. 102, 1217 (1956).
  * Maleev [1958] S. Maleev, Multimagnon processes in the scattering of slow neutrons in ferro-magnets, Soviet Physics JETP 34 (1958).
  * Schmidt and Uhrig [2006] K. P. Schmidt and G. S. Uhrig, Hard-core magnons in the $S=1/2$ Heisenberg model on the square lattice, Phys. Rev. B 73, 172407 (2006).
  * Uhrig and Majumdar [2013] G. S. Uhrig and K. Majumdar, Varied perturbation theory for the dispersion dip in the two-dimensional Heisenberg quantum antiferromagnet, Eur. Phys. J. B 86, 1 (2013).
  * Takahashi [1989] M. Takahashi, Modified spin-wave theory of a square-lattice antiferromagnet, Phys. Rev. B 40, 2494 (1989).
  * Uhrig and Schulz [1996] G. S. Uhrig and H. J. Schulz, Magnetic excitations spectrum of dimerized antiferromagnetic chains, Phys. Rev. B 54, 9624(R) (1996).
  * Uhrig and Schulz [1998] G. S. Uhrig and H. J. Schulz, Erratum:“Magnetic excitations spectrum of dimerized antiferromagnetic chains”, Phys. Rev. B 58, 2900 (1998).
  * Zhitomirsky and Chernyshev [2013] M. E. Zhitomirsky and A. L. Chernyshev, Spontaneous magnon decays, Rev. Mod. Phys. 85, 219 (2013).
  * Ahnert and Mulansky [2011] K. Ahnert and M. Mulansky, Odeint-solving ordinary differential equations in C++ (2011).
  * Okunishi  _et al._ [2001] K. Okunishi, Y. Akutsu, N. Akutsu, and T. Yamamoto, Universal relation between the dispersion curve and the ground-state correlation length in one-dimensional antiferromagnetic quantum spin systems, Phys. Rev. B 64, 104432 (2001).
  * Bishop  _et al._ [2017] R. Bishop, P. Li, R. Zinke, R. Darradi, J. Richter, D. Farnell, and J. Schulenburg, The spin-half XXZ antiferromagnet on the square lattice revisited: A high-order coupled cluster treatment, J. Magn. Magn. Mater. 428, 178 (2017).
  * Schmiedinghoff  _et al._ [2022] G. Schmiedinghoff, L. Müller, G. S. Uhrig, U. Kumar, and B. Fauseweh, Three-body bound states in antiferromagnetic spin ladders, Comm. Phys. 5, 218 (2022).
  * Wang  _et al._ [2006] L. Wang, K. Beach, and A. Sandvik, High-precision finite-size scaling analysis of the quantum-critical point of $S=1/2$ Heisenberg antiferromagnetic bilayers, Phys. Rev. B 73, 014431 (2006).
  * Lohöfer  _et al._ [2015] M. Lohöfer, T. Coletta, D. Joshi, F. Assaad, M. Vojta, S. Wessel, and F. Mila, Dynamical structure factors and excitation modes of the bilayer Heisenberg model, Phys. Rev. B 92, 245137 (2015).
  * Uhrig  _et al._ [2009] G. S. Uhrig, M. Holt, J. Oitmaa, O. P. Sushkov, and R. Singh, Pnictides as frustrated quantum antiferromagnet close to a quantum phase transition, Phys. Rev. B 79, 092416 (2009).
  * Jiang  _et al._ [2012] H.-C. Jiang, H. Yao, and L. Balents, Spin liquid ground state of the spin-$\frac{1}{2}$ square $J_{1}$-$J_{2}$ Heisenberg model, Phys. Rev. B 86, 024424 (2012).
  * Gong  _et al._ [2014] S.-S. Gong, W. Zhu, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Plaquette Ordered Phase and Quantum Phase Diagram in the Spin-$\frac{1}{2}$ $J_{1}$-$J_{2}$ Square Heisenberg Model, Phys. Rev. Lett. 113, 027201 (2014).
  * Richter  _et al._ [2015] J. Richter, R. Zinke, and D. J. Farnell, The spin-1/2 square-lattice $J_{1}$-$J_{2}$ model: the spin-gap issue, Eur. Phys. J. B 88, 2 (2015).
  * Morita  _et al._ [2015] S. Morita, R. Kaneko, and M. Imada, Quantum Spin Liquid in Spin $1/2$ $J_{1}$–$J_{2}$ Heisenberg Model on Square Lattice: Many-Variable Variational Monte Carlo Study Combined with Quantum-Number Projections, J. Phys. Soc. Japan 84, 024720 (2015).

##  Appendix A Mean-field decoupling in real space

Here we provide the explicit expressions for each step of the derivation of
the starting point of the CST, i.e., for determining the self-consistent mean-
field solution. Apart from the spin anisotropy $\lambda$, the calculations are
analogous to the ones in Refs. [4, 5, 35]. The Hamiltonian in Dyson-Maleev
representation reads

$\displaystyle\mathcal{H}$
$\displaystyle=J\sum_{i\in\Gamma_{A},\delta}\left[-S^{2}\right.$
$\displaystyle+S\left(a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}+b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}+\lambda
a^{\phantom{\dagger}}_{i}b^{\phantom{\dagger}}_{i+\delta}+\lambda
a^{\dagger}_{i}b^{\dagger}_{i+\delta}\right)$ (9)
$\displaystyle-\left.{a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}}-\frac{\lambda}{2}{a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}a^{\phantom{\dagger}}_{i}b^{\phantom{\dagger}}_{i+\delta}}-\frac{\lambda}{2}{a^{\dagger}_{i}b^{\dagger}_{i+\delta}b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}}\right],$

where we switched the sum over all adjacent pairs$\langle i,j\rangle$ in Eq.
(1) to a sum over all sites in sublattice A, i.e., $\Gamma_{A}$, and their
nearest neighbors on the square lattice reached via the distances $\delta$.
The first line represents the classical energy, the second the bilinear one-
magnon terms and the last line the quadrilinear two-magnon terms. Next, we
perform the standard mean-field decoupling [36] in real space

$\displaystyle
b^{(\dagger)}_{i}\tilde{b}^{(\dagger)}_{j}=:{b^{(\dagger)}_{i}\tilde{b}^{(\dagger)}_{j}}:+\Braket{b^{(\dagger)}_{i}\tilde{b}^{(\dagger)}_{j}}_{0},$
(10)

where $:{...}:$ indicates normal-ordered operators and $\Braket{...}_{0}$ are
the expectation values in the mean-field ground-state. For quadrilinear terms
we apply the corresponding relations based on Wick’s theorem.

To make the resulting expression more concise, we use exact relations and
conservation laws. First, the Hamiltonain is invariant under the exchange of
the A and B sublattice. For the Dyson-Maleev bosons this translates into a
complex conjugation of the prefactors and swapping
$a_{i}^{\dagger}\leftrightarrow b_{i+\delta}$ and $a_{i}\leftrightarrow
b_{i+\delta}^{\dagger}$. We define the expectation values

$\displaystyle n$
$\displaystyle\coloneqq\Braket{a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}}_{0}=\Braket{b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}}_{0}$
(11a) $\displaystyle\Delta$
$\displaystyle\coloneqq\Braket{a^{\phantom{\dagger}}_{i}b^{\phantom{\dagger}}_{i+\delta}}_{0}=\Braket{a^{\dagger}_{i}b^{\dagger}_{i+\delta}}_{0},$
(11b)

which are determined self-consistently below.

Second, the Hamiltonian (1) conserves the total spin component

$S^{z}_{\rm
tot}=\sum_{i\in\Gamma_{A}}(S^{z}_{i}+S^{z}_{i+\delta})=\sum_{i\in\Gamma_{A}}(a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}-b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}).$
(12)

From this conservation we can infer

$\displaystyle\Braket{a^{\phantom{\dagger}}_{i}a^{\phantom{\dagger}}_{i}}_{0}=\Braket{b^{\dagger}_{i+\delta}b^{\dagger}_{i+\delta}}_{0}=\Braket{a^{\phantom{\dagger}}_{i}b^{\dagger}_{i+\delta}}_{0}=\Braket{a^{\dagger}_{i}b^{\phantom{\dagger}}_{i+\delta}}_{0}=0.$
(13)

In total, the mean-field decoupled Hamiltonian reads

$\displaystyle\mathcal{H}$
$\displaystyle=J\sum_{i\in\Gamma_{A},\delta}\Big{[}-S^{2}+2Sn+2S\lambda\Delta-n^{2}-\Delta^{2}-2\lambda
n\Delta$
$\displaystyle\qquad+\left(S-n-\lambda\Delta\right)\big{(}:{a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}}:+:{b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}}:\big{)}$
$\displaystyle\qquad+\left(\lambda S-\lambda
n-\Delta\right)\big{(}:{a^{\phantom{\dagger}}_{i}b^{\phantom{\dagger}}_{i+\delta}}:+:{a^{\dagger}_{i}b^{\dagger}_{i+\delta}}:\big{)}$
$\displaystyle\qquad-:{a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}}:$
(14)
$\displaystyle\qquad-\frac{\lambda}{2}:{a^{\dagger}_{i}a^{\phantom{\dagger}}_{i}a^{\phantom{\dagger}}_{i}b^{\phantom{\dagger}}_{i+\delta}}:-\frac{\lambda}{2}:{a^{\dagger}_{i}b^{\dagger}_{i+\delta}b^{\dagger}_{i+\delta}b^{\phantom{\dagger}}_{i+\delta}}:\Big{]}.$

The first line is the zero-magnon part, the second and third line the one-
magnon part and the last two lines comprise the two-magnon terms. We stress
that no approximation is performed. Next, we apply a Fourier transformation
yielding

$\displaystyle\mathcal{H}$
$\displaystyle=J\Big{[}N\left(-S^{2}+2Sn+2S\lambda\Delta-n^{2}-\Delta^{2}-2\lambda
n\Delta\right)$
$\displaystyle+\sum_{\bm{k}\in\text{MBZ}}A_{\bm{k}}\big{(}:{a^{\dagger}_{\bm{k}}a^{\phantom{\dagger}}_{\bm{k}}}:+:{b^{\dagger}_{\bm{k}}b^{\phantom{\dagger}}_{\bm{k}}}:\big{)}$
$\displaystyle+\sum_{\bm{k}\in\text{MBZ}}B_{\bm{k}}\big{(}:{a^{\phantom{\dagger}}_{\bm{k}}b^{\phantom{\dagger}}_{-\bm{k}}}:+:{a^{\dagger}_{\bm{k}}b^{\dagger}_{-\bm{k}}}:\big{)}$
(15)
$\displaystyle-\frac{Z}{N}\sum_{\bm{1},\bm{2},\bm{3},\bm{4}\in\text{MBZ}}\delta(\bm{1}\\!-\\!\bm{2}\\!+\\!\bm{3}\\!-\\!\bm{4})\gamma(\bm{3}\\!-\\!\bm{4}):{a^{\dagger}_{\bm{1}}a^{\phantom{\dagger}}_{\bm{2}}b^{\dagger}_{\bm{3}}b^{\phantom{\dagger}}_{\bm{4}}}:$
$\displaystyle-\frac{\lambda
Z}{2N}\sum_{\bm{1},\bm{2},\bm{3},\bm{4}\in\text{MBZ}}\delta(\bm{1}\\!-\\!\bm{2}\\!-\\!\bm{3}\\!-\\!\bm{4})\gamma(\bm{4}):{a^{\dagger}_{\bm{1}}a^{\phantom{\dagger}}_{\bm{2}}a^{\phantom{\dagger}}_{\bm{3}}b^{\phantom{\dagger}}_{\bm{4}}}:$
$\displaystyle-\frac{\lambda
Z}{2N}\\!\\!\\!\sum_{\bm{1},\bm{2},\bm{3},\bm{4}\in\text{MBZ}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\delta(\bm{1}\\!+\\!\bm{2}\\!+\\!\bm{3}\\!-\\!\bm{4})\gamma(\bm{4}\\!-\\!\bm{2}\\!-\\!\bm{3}):{a^{\dagger}_{\bm{1}}b^{\dagger}_{\bm{2}}b^{\dagger}_{\bm{3}}b^{\phantom{\dagger}}_{\bm{4}}}:\\!\Big{]},$

where we use the shorthand $\bm{j}$ for $\bm{k}_{j}$. The Kronecker $\delta$
symbols hold modulo reciprocal lattice vectors. If a term describes an Umklapp
process, i.e., a reciprocal lattice vector $\bm{G}$ is needed to fulfill
momentum conservation, the Kronecker $\delta$ takes the value
$\gamma(\bm{G})=\pm 1$. In the above equation, we used the definitions

$\displaystyle\gamma(\bm{k})$
$\displaystyle\coloneqq\frac{1}{Z}\sum_{\delta}{\text{e}}^{-\text{i}\delta\bm{k}}=\frac{2}{Z}\left[\cos(k_{x})+\cos(k_{y})\right]$
(16a) $\displaystyle A_{\bm{k}}$ $\displaystyle\coloneqq
Z\left(S-n-\lambda\Delta\right)$ (16b) $\displaystyle B_{\bm{k}}$
$\displaystyle\coloneqq Z\gamma(\bm{k})\left(\lambda S-\lambda
n-\Delta\right),$ (16c)

where we introduced the coordination number $Z=4$. Next, we diagonalize the
bilinear part in Eq. (A) by the Bogoliubov transformation

$\displaystyle a^{\dagger}_{\bm{k}}$
$\displaystyle=l_{\bm{k}}\alpha^{\dagger}_{\bm{k}}+m_{\bm{k}}\beta^{\phantom{\dagger}}_{-\bm{k}}$
(17a) $\displaystyle b^{\dagger}_{\bm{k}}$
$\displaystyle=m_{\bm{k}}\alpha^{\phantom{\dagger}}_{\bm{-k}}+l_{\bm{k}}\beta^{\dagger}_{\bm{k}}$
(17b)

with

$l_{\bm{k}}^{2}-m_{\bm{k}}^{2}=1.$ (18)

Explicitly, we parametrize as in Refs. [4, 5, 35]

$\displaystyle\mu_{\bm{k}}$
$\displaystyle\coloneqq\sqrt{1-\left(\frac{B_{\bm{k}}}{A_{\bm{k}}}\right)^{2}}$
(19a) $\displaystyle l_{\bm{k}}$
$\displaystyle\coloneqq\sqrt{\frac{1-\mu_{\bm{k}}}{2\mu_{\bm{k}}}}$ (19b)
$\displaystyle m_{\bm{k}}$
$\displaystyle\coloneqq-\operatorname{sign}(\gamma(\bm{k}))\sqrt{\frac{1+\mu_{\bm{k}}}{2\mu_{\bm{k}}}}=-\operatorname{sign}(\gamma(\bm{k}))x_{\bm{k}}l_{\bm{k}}$
(19c) $\displaystyle x_{\bm{k}}$
$\displaystyle\coloneqq\sqrt{\frac{1+\mu_{\bm{k}}}{1-\mu_{\bm{k}}}}.$ (19d)

With these steps we rewrite the Hamiltonian (2) in the form

$\displaystyle\mathcal{H}=E_{0}+\sum_{\bm{k}}\omega_{0}(\bm{k})\big{(}:{\alpha^{\dagger}_{\bm{k}}\alpha^{\phantom{\dagger}}_{\bm{k}}}:+:{\beta^{\dagger}_{\bm{k}}\beta^{\phantom{\dagger}}_{\bm{k}}}:\big{)}+\Gamma_{0}+\mathcal{V}_{0},$
(20)

with

$\displaystyle E_{0}$
$\displaystyle=JN\left(-S^{2}+2Sn+2S\lambda\Delta-n^{2}-\Delta^{2}-2\lambda
n\Delta\right)$ (21a) $\displaystyle\omega_{0}(\bm{k})$
$\displaystyle=A_{\bm{k}}\mu_{\bm{k}}$ (21b) $\displaystyle\Gamma_{0}$
$\displaystyle=\sum_{\bm{k}\in\text{MBZ}}\left(2A_{\bm{k}}l_{\bm{k}}m_{\bm{k}}+B_{\bm{k}}(m_{\bm{k}}^{2}+l_{\bm{k}}^{2})\right)$
$\displaystyle\qquad\qquad\times\big{(}:{\alpha^{\phantom{\dagger}}_{\bm{k}}\beta^{\phantom{\dagger}}_{-\bm{k}}}:+:{\alpha^{\dagger}_{\bm{k}}\beta^{\dagger}_{-\bm{k}}}:\big{)}.$
(21c)

The term $\mathcal{V}_{0}$ comprises all the quadrilinear two-magnon terms. We
do not provide here explicitly because is already published, for instance in
Ref. [35], except that we denote the normal-ordering of the quadrilinear terms
by $:{...}:$, and add a factor $\lambda$ in front of $\gamma$ with one or a
sum of three momenta as arguments. For example the term $x_{1}\gamma(1)$
becomes $\lambda x_{1}\gamma(1)$, the term $x_{1}x_{3}x_{4}\gamma(1-3-4)$
becomes $\lambda x_{1}x_{3}x_{4}\gamma(1-3-4)$, while $x_{1}x_{3}\gamma(1-3)$
stays $x_{1}x_{3}\gamma(1-3)$. This can be directly traced back to the
addtional $\lambda$ factors in the quartic part of Eq. A compared to the
Heiseberg model.

In order to retrieve the values of $n$ and $\Delta$, we exploit translational
invariance and finally perform the Fourier and Bogoliubov transformation on
the right hand sides of Eqs. (11), evaluate the vacuum expectation value for
the Bogoliubov bosons $\alpha^{\phantom{\dagger}}_{\bm{k}}$ and
$\beta^{\phantom{\dagger}}_{\bm{k}}$ to obtain the self-consistency conditions

$\displaystyle n$
$\displaystyle=\frac{1}{N}\sum_{\bm{k}}\Braket{a^{\dagger}_{\bm{k}}a^{\phantom{\dagger}}_{\bm{k}}}_{0}=\frac{1}{N}\sum_{\bm{k}}l_{\bm{k}}^{2}$
(22a) $\displaystyle\Delta$
$\displaystyle=\frac{1}{N}\sum_{\bm{k}}\gamma(\bm{k})\Braket{a^{\dagger}_{\bm{k}}b^{\phantom{\dagger}}_{\bm{k}}}_{0}=\frac{1}{N}\sum_{\bm{k}}\gamma(\bm{k})l_{\bm{k}}m_{\bm{k}}.$
(22b)

They are solved numerically by iteration to convergence for $N_{\emptyset}$
mesh with $L=30001$ sites with a tolerance $10^{-13}$.