From deconfined spinons to coherent magnons in an antiferromagnetic Heisenberg chain with long range interactions

We study the spectrum and the nature of the excitations of an antiferromagnetic (AFM) Heisenberg chain with staggered long range interactions, both numerically using the time-dependent density matrix renormalization group (tDMRG) method and by means of a multi-spinon approximation that qualitatively explains its main features. The unfrustrated long-range nature of the exchange effectively increases the dimensionality of the system and the chain is able to undergo true symmetry breaking and develop long range order, transitioning from a gapless spin liquid to a gapless ordered AFM phase. We calculated the momentum resolved spin dynamical structure factor and found that for weakly decaying interactions the emergence of N\'eel order can be associated to the formation of bound states of spinons that become coherent magnons. The quasiparticle band leaks out from the two-spinon continuum that is pushed up to higher energies. Our physical picture is also supported by an analysis of the behavior of the excitations in real-time.

Spin-chains are not well ordered antiferromagnets: their correlations decay algebraically and they do not develop true long-range order. Higher dimensional magnets may develop long-range order and, in such a case, excitations are gapless magnons with well-defined Goldstone modes. Spinons, or fractionalized excitations, are not only a feature of 1D spin chains but are also expected to emerge in 2D spin liquids -states that do not break continuous symmetries -as postulated by the theory of deconfined quantum criticality [19][20][21] . To reconcile these two pictures we interpret magnons (which carry spin 1) as bound states of spinons (that carry spin 1/2). Notice , however, that in 1D and even in 2D spin liquids 22 , it is possible to have bound states of spinons without long range order. In such a case, these excitations are referred to as "triplons" (the simplest example is a triplet excitation on top of a dimerized gapped valence bond solid) 23 .
On the 2D square lattice, a prototypical antiferromagnet, the spin wave dispersion agrees with numerical results with high accuracy in the entire Brillouin zone, with the only deviations along the (π, 0) − (π/2, π/2) path 25,26 . Along this segment, the spin-wave theory dispersion is essentially flat, while numerical results indicate a dip. Recent experiments are in excellent agreement with numerics [27][28][29] which has prompted the speculation of physics beyond magnons. In particular, in recent lowtemperature polarized neutron scattering experiments 29 , the broad and spin-isotropic continuum in S z (k, ω) at q = (π, 0) was interpreted as a sign of deconfinement of spinons in a region of momentum space. Recent numerical Monte Carlo results show, however, that magnons are still present and, even though their spectral weight may become very small, it never vanishes 30 .
Similarly, neutron scattering experiments on the 2D triangular antiferromagnet Ba 3 CoSb 2 O 9 31,32 indicate that the spectrum consists not only of low energy magnon branches, but also high energy continua with a separation of the order of the exchange interaction J. The disagreement with expectations from spin-wave theory stimulated further theoretical work 33,34 that pointed at deconfined spinons being the culprits of the high-energy features, with magnons consisting of bound states of spinons.
The above considerations beg the questions: how do gapless spinons evolve into gapless magnons and singularities in the spectrum into coherent Goldstone modes as we transition from a gapless spin liquid into a gapless antiferromagnet with long range order? 35,36 In order to address these issues we resort to one dimensional spin-1/2 chains with staggered SU (2)-symmetric long-range interactions that allow us to realize actual spontaneous symmetry breaking and true antiferromagnetic order: The general problem can be formulated by means of the following Hamiltonian: The long range nature of the interactions artificially increases the dimensionality of the problem and circumvents the restrictions imposed by the Mermin-Wagner theorem. At the same time, they introduce volume-law entanglement, making the calculations more challenging. However, the staggered phase enhances antiferromagnetism and avoids frustration, also making it amenable to quantum Monte Carlo calculations 24,37,38 .
The phase diagram of the extended Heisenberg chain as a function of the coupling λ and exponent α was obtained by Laflorencie et al in Ref. 24 using quantum Monte Carlo, who found a critical line separating Néel ordered and disordered phases with dynamical exponent z < 1 (see Fig.1(a)). This indicates that the system generally does not admit a description in terms of conformal field theory (which is not surprising), which should be manifested in the finite-size scaling of the spin gap and the curvature of the dispersion (which is sublinear), as well as the behavior of the entanglement entropy. Therefore, the chain undergoes a transition from a gapless ordered phase with strong AFM correlations, to a gapless disordered one with fractionalized excitations. The two regimes are characterized by an order parameter, the staggered magnetization m = S z (k = π), and by the nature of the excitations that should be reflected in the spectrum, given by its dynamical spin structure factor S z (k, ω).
In this work, we focus on the particular case of λ = 1. In this limit, the Hamiltonian becomes: (2) The manuscript is organized as follows: in section II we discuss the spectral function obtained by means of the time-dependent density matrix renormalization group method (tDMRG) [39][40][41][42] . In section III, we describe a multi-spinon analytical approach and analyze the results from the point of view of confined spinon excitations. We provide an intuitive picture of the nature of Momentum resolved dynamical structure factor S z (k, ω) of the Heisenberg chain with long range interactions, λ = 1, and different values of α across the phase transition.
FIG. 3. Spin dynamical structure factor across k = π/2 cuts for λ = 1 and different values of α. The emergence of a sharp isolated peak leaking out of the continuum is clearly observed as α decreases and long range order is developed. The curves are shifted for clarity. Negative values are artifacts of the Fourier transform, as described in the text. the excitations by studying their behavior in real time in section IV and we finally conclude with a summary and discussion.

II. SPIN DYNAMICS
In the disordered phase of the model Eq.
(2), excitations are described in terms of deconfined spinons. Assuming a spinon dispersion (k), the two-spinon continuum is constructed by all possible energies 2 (k) = (k 1 ) + (k 2 ), with k = k 1 + k 2 . For the conventional nearest-neighbor Heisenberg chain, the resulting spectrum is bounded from below by the des Cloizeaux-Pearson dispersion πJ/2| sin k| 43 , and the upper boundary of the continuum is πJ| sin (k/2)| 44 . Therefore, it will be characterized by singularities and will not realize coherent quasiparticles, that in the spectrum would appear as δ-like peaks, accompanied by and incoherent background at high energies (that can correspond to spinons or a two-magnon continuum). Magnons are associated to symmetry breaking and the emergence of gapless Goldstone modes after some gapped mode condenses. However, it is possible to transition from a gapless phase without long range order (a gapless spin liquid) to an ordered one with a well defined order parameter. In this case, it is expected that the gapless deconfined excitations of the spin liquid will form bound states in the ordered phase. For this to happen, an attractive confining potential should be strong enough to overcome the kinetic energy of the free spinons. This is precisely what occurs in our model.
In the conventional 1D Heisenberg model, flipping a spin would create two domain walls. The energy cost of such excitation does not scale with the separation between the particles. However, in the case of long range staggered interactions, all spins interact with each other and the local disturbance is felt by the bulk of the chain. Separating the domain walls costs an extensive amount of energy, as depicted in Fig. 1(b), where we show the confining potential for different values of α and λ = 1.
We are interested in determining signatures of confinement in the excitation spectrum of the model. In order to obtain the spin dynamical structure factor we used the time-dependent DMRG method (tDMRG) 39-42 following the prescription detailed in the original work Ref. 39. The idea consists of calculating the two-time spin-spin correlator: where S z 0 here is defined at the center of the chain, and r is the distance from center. Fourier transforming to momentum space and frequency, we reconstruct the momentum resolved spectral function. The Fourier transform is carried out over a finite time range (in our case t max = 20, which requires the use of a windowing technique to attenuate artificial ringing (satellite oscillations associated to the natural frequencies that lead to artifacts, such as negative values). The spectrum will exhibit an artificial broadening that is inversely proportional to the width of our time window. In addition, good resolution at high frequencies can be improved by using a small time-step, while a long time-window is necessary to improve resolution at low frequencies. In order to timeevolve the wave function, we use a time-step targeting procedure with a Krylov expansion of the time-evolution operator 45 and a time step δt = 0.05 (time is measured in units of J −1 and J is our unit of energy).
We study chains of length L = 48 using 400 DMRG states that guarantees that the truncation error remains smaller than 10 −6 over the time window. Results for the dynamical structure factor are displayed in Fig.2 for λ = 1 and different values of α across the phase transition. The spectrum is bounded from below by a sharp peak, which for large α > 2.2 corresponds to the edge of the two spinon continuum. For smaller α, as we cross over to the ordered phase, the peak splits out from the continuum, that moves to higher energies. This can be seen more clearly in Fig.3, where we show cuts along the k = π/2 direction. The splitting of the peak and the shifting of the continuum to higher energies are signatures of the formation of bound states, which become coherent quasiparticles in the symmetry broken phase.

III. TWO-SPINON BOUND STATES
In order to develop intuition on the nature of the excitations and the confining mechanism, we study a related toy problem that will serve as a close approximation to the present situation. We will follow Villain 46 and assume that the ground state of the system has uni-axial FIG. 5. Two-spinon spectra for λ = 1 and different values of α obtained by using Villain's approach for L = 60, as described in the text. symmetry (the Ising limit), and consider the dynamics of mobile domain walls. This is done by considering the motion of the spinons by means of spin flips, ignoring the action of these terms on pairs of spins that do not involve domain walls. The procedure was clearly outlined in Ref. 47 but, in our case, we need to consider the longrange nature of the interactions (a similar procedure was carried out for the ferromagnetic case in Refs.48 and 49). We firstly define the space spanned by all the possible configurations with two spinons in an antiferromagnetic background |i, j , where i,j are the positions of the domain walls, illustrated by Fig.4(a), (b) and (c). We then exploit translational symmetry and define wave functions in a sector with momentum k = 2πn/L (n = 0, · · · L − 1): where the translation operator acts as T d |i, j = |i+d, j + d (periodic boundary conditions enforce position to be defined mod (L)). The Hamiltonian matrix elements are easily calculated and each momentum sector is diagonalized independently. Results for the two spinon excitation spectrum for chains of length L = 60 are shown in Fig.5 for λ = 1 and several values of α. We plot the energy of the two spinon state E(k). For small α we see several bound states leaking out of the two spinon continuum. As α increases above the expected value for the transition in the isotropic case α c 2.2, the continuum tends to collapse and merge with the two-spinon bound states. This scenario agrees qualitatively with the observed behavior in the SU (2) symmetric case. For small α, the two spinon continuum is pushed to higher energies, and is clearly separated from the "magnon" band. So far, our approach has ignored other possibilities that can be realized through long range spin flips. In fact, it is easy to convince oneself that two-spinon configurations can only propagate via nearest-neighbor spin-flips. In order to account for the long-range off-diagonal processes we have to consider states with three spinons and one anti-spinon, as illustrated in Fig.4(d) and (e). As a matter of fact, long range spin flips can create a proliferation of spinons and anti-spinons throughout the chain, but this is prevented by energetic considerations. For this reason and due to the fast growth of the number of possible configurations (scaling as ∼ L 3 ), we only preserve those with one anti-spinon, and even so we can solve for sizes up to L = 36. The spectrum is now modified as seen in Fig.6, but the low energy features remain qualitatively similar. However, we see a continuum of high energy states corresponding to the new sector that gets mixed with the two-spinon continuum. A significant different that we observe between these results and those in Fig.5 is that the number of bands leaking out of the continuum is suppressed. It is reasonable to expect that as more spinon and anti-spinon states are included, only one band will finally survive. Clearly, the low-energy physics is well described as consisting of bound states of spinons but, as we shall see, accounting for the additional spin fluctuations becomes important when it comes to understanding the real-time evolution of the system.

IV. REAL-TIME EVOLUTION
It is natural to ask whether the spinon confinement can be identified in a numerical "time-of-flight" experiment, in which a spinon is created at the center of a chain by FIG. 7. Domain wall expansion for λ = 1 and different values of α, obtained with tDMRG for a chain with L = 48 spins. Results for the nearest neighbor case are also included. We show the "spinons density": N ↑ (r, t)N ↑ (r + 1, t) . Color density is in a log scale.
the application of the S + operator, and left to evolve under the action of the Hamiltonian. Results obtained with tDMRG as shown in Fig.7 where we plot the correlations N ↑ (r, t)N ↑ (r + 1, t) , where N ↑ = (2S z + 1). Notice that bound states are not necessarily localized in nearest neighbors, but actually are extended objects that have a characteristic size 37 that gets smaller with decreasing α. However, to a certain extent, these correlations can help us develop intuition and, for this purpose, we shall refer to them as the "spinon density". Without long-range interactions, spinons propagate ballistically 50,51 and this is seen in Fig.7 as a "lightcone" with a velocity determined by the maximum slope of the spinon dispersion v = π/2. As the value of alpha decreases, spinons become more and more confined and, consequently, "heavier" (or slower), since bound states of spinons move coherently through second-order processes by means of two consecutive spin flips. However, a side effect of the long-range interactions is that spinons can "hop" longer distances. As a consequence, the originally sharp edges of the lightcone become more diffuse and, moreover, they acquire an apparent curvature that gives the misleading impression of an underlying "accelleration". As it turns out, this illusion is due to the superposition of two characteristic velocities, as we discuss next.
Our first attempt to explain the observed behavior is to use the Villain approximation with two spinons, presented in Fig.8. As we mentioned earlier, the long range interactions in this case enter only as a diagonal contribution, since long range spin flips produce a proliferation of spinons that take us outside of the two-spinon sector. However, we can already identify very interesting features as we change the exponent from basically α → ∞ (nearest neighbor interactions, only), to α = 1.8. In the former case, we see a coherent ballistic propagation of the spinons with a well defined characteristic velocity, as expected, since the problem is equivalent to two non-interacting particles. For α = 2.2 we identify two coexisting lightcones: a fainter one preserves the same slope as the free deconfined spinons, while the second one, with larger weight, describes coherent particles moving with roughly half the spinon velocity. It only makes sense to attribute these features to a bound state of spinons, a "magnon". As we reduce α even further, the free spinon lightcone loses weight, which is transferred to the magnons. For α = 1.8 we can clearly identify a single magnon lightcone, as free spinons move to higher energies. In order to support these observations we also calculate N ↑ (r, t)N ↑ (r + 1, t)N ↑ (r + 2, t) , or "magnon density", in Fig.9. In the nearest-neighbor limit we only see a very faint feature that loses weight as time evolves: the original flipped spin creates a state like the one depicted in Fig.4(b), but it is short lived and breaks into two spinons. As α decreases, the magnon lightcone becomes more and more coherent and correlates exactly with the features observed in Fig.8.
Having determined the coexistence of deconfined spinons and magnons in the lightcone, it rests to explain the apparent curvature in the DMRG results. For this, we need to extend our treatment by considering the possibility of three spinons and one anti-spinon to account for the long-range spin flips. The results for the spinon and magnon densities are presented in Fig.10 and Fig.11, respectively. We basically observe a "fan" or excitations covering the region between the magnon and free spinon wavefronts. Moreover, attempting to identify a characteristic velocity is an ill-defined problem, since spins are allowed to hop to all distances. This is also reflected in the magnon dispersion no longer having a linear dispersion, as previously observed.

V. SUMMARY AND CONCLUSIONS
We have analyzed the spectrum and studied the nature of the excitations of a Heisenberg chain with staggered long range interactions. The unfrustrated long-range nature of the exchange effectively increases the dimensionality of the system and the chain is able to undergo true symmetry breaking and develop long range order. For FIG. 11. Magnon density N ↑ (r, t)N ↑ (r + 1, t)N ↑ (r + 2, t) for λ = 1 and different values of α, obtained using Villain's approximation including three spinons and one anti-spinon.
weakly decaying interactions, our tDMRG calculations show that the emergence of Néel order can be associated to the formation of bound states of spinons that become coherent quasiparticles (magnons). At the same time, the two-spinon continuum is pushed to higher energies. This is supported by two-spinon and three-spinon approximations that reproduce the main features and explain the formation of bound states due to a confining potential that grows logarithmicaly. The observed behavior bears very close resemblance to the one found in actual neutron experiments in higher dimensional materials.