Quantum dynamics of the square lattice Heisenberg model
Abstract
Despite nearly a century of study of the Heisenberg model on the square lattice, there is still disagreement on the nature of its highenergy excitations. By tuning toward the Heisenberg model from the exactly soluble Ising limit, we find that the strongly attractive magnon interactions of the latter naturally account for a number of spectral features of the Heisenberg model. This claim is backed up both numerically and analytically. Using the density matrix renormalization group method, we obtain the dynamical structure factor for a cylindrical geometry, allowing us to continuously connect both limits. Remarkably, a semiquantitative description of certain observed features arises already at the lowest nontrivial order in perturbation theory around the Ising limit. Moreover, our analysis uncovers that highenergy magnons are localized on a single sublattice, which is related to the entanglement properties of the ground state.
At its ripe age of 90 yearsHeisenberg (1928), the square lattice antiferromagnetic Heisenberg model has had its dynamical properties studied intensely. Spin wave theory was in large part developed to investigate this model’s lowenergy propertiesHulthén (1936); Anderson (1952). Its anomalous terms lead to a ‘relativistic’ (linear) lowenergy dispersion, related to coupling the two sublattices of the ground state’s spontaneous Néel orderingKeffer et al. (1953). While evading an exact treatment, its dynamics has been studied numerically via quantum Monte Carlo simulationsSandvik and Singh (2001); Shao et al. (2017) and exact diagonalizationLüscher and Läuchli (2009); also several highorder perturbative expansions have been devised, such as around the Ising limitSingh (1989); Gelfand (1996); Singh and Gelfand (1995); Zheng et al. (2005) or via a continuous unitary transformations (CUT)Powalski et al. (2015, 2018). In addition, a number of adhoc approaches have been motivated by a range of different physical pictures. Moreover, this model is related to experimental systems, not least to the parent states of the cupratesAnderson (1987); Rice (1987); Chakravarty et al. (1988); Ho et al. (2001), the study of which has led to much of the modern theory of quantum magnetism.
It may thus seem all the more surprising that there is still no consensus on the appropriate physical picture for certain regions in momentum space. Proposals include stronglyinteracting magnonsPowalski et al. (2015, 2018), deconfined spinonsHsu (1990); Ho et al. (2001); Syljuåsen and Lee (2002); Dalla Piazza et al. (2015); Ferrari and Becca (2018), all the way to a connection to deconfined quantum criticality Shao et al. (2017). The disagreement is not limited to the underlying physical mechanism, but also pertains to quantitative aspects of spectral properties. One uniting factor, at least, is a need to go beyond a perturbativelydressed singlemagnon picture.
Our intention is not to propose yet another scenario. Rather, we adopt the perspective that away from the unassailable hydrodynamic limit—accounting for the lowenergy Goldstone modesKeffer et al. (1953); Nielsen and Chadha (1976)—other features which have caught the attention of the community may not even be uniquely described by one picture as opposed to another. Instead, we seek to provide a simple account of salient features of the intermediate and highenergy part of the spectrum.
Perhaps the most controversial region concerns magnons with momenta . Spin wave theory predicts that these are dispersionless, in disagreement with both experimentRønnow et al. (2001); Christensen et al. (2007); Dalla Piazza et al. (2015) and theoretical methods. What is instead observed, is a local mininum at —commonly referred as the roton mode, in analogy with the quasiparticle dispersion in liquid HeliumLandau (1941). Moreover, spin wave theory is unable to account for the high spectral weight in the continuum just above this mode.
In this paper, we advance along two complementary tracks. First, we determine the dynamical structure factor using a method based on the density matrix renormalization group (DMRG)White (1992, 1993); Schollwöck (2011); Stoudenmire and White (2013); Zaletel et al. (2015); Gohlke et al. (2017), with systematic errors distinct from those of previous approaches. This novel method has at least two welcome features: it confirms that the phenomenology of the roton mode is indeed beyond the dressed singlemagnon picture, and it uncovers a hithertounrecognized property of of magnons with , which we refer to as sublatticelocalization. We also clarify how the latter is related to the entanglement of the ground state.
Second, we use existing data from an Ising expansion developed by Singh and GelfandSingh (1989); Gelfand (1996); Singh and Gelfand (1995)—and pushed further by Zheng, Oitmaa and HamerZheng et al. (2005)—to point out that some known results at the isotropic point are already semiquantitatively accounted for by the lowest nontrivial order. Moreover, our numerical method allows us to study such an XXZ model (with dominant easyaxis anisotropy) without any perturbative approximation.
The main message of our paper is that aspects of the attractive magnon interactions, i.e. the physics beyond spin wave theory, arise naturally from domainwallcounting in the Ising limit, sometimes connecting all the way to the isotropic Heisenberg point. In particular, the numerics shows this at a phenomenological level, but a simple perturbative calculation also sheds light on, e.g., the smallyetnonzero magnitude and shape of the roton mode’s dispersion. Moreover, even the aforementioned phenomenon of sublatticelocalization can be accounted for within a loworder perturbative picture. In addition, we provide a quantitative analysis of the roton mode.
The remainder of this paper is structured as follows. In section I we give a brief overview of the model’s salient features, relating them to previous literature whenever possible. The spectral functions obtained using DMRG are shown in section II: first for the Heisenberg model, which is then connected to the Ising limit. Section III supplements this by showing how various features, such as the roton minimum or sublatticelocalization, naturally arise within a loworder perturbative picture. The apparently hithertounexplored phenomenon of sublatticelocalization is studied numerically in section IV, emphasizing its link to entanglement (or absence thereof). Section V contains a quantitative analysis of the roton mode with comparison to results from the literature.
I Square lattice Heisenberg model
We study the spin antiferromagnetic Heisenberg model (AFH) on the square lattice, allowing for easyaxis anisotropy:
(1) 
where . We are principally interested in the isotropic point , where the Néel order of the ground state spontaneously breaks the symmetry down to a group generated by (where we define the ordering direction to be along the spin axis). As we will argue, it is also useful to consider , where the model is in a gapped Ising phase which spontaneously breaks the symmetry .
i.1 Dynamical structure factor and quantum numbers
Spectral functions give direct insight into the properties of excitations. In this work we focus on the dynamical structure factor, which is experimentally accessible through, for example, inelastic neutron scattering. It can be expressed in terms of the dynamical correlation functions :
(2) 
which is normalized as . We focus on the transverse spectral function:
(3) 
This object gives direct insight into the excitations above the ground state. If , one can show^{1}^{1}1To argue that , one can use that has a welldefined quantum number with respect to the modified translation operator . that
(4) 
where . It is natural to choose a basis , where is the momentum with respect to the translation symmetry of the twosite unit cell. Eq. (4) tells us that the spectral function gives information about the existence of energy eigenstates with momentum and and .
Note that when labeling states, lives in the reduced (magnetic) Brillouin zone, , but the spectral function itself is periodic only with respect to the original (lattice) Brillouin zone, (taking the lattice constant to be unity).
i.2 Spin wave theory
In terms of the above quantum numbers, spin wave theory predicts two bands^{2}^{2}2One could instead work in a local frame where there is only one band, but this obscures some of the physics. In particular, then is no longer a good quantum number.. These exactly coincide and are distinguished by . The dispersion relation to order^{3}^{3}3To wit, the HolsteinPrimakoff transformation is . The dominant term of the square root is of order , which leads to LSWT. , i.e. linear spin wave theory (LSWT), isAnderson (1952); Kubo (1952)
(5) 
Hence for , there are two linearlydispersing Goldstone modes at the zone center (in sectors ), consistent with two of three generators of being spontaneously brokenNielsen and Chadha (1976). However, based on general sum rulesStringari (1994), it is known that as , such that the Goldstone modes will have vanishing intensity in the transverse spectral function at the zone center. Instead, they show up near the ordering wavevector M , since the same sum rules imply an (integrable) divergence as .
The first order corrections to the dispersion within spin wave theory areOguchi (1960); Weihong et al. (1991)
(6) 
where and are (dependent) constants^{4}^{4}4 and , where the integration is over .. At the isotropic point, this correction is only a momentumindependent rescaling. Higherorder corrections ( and ) are also knownHamer et al. (1992); Weihong and Hamer (1993); Syromyatnikov (2010), which we discuss in section V.
i.3 Phenomenology of diagonal magnons: a short review
The purpose of this section is to give a brief (and, unavoidably, partial) overview of some of the salient features which have been the focus of much of the previous work on the excitations of this model.
There is a peculiar property of the LSWT prediction: is constant along . For convenience, we refer to magnons with these momenta as being diagonal. This dispersionless feature is a consequence of the more basic fact that at these momenta, the loworder spin wave Hamiltonian vanishes. At this order in spin wave theory, the Hamiltonian consists purely of anomalous terms, creating pairs of magnons—one on each sublattice. Hence, its vanishing has the even more peculiar consequence that, within LSWT, the diagonal magnons are purely localized on a single sublattice. This seems to have gone hitherto unnoticed.
In fact, this onedimensional flatness in the spectrum means one has a freedom in choosing a basis of energy eigenstates. By Fourier transforming the momentum eigenstates along one direction, one can thus construct eigenstates which are spatially localized onto a single diagonal (of a given sublattice), with alternating signs along this diagonal. In summary, at low order in SWT, diagonal magnons are both localized on a sublattice and a diagonal. (In section III, we show that the same features arise naturally at low order in the Ising expansion.)
Despite being flat in LSWT and LSWT+, the diagonal magnons acquire a finite but very small dispersion at higher orderHamer et al. (1992); Weihong and Hamer (1993); Syromyatnikov (2010); Uhrig and Majumdar (2013). Equivalently, this means that magnons can no longer be confined onto a single diagonal. However, it has not yet been investigated whether the aforementioned sublatticelocalization persists. We study this both numerically (section IV) and perturbatively (section III.5).
The SWT predictions at diagonal momenta, , do not agree well with other methods or experimentsRønnow et al. (2001); Christensen et al. (2007); Dalla Piazza et al. (2015)—both with respect to single and multimagnon features. Examples of previous studies include methods based on quantum Monte Carlo (QMC) combined with analytic continuationSandvik and Singh (2001); Shao et al. (2017), series expansions in (up to 14th order)Singh and Gelfand (1995); Zheng et al. (2005), exact diagonalization (ED)Lüscher and Läuchli (2009) and the continuous unitary transform (CUT) Powalski et al. (2015, 2018). All these methods predict a more pronounced local minimum of the magnon at , referred to as the roton mode, although they do not agree on its exact magnitude or shape (a quantitative discussion is deferred to section V). More strikingly, they also predict an unusually large weight in the continuum above this local minimum.
The latter phenomenology is also observed in experimentChristensen et al. (2007); Headings et al. (2010); Dalla Piazza et al. (2015), and exotic scenarios have been given to explain it. For example, it has been argued that near , the magnon can be seen as two (nearly) deconfined spinons Ho et al. (2001); Dalla Piazza et al. (2015); Shao et al. (2017); Ferrari and Becca (2018). This interpretation has subsequently been challenged by the CUT method Powalski et al. (2015, 2018), which reproduces various salient features based on a picture of stronglyinteracting magnons. The intuitive nature of said strong interactions, however, has not yet been clarified. We will argue that an Isinglike domainwall interaction naturally accounts for it.
Ii Spectral functions
In this section we discuss the transverse spectral function as defined in Eq. (3). For this, we use the numerical method introduced in Ref. Gohlke et al., 2017, which we briefly outline here. Firstly, the model (1) is put on an infinitely long cylinder whose finite, periodic direction is along a zigzag/staircase path. We define the circumference in Manhattan distance, i.e. the minimal number of bonds needed to wrap around the cylinder. In this work, . The infinite density matrix renormalization group (iDMRG) method White (1992, 1993); Kjäll et al. (2013) is used to obtain the ground stateStoudenmire and White (2013). The dynamical spinspin correlations can then be calculated by using a matrixproductoperatorbased time evolutionZaletel et al. (2015). The spectral function follows directly from Eq. (2) Schollwöck (2011).
Let us mention a few technical details before discussing the results. To minimize the effects of Fourier transforming a finitetime window, we use linear predictionWhite and Affleck (2008) to increase the time window, after which we multiply the data with a Gaussian envelope^{5}^{5}5The width of the Gaussian is chosen such that there is only little weight on the data generated by linear prediction.. This effectively introduces an artificial broadening of the spectral function with fullwidthathalfmaximum . For a given circumference, we confirm our results are converged in both bond dimension and inverse timestep by increasing both until the results no longer change. Due to the expensive nature of timeevolving large cylinders, in this work we are limited to bond dimension for the largest circumference considered (). A typical size that we used for the timestep is . The conservation of was implemented explicitly.
ii.1 Isotropic/Heisenberg model
Fig. 1 shows the transverse spectral function at the isotropic point (). Because of the cylindrical geometry, momentum is discrete along one direction and continuous along the other. This is indicated by the red lines in the Brillouin zone in Fig. 1. Since the periodic direction is along a zigzag/staircase path, the momentum cuts are lines of constant . This means we can directly access a line of diagonal magnons (as defined in section I), in the figure denoted by the line segment Y–X, where and (by symmetry). In fact, while it is true that and are symmetryequivalent in 2D, this is not strictly true on the cylinder geometry. However, such finitesize effects turn out to be small, as discussed in Appendix B. The same line of diagonal magnons, X–Y, can be accessed for if there are antiperiodic boundary conditions along the finite direction, shifting the momentum cuts as shown.
We numerically observe the Goldstone modes at the zone center and the ordering wavevector . Moreover, the intensity vanishes at the zone center, and diverges at M, consistent with the sum rules discussed in section I. This agrees with the Goldstone modes predicted by LSWT+ (solid blue line), whereas the naive evaluation of the series expansion data (up to ) does not reproduce this^{6}^{6}6More sophisticated extrapolation techniques could be and have been used for the lowenergy modes near the isotropic point Singh and Gelfand (1995); Zheng et al. (2005), however in this work we focus on the highenergy modes. (dashed black line).
On the other hand, along the Y–X line, the series expansion data fares better at reproducing the local minimum at . As discussed in section I, SWT predicts a flat dispersion along YX. Moreover, even in this linear color scale, we can see spectral weight above the singlemagnon curve near . These single and multimagnon features at the isotropic point are analyzed in more detail in section V. Here, we limit ourselves to a few general, conceptual remarks.
For , the system is gapless at the zone center, hence the multimagnon continuum starts at the onemagnon branch (up to the miniscule gap introduced by using a finite bond dimension). For the geometry, however, the antiperiodic boundary condition imply that we do not pass through the Goldstone mode, such that the multimagnon continuum is separated from the onemagnon branch. Nevertheless, these antiperiodic boundary conditions have some useful sideeffects. Due to now simulating a gapped system, it is easier to converge the numerics in the bond dimension parameter of the matrix product state describing the ground state. Moreover, it allows the ground state to spontaneously break the symmetry, even at the isotropic point. This is nontrivial given our setup, since the cylinder is effectively a onedimensional system (with a large unit cell), such that the MerminWagnerColeman theoremMermin and Wagner (1966); Coleman (1973) should prevent ordering. The catch is that the antiperiodic boundary conditions explicitly break the symmetry, although this is not locally noticeable.
This effective gap for can give us further insight into the physics beyond that of a single magnon. In Fig. 2, we show the same data in logscale. We see a continuum right above the singlemagnon branch near . However, this continuum does not fall within the frequency region of the kinematic (noninteracting) threemagnon continuum^{7}^{7}7Note that due to parity symmetry, the transverse spectral function does not pick up the twomagnon continuum.. To emphasize this, we have plotted the threemagnon continuum for this cylinder geometry in the grey shaded region. Hence, the continuum above the roton mode is instead related to (quasi)bound states. More precisely, using the insights from the upcoming section II.2, this continuum is a combination of closely packed threemagnon (quasi)bound states and a continuum made out of a single magnon and a twomagnon (quasi)bound state. This is strongly suggestive that the roton mode arises by being repelled from these stronglyinteracting states. This agrees with the conclusions of the CUT approachPowalski et al. (2015, 2018).
ii.2 Interpolating between the Ising and Heisenberg limits
Fig. 1 confirms that the DMRG method can reproduce the roton mode and the strong presence of multimagnon features near . However, to gain insight into how and from where these features appear, it is useful to interpolate from the Ising limit () to the symmetric point (). This is illustrated in Fig. 3 for . In particular, we demonstrate that the features at the isotropic point can be traced back to those in the Ising limit.
To this end, let us first identify the spectral features close to the Ising limit (). We can do this by using simple energetic arguments. Note that when , a ‘magnon’ corresponds to a single localized spin flip with energy cost per bond, totaling . In Fig. 3, we see that for for the magnon has gained some dispersion, but its energy is still roughly . Since the transverse spectral function picks up states with , the next excitation contains three magnons. Energetically, these magnons prefer to form a bound state whose domain wall counts eight bonds. Indeed, we observe bound states at energy . There are several such states at this energy due to the internal degree of freedom corresponding to orientation and shape. At even higher energies, there is the kinematic continuum made out of a twomagnon bound state () and a free magnon () with total energy around .
Having identified all spectral features for , we track their evolution as we tune in Fig. 3. At , some of the threemagnon bound states have merged. When , several of the threemagnon bound states have already been absorbed into the threemagnon continuum. Closer to the isotropic point, , the threemagnon continuum continues to come down in energy. This trend gradually continues up to .
We see that the spectral features vary continuously as a function of . In particular, we see that there is no restructuring of the magnon near for any . This was a priori not a given. Read backwards, this means that the features near the isotropic point can be continuously traced back to those in the Ising limit. Relatedly, it is worth pointing that even at the isotropic point, the multimagnon continuum is not featureless. We discuss this substructure more quantitatively in section V.
In Fig. 2, we saw how for , there is a continuum above the roton mode which is not made out of kinematic combinations of magnons. In section II.1, we claimed that it is instead a continuum made up out of (quasi)bound states. The justification for this claim is that by smoothly decreasing , the observed continuum indeed splits up into threemagnon bound states and a continuum made up out of a magnon and a twomagnon bound state.
Iii Perturbative understanding from the Ising limit
In section II we saw that we could connect spectral features of the isotropic model to those near the Ising limit. The purpose of this section is to complement this by gaining insights from loworder perturbation theory in . The point is not to see how well the isotropic point can be described quantitatively by a series expansion in Singh and Gelfand (1995); Zheng et al. (2005). Instead, we ask what the lowest order processes are that qualitatively explain certain features at the isotropic point. Intriguingly, this already naturally leads to a semiquantitative description.
We rewrite Hamiltonian (1) as with
(7) 
The Ising limit is exactly solvable: the ground state is a product Néel state, and the singlemagnon excitations consist of localized spin flips. The perturbation introduces dynamics to these static excitations. Before going through this in detail, let us give the broad brush strokes to emphasize the simplicity of both the ingredients and results.
iii.1 Overview and summary of the perturbative picture
As we will argue, the effective Hamiltonian has contributions at even order in only. The Ising magnons start to hop at order by going through a virtual threemagnon bound state (see Fig. 4(a)). Nevertheless, as we explain in section III.3, magnons with diagonal momentum are still dispersionless at this order. This is equivalent to the statement that if one builds a onemagnon state which is entirely localized on a single diagonal and has momentum along it (see Fig. 4(b)), then it cannot hop off due to destructive interference (see Fig. 4(c)).
The key to the destructive interference traces back to the fact that all threemagnon bound states have the same energy, and hence the virtual paths—half of which come with opposite signs due to the momentum—can cancel exactly. Thus from the viewpoint of the Ising expansion, such destructive interference and the resulting flatness of the diagonal magnons seems accidental. It is hence not surprising that if one goes to nexttoleading order, i.e. , the diagonal magnons acquire a dispersion. Indeed, now virtual fivemagnon bound states appear, which can have differing energies (see Fig. 4(d)).
Since the emergence of the roton mode is due to the physics of (virtual) bound states, one can indeed say that this phenomenology is due to the attractive interactions between magnons. Since this interaction is so natural in the Ising language, the qualitatively correct physics arises rather easily. Indeed, the resulting dispersion at order does not just correctly reproduce the qualitative features of having a local minimum at and maximum at , but evaluating it at even gives a semiquantitative description for the isotropic model, as we discuss in section III.4. It is moreover in remarkable proximity to the CUT predictionPowalski et al. (2015, 2018), which is a sophisticated framework for stronglyinteracting magnons. This success at relatively low order is in contrast with higherorder SWT.
That yethigherorder corrections don’t radically change the physics at hand can be confirmed by repurposing results from previous studiesSingh and Gelfand (1995); Zheng et al. (2005). In Fig. 5, we show how the dispersion along the line of diagonal magnons has certain ‘harmonics’ generated at distinct orders in (the first nontrivial harmonic appearing at ). We see that the higher harmonics die off exponentially fast, justifying a loworder picture. Note that such an exponential decay is a priori not obvious, considering that perturbation theory generically leads to an exponential proliferation of the number of terms.
Lastly, we also consider the perturbed wavefunctions at leading order. In particular, the ground state is dressed by pairs of correlated spin flips, introducing entanglement. This would usually allow one to create magnons associated to one sublattice by acting on the other sublattice. Surprisingly, our perturbative analysis implies that this is not possible for diagonal magnons. In other words, they seem to be localized on a given sublattice. This is discussed in section III.5. (We also revisit and confirm such sublatticelocalization numerically in section IV.)
We now provide the details of the story as just described.
iii.2 Ising limit and defining magnons
In the Ising limit , the ground state is a product Néel state. Let us define the A(B)sublattice to be where spins point down (up) in the ground state. For any product state in the spin basis, we can count the number of flipped spins on the Asublattice, relative to the ground state, which we denote by (similarly for ). Hence, the ground state corresponds to , whereas the singlemagnon states have , (called Amagnons) or , (Bmagnons).
The perturbation will mix states with different and , however remains a welldefined quantum number. If we perturb the system such that, for a given momentum, the onemagnon energy scale does not overlap with multimagnon states, we can nonperturbatively label the singlemagnon states by or , respectively referred to as A and Bmagnons. Since can thus not connect Amagnons to Bmagnons, we can limit our study to Amagnons. More precisely, we are interested in the effective Hamiltonian on , the Hilbert space of states satisfying and .
iii.3 Dispersionless diagonal magnons at leadingorder
The singlemagnon states are completely static and localized in the Ising limit. They moreover stay immobile at first order in . More precisely, denoting the projector onto as , then at first order we have . This is because creates a pair of A and Bmagnons out of the vacuum: it, e.g., maps , whereas it annihilates ferromagnetic bonds. More generally, flips the parity of , hence the conservation of shows that there are no contributions to at any odd order .
Thus, by standard perturbation theory, the lowestorder effective Hamiltonian on is
(8) 
where . This indeed introduces hopping, as shown in Fig. 4(a): the Amagnon can hop to any of the eight nearest sites (on the same sublattice) by going through a virtual threemagnon bound state. As discussed in section II, such a bound state involves eight ferromagnetic bonds, with total energy cost —whereas the cost of a single magnon is . Thus, the path shown in Fig. 4(a) carries a weight .
Magnons can thus hop at order . However, certain superpositions are immobile due to destructive interference. Consider, for example, the state shown in Fig. 4(b): it is localized on a single Adiagonal, with an alternating sign structure. (We can say its momentum along the diagonal is .) The magnon is unable to hop off the diagonal at order . This is illustrated in Fig. 4(c), showing two destructively interfering paths. It is important that both paths go through a virtual threemagnon bound state (both with energy ), such that the two weights cancel exactly. Equivalently, all Amagnon momentum eigenstates with —which we referred to as diagonal magnons in the previous sections—have constant energy. In summary, the diagonal magnons are dispersionless in the Ising expansion up to order . It is interesting to note that this coincides with the LSWT predictions in section I.
iii.4 Roton mode at nexttoleading order
It is hence important to see what happens at subleading order in . Here we follow the perturbative scheme by KatoKato (1949) and TakahashiTakahashi (1977). In terms of and , they constructed a generalpurpose unitary mapping which embeds the unperturbed states into the space spanned by the true eigenstates (the latter being dependent). This object gives us access to the effective Hamiltonian
(9) 
As argued before, . Moreover, and . From knowing the aforementioned object (which, for completeness, we reproduce as a function of and in Appendix A), one can derive that
(10) 
where and .
From this, one can calculate that the diagonal magnons acquire a dispersion at this order. The reason for this is in fact simple, as hinted at in section III.1. Fig. 4(d) shows two possible virtual fivemagnon bound states that can appear as intermediate states. These two states have domain walls extending over, respectively, twelve and ten bonds. Their energy is thus different, and one should not expect perfect destructive interference.
More precisely, the resulting dispersion at order is described by a simple cosinelike dispersion for the diagonal magnons: with . This has a local (roton) minimum at and a maximum at . Moreover, evaluating this at already gives a semiquantitative description of the isotropic model. We refer the reader interested in quantitative details to section V, which is devoted to an indepth comparison between various different methods.
In summary, the roton mode naturally appears in the Ising expansion. Through the property of all threemagnon bound states having the same energy, the local minimum at is absent at leading order, already indicating that it is less pronounced at the isotropic point. At the same time, its salient features readily appear at nexttoleading order, leading to a semiquantitatively correct description. From this point of view, it is indeed an interactingmagnon effect, where the interaction is based on a simple domainwall counting in the Ising limit.
iii.5 Sublatticelocalization of diagonal magnons
Aside from looking at the effective Hamiltonian, it can be instructive to consider how the eigenstates evolve with . This is exactly the information encoded in . In Fig. 6 we show the leadingorder results, both for the ground state as well as a localized Amagnon. (One would have to Fourier transform the latter to obtain an energy eigenstate.) We see that these states are dressed with ‘pair fluctuations’ whilst staying within a welldefined sector.
Having access to the perturbed states, we can ask what excitations are created upon acting with a local operator on the ground state. Fig. 6 shows that at this order, a operator does not just create a single magnon, but also threemagnon bound states. This is to be expected and is directly in line with the spectral weight observed in Fig. 3.
It is more interesting to consider what happens when applying on the Asublattice. This brings us into the sector . In other words, by acting with this operator on the Asublattice, we create a Bmagnon. This is not possible in the product state Ising limit , where acting with on the Asublattice annihilates the ground state. But as shown in Fig. 6, the perturbation introduces entanglement, such that is nonzero and has a Bmagnon on the four Bsites adjacent to the original site .
However, something unusual happens for diagonal momenta. Note that for any given Bsite, there are four adjacent Asites. If the operator we acted with on the Asublattice has momentum , then half of these four sites carry a positive sign, and half a negative sign, so that there would be perfect cancellation. In other words: we are not able to create a Bmagnon with by acting on the Asublattice.
The above used an explicit calculation, but the essential mechanism at play should hold at all orders. If we act on a given site of the Asublattice, then by the symmetry of the model, the signs and weights will be the same in all four directions. However, this is incompatible with the alternating sign structure of diagonal momenta. This argument suggests that as long as we act on a single site of the Asublattice, we cannot create a Bmagnon with a diagonal momentum. We investigate this claim of sublatticelocalization nonperturbatively in the following section, detailing its relationship to entanglement.
Iv Entanglement and sublatticelocalization of diagonal magnons
In this section we discuss a peculiar property of the entanglement in this model. As before, we define the Asublattice where the spins point down in the ground state (opposite for Bsublattice). The ground state is in the sector , and a magnon associated to the Bsublattice, a Bmagnon, is in the sector .
If the ground state had no entanglement between the two sublattices, then by acting on the Asublattice, one could not create a Bmagnon. The intuitive idea is sketched in Fig. 7(a), but the more precise wording is as follows: if (where lives on the sublattice), then, e.g., —which puts us in the sector of a Bmagnon—annihilates the ground state if . To argue this, note that since the ground state is an eigenstate of , then the factorization implies that must be an eigenstate of . Moreover, since the product Néel state has a finite^{8}^{8}8For overlaps to be welldefined, one can apply the argument to finite systems. overlap with , this fixes the eigenvalue of to be as (algebraically) small as possible. Hence, acting with the lowering operator on A must annihilate the state.
The actual ground state will of course have entanglement; for example, at leading order in , there are twosite spinflips (‘pair fluctuations’) which entangle the two sublattices, see e.g. Fig. 7(a) or Fig. 6. Thus, it is indeed generically possible to create a Bmagnon by acting on the Asublattice. However, this does not seem to be true for diagonal magnons, i.e. when . To make this precise, we introduce what we call the sublattice spectral function,
(11) 
where . There are two crucial differences that distinguish it from the usual transverse spectral function as in Eq. (2). Firstly, we only act on the Asublattice, where spins point down in the symmetrybroken ground state. Secondly, we act with the lowering operator , putting us in the sector of Bmagnons.
In Fig. 7, we show this sublattice spectral function . For convenience, we consider instead of the isotropic point, as this allows to tell the onemagnon branch straightforwardly apart from the multimagnon sector. We see the response is nonzero on almost the whole singlemagnon branch (dashed lines). However, the spectral weight is exactly zero for any of the diagonal magnons (i.e. along the border of the shaded region in the Brillouin zone). We conclude that the diagonal magnons appear to be localized on their respective sublattices.
In section III.5 we gave a symmetrybased argument for the sublatticelocalization within the perturbative framework, using the fact that we act with a singlesite operator. However, we have also numerically confirmed that the same absence of spectral weight occurs even if we act with multisite operators localized on the Asublattice (not shown). We have not found an explanation for this and it would be interesting to study this in more detail. It is an open question whether there is a probe that could directly access in an experimental setup.
V Quantitative analysis at diagonal momenta
In this section, we analyze the roton mode at and its associated multimagnon features in more quantitative detail, including a comparison to previous work.
v.1 Depth of the anomalous mode at
In Fig. 8, we consider the depth of the roton mode, i.e. the maximum of the dispersion at relative to the local minimum at . This is shown as a function of the parameter . As derived in section III, we expect the dispersion to scale as for small . For that reason, we scale our axis accordingly. The numerical results obtained with our DMRGbased method for (up to ) are shown as yellow dots. At the isotropic point (), we plot the predictions of QMCSandvik and Singh (2001); Shao et al. (2017) (extrapolated from up to ), CUTPowalski et al. (2015, 2018), EDLüscher and Läuchli (2009) (extrapolated from up to ) and SWTAnderson (1952); Kubo (1952); Oguchi (1960); Weihong et al. (1991); Hamer et al. (1992); Weihong and Hamer (1993); Syromyatnikov (2010).
This quantity is extracted from the spectral function by fitting the singlemagnon response with a Gaussian^{9}^{9}9Note that the width of the Gaussian is known; see section II.. Near the isotropic point, this fitting is somewhat subtle, as the onemagnon response is not easily separated from the multimagnon weight. Fitting the (quasi)bound states just above the singlemagnon branch—which are highly relevant near —as well, we obtain results which are stable with respect to the numerical parameters.
At the isotropic point, the method that DMRG is closest to is QMC. For , we obtain and . This can be compared to the QMCShao et al. (2017) extrapolations and . We are unable to perform a finitecircumference analysis, since for there are domainwall excitations (wrapping around the circumference), while for the system is gapless such that we expect stronger finitecircumference effects. However, one can compare our results to the finitesize QMC data with linear dimension , corresponding to our largest cylinder circumference. In that case, QMC obtainsShao et al. (2017) . Taking this to give a rough finitesize estimate, we note that we are within the same distance to the extrapolated QMC data (although at the opposite extreme).
It is illuminating to not just focus on the isotropic case, but to track the evolution as a function of as shown in Fig. 8. The dashed lines are from series expansionsSingh and Gelfand (1995); Zheng et al. (2005) to different orders in . If we track the lowestorder result toward the isotropic point , we already obtain the correct order of magnitude. This is moreover in striking proximity to the CUT prediction. As we include higher order terms, we see that the dispersion gradually creeps up, showing no real sign of convergence. However, since SWT results are analytic in the modified parameter , it is suggestive to rewrite the series expansion in terms of . Doing so, and building a Padé approximant out of it, we obtain the solid line in Fig. 8.
We find that the Padé approximant is remarkably robust: the approximants , , , , all give virtually indistinguishable results! This stability suggests that the solid line could be a reasonable prediction for the true evolution of the dispersion as a function of . Exactly at the isotropic point, there is a small caveat: any power series in , when rewritten in terms of , generically predicts a diverging slope at . Hence, also in this case, we find that the dashed curve is finite at but its slope is vertical. It is not clear whether this particular feature is physical or not. Other than that, we expect that the Padé approximant should be reliable and we are encouraged by the fact that our numerical results agree so well with the Padé curve, indicating that finitesize corrections for are already rather small.
It would be interesting to investigate to what extent the Padé approximant gives the correct prediction. In particular, it could be worthwhile to test and compare the other methods^{10}^{10}10We note that applying the CUT method developed in Refs. Powalski et al., 2015, 2018 to would be distinct from the CUT method that was applied to the XXZ model in Ref. Dusuel et al., 2010: in the latter case, the CUT method was perturbative in the parameter , hence being an alternative way of calculating the series expansion coefficients. at .
v.2 Dispersion relation
Aside from studying the depth of the roton mode, we also consider its shape. Our numerical result is shown in Fig. 9 (solid black line). We compare it to the functional forms obtained by CUT, QMC and series expansion. As discussed in section III, the lowestorder nontrivial prediction from the Ising expansion is a simple cosinelike dispersion. This is not significantly altered at higher orders in (or at the very least, only very slowly so), as was shown in Fig. 5. The purpose of the inset is to show a comparison with this simple cosine.
Remarkably, the CUT result is perfectly fit by a single harmonic. In conjunction with section V.1, we thus conclude that the fourth order extrapolation from the Ising expansion seems to be in striking proximity to the CUT prediction of the roton dispersion. Our result, on the other hand, while being dominated by the same cosine, also contains the higher harmonics (which are qualitatively generated in the Ising expansion). We point out that the QMC curve has more structure near , with a possible small subsidiary maximum at the X point itself; it would be interesting to investigate its origin.
v.3 Multimagnon features
Lastly, in Fig. 10 we show a more detailed slice of the transverse spectral function first shown in in Fig. 1. At two values of the momentum, and , we show the spectral weight as a function of . This numerical data is for the smaller circumference , since as discussed in section II, in that case the system is gapless. Being gapless, one expects there to be more significant finitesize effects on a quantitative level, but the qualitative shape should look more like the 2D limit than the data would.
For either momentum, we clearly see the singlemagnon peak (broadened due to our finitetime window, as explained in section I) and a broad threemagnon continuum. Moreover, for , we recognize a second, smaller peak. This is a threemagnon (quasi)bound state.
We would like to comment on the following two features of Fig. 10. Firstly, there is considerable weight in the multimagnon sector at , and not at . Due to not having a tight grasp on finitesize effects for , we do not believe there is great value in quoting precise numbers, but at , only roughly half the weight is on the single magnon. Secondly, there is certainly a substructure to the multimagnon weight. This seems to be in contrast to the featureless spectral function observed in a recent Monte Carlo studyShao et al. (2017). Due to the absence of a finitesize analysis, this substructure is not conclusive and it would be interesting to investigate this further.
Vi Conclusion
We have studied the spectral properties of the Heisenberg model, allowing for Ising anisotropy, using two complementary methods: a DMRGbased approach to obtain the unbiased dynamical structure factor (for certain circumferences), and a loworder perturbative expansion around the Ising limit to give insights into the physical mechanisms at play.
One of our key messages is that some of its salient dynamical features are naturally accounted for starting from the Ising limit. The exactly soluble Ising limit itself has stronglyattractive magnon interactions based on simple domainwallcounting. One does not expect such an Isingbased picture to be applicable to the lowenergy hydrodynamic Goldstone modes as , but our work shows that it does remain relevant for the highenergy magnons.
In particular, the lowest nontrivial order in the Ising expansion already captures the physics of the roton mode—i.e. the dispersion along —at a semiquantitative level. Furthermore, we clarified its physical origin in terms of the properties of virtual bound states which mediate the magnon’s hopping. On a more phenomenological level, the spectral function for the Heisenberg model obtained with DMRG shows that the anomalous local minimum at grows monotonically when coming from the Ising limit. Moreover, even the strong continuum above this mode is continuously connected to spectral features near the Ising limit.
The spectral function on the geometry with directly supports the point of view that the physics near is beyond that of a perturbativelydressed magnon: in Fig. 2 we saw that there is a continuum directly above the magnon which is not a standard kinematic threemagnon continuum. Instead, the relevant physics is due to (quasi)bound states arising from attractive interactions of magnons ‘sharing’ their domain walls. This agrees with the insights of the CUTbased analysisPowalski et al. (2015, 2018). It would be interesting to further explore the potential link between the interactions arising in the CUT framework and that of the Isingbased picture. Remarkably, as far as the roton mode is concerned, a loworder Ising expansion gives predictions close to that of the CUT approach, but at this point it is not clear whether this is accidental or not.
Our study has also uncovered a curious spectral property. It turns out that magnons with are localized on their respective sublattices. This means that any operator localized on the Asublattice cannot create such a magnon associated to the Bsublattice. We have emphasized that this is rather unusual, since entanglement in the ground state would generically allow for it. Interestingly, such sublatticelocalization is also predicted at low order in spin wave theory. It is not yet clear to what extent it is compatible with higherorder corrections in .
Having established a link between the spectral properties of the Heisenberg model and an Isingbased picture, several questions can be raised. Firstly, does the latter simple picture also give an intuitive explanation for why the (quasi)bound states bunch up near ? Secondly, as already mentioned in the introduction, the spectral features under discussion may not be uniquely described by one picture as opposed to another. Hence, we are not proposing the Isingbased picture as a complete framework, but rather as a very simple account of various features. Hence, it leaves open a link between the Ising limit and the other approaches that have been explored thus far. This could be interesting to investigate further. For example, a recent work interpreted the properties of the Heisenberg model in the context of the larger modelShao et al. (2017), and it could be worthwhile to explore its interplay with an Ising anisotropy.
Vii Acknowledgements
The authors would like to thank Ehud Altman, Sylvain Capponi, Alexander Chernyshev, Fabian Essler, Efstratios Manousakis, Kai Schmidt and Götz Uhrig for helpful discussions. In particular, we thank Jan Oitmaa for sharing the series expansion data of Ref. Zheng et al., 2005 up to twelfth order. RV was supported by the German Research Foundation (DFG) through the Collaborative Research Centre SFB 1143 and FP acknowledges the support of the DFG Research Unit FOR 1807 through grants no. PO 1370/21, TRR80, the Nanosystems Initiative Munich (NIM) by the German Excellence Initiative, and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 771537).
References
 Heisenberg (1928) W. Heisenberg, Zeitschrift für Physik 49, 619 (1928).
 Hulthén (1936) L. Hulthén, Proc. Roy. Acad. Sci. Amsterdam 39, 190 (1936).
 Anderson (1952) P. W. Anderson, Phys. Rev. 86, 694 (1952).
 Keffer et al. (1953) F. Keffer, H. Kaplan, and Y. Yafet, American Journal of Physics 21, 250 (1953).
 Sandvik and Singh (2001) A. W. Sandvik and R. R. P. Singh, Phys. Rev. Lett. 86, 528 (2001).
 Shao et al. (2017) H. Shao, Y. Q. Qin, S. Capponi, S. Chesi, Z. Y. Meng, and A. W. Sandvik, Phys. Rev. X 7, 041072 (2017).
 Lüscher and Läuchli (2009) A. Lüscher and A. M. Läuchli, Phys. Rev. B 79, 195102 (2009).
 Singh (1989) R. R. P. Singh, Phys. Rev. B 39, 9760 (1989).
 Gelfand (1996) M. P. Gelfand, Solid State Communications 98, 11 (1996).
 Singh and Gelfand (1995) R. R. P. Singh and M. P. Gelfand, Phys. Rev. B 52, R15695 (1995).
 Zheng et al. (2005) W. Zheng, J. Oitmaa, and C. J. Hamer, Phys. Rev. B 71, 184440 (2005).
 Powalski et al. (2015) M. Powalski, G. S. Uhrig, and K. P. Schmidt, Phys. Rev. Lett. 115, 207202 (2015).
 Powalski et al. (2018) M. Powalski, K. P. Schmidt, and G. S. Uhrig, SciPost Phys. 4, 001 (2018).
 Anderson (1987) P. W. Anderson, Science 235, 1196 (1987).
 Rice (1987) T. M. Rice, Zeitschrift für Physik B Condensed Matter 67, 141 (1987).
 Chakravarty et al. (1988) S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. Lett. 60, 1057 (1988).
 Ho et al. (2001) C.M. Ho, V. N. Muthukumar, M. Ogata, and P. W. Anderson, Phys. Rev. Lett. 86, 1626 (2001).
 Hsu (1990) T. C. Hsu, Phys. Rev. B 41, 11379 (1990).
 Syljuåsen and Lee (2002) O. F. Syljuåsen and P. A. Lee, Phys. Rev. Lett. 88, 207207 (2002).
 Dalla Piazza et al. (2015) B. Dalla Piazza, M. Mourigal, N. B. Christensen, G. J. Nilsen, P. TregennaPiggott, T. G. Perring, M. Enderle, D. F. McMorrow, D. A. Ivanov, and H. M. Rønnow, Nature Physics 11, 62 (2015).
 Ferrari and Becca (2018) F. Ferrari and F. Becca, ArXiv eprints (2018), arXiv:1805.09287 [condmat.strel] .
 Nielsen and Chadha (1976) H. Nielsen and S. Chadha, Nuclear Physics B 105, 445 (1976).
 Rønnow et al. (2001) H. M. Rønnow, D. F. McMorrow, R. Coldea, A. Harrison, I. D. Youngson, T. G. Perring, G. Aeppli, O. Syljuåsen, K. Lefmann, and C. Rischel, Phys. Rev. Lett. 87, 037202 (2001).
 Christensen et al. (2007) N. B. Christensen, H. M. Rønnow, D. F. McMorrow, A. Harrison, T. G. Perring, M. Enderle, R. Coldea, L. P. Regnault, and G. Aeppli, Proceedings of the National Academy of Science 104, 15264 (2007).
 Landau (1941) L. Landau, Phys. Rev. 60, 356 (1941).
 White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
 White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
 Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011).
 Stoudenmire and White (2013) E. M. Stoudenmire and S. R. White, Phys. Rev. B 87, 155137 (2013).
 Zaletel et al. (2015) M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B 91, 165112 (2015).
 Gohlke et al. (2017) M. Gohlke, R. Verresen, R. Moessner, and F. Pollmann, Phys. Rev. Lett. 119, 157203 (2017).
 (32) To argue that , one can use that has a welldefined quantum number with respect to the modified translation operator .
 (33) One could instead work in a local frame where there is only one band, but this obscures some of the physics. In particular, then is no longer a good quantum number.
 (34) To wit, the HolsteinPrimakoff transformation is . The dominant term of the square root is of order , which leads to LSWT.
 Kubo (1952) R. Kubo, Phys. Rev. 87, 568 (1952).
 Stringari (1994) S. Stringari, Phys. Rev. B 49, 6710 (1994).
 Oguchi (1960) T. Oguchi, Phys. Rev. 117, 117 (1960).
 Weihong et al. (1991) Z. Weihong, J. Oitmaa, and C. J. Hamer, Phys. Rev. B 43, 8321 (1991).
 (39) and , where the integration is over .
 Hamer et al. (1992) C. J. Hamer, Z. Weihong, and P. Arndt, Phys. Rev. B 46, 6276 (1992).
 Weihong and Hamer (1993) Z. Weihong and C. J. Hamer, Phys. Rev. B 47, 7961 (1993).
 Syromyatnikov (2010) A. V. Syromyatnikov, Journal of Physics: Condensed Matter 22, 216003 (2010).
 Uhrig and Majumdar (2013) G. S. Uhrig and K. Majumdar, The European Physical Journal B 86, 282 (2013).
 Headings et al. (2010) N. S. Headings, S. M. Hayden, R. Coldea, and T. G. Perring, Phys. Rev. Lett. 105, 247001 (2010).
 Kjäll et al. (2013) J. A. Kjäll, M. P. Zaletel, R. S. K. Mong, J. H. Bardarson, and F. Pollmann, Phys. Rev. B 87, 235106 (2013).
 White and Affleck (2008) S. R. White and I. Affleck, Phys. Rev. B 77, 134437 (2008).
 (47) The width of the Gaussian is chosen such that there is only little weight on the data generated by linear prediction.
 Oitmaa (2018) J. Oitmaa, private communications (2018).
 (49) More sophisticated extrapolation techniques could be and have been used for the lowenergy modes near the isotropic point Singh and Gelfand (1995); Zheng et al. (2005), however in this work we focus on the highenergy modes.
 Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
 Coleman (1973) S. Coleman, Communications in Mathematical Physics 31, 259 (1973).
 (52) Note that due to parity symmetry, the transverse spectral function does not pick up the twomagnon continuum.
 Kato (1949) T. Kato, Progress of Theoretical Physics 4, 514 (1949).
 Takahashi (1977) M. Takahashi, Journal of Physics C: Solid State Physics 10, 1289 (1977).
 (55) For overlaps to be welldefined, one can apply the argument to finite systems.
 (56) Note that the width of the Gaussian is known; see section II.
 (57) We note that applying the CUT method developed in Refs. \rev@citealpnumPowalski15,Powalski17 to would be distinct from the CUT method that was applied to the XXZ model in Ref. \rev@citealpnumDusuel10: in the latter case, the CUT method was perturbative in the parameter , hence being an alternative way of calculating the series expansion coefficients.
 Dusuel et al. (2010) S. Dusuel, M. Kamfor, K. P. Schmidt, R. Thomale, and J. Vidal, Phys. Rev. B 81, 064412 (2010).
Appendix A effective Hamiltonians to arbitrary order
For completeness, we reproduce the general perturbative scheme that allows to obtain a welldefined effective Hamiltonian to any order. As in the main text, we consider a Hamiltonian of the form . Let be the Hilbert space associated to the degenerate eigenvalue of . Moreover, let be the projector onto , and . Note that we can decompose the total Hilbert space as .
Suppose that for , we can decompose the Hilbert space into in such a way that

;

is a smooth function of ;

the Hamiltonian respects the decomposition.
Physically speaking, this formalizes the idea that we want the the energy levels of the sector we are interested to stay separated from the remaining levels; otherwise the idea of an effective Hamiltonian is misguided.
If those conditions hold, the work by KatoKato (1949) and TakashiTakahashi (1977) showed that for the same range , one can explicitly construct a smooth unitary mapping which maps the unperturbed eigenstates into the perturbed ones. Hence, the desired effective Hamiltonian on is then simply .
To perturbatively express as a function of the known quantities , and , it is useful to define a few other quantities. Firstly, let be the projector onto ; we will derive a perturbative expression for this object as well. Secondly, define
(12)  
(13) 
Note that is expressed in terms of known quantities.
Moreover, it can be shown that the following function then has the desired properties:
(15) 
It can be proven that as thus defined indeed satisfies .
In terms of the above quantities, we thus have that
(16) 
The result in Eq. (16), namely that , is a direct consequence of and the fact that .
Appendix B finitecircumference effects for symmetryequivalent points
The twodimensional models enjoys symmetries which are broken when putting the model on a cylinder. Rotating the square lattice by gives an example. One can turn this curse into a blessing, since it gives us a direct probe of the finitecircumference effects. More precisely, suppose and are two distinct momenta which are symmetryequivalent in the twodimensional limit, but which are not symmetryequivalent on a cylinder with circumference . Hence, gives us a rough sense of how strong the finitecircumference effects are. Fig. 11 plots the maximum of this over all pairs of symmetryequivalent points. We see that it goes down as , as expected.