Dynamical structure factor of the $J_1-J_2$ Heisenberg model on the triangular lattice: magnons, spinons, and gauge fields

Understanding the nature of the excitation spectrum in quantum spin liquids is of fundamental importance, in particular for the experimental detection of candidate materials. However, current theoretical and numerical techniques have limited capabilities, especially in obtaining the dynamical structure factor, which gives a crucial characterization of the ultimate nature of the quantum state and may be directly assessed by inelastic neutron scattering. In this work, we investigate the low-energy properties of the $S=1/2$ Heisenberg model on the triangular lattice, including both nearest-neighbor $J_1$ and next-nearest-neighbor $J_2$ super-exchanges, by a dynamical variational Monte Carlo approach that allows accurate results on spin models. For $J_2=0$, our calculations are compatible with the existence of a well-defined magnon in the whole Brillouin zone, with gapless excitations at $K$ points (i.e., at the corners of the Brillouin zone). The strong renormalization of the magnon branch (also including roton-like minima around the $M$ points, i.e., midpoints of the border zone) is described by our Gutzwiller-projected state, where Abrikosov fermions are subject to a non-trivial magnetic $\pi$-flux threading half of the triangular plaquettes. When increasing the frustrating ratio $J_2/J_1$, we detect a progessive softening of the magnon branch at $M$, which eventually becomes gapless within the spin-liquid phase. This feature is captured by the band structure of the unprojected wave function (with $2$ Dirac points for each spin component). In addition, we observe an intense signal at low energies around the $K$ points, which cannot be understood within the unprojected picture and emerges only when the Gutzwiller projection is considered, suggesting the relevance of gauge fields for the low-energy physics of spin liquids.


I. INTRODUCTION
The antiferromagnetic Heisenberg model for S ¼ 1=2 spins interacting on the triangular lattice represents the simplest example in which quantum fluctuations give rise to strong modifications of the classical picture, where the minimum energy configuration shows 120°order. Indeed, this was the first microscopic model that has been proposed for the realization of the so-called resonating valence-bond state [1,2]. Within this approach, the ground state is described by a superposition of an exponentially large number of singlet coverings of the lattice, generalizing the concept of resonance introduced and developed by Rumer [3] and Pauling [4] to describe the chemical bond. Even though recent numerical investigations [5,6] have shown that the ground state possesses a finite magnetization in the thermodynamic limit, the results confirmed large deviations from classical and semiclassical limits. In addition, small perturbations on top of the nearest-neighbor Heisenberg model have been shown to drive the system into magnetically disordered phases [7,8]. By keeping the spin SU(2) symmetry, a natural way to induce further magnetic frustration is to include a next-nearest-neighbor superexchange coupling, leading to the following Hamiltonian: where hÁ Á Ái and ⟪ Á Á Á ⟫ indicate nearest-neighbor and next-nearest-neighbor sites in the triangular lattice, S i ¼ ðS x i ; S y i ; S z i Þ is the spin-1=2 operator at the site i, and, finally, J 1 and J 2 are the antiferromagnetic coupling constants. This model has been intensively investigated in the past, from the semiclassical approaches of the early days [9,10] to the recent numerical approaches [11][12][13]. The latter ones indicated a rather fragile 120°magnetic order, which is melted for J 2 =J 1 ≈ 0.07ð1Þ (a value that is in very good agreement among these calculations). For larger values of the frustrating ratio J 2 =J 1 , the nature of the nonmagnetic phase is not settled, with evidences for either a gapped [11,12] or a gapless [13] spin liquid.
An important information about the physical properties is given by the features of the low-energy spectrum. In particular, the dynamical structure factor Sðq; ωÞ gives a direct probe to assess the nature of the relevant excitations. These can be divided into two broad classes: standard gapless magnons (or gapped triplons), which exist in magnetically ordered phases (or valence-bond solids), and more exotic (gapped or gapless) spinons, which exist in deconfined spin liquids. In addition to spinons, another kind of excitation is present, due to the emergence of gauge fluctuations in the low-energy effective theory of spin liquids [14].
For the Heisenberg model with only nearest-neighbor couplings on the triangular lattice, semiclassical approaches, based upon the large-S expansion, suggested that the excitation spectrum obtained within the leading order (i.e., within the linear spin-wave approximation) is subjected to significant corrections when interactions between spin waves are taken into account [15]. This fact is mainly due to the noncollinearity of the magnetization, which allows for three-magnon interactions. Then, despite the presence of long-range order, the Goldstone modes are not stable, but they may decay in a large part of the Brillouin zone (see Fig. 1); in particular, the existence of more than one Goldstone mode, with different velocities, immediately implies that magnons may be kinematically unstable, decaying into two magnons with lower energy [16,17]. A detailed analysis, which includes interactions among spin waves, corroborated this outcome, also showing rotonlike minima at M ¼ ð0; 2π= ffiffi ffi 3 p Þ and symmetry-related points (i.e., midpoints of the edges of the Brillouin zone) [16][17][18]. The latter aspect shares similarities with the Heisenberg model on the square lattice, where minima of the magnon dispersion are present around ðπ; 0Þ and ð0; πÞ [19,20]. As far as the triangular lattice is concerned, aspects of the strong renormalization of the magnon dispersion at high energies have been confirmed by series expansions [21]. Moreover, within these numerical calculations, a huge downward renormalization of the onemagnon excitations is recovered, leading to a relatively dispersionless mode.
While there are a number of materials whose low-energy behavior can be well described by the S ¼ 1=2 Heisenberg model on the square lattice (among them, we just mention La 2 CuO 4 for its relevance to cuprate superconductors [22]), until very recently there were no compounds that could be well approximated by the same model on the equilateral triangular lattice. For example, in Cs 2 CuCl 4 the superexchange couplings are not isotropic in the nearestneighbor bonds, one out of the three being much stronger than the other ones (thus defining weakly coupled zigzag chains) [23]. Here, inelastic neutron scattering measurements have shown the existence of a very broad continuum, which has been associated to spin fractionalization and spin-liquid behavior [23].
Recently, measurements on Ba 3 CoSb 2 O 9 have been reported, providing evidence that it can be described by a S ¼ 1=2 Heisenberg model on the undistorted triangular lattice with predominant nearest-neighbor superexchange couplings (a small easy-plane anisotropy is present, in addition to a small interlayer coupling) [24]. The initial interest was aimed at the study of the magnetization curve and the stabilization of magnetization plateaus [24,25], and the proximity to a spin-liquid phase [26]. Later, inelastic neutron scattering measurements were performed, in order to clarify the nature of the magnetic excitations on top of the ground state [27,28]. Even though Ba 3 CoSb 2 O 9 possesses long-range magnetic order (with 120°ordering), several aspects of the magnon dispersion and the multimagnon continuum reveal an unconventional behavior, which can only be partially explained within semiclassical approaches. First of all, at low energies, the magnon dispersion is strongly renormalized with respect to the linear spin-wave approximation; an anomalous line broadening has also been detected, leading to the conclusion that magnon decay may be plausible; finally, the continuum  6. In both lower panels the orange shaded area corresponds to the region of the Brillouin zone in which magnon decay is predicted by the spinwave approximation [16,17], and the dashed line delimits the magnetic Brillouin zone. presents unexpected dispersive features at high energies. It should be noticed that, since neutron scattering data are sensitive to the full dynamical spin structure factor, three copies of the magnon dispersion (translated by the ordering vectors) are visible in the spectrum. Experimental investigations have also been performed to infer the nature of the magnon excitations on top of the gapped phase that is stabilized at the one-third magnetization plateau [29]. In this case, the situation seems to be more conventional, with the experimental results in relatively good agreement with theoretical predictions.
Motivated by these experimental findings, there have been a few attempts to investigate the Heisenberg model (also including small perturbations) with both analytical and numerical tools [30][31][32][33]. In particular, by using density-matrix renormalization group (DMRG) calculations, Verresen et al. [32] claimed that the magnon decay does not take place, because of the strong coupling interactions between quasiparticles (i.e., magnons) in the Heisenberg model [34]. As a result of the avoided decay, the midpoint of the edge of the magnetic Brillouin zone (dubbed Y 1 ) displays a minimum of the magnon dispersion, possibly explaining the high-energy features seen around the M point in Ref. [28].
Within this context, the discovery of YbMgGaO 4 [35] and, more recently, NaYbO 2 [36] will also give a further impetus to study (generalized) spin models on the triangular lattice. In both cases, no signatures of magnetic order appear down to very low temperatures, suggesting the existence of a quantum spin liquid. While both materials host effective J ¼ 1=2 spin degrees of freedom (d.o.f.), the actual low-energy Hamiltonian may be more complicated than the SU(2)-invariant one of Eq. (1); still, the physical properties can share many similarities with the ground state of the J 1 − J 2 model, as suggested in Ref. [7].
In this work, we employ a dynamical variational Monte Carlo approach [37] to compute the out-of-plane dynamical spin structure factor for the Heisenberg model on the triangular lattice, also in the presence of a nextnearest-neighbor coupling J 2 . First, we focus our attention on the model with J 2 ¼ 0 for which we confirm huge corrections from the linear spin-wave calculations. Our results support the idea that the magnon excitations are stable in the whole Brillouin zone; indeed, even though a discrete set of excitations is obtained within our numerical method, the lowest-energy state for each momentum q appears to be rather well separated from the rest of the spectrum at higher energies, suggesting the existence of a faint continuum just above the magnon branch. The second part of this work deals with the J 1 − J 2 model, to highlight the modifications in the dynamical structure factor that take place when entering the spin-liquid phase (which, according to our variational approach, is gapless [13]). Here, the spectrum shows gapless excitations at M points; in addition, a strong signal at low energies is present at the corners of the Brillouin zone, i.e., K ¼ ð2π=3; 2π= ffiffi ffi 3 p Þ and K 0 ¼ ð4π=3; 0Þ. While the former aspect can be easily understood by inspecting the noninteracting spinon band structure, the latter one is a genuine feature that emerges from the Gutzwiller projector, which includes interactions between spinons and gauge fields. Indeed, while the noninteracting wave function corresponds to a mean-field approximation, in which gauge fields are completely frozen, the Gutzwiller projection has the effect of inserting back the temporal fluctuations of those fields [38]. In this respect, it is worth mentioning that a recent field-theoretical analysis indicated the existence of low-energy (triplet) monopole excitations at the zone corners, which are expected to contribute to the dynamical structure factor [39].

II. DYNAMICAL VARIATIONAL MONTE CARLO METHOD
The dynamical structure factor, which is directly measured within inelastic neutron scattering experiments, can be used to unveil the nature of the elementary excitations of the models and materials under investigation. In its spectral form, this quantity reads as where jϒ 0 i and fjϒ q α ig α are the ground state and the set of all excited states with momentum q, whose corresponding energies are E 0 and fE q α g α , respectively. In this work, we evaluate the dynamical structure factor of the spin model Eq. (1) by directly constructing accurate variational Ansätze for its ground state and a few low-energy excited states. Our variational approach is based on the so-called parton construction, in which the spin d.o.f. of the model are rewritten in terms of auxiliary fermionic operators [14,40]. The fermionic language constitutes a versatile framework to define variational wave functions for both magnetically ordered and disordered phases of matter. This section is dedicated to the introduction of the fermionic wave functions for spin models and to the description of the variational Monte Carlo method employed for the calculation of the dynamical structure factor.
A. Gutzwiller-projected fermionic wave functions for the ground state Here, for the sake of generality, we consider a generic SU(2) model for frustrated spin systems, which consists of a set of spin-1=2 d.o.f. sitting on the sites of a lattice and interacting through the Heisenberg exchange couplings J i;j : The interplay of the different interactions can lead to the stabilization of different phases of matter. In the absence of frustration, i.e., when no competing couplings are present, the ground state may develop some kind of magnetic order, which minimizes the classical energy of the model. On the contrary, when different interactions compete with each other, magnetically disordered phases can arise, such as spin liquids.
The first attempt to describe spin-liquid states dates back to the resonating valence-bond approach, where a variational wave function is defined in terms of a linear superposition of singlet coverings of the lattice [1]. More recently, Wen [40] developed a general approach to classify and construct spin-liquid states, which satisfy all the symmetries of a given lattice model. This method is built upon the introduction of auxiliary Abrikosov fermions, which form a projective representation of S ¼ 1=2 spin operators: Here, c i;α (c † i;α ) destroys (creates) a fermion with spin α ¼ ↑; ↓ on site i, and the vector σ ¼ ðσ x ; σ y ; σ z Þ is the set of Pauli matrices. The anticommutation relations among fermions ensure that the Abrikosov representation yields the correct commutation relations among different spin components. Still, in order to faithfully reproduce the Hilbert space of the original spin model, only configurations with one fermion per site must be considered, which implies that the Abrikosov fermions must satisfy the constraint or, equivalently, Besides constant terms, the Hamiltonian of Eq. (3) can be rewritten in terms of Abrikosov fermions as follows: At this stage, the Hamiltonian Eq. (7) with the constraints of Eqs. (5) and (6) gives an exact representation of the original model. In order to tackle the above interacting fermionic system, one possibility is to perform a mean-field decoupling [40]. For the purpose of studying spin-liquid phases, we keep only the mean-field terms that do not break the SU(2) symmetry of the original spins. The result is a quadratic Hamiltonian: which contains a hopping term t i;j and a singlet pairing term Δ i;j , which are related to the expectation values hc † j;σ c i;σ i and hc i;σ c j;−σ i, respectively. In addition, the one-fermion-per-site constraint of the parton construction is enforced in a global fashion by including a chemical potential μ i and an on-site pairing ζ i as Lagrange multipliers in H 0 [40]. Within the mere mean-field approach, the parameters of H 0 are computed self-consistently and define a low-energy effective theory for the spin model under investigation. However, the ground state of H 0 , named jΦ 0 i, satisfies the constraints of Eqs. (5) and (6) only on average and, therefore, does not represent a valid wave function for spins. Within this approach, a full treatment of the original spin model requires the inclusion of all fluctuations of the parameters around the mean-field solution. Since this task is in general unfeasible, an alternative approach can be pursued, in which the Hamiltonian H 0 is exploited as a starting point for the definition of a variational wave function for the initial spin model. Indeed, the one-fermion-per-site constraint can be enforced exactly by applying the Gutzwiller projector, to the ground-state wave function of H 0 . We emphasize that in general the Gutzwiller projection cannot be treated analytically, due to its intrinsic many-body character; however, it can be considered within Monte Carlo sampling. At variance with the mean-field treatment, in the variational approach the parameters of H 0 are not computed selfconsistently, but are optimized in order to minimize the energy of the Gutzwiller-projected Ansatz P G jΦ 0 i. The artificial enlargement of the Hilbert space introduced by the parton construction gives rise to a gauge redundancy in the representation of the spin d.o.f. Specifically, the mapping Eq. (4) is invariant under local SU(2) transformations of the Abrikosov fermions operators [40]. As a consequence, all physical properties of the spins are independent of the gauge choice for fermions. For example, whenever we perform SU(2) transformations to the unprojected Hamiltonian H 0 , the variational wave function with the Gutzwiller projector remains invariant. Exploiting this gauge redundancy, it is possible to classify all the quadratic Hamiltonians H 0 whose Gutzwiller-projected ground states fulfill the symmetries of the lattice model. This procedure, known as projective symmetry group analysis [40], provides a recipe to construct all the distinct spin-liquid Ansätze for a given spin model. From a variational point of view, the spin-liquid wave function with the lowest variational energy is the one which better describes the true ground state of the model.
In general, the variational Ansätze defined by Gutzwiller projecting the ground state of Eq. (8) do not display any magnetic order [41]. For the purpose of defining suitable wave functions for magnetically ordered phases, an additional term can be added to H 0 : Here, h is a fictitious magnetic field which lies in the XY plane and displays a periodic pattern defined by the pitch vector Q. Since the ground-state wave function of the Hamiltonian Eq. (10) tends to overestimate the magnetic order [42], further transverse quantum fluctuations are added through the application of a spin-spin Jastrow factor, to the Gutzwiller-projected state. Specifically, the complete form of the variational wave functions employed in this work is where, in addition to the Gutzwiller projection and the Jastrow factor, we apply a projector enforcing zero value for the z component of the total spin (P S z ). By using this approach, the variational phase diagram for the J 1 − J 2 model on the triangular lattice has been obtained in Ref. [13]: the system undergoes a phase transition between a magnetically ordered phase to a gapless spin liquid at J 2 =J 1 ≈ 0.08. For this model, the optimal variational wave functions are obtained by considering only a hopping term (no pairing) and the fictitious magnetic field in the quadratic Hamiltonian: Here, t is a first-neighbor hopping with a nontrivial sign structure (s i;j ¼ AE1) which generates a pattern of alternating 0 and π fluxes through the triangular plaquettes of the lattice, see Fig. 1, and h is a fictitious magnetic field which displays the classical 120°order with Q ¼ ð2π=3; 2π= ffiffi ffi 3 p Þ, see Fig. 1 [considering Q ¼ ð4π=3; 0Þ would not change the physical content of the ground state wave function]. All the parameters included in H 0 and the pseudopotential v i;j (one parameter for each distance jR i − R j j in the translational invariant lattice) entering the Jastrow factor can be optimized to minimize the variational energy. While in the magnetic phase of the system the optimal value for the ratio h=t is finite, for J 2 =J 1 ≳ 0.08 the system enters the spin-liquid phase and the magnetic field parameter vanishes in the thermodynamic limit [13]. The values of the fictitious magnetic field as a function of J 2 =J 1 can be found in Ref. [13].
In this work, we compute the dynamical structure factor for the J 1 − J 2 model on the 30 × 30 triangular lattice. For J 2 ¼ 0, we first consider the crudest approximation for the ground state, which consists in setting the hopping term t to zero. The resulting wave function is equivalent to the state of Ref. [43] with only a two-body Jastrow factor. Much more accurate results are then obtained by restoring the hopping term in the Hamiltonian and optimizing all the variational parameters for the cases J 2 ¼ 0 and J 2 =J 1 ¼ 0.07. On the other hand, when the system is in the spin-liquid regime (J 2 =J 1 ¼ 0.09 and J 2 =J 1 ¼ 0.125), the fictitious magnetic field is vanishing and the Jastrow factor is not considered, because of its negligible effects on the variational results. According to the projective symmetry group classification, the wave function obtained by considering only the hopping term in H 0 is a fully symmetric U(1) spin liquid [44].

B. Dynamical structure factor
As already mentioned, the dynamical structure factor of the J 1 − J 2 model is computed by constructing variational Ansätze to approximate the low-energy excited states of the system. Here we limit ourselves to the calculation of the out-of-plane component S z ðq; ωÞ, and we employ the technique outlined in Refs. [37,45,46], which is briefly summarized in the following.
First, we find the optimal variational Ansatz for the ground state of the model, which has the form of Eq. (12), by minimizing the variational energy. The resulting wave function is employed as a reference state to construct a set of projected particle-hole excitations with a given momentum q: These states are labeled by R, which runs over all lattice vectors. We approximate the low-energy excited states of the model by using linear combinations of the elements of the basis set fjq; Rig R : For a certain momentum q, we consider the Schrödinger equation for the J 1 − J 2 Hamiltonian restricting the form of its eigenvectors to the one of Eq. (15), i.e., HjΨ q n i ¼ E q n jΨ q n i. Expanding everything in terms of fjq; Rig R , DYNAMICAL STRUCTURE FACTOR OF THE J 1 − J 2 HEISENBERG … PHYS. REV. X 9, 031026 (2019) 031026-5 we arrive at the following generalized eigenvalue problem, X R 0 hq; RjHjq; R 0 iA n;q R 0 ¼ E q n X R 0 hq; Rjq; R 0 iA n;q R 0 ; ð16Þ which is solved to find the expansion coefficients A n;q R and the energies E q n of the excitations. All the matrix elements, hq; RjHjq; R 0 i and hq; Rjq; R 0 i, are evaluated within the Monte Carlo procedure, by sampling according to the variational ground-state wave function. Finally, the dynamical structure factor is computed by where E var 0 is the variational energy of jΨ 0 i.

III. RESULTS
In this section, we present the numerical calculations for the dynamical structure factor Sðq; ωÞ obtained by the variational approach described in the previous section. First, we discuss the case of the Heisenberg model with only nearest-neighbor superexchange J 1 , also comparing our results with recent DMRG calculations [32]. Then, we include the next-nearest-neighbor coupling J 2 to increase frustration and melt the magnetic order. In this way, a gapless spin-liquid regime is reached for J 2 =J 1 ≈ 0.08 [13].
A. Nearest-neighbor model with J 2 = 0 Let us start our analysis by considering the case in which the ground-state wave function contains only the fictitious magnetic field, i.e., t ¼ 0. In this case, the Abrikosov fermions are completely localized (e.g., the eigenvalues of the auxiliary Hamiltonian define flatbands) and the wave function corresponds to the Jastrow state of Ref. [43] with only a two-body Jastrow factor. The results for the dynamical structure factor on the 30 × 30 cluster are shown in Fig. 2. Here, the spectrum consists of a single mode, which is identified as the magnon excitation (no continuum is visible). Notice that only one magnon branch is visible, related to the magnon dispersion ϵ q , since we consider the out-of-plane dynamical structure factor (the folded branches ϵ qAEK do not contribute to the signal). Remarkably, the dispersion of the magnon branch is possible thanks to the Jastrow factor, since the wave function without it would give rise to a trivially flat (gapped) excitation spectrum, reflecting the noninteracting band structure of fermions. By contrast, the long-range Jastrow term is able to produce a reasonable magnon mode, which agrees fairly well with the spinwave calculations. In particular, the spectrum is gapless at Γ ¼ ð0; 0Þ (with a vanishingly small weight). Instead, in contrast to spin waves, which correctly predict gapless magnons at K and K 0 due to the coplanar 120°order, this simple wave function leads to a gapped spectrum at the corners of the Brillouin zone. In connection to that, the outof-plane static structure factor S z ðqÞ ¼ R dωS z ðq; ωÞ does not diverge at K or K 0 when L → ∞, showing only a maximum.
A much more realistic spectrum is obtained when considering a finite fermion hopping t (with the π-flux pattern shown in Fig. 1), as well as the optimized value of the fictitious magnetic field h (and the Jastrow factor). The results for the 30 × 30 lattice are reported in Fig. 3. In this case, there are several excitations with a finite weight for each momentum, thus reproducing the existence of a broad continuum, which extends up to relatively large energies. We mention that, with respect to the square lattice [46][47][48], here many more excitations for each momentum possess a visible spectral weight. Within this calculation, we identify the lowest-energy excitation E q 0 as the magnon peak. This assumption is corroborated by the results shown in Fig. 4, where the variational energies E q 0 closely follow the magnon branch obtained by series expansions. Instead, identifying the lowest-energy peak as the bottom of the continuum is not very plausible, since a much broader signal should be present in this case. In this regard, the basis set that is used here for the excited states is made of particle-hole spinon excitations on top of the ground state of the auxiliary Hamiltonian of Eq. (13), before Gutzwiller projection. For this reason, we argue that, in general, our approach is particularly suited to capture (i) two-spinon excitations or (ii) bound states of spinons, e.g., magnons. Multimagnon excitations are expected to show up with a reduced intensity. In order to discuss the issue of magnon decay, we apply a kinematic argument (as done both in the FIG. 2. Dynamical structure factor of the nearest-neighbor Heisenberg model on the triangular lattice obtained by using the variational wave function of Eqs. (12) and (13) with t ¼ 0 on the 30 × 30 cluster. The path along the Brillouin zone is shown in Fig. 1. A Gaussian broadening of the spectrum has been applied (σ ¼ 0.02J 1 ). The spin-wave energies of the magnon branch (ϵ q ), on the same cluster size, are represented by the white dots connected with a solid line. The dashed line corresponds to the bottom of the continuum within linear spin waves, i.e., E q ¼ min k fϵ q−k þ ϵ k g. Notice that E q < ϵ q in most of the Brillouin zone, as obtained in Refs. [16,17].
linear spin-wave approach [16,17] and within DMRG [32]) and we consider all the possible two-magnons decays, which fulfill the conservation of momenta, i.e., E q ¼ min k fE q−k 0 þ E k 0 g. For this purpose, we computed the spectrum E k 0 for all the k vectors in the Brillouin zone on the 30 × 30 lattice. The outcome is that the bottom of the two-magnon continuum, defined by the kinematic analysis, lies above the magnon branch. These results clearly indicate an avoided decay in a large part of the Brilloiun zone, as suggested by DMRG calculations, which considered certain (high-energy) parts of the magnon dispersion [32]. Still, we cannot exclude the existence of small regions where the magnon decay may persist, especially close to the gapless points. In this respect, within the linear spin-wave approach, the different velocities of the excitation spectrum at Γ and K immediately lead to an unstable magnon branch close to the Γ point [16,17]. Should this aspect be a genuine feature of the model, the magnon would be unstable in a small part around the center of the Brillouin zone. Unfortunately, given the finiteness of the cluster used in our numerical calculations, we cannot reliably estimate the slope of the magnon spectrum at Γ and K and, therefore, make definitive statements for this issue.
Here, we would like to notice the strong renormalization of the magnon branch with respect to spin-wave calculations; see Fig. 4. Most importantly, we emphasize that, within this most accurate calculation, the magnon branch shows a rotonlike minimum not only at M but also at Y 1 , i.e., the midpoint of the edge of the magnetic Brillouin zone (see also Fig. 5), as already detected by neutron scattering measurements in Ba 3 CoSb 2 O 9 [28]. This feature was not captured by the previous series expansion calculations [21] but, instead, has also been observed by recent DMRG calculations on an infinitely long cylinder (with a small circumference L ¼ 6) [32] and has been interpreted as the hallmark for the absence of magnon decay. In order to make a closer comparison with DMRG data, we perform the variational calculations on a long cylinder (84 × 6) along the same path in the Brillouin zone as the one that has been considered in Ref. [32]. The results are shown in Fig. 6. Here, the large number of lattice points along the cylinder FIG. 3. The same as Fig. 2 but for the optimal variational wave function with both hopping t and fictitious magnetic field h. The path along the Brillouin zone is shown in Fig. 1. The dotted line denotes the bottom of the continuum E q ¼ min k fE q−k 0 þ E k 0 g, where E q 0 is the lowest energy for a given momentum q obtained within our variational approach. Since the spectrum is gapless at the Γ point, we exclude the cases k ¼ ð0; 0Þ and k ¼ q in the search of the minimum, because the resulting E q would simply coincide with the energy of the magnon branch E q 0 all over the Brillouin zone. The purpose of this kinematic analysis is to show that no magnon decay can yield an energy E q which is lower than the one of the magnon branch E q 0 (in contrast to spinwave results).

FIG. 4. Energies of the magnon branch for the nearest-neighbor
Heisenberg model on the triangular lattice obtained with different methods. The path along the Brillouin zone is shown in Fig. 1. The black line corresponds to linear spin wave, the blue squares to series expansion [21], and the orange circles to our variational results (on the 30 × 30 cluster). allows us to have a detailed resolution of the magnon branch, which closely follows the one obtained by DMRG. In particular, we can estimate the bottom of the continuum by evaluating where we consider the possible decays involving a magnon at K and −K. In doing so, we find that the lowest-energy excitation E q 0 is always below E q , indicating that a welldefined branch exists and magnon decay is avoided. We finally remark that a roton minimum is detected along the same path as the one studied by Verresen et al. [32], strongly suggesting that this is a genuine feature of the Heisenberg model.
We now move to the case where also a next-nearestneighbor coupling J 2 is present. Within our variational approach, a gapless spin-liquid phase is stabilized for 0.08 ≲ J 2 =J 1 ≲ 0.16; here, the fictitious magnetic field vanishes in the thermodynamic limit and the best wave function contains only fermionic hopping (with π flux threading half of the triangular plaquettes) [13]. On a finite size, a small value of h can be stabilized, as well as a tiny Jastrow pseudopotential. Still, we verified that these ingredients do not cause significant differences in the dynamical structure factor. In Fig. 7, we show the results for the 30 × 30 cluster and for two values of J 2 =J 1 , which are very close to the transition point, one still inside the magnetic phase (J 2 =J 1 ¼ 0.07) and the other one in the spin-liquid region (J 2 =J 1 ¼ 0.09). By approaching the quantum phase transition, the major modification of the spectrum comes from the softening of the magnon excitation at the M points. This feature closely resembles the case of the frustrated J 1 − J 2 model on the square lattice, previously studied with the same numerical technique [46], where a softening is clearly detected for q ¼ ðπ; 0Þ [and ð0; πÞ]. In this latter case, this fact has been connected to the progressive deconfinement of spinons that have gapless (Dirac) points at q ¼ ðAEπ=2; AEπ=2Þ. We would like to mention that the possibility to have (gapped) almostdeconfined spinons in the unfrustrated Heisenberg model has been suggested by a recent quantum Monte Carlo calculation [49]; moreover, clear signatures for deconfined spinons at the transition between an antiferromagnetically ordered phase and a valence-bond crystal have been reported in the so-called J − Q model [50]. On the triangular lattice, the softening of the spectrum at the M points is a direct consequence of the Dirac points at q ¼ ð0; AEπ= ffiffi ffi 3 p Þ in the spinon band structure. Therefore, we expect both M and K points to be gapless at the transition (as well as Y 1 , which can be obtained by combining M and K vectors). Indeed, this is necessary for a continuous phase transition, as the one that appears in the J 1 − J 2 Heisenberg model, according to ground-state calculations [13].
FIG. 6. The dynamical structure factor for the nearest-neighbor Heisenberg model on a cylindrical geometry (84 × 6), to make a close comparison with DMRG calculations by Verresen et al. [32]. We apply a Gaussian broadening to the spectrum, which is equivalent to the one of the aforementioned DMRG result (σ ¼ 0.077J 1 ). The path in the Brillouin zone is shown in the inset and in Fig. 1 0 is the lowest energy for a given momentum q obtained within our variational approach and K ¼ ð2π=3; 2π= ffiffi ffi 3 p Þ.
FIG. 7. The dynamical structure factor for the J 1 − J 2 Heisenberg model on the 30 × 30 cluster with J 2 =J 1 ¼ 0.07 (top) and J 2 =J 1 ¼ 0.09 (bottom). The path along the Brillouin zone is shown in Fig. 1, and a Gaussian broadening of the spectrum has been applied (σ ¼ 0.02J 1 ).
In Fig. 8, we report the dynamical structure factor for J 2 =J 1 ¼ 0.125. The spin-liquid state is characterized by a broad continuum that extends up to relatively large energies. In particular, around the M points, the magnon rotonlike minima of the ordered phase fractionalize into an incoherent set of excitations at low energies. This feature is compatible with the existence of Dirac points in the unprojected spectrum of the auxiliary Hamiltonian H 0 ; see Fig. 8. By contrast, a strong signal in the lowestenergy part of the spectrum is detected around the K points, where the unprojected spinon spectrum is instead gapped. In this respect, the Gutzwiller projection is fundamental to include interaction among spinons in a nonperturbative way and give a drastic modification of the low-energy features. This is a distinctive aspect of the triangular lattice, since, on the square lattice, all the low-energy (gapless) points observed in the presence of the Gutzwiller projector [i.e., q ¼ ð0; 0Þ, ðπ; πÞ, ðπ; 0Þ, and ð0; πÞ] already exist in the noninteracting picture [51]; see Fig. 9. We would like to emphasize that, in contrast to the magnetically ordered phase, where no visible spectral weight is present right above the magnon branch (see Fig. 3), in the spin-liquid phase the continuum is not separated from the lowest-energy excitation. This outcome corroborates the fact of having deconfined spinons in the magnetically disordered phase. The intense signal at K points immediately implies strong (but short-range) antiferromagnetic correlations in the variational wave function, which are absent in the unprojected π-flux state (by contrast, on the square lattice, the π-flux state already has significant antiferromagnetic correlations built in it).  Fig. 1. We applied a Gaussian broadening of σ ¼ 0.02J 1 to the variational results. Notice that, for the unprojected data, the energy scale is given by the hopping amplitude t of the unprojected Hamiltonian Eq. (13), instead of J 1 . In addition, the broadening has been rescaled in order to account for the larger bandwidth of the spectrum.  [46] for details). We applied a Gaussian broadening of σ ¼ 0.02J 1 to the variational results. Notice that, for the unprojected data, the energy scale is given by the hopping amplitude t of the unprojected Hamiltonian of Ref. [46], instead of J 1 . In addition, the broadening has been rescaled in order to account for the larger bandwidth of the spectrum. The presence of low-energy spectral weight at the corners of the Brillouin zone could be ascribed to the existence of critical monopole excitations, as suggested by the analysis of Ref. [39]. In fact, the Gutzwiller projector, which imposes single occupancy on each lattice site, introduces temporal fluctuations of the gauge fields that are completely frozen within the noninteracting picture (i.e., within the unprojected wave function). Even though we cannot exclude a more conventional picture where a bound state of spinons is responsible for the intense signal around K, it is plausible that this feature originates from the existence of gauge fields, which emerge in the fieldtheoretical description of spin liquids [14]. While gauge fields are known to predominantly contribute to spectral functions of specific Kitaev spin liquids with Z 2 magnetic fluxes [52], our calculations suggest that monopole excitations may give some relevant signature in the spin-liquid phase of the J 1 − J 2 Heisenberg model on the triangular lattice. Remarkably, on the 30 × 30 cluster, the lowestenergy excitation at K is slightly higher inside the spinliquid phase (i.e., for J 2 =J 1 ¼ 0.125) than close to the critical point (i.e., for J 2 =J 1 ≈ 0.08); see Figs. 7 and 8. This fact may suggest the possibility that this kind of excitation may be slightly gapped in the spin-liquid region, while being gapless at the critical point. We finally highlight the existence of an unexpected high-energy dispersing mode, which bends from the Γ point down into the continuum, being seemingly connected to the low-energy excitation at K. A comparison with other numerical techniques will be needed to clarify whether this feature is a genuine aspect of the model or an artifact of the present variational approach.

IV. CONCLUSIONS
In this work, we performed variational Monte Carlo calculations to estimate the dynamical structure factor of the J 1 − J 2 Heisenberg model on the triangular lattice. The results for J 2 ¼ 0 are consistent with the existence of a well-defined magnon branch in the whole Brillouin zone, in agreement with recent DMRG calculations [32]. This outcome contrasts the semiclassical predictions [16,17], which suggested the presence of magnon decay in a large portion of the Brillouin zone. When a finite J 2 superexchange is included and the spin-liquid phase is approached, a clear softening of the spectrum is detected around the M points, in close similarity to what happens on the square lattice [46]. Remarkably, the low-energy physics of the spin-liquid phase cannot be fully described by the unprojected spinon picture, since, besides gapless excitations at M and M 0 , there are anomalously low-energy states appearing around the K points. Our numerical calculations provide indisputable evidence of the fact that the noninteracting (i.e., unprojected) spinon spectrum is not sufficient to fully explain the low-energy spectrum detected by the dynamical structure factor. In light of the recent field-theoretical analysis [39], the natural interpretation of the spectral features around the corners of the Brillouin zone comes from the existence of low-energy monopole excitations. This outcome is particularly important, since it would give a direct signature of the fact that these theoretical approaches correctly capture the nature of the spin-liquid phase. We hope that the present results will motivate future investigations in this direction.