Quantum dynamics of the square lattice Heisenberg model

# Quantum dynamics of the square lattice Heisenberg model

Ruben Verresen Department of Physics, T42, Technische Universität München, 85748 Garching, Germany Max-Planck-Institute for the Physics of Complex Systems, 01187 Dresden, Germany    Frank Pollmann Department of Physics, T42, Technische Universität München, 85748 Garching, Germany    Roderich Moessner Max-Planck-Institute for the Physics of Complex Systems, 01187 Dresden, Germany
July 16, 2019
###### Abstract

Despite nearly a century of study of the Heisenberg model on the square lattice, there is still disagreement on the nature of its high-energy 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 semi-quantitative description of certain observed features arises already at the lowest non-trivial order in perturbation theory around the Ising limit. Moreover, our analysis uncovers that high-energy 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 low-energy propertiesHulthén (1936); Anderson (1952). Its anomalous terms lead to a ‘relativistic’ (linear) low-energy 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 high-order 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 ad-hoc 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 strongly-interacting 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 perturbatively-dressed single-magnon 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 low-energy 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 high-energy 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 quasi-particle 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 single-magnon picture, and it uncovers a hitherto-unrecognized property of of magnons with , which we refer to as sublattice-localization. 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 semi-quantitatively accounted for by the lowest non-trivial order. Moreover, our numerical method allows us to study such an XXZ model (with dominant easy-axis 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 domain-wall-counting 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 small-yet-nonzero magnitude and shape of the roton mode’s dispersion. Moreover, even the aforementioned phenomenon of sublattice-localization can be accounted for within a low-order 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 sublattice-localization, naturally arise within a low-order perturbative picture. The apparently hitherto-unexplored phenomenon of sublattice-localization 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 easy-axis anisotropy:

 H=J∑⟨n,m⟩(SznSzm+λ[SxnSxm+SynSym]) (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 :

 Sγγ(k,ω)=12π∑r∫∞−∞ei(ωt−k⋅r)Cγγ(r,t)dt, (2)

which is normalized as . We focus on the transverse spectral function:

 St(k,ω)=Sxx(k,ω)+Syy(k,ω). (3)

This object gives direct insight into the excitations above the ground state. If , one can show111To argue that , one can use that has a well-defined quantum number with respect to the modified translation operator . that

 Sγγ(k,ω)=∑αδ(ω−(ωα−ω0))|⟨α|~Sγk|0⟩|2 (4)

where . It is natural to choose a basis , where is the momentum with respect to the translation symmetry of the two-site 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 bands222One 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 order333To wit, the Holstein-Primakoff 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)

 εLSWTk=√4−λ2(cos(kx)+cos(ky))2. (5)

Hence for , there are two linearly-dispersing 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)

 εLSWT+1/Sk=aεLSWTk−(1−λ2)b⎛⎝2εLSWTk−εLSWTk2⎞⎠ (6)

where and are (-dependent) constants444 and , where the integration is over .. At the isotropic point, this correction is only a momentum-independent rescaling. Higher-order 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 low-order 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 one-dimensional 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 sublattice-localization 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 multi-magnon 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 strongly-interacting magnons. The intuitive nature of said strong interactions, however, has not yet been clarified. We will argue that an Ising-like domain-wall 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 spin-spin correlations can then be calculated by using a matrix-product-operator-based 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 finite-time window, we use linear predictionWhite and Affleck (2008) to increase the time window, after which we multiply the data with a Gaussian envelope555The 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 full-width-at-half-maximum . For a given circumference, we confirm our results are converged in both bond dimension and inverse time-step by increasing both until the results no longer change. Due to the expensive nature of time-evolving large cylinders, in this work we are limited to bond dimension for the largest circumference considered (). A typical size that we used for the time-step 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 symmetry-equivalent in 2D, this is not strictly true on the cylinder geometry. However, such finite-size 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 this666More sophisticated extrapolation techniques could be and have been used for the low-energy modes near the isotropic point Singh and Gelfand (1995); Zheng et al. (2005), however in this work we focus on the high-energy 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 Y-X. Moreover, even in this linear color scale, we can see spectral weight above the single-magnon curve near . These single- and multi-magnon 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 multi-magnon continuum starts at the one-magnon 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 multi-magnon continuum is separated from the one-magnon branch. Nevertheless, these antiperiodic boundary conditions have some useful side-effects. 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 non-trivial given our set-up, since the cylinder is effectively a one-dimensional system (with a large unit cell), such that the Mermin-Wagner-Coleman 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 log-scale. We see a continuum right above the single-magnon branch near . However, this continuum does not fall within the frequency region of the kinematic (non-interacting) three-magnon continuum777Note that due to parity symmetry, the transverse spectral function does not pick up the two-magnon continuum.. To emphasize this, we have plotted the three-magnon 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 three-magnon (quasi-)bound states and a continuum made out of a single magnon and a two-magnon (quasi-)bound state. This is strongly suggestive that the roton mode arises by being repelled from these strongly-interacting 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 multi-magnon 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 two-magnon 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 three-magnon bound states have merged. When , several of the three-magnon bound states have already been absorbed into the three-magnon continuum. Closer to the isotropic point, , the three-magnon 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 multi-magnon 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 three-magnon bound states and a continuum made up out of a magnon and a two-magnon 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 low-order 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 semi-quantitative description.

We rewrite Hamiltonian (1) as with

 H0=J∑⟨n,m⟩SznSzm,V=J2∑⟨n,m⟩(S+nS−m+S−nS+m). (7)

The Ising limit is exactly solvable: the ground state is a product Néel state, and the single-magnon 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 three-magnon 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 one-magnon 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 three-magnon 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 next-to-leading order, i.e. , the diagonal magnons acquire a dispersion. Indeed, now virtual five-magnon 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 semi-quantitative 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 strongly-interacting magnons. This success at relatively low order is in contrast with higher-order SWT.

That yet-higher-order 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 non-trivial harmonic appearing at ). We see that the higher harmonics die off exponentially fast, justifying a low-order 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 sublattice-localization 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 A-sublattice, relative to the ground state, which we denote by (similarly for ). Hence, the ground state corresponds to , whereas the single-magnon states have , (called A-magnons) or , (B-magnons).

The perturbation will mix states with different and , however remains a well-defined quantum number. If we perturb the system such that, for a given momentum, the one-magnon energy scale does not overlap with multi-magnon states, we can non-perturbatively label the single-magnon states by or , respectively referred to as A- and B-magnons. Since can thus not connect A-magnons to B-magnons, we can limit our study to A-magnons. More precisely, we are interested in the effective Hamiltonian on , the Hilbert space of states satisfying and .

Note that these Ising magnons are exactly those encountered in Fig. 3. As discussed in section II, a single magnon has a domain wall crossing four bonds and hence has energy relative to the ground state.

### iii.3 Dispersionless diagonal magnons at leading-order

The single-magnon 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 B-magnons 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 lowest-order effective Hamiltonian on is

 Heff=E0P0+λ2P0VG0VP0+O(λ4), (8)

where . This indeed introduces hopping, as shown in Fig. 4(a): the A-magnon can hop to any of the eight nearest sites (on the same sublattice) by going through a virtual three-magnon 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 A-diagonal, 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 three-magnon bound state (both with energy ), such that the two weights cancel exactly. Equivalently, all A-magnon 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 next-to-leading order

It is hence important to see what happens at sub-leading order in . Here we follow the perturbative scheme by KatoKato (1949) and TakahashiTakahashi (1977). In terms of and , they constructed a general-purpose 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

 Heff:=Γ†λHΓλ=:∞∑n=0λnH(n)eff. (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

 H(4)eff=P0V~G0V~G0V~G0VP0−12{H(2)eff,P0V~G20VP0} (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 five-magnon 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 cosine-like dispersion for the diagonal magnons: with . This has a local (roton) minimum at and a maximum at . Moreover, evaluating this at already gives a semi-quantitative description of the isotropic model. We refer the reader interested in quantitative details to section V, which is devoted to an in-depth comparison between various different methods.

In summary, the roton mode naturally appears in the Ising expansion. Through the property of all three-magnon 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 next-to-leading order, leading to a semi-quantitatively correct description. From this point of view, it is indeed an interacting-magnon effect, where the interaction is based on a simple domain-wall counting in the Ising limit.

### iii.5 Sublattice-localization 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 leading-order results, both for the ground state as well as a localized A-magnon. (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 well-defined 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 three-magnon 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 A-sublattice. This brings us into the sector . In other words, by acting with this operator on the A-sublattice, we create a B-magnon. This is not possible in the product state Ising limit , where acting with on the A-sublattice annihilates the ground state. But as shown in Fig. 6, the perturbation introduces entanglement, such that is nonzero and has a B-magnon on the four B-sites adjacent to the original site .

However, something unusual happens for diagonal momenta. Note that for any given B-site, there are four adjacent A-sites. If the operator we acted with on the A-sublattice 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 B-magnon with by acting on the A-sublattice.

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 A-sublattice, 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 A-sublattice, we cannot create a B-magnon with a diagonal momentum. We investigate this claim of sublattice-localization non-perturbatively in the following section, detailing its relationship to entanglement.

## Iv Entanglement and sublattice-localization of diagonal magnons

In this section we discuss a peculiar property of the entanglement in this model. As before, we define the A-sublattice where the spins point down in the ground state (opposite for B-sublattice). The ground state is in the sector , and a magnon associated to the B-sublattice, a B-magnon, is in the sector .

If the ground state had no entanglement between the two sublattices, then by acting on the A-sublattice, one could not create a B-magnon. 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 B-magnon—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 finite888For overlaps to be well-defined, 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 two-site spin-flips (‘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 B-magnon by acting on the A-sublattice. 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,

 SA→B(k,ω)=∑αδ(ω−(ωα−ω0))|⟨α|~S−A,k|0⟩|2 (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 A-sublattice, where spins point down in the symmetry-broken ground state. Secondly, we act with the lowering operator , putting us in the sector of B-magnons.

In Fig. 7, we show this sublattice spectral function . For convenience, we consider instead of the isotropic point, as this allows to tell the one-magnon branch straightforwardly apart from the multi-magnon sector. We see the response is non-zero on almost the whole single-magnon 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 symmetry-based argument for the sublattice-localization within the perturbative framework, using the fact that we act with a single-site operator. However, we have also numerically confirmed that the same absence of spectral weight occurs even if we act with multi-site operators localized on the A-sublattice (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 set-up.

## V Quantitative analysis at diagonal momenta

In this section, we analyze the roton mode at and its associated multi-magnon features in more quantitative detail, including a comparison to previous work.

### v.1 Depth of the anomalous mode at k=(π,0)

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 DMRG-based 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 single-magnon response with a Gaussian999Note that the width of the Gaussian is known; see section II.. Near the isotropic point, this fitting is somewhat subtle, as the one-magnon response is not easily separated from the multi-magnon weight. Fitting the (quasi-)bound states just above the single-magnon 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 finite-circumference analysis, since for there are domain-wall excitations (wrapping around the circumference), while for the system is gapless such that we expect stronger finite-circumference effects. However, one can compare our results to the finite-size 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 finite-size 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 lowest-order 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 finite-size 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 methods101010We 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 lowest-order non-trivial prediction from the Ising expansion is a simple cosine-like 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 Multi-magnon 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 finite-size 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 single-magnon peak (broadened due to our finite-time window, as explained in section I) and a broad three-magnon continuum. Moreover, for , we recognize a second, smaller peak. This is a three-magnon (quasi-)bound state.

We would like to comment on the following two features of Fig. 10. Firstly, there is considerable weight in the multi-magnon sector at , and not at . Due to not having a tight grasp on finite-size 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 multi-magnon 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 finite-size 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 DMRG-based approach to obtain the unbiased dynamical structure factor (for certain circumferences), and a low-order 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 strongly-attractive magnon interactions based on simple domain-wall-counting. One does not expect such an Ising-based picture to be applicable to the low-energy hydrodynamic Goldstone modes as , but our work shows that it does remain relevant for the high-energy magnons.

In particular, the lowest non-trivial order in the Ising expansion already captures the physics of the roton mode—i.e. the dispersion along —at a semi-quantitative 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 perturbatively-dressed magnon: in Fig. 2 we saw that there is a continuum directly above the magnon which is not a standard kinematic three-magnon 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 CUT-based 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 Ising-based picture. Remarkably, as far as the roton mode is concerned, a low-order 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 A-sublattice cannot create such a magnon associated to the B-sublattice. We have emphasized that this is rather unusual, since entanglement in the ground state would generically allow for it. Interestingly, such sublattice-localization is also predicted at low order in spin wave theory. It is not yet clear to what extent it is compatible with higher-order corrections in .

Having established a link between the spectral properties of the Heisenberg model and an Ising-based 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 Ising-based 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/2-1, 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).

## Appendix A effective Hamiltonians to arbitrary order

For completeness, we reproduce the general perturbative scheme that allows to obtain a well-defined 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

1. ;

2. is a smooth function of ;

3. 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

 S0 :=−P0, (12) Sk :=~G0(E0)k:=(Q01E0−H0Q0)k(k≠0). (13)

Note that is expressed in terms of known quantities.

One can then deriveKato (1949); Takahashi (1977) that

 PλP0=∞∑n=0λn∑k1+k2+⋯+kn=n,ki≥0Sk1VSk2V⋯SknVP0. (14)

Moreover, it can be shown that the following function then has the desired properties:

 Γλ:=PλP0(P0+∞∑n=1(2n−1)!!(2n)!![P0−P0PλP0]n). (15)

It can be proven that as thus defined indeed satisfies .

In terms of the above quantities, we thus have that

 Heff :=Γ†λHΓλ=E0P0+λΓ†λVΓλ−Γ†λS−1Γλ. (16)

The result in Eq. (16), namely that , is a direct consequence of and the fact that .

## Appendix B finite-circumference effects for symmetry-equivalent points

The two-dimensional 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 finite-circumference effects. More precisely, suppose and are two distinct momenta which are symmetry-equivalent in the two-dimensional limit, but which are not symmetry-equivalent on a cylinder with circumference . Hence, gives us a rough sense of how strong the finite-circumference effects are. Fig. 11 plots the maximum of this over all pairs of symmetry-equivalent points. We see that it goes down as , as expected.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters