# Dynamical structure factor of the Heisenberg model in one dimension: the variational Monte Carlo approach

## Abstract

The dynamical spin structure factor is computed within a variational framework to study the one-dimensional Heisenberg model. Starting from Gutzwiller-projected fermionic wave functions, the low-energy spectrum is constructed from two-spinon excitations. The direct comparison with Lanczos calculations on small clusters demonstrates the excellent description of both gapless and gapped (dimerized) phases, also including incommensurate structures for . Calculations on large clusters show how the intensity evolves when increasing the frustrating ratio and give an unprecedented accurate characterization of the dynamical properties of (non-integrable) frustrated spin models.

## I Introduction

Quantum spin liquids are unconventional phases of matter in which quantum fluctuations endure any tendency to develop local (e.g., magnetic) order, down to zero temperature. Most importantly, the lack of any symmetry breaking mechanism is accompanied by topological order and the presence of excitations with fractional quantum numbers balents2010 (). For concreteness, let us consider the spin- Heisenberg model on a given two-dimensional lattice. Whenever magnetic order sets in (e.g., on the square lattice with nearest-neighbor antiferromagnetic interactions), the elementary excitations are magnons or spin waves, which can be pictured as coherent Bloch waves made from localized spin flip excitations. Instead, in spin liquids, the elementary objects are spinons, while excitations decay in two spinons that are asymptotically free at long distances balents2010 (). Since, for any lattice with fixed number of sites, the minimal change of the total spin is , the existence of objects with implies a fractionalization of the spin quantum number. Spinons do exist in the one-dimensional Heisenberg model faddeev1981 (), where the magnetic long-range order is hampered in agreement with the generalized Mermin-Wagner theorem pitaevskii1991 (). In this case, the spectrum is gapless and characterized by the presence of a broad continuum of excitations. A similar feature is present in the one-dimensional Heisenberg model with inverse-square superexchange haldane1988 (); in this case, spinons are non interacting and the whole excitation spectrum can be found explicitly in a closed form haldane1991 (). Moreover, objects are elementary excitations also in gapped systems, as in the case of the Majumdar-Ghosh point of the frustrated Heisenberg model (where the ratio between the first-neighbor coupling and the second-neighbor one is equal to ). Here, the ground state is doubly degenerate (with long-range dimer order) majumdar1969 () and elementary excitations can be seen as propagating defect boundaries between the two ground states (they are analogous to solitons) shastry1981 (). Anyhow, in one-dimensional systems, fractional excitations are ubiquitous and represent the general feature of spin models. By contrast, in two spatial dimensions, neat examples of spin liquids are rarer and the possibility to have free (i.e., deconfined) spinons when magnetic order is destroyed is not taken for granted. A beautiful example in which fractional excitations are present is given by the Kitaev (compass) model on the honeycomb lattice kitaev2006 (), where elementary excitations are Majorana fermions and gauge fluxes.

In the last 20 years, a huge effort has been devoted to assess the low-energy behavior of various frustrated spin models, in order to unveil possible spin-liquid ground states. Many different lattice structures and various kind of interactions, including long-range superexchange or ring-exchange terms, have been considered by employing a large variety of analytical and numerical techniques lacroixbook (). Nevertheless, most of these studies focused on the ground-state properties, by computing correlation functions or different quantities that may give information about the presence/absence of a spin gap. The direct computation of the spin gap has been performed in few cases, notably for the the Heisenberg model on the kagome lattice where this issue has been addressed by density-matrix renormalization group yan2011 (); depenbrock2012 (), variational Monte Carlo iqbal2014 (), and Lanczos lauchli2016 () approaches. By contrast, only few works, mainly based upon analytical or semi-analytical approaches, have focused on the whole dynamical properties of frustrated spin models fuhrman2012 (); mourigal2013 (); zhitomirsky2013 (); ghioldi2018 (). In this respect, an important quantity that gives direct access to the nature of the excitation spectrum is the dynamical structure factor:

(1) |

where is the ground state of the system with energy , are the excited states with momentum (relative to the ground state) and energy , and

(2) |

is the Fourier transformed spin operator for the components , being the number of sites of the cluster. In this regard, inelastic neutron scattering provides a direct measurement of the dynamical structure factor, as a function of momentum and energy . Therefore, the theoretical computation of has an immediate connection to experiments, thus validating/disproving the modellization of a given real material coldea2001 (); ronnow2001 (); lake2013 ().

Unfortunately, the theoretical evaluation of the dynamical structure factor is possible only for very limited cases. In one dimension, an exact calculation is possible for the Heisenberg model with inverse-square superexchange haldane1993 (), while very accurate approximations are now available for the simple Heisenberg model with nearest-neighbor interaction karbach1997 (); caux2005 (); caux2006 (). In two spatial dimensions, an exact computation of the dynamical structure factor is possible for the Kitaev model in both the gapless and gapped regimes knolle2014 (). Besides these fortunate cases in which the model is integrable, numerical techniques have been employed to study generic models, especially in one dimension. Here, exact diagonalizations can be performed on relatively small clusters yokoyama1997 () and their results can be compared with semi-analytical calculations lavarelo2014 (). Moreover, density-matrix renormalization group barthel2009 () or matrix-product states xie2018 () can be used. Alternatively, a variational technique, implemented within a quantum Monte Carlo method, has been suggested to approximate the exact sprectrum with states for each momentum li2010 (). The important advantage of this approach is that the dynamic structure factor is directly accessible, without any computation requiring (unstable) transformations from real/imaginary times to frequencies. Curiously, this approach has been applied in very few cases dallapiazza2015 (); mei2015 (), without systematic benchmarks and comparisons with other methods.

In this paper, we employ the variational approach that has been introduced in Ref. li2010 () to compute the dynamical structure factor of Eq. (1) for momentum and energy in the spin- Heisenberg model in one dimension:

(3) |

where are the (integer) coordinates of the sites and is the spin operator on the site . Periodic boundary conditions are considered in the spin Hamiltonian. In the following, we will consider the case with in Eq. (1), since, given the symmetry of the Heisenberg model, any component of the structure factor gives the same result. The variational Monte Carlo results are compared with Lanczos diagonalizations on a small cluster, in order to show the accuracy of the method for different values of the frustrating ratio . Then, calculations are reported for large systems, illustrating how the various features of the dynamical structure factor evolve from the gapless to the gapped phase, also entering in the incommensurate region with .

## Ii Variational method

The variational approach is based on a Gutzwiller projected fermionic wave function, which is constructed from an auxiliary superconducting (BCS) Hamiltonian:

(4) | |||||

here, () creates (destroys) an electron with spin on the site ; and are hopping and singlet pairing terms, respectively. This Hamiltonian is quadratic in the fermionic operators and, therefore, can be easily diagonalized. Its ground state is denoted by . Of course, this quantum state is not suitable to describe a spin system, since it is defined in the “enlarged” Hilbert space with also empty and doubly occupied sites. Then, a suitable variational wave function for the Heisenberg model of Eq. (3) can be obtained by projecting out all the configurations with at least one empty or doubly occupied site:

(5) |

where ( being the local electron density per spin on site ) is the Gutzwiller projector, which enforces single fermionic occupation on each site. It has been shown that Gutzwiller-projected fermionic wave functions are very accurate to describe the exact ground state of the one-dimensional Heisenberg model with , as well as the lowest-energy spinon excitations yunoki2006 ().

The parameters of , i.e., hopping and pairing amplitudes, are taken to be real and fully optimized by mean of the stochastic reconfiguration technique, in order to minimize the variational energy of sorella2005 (). In most of the calculations, we will impose translational symmetry in the quadratic Hamiltonian , which strongly reduces the number of independent parameters to be treated. We must emphasize the fact that both periodic boundary conditions (PBC) and anti-periodic boundary conditions (APBC) are allowed within the auxiliary BCS Hamiltonian (leading to a real wave function). However, while in presence of a gapped fermionic spectrum either options will lead to a unique ground state, the same may not be true for a gapless spectrum. For example, if there are gapless points at , the ground state is unique if PBC (APBC) are considered for (), where is an integer.

In order to tackle the problem of computing , we need to devise a way to construct excited states. Following the procedure of Ref. li2010 (), we first introduce a set of two-spinon triplet excitations with momentum , which are obtained by the Gutzwiller projection of particle-hole fermionic excitations:

(6) |

Notice that, since the Gutzwiller projector commutes with , we have that . Then, for each momentum , defines a (non-orthogonal) basis set that can be used to approximate the exact low-energy excitations. In other words, we can define a set of states note1 () (for each momentum ), which are labeled by :

(7) |

The coefficients are obtained by diagonalizing the Heisenberg Hamiltonian within the subspace generated by for each , namely by solving the generalized eigenvalue problem:

(8) |

where we have introduced two matrices, (Hamitonian) and (overlap). Finally, the dynamical structure factor of Eq. (1) is approximated by taking:

(9) |

where, compared to the exact form of Eq. (1), the variational states and are considered (instead of the exact eigenstates), and the variational energies (corresponding to ) and are taken (instead of the exact ones). Most importantly, the sum over excited states runs over at most states (instead of an exponentially large number). By using Eq. (7), we have:

(10) |

Remarkably, all the quantities that define the dynamical structure factor of Eq. (10) can be computed within a variational Monte Carlo scheme (i.e., without any sign problem). In fact, the entries of the Hamiltonian and overlap matrices are given by:

(11) | |||||

(12) |

where is a set of normalized and orthogonal states, which can be sampled by using the variational wave function as the probability distribution:

(13) | |||||

(14) |

At this stage, it is worth making two remarks. First of all, our sampling procedure is possible because both the ground-state wave function and the particle-hole excitations of Eq. (6) have and, therefore, the set of configurations can be chosen to also have . This would not be possible whenever considering excitations involving a spin flip. For this case, a different sampling procedure has been proposed in Ref. li2010 (). The most important advantage of our approach is that all the values of the momentum can be computed with a single Monte Carlo simulation, at variance with the previous technique, where each needs a separate calculation. Second, within our formulation, the sampling is correct whenever the ground-state wave function is nonzero for all the configurations ; otherwise, the sampling procedure neglects the contributions from these “vanishing” configurations. We checked that, for most of the cases we considered, the number of these configurations is negligible and, therefore, they do not affect the final results.

We emphasize that, within this procedure, once the ground-state wave function is optimized, the only remaining parameters are the coefficients , which are completely determined by solving Eq. (8). In other words, the particle-hole excitations are applied to a fixed reference state, i.e., , which is optimized, once for all, to minimize the ground-state variational energy.

The evaluation of and essentially boils down to computing the following quantities:

(15) | |||||

(16) |

Once and have been evaluated, the generalized eigenvalue problem (8) for the excited states of momentum can be solved. In doing so, we need to get rid of the linear dependence, which may affect the set . This is achieved by restricting Eq. (8) to the subspace of eigenvectors of the overlap matrix that have non-zero eigenvalues. In practice, since the entries are computed stochastically, we need to discard the eigenvectors whose eigenvalues are smaller than some given tolerance.

## Iii The ground-state variational wave function

The phase diagram of the one-dimensional Heisenberg model is well-known white1996 (): for small values of the frustrating ratio, the system is gapless (i.e., a Luttinger fluid) with power-law spin-spin correlations, while for large values of , the system is in a gapped phase characterized by long-range dimer order. In addition, (short-ranged) spin-spin correlations show an incommensurate periodicity for . The critical (Kosteritz-Thouless) point that separates gapless and gapped phases has been estimated with high accuracy, eggert1996 ().

The simplest Ansatz that can be used to describe both phases is the one obtained from a pure hopping Hamiltonian with broken translational symmetry. This can be achieved by doubling the unit cell and taking different intra-cell () and inter-cell () hoppings. When , recovers translational invariance and reduces to the case of free fermions in one-dimension, which have a Fermi sea ground state and gapless excitations. Instead, when , there are two fermionic bands separated by a finite gap. The uniform and dimerized states are dubbed UFS and DFS, respectively. The accuracy for a cluster with sites is shown in the top panel of Fig. 1. For the optimal wave function does not break the translational symmetry (i.e., ); by contrast, for larger values of the frustrating ratio, . At the Majumdar-Ghosh point (), one of the two hopping parameters is equal to zero, indicating that the wave function is a product of nearest-neighbor singlets. Here, the variational state becomes exact. Actually, the fully dimerized wave function remains the optimal solution for , but its accuracy quickly worsens, since its energy is independent on .

More accurate wave functions can be built from translationally invariant Ansätze, which include both hopping and pairing terms (with and ). Nonetheless, even by considering translational symmetry, a “sponteneous symmetry breaking” mechanism is possible after Gutzwiller projection is included, leading for example to dimer order becca2011 (); kaneko2016 (). Within a gapless regime, an extremely accurate state, dubbed WFA is constructed from a fermionic Hamiltonian that contains first- and third-neighbor hoppings ( and ), as well as first-neighbor pairing (). This choice gives a gapless fermionic band at and can be stabilized up to . For larger values of the frustrating ratio, the pairing term goes to zero and the wave function coincides with the UFS state. A different possibility, which allows the existence of a gap in the fermionic spectrum, is given by taking first-neighbor hopping and both onsite () and second-neighbor () pairings. This Ansatz, which is dubbed WFB, is gapped unless . Optimizing the parameters of this wave function for sites, we find that it reduces to the simple UFS state (i.e. ) for . Then, the optimal pairing terms become non-zero and the wave function proves to be more accurate than the DFS state and stable for all the values of the frustrating ratios, see Fig. 1.

We expect that, in the thermodynamic limit, the gap in the fermionic spectrum opens in the vicinity of the exact transition point and follows the behavior of the spin gap. However, it is extremely hard to locate this point by performing a finite size-scaling analysis, since the gap is exponentially small in an extended region after the critical point.

## Iv Results

Here, we present the numerical results for the spin structure factor of a one-dimensional chain. Let us start by considering a small cluster with sites, where exact diagonalizations are possible by using the Lanczos method. First of all, we demonstrate that the variational results do not change appreciably when considering the wave function WFA or WFB to compute the dynamical structure factor, see Fig. 2. Indeed, even though the latter state is about five times less accurate than the former one for (see Fig. 1), the actual differences between the two dynamical calculations are negligible (and either option gives an excellent description of the exact results). Therefore, in the following, we will consider only the WFB wave function to compute the dynamical structrure factor.

In order to best quantify the agreement between the variational and the exact calculations, we directly report for several momenta as a function of the frequency for two values of the frustrating ratio, and , see Figs. 3 and 4. Here, the delta-functions related to the exact and variational energies entering the Eqs. (1) and (9) have been replaced by normalized Gaussians with . The agreement is very good, not only for the unfrustrated case with (see Fig. 3), but also in presence of a sizable frustration, (see Fig. 4). Similar results are also obtained for larger values of the ratio (see below). Therefore, it is expected that, within this approach, both gapless and gapped regimes are correctly described. The accuracy of the variational method is highlighted in Fig. 5, where we report the overlaps between the variational excited states of Eq. (7) and the exact ones. In particular, whenever the excited states are well separated in energy, it is easy to match each exact excitation with a corresponding variational one (in this case, the overlap is very large, as for ). Instead, when two or even more excitations are close in energy, this correspondence is not easy to resolve and, for each variational state, we computed the overlap with all the exact states and plotted the maximum value. In this case, a reduction of the overlap is observed, as for a number of cases at . Nevertheless, in all cases, the relevant excitations, which carry sizable spectral weigth, are well reproduced by the variational approach and a reduced overlap is detected for states which do not contribute much to the whole intensity of the dynamical structure factor.

Finally, the results obtained with the exact and the variational approaches are compared in Fig. 6, where of a chain of sites is represented using color maps for , , and . In all the cases, the variational results follow the exact ones, including the development of incommensurate features when . In fact, by increasing the frustrating ratio, the intensity progressively shifts from (low energies) to (high energies). At even larger values of , the modes at soften and eventually become gapless for (in this limit, the spin Hamiltonian consists of two decoupled Heisenberg models, one for each sublattice, and represents a nearest-neighbor superexchange on each sublattice). Remarkably, the variational approach is able to perfectly reproduce all the relevant features of the dynamical structure factor. We mention the fact that the only case where our sampling technique fails is at the Majumdar-Ghosh point , where the number of vanishing configurations in the ground-state wave function (exactly reproduced by our Gutzwiller-projected fermionic state) is exponentially large.

The results for a large cluster with sites are reported in Fig. 7 for , , , , , and . In the unfrustrated case, it is known karbach1997 (), that most of the total intensity of the dynamical structure factor is carried by the two-spinon contributions. For these excitations the lower and upper energy limits are given by:

(17) | |||

(18) |

Indeed, we find that our dynamical structure factor is bounded by these limits and closely resembles the one that has been recently obtained with a Bethe Ansatz approach lake2013 (); caux2006 ().

It should be stressed that, for a relatively large region within the gapped phase, the value of the spin gap remains very small, since the transition from the gapless to the dimerized phase belongs to the Kosterlitz-Thouless universality class. Therefore, even for a relatively large system size, it is very hard to detect the presence of a finite gap in the excitation spectrum: for example, the dynamical structure factors at and (see Fig. 7) look very similar, even though the former case corresponds to a gapless phase and the latter one to a gapped spectrum. On this large cluster, the gradual shift of the intensity from to is evident, as well as the presence of a “rounding” around within the gapped phase for . Within such a large size, incommensurate features appear clearly for , namely, the excitations with lowest energy move from to , giving rise to a non-trivial form of the spectral function.

## V Conclusions

In conclusion, by using a variational approach li2010 () that does not involve any possible instability (which can be due to sign-problem sampling and/or analytic continuation from imaginary/real times to frequencies), we have studied the dynamical spin structure factor of the frustrated model in one spatial dimension. We have reported the unprecedent accuracy of this method, not only at or close to the integrable point with , but also for generic values of the frustrating ratio. The remarkable advantage of this variational procedure is given by the fact that the relevant part of the low-energy spectrum can be described by considering particle-hole excitations on top of a fixed “reference” state. Indeed, once the variational wave function has been optimized for the ground state, only parameters for each are used to reproduce the low-energy part of the spectrum. This fact suggests that Gutzwiller-projected fermionic wave functions not only may accurately reproduce the ground-state properties of frustrated spin models, but also constitute a good framework to generate low-energy excitations.

### References

- L. Balents, Nature (London) 464, 199 (2010).
- L.D.Faddeev and L.A.Takhtajan, Phys. Lett. A 85, 375 (1981).
- L. Pitaevskii and S. Stringari, J. Low Temp. Phys. 85, 377 (1991).
- F.D.M. Haldane, Phys. Rev. Lett. 60, 635 (1988); B.S. Shastry, Phys. Rev. Lett. 60, 639 (1988).
- F.D.M. Haldane, Phys. Rev. Lett. 66, 1529 (1991).
- C.K. Majumdar and D. Ghosh, J. Math. Phys. 10, 1388 (1969).
- B.S. Shastry and B. Sutherland, Phys. Rev. Lett. 47, 964 (1981).
- A. Kitaev, Ann. Phys. (N.Y.) 321, 2 (2006).
- See for example, C. Lacroix, P. Mendels, and F. Mila, Introduction to Frustrated Magnetism Materials, Experiments, Theory (Springer, 2013).
- S. Yan, D.A. Huse, and S.R. White, Science 332, 1173 (2011).
- S. Depenbrock, I.P. McCulloch, and U. Schollwöck, Phys. Rev. Lett. 109, 067201 (2012).
- Y. Iqbal, D. Poilblanc, and F. Becca, Phys. Rev. B89, 020407(R) (2014).
- A.M. Läuchli, J. Sudan, R. Moessner, arXiv:1611.06990 (unpublished).
- W.T. Fuhrman, M. Mourigal, M.E. Zhitomirsky, and A.L. Chernyshev, Phys. Rev. B85, 184405 (2012).
- M. Mourigal, W. T. Fuhrman, A. L. Chernyshev, and M.E. Zhitomirsky, Phys. Rev. B88, 094407 (2013).
- E.A. Ghioldi, M.G. Gonzalez, S. Zhang, Y. Kamiya, L.O. Manuel, A.E. Trumper, and C.D. Batista, arXiv:1802.06878 (unpublished).
- M.E. Zhitomirsky and A.L. Chernyshev, Rev. Mod. Phys. 85, 219 (2013).
- R. Coldea, S.M. Hayden, G. Aeppli, T.G. Perring, C.D. Frost, T.E. Mason, S.-W. Cheong, and Z. Fisk, Phys. Rev. Lett. 86, 5377 (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).
- B. Lake, D.A. Tennant, J.-S. Caux, T. Barthel, U. Schollwöck, S.E. Nagler, and C.D. Frost, Phys. Rev. Lett. 111, 137205 (2013).
- F.D.M. Haldane and M.R. Zirnbauer, Phys. Rev. Lett. 71, 4055 (1993).
- M. Karbach, G. Müller, A.H. Bougourzi, A. Fledderjohann, and K.-H. Mütter, Phys. Rev. B55, 12510 (1997).
- J.-S. Caux and J.M. Maillet, Phys. Rev. Lett. 95, 077201 (2005).
- J.-S. Caux and R. Hagemans, J. Stat. Mech. (2006) P12013.
- J. Knolle, D.L. Kovrizhin, J.T. Chalker, and R. Moessner, Phys. Rev. Lett. 112, 207203 (2014); Phys. Rev. B92, 115127 (2015).
- H. Yokoyama and Y. Saiga, J. Phys. Soc. Jpn 66, 3617 (1997).
- A. Lavarélo and G. Roux, Eur. Phys. J. B 87, 229 (2014).
- T. Barthel, U. Schollwöck, and S.R. White, Phys. Rev. B79, 245101 (2009).
- H.D. Xie, R.Z. Huang, X.J. Han, X. Yan, H.H. Zhao, Z.Y. Xie, H.J. Liao, and T. Xiang, Phys. Rev. B97, 075111 (2018).
- T. Li and F. Yang, Phys. Rev. B81, 214509 (2010).
- B. Dalla Piazza, M. Mourigal, N.B. Christensen, G.J. Nilsen, P. Tregenna-Piggott, T.G. Perring, M. Enderle, D.F. McMorrow, D.A. Ivanov, and H.M. Rønnow, Nature Physics 11, 62 (2015).
- J.-W. Mei and X.-G. Wen, arXiv:1507.03007 (unpublished).
- S. Yunoki and S. Sorella, Phys. Rev. B74, 014408 (2006).
- S. Sorella, Phys. Rev. B71, 241103(R) (2005).
- In many cases, the set of states is linearly dependent. This reduces the effective number of excitations which can be constructed.
- S.R. White and I. Affleck, Phys. Rev. B54, 9862 (1996).
- S. Eggert, Phys. Rev. B54, 9612 (1996).
- F. Becca, L. Capriotti, A. Parola, and S. Sorella, Variational Wave Functions for Frustrated Magnetic Models, Springer Series in Solid-State Sciences edited by C. Lacroix, P. Mendels, and F. Mila (Springer, 2011), p. 379.
- R. Kaneko, L.F. Tocchio, R. Valenti, F. Becca, and C. Gros, Phys. Rev. B93, 125127 (2016).