Spin Supersolid Phase and Double Magnon-Roton Excitations in a Cobalt-based Triangular Lattice
Yuan Gao Chuandi Zhang Junsen Xiang Dehong Yu Xingye Lu Peijie Sun Wentao Jin Gang Su Wei Li




These authors contributed equally to this work.
School of Physics, Beihang University, Beijing 100191, China
Institute of Theoretical Physics, Chinese Academy of Sciences, 
Beijing 100190, China
		

These authors contributed equally to this work.
School of Physics, Beihang University, Beijing 100191, China


Beijing National Laboratory for Condensed Matter Physics, 
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China


Australian Nuclear Science and Technology Organisation, 
Lucas Heights NSW 2234, Australia


Center for Advanced Quantum Studies and Department of Physics,
Beijing Normal University, Beijing 100875, China


Beijing National Laboratory for Condensed Matter Physics, 
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China


wtjin@buaa.edu.cn
School of Physics, Beihang University, Beijing 100191, China


gsu@ucas.ac.cn
Kavli Institute for Theoretical Sciences, University of Chinese 
Academy of Sciences, Beijing 100190, China


w.li@itp.ac.cn
Institute of Theoretical Physics, Chinese Academy of Sciences, 
Beijing 100190, China
Peng Huanwu Collaborative Center for Research and Education,
Beihang University, Beijing 100191, China




Supersolid is an exotic quantum 
state of matter that hosts spontaneously the features of both solid and superfluidity,
which breaks the lattice translational symmetry and U(1) gauge symmetry. Here we conduct inelastic 
neutron scattering (INS) measurements and tensor-network calculations on the 
triangular-lattice cobaltate , which is proposed in [Xiang et al., 
Nature 625, 270-275 (2024)] as a quantum magnetic analog of supersolid.
We uncover characteristic dynamical signatures, which include distinct magnetic 
Bragg peaks indicating out-of-plane spin solidity and gapless Goldstone modes 
corresponding to the in-plane spin superfluidity, offering comprehensive spectroscopic 
evidence for spin supersolid in . We also compute spin dynamics of the 
easy-axis triangular-lattice model, and reveal magnon-roton excitations containing 
U(1) Goldstone and roton modes associated with the in-plane spin superfluidity, 
as well as pseudo-Goldstone and roton modes related to the out-of-plane spin 
solidity, rendering double magnon-roton dispersions in the spin supersolid. Akin 
to the role of phonon-roton dispersion in shaping the helium thermodynamics, 
the intriguing magnetic excitations also strongly influence the low-temperature 
thermodynamics of spin supersolid down to sub-Kelvin regime, explaining 
the recently observed giant magnetocaloric effect in .


Spin Supersolid Phase and Double Magnon-Roton Excitations in a Cobalt-based Triangular Lattice
    Wei Li
    April 29, 2024
==============================================================================================




Introduction.— 
As a paradigmatic example of frustrated quantum magnet, 
the triangular-lattice antiferromagnets (TLAF) have garnered significant 
attention in the past <cit.>. 
The perfect isosceles triangular-lattice compounds, including the cobaltate 
Ba_3CoSb_2O_9 <cit.>, rare-earth triangular 
compounds such as REMgGaO_4 (with RE = rare earth)
<cit.> and structurally similar compounds 
ARECh_2 with A = Na, K, Cs, and Ch = O, S, Se <cit.>, have been synthesized recently. There are 
investigations into possible quantum spin liquid (QSL) <cit.>, fractional magnetization plateau <cit.>, 
anomalous low-temperature thermodynamics <cit.>, 
topological phase transitions <cit.>, and exotic 
spin excitations <cit.>, 
etc, making the TLAF systems a very intriguing, fertile ground for studying 
emergent quantum phenomena. 





Recently, an easy-axis Co-based triangular-lattice antiferromagnet 
 (NBCP) has raised great research interests <cit.>. 
The Co^2+ ions form stacked triangular lattices [c.f., Fig. <ref>(a)], 
which carry effective spin S=1/2 under the effects of spin-orbit coupling 
and crystal electric field. The spin-spin couplings are highly two-dimensional, 
i.e., the intra-layer spin exchanges are dominating over those between 
the layers <cit.>. 
The highly frustrated quantum magnet was first proposed to host a quantum 
spin liquid, where magnetic ordering was absent down to 300 mK
<cit.>. A later specific heat measurement finds an 
anomaly at about T_ N≃150 mK <cit.>, which 
instead suggests the formation of certain spin order at low temperature. 
The existence of residual thermal conductivity has also been reported in 
this compound <cit.>, although a different conclusion was drawn 
based on an independent measurement <cit.>. 

An easy-axis TLAF model has been put forward to elucidate the diverse experimental findings in NBCP. The model Hamiltonian reads H = 
∑_⟨ i,j⟩ J_xy (S_i^xS_j^x+S_i^yS_j^y) + J_z S_i^zS_j^z, 
where the nearest-neighbor interactions between i, j are J_xy=0.88
 K and J_z=1.48  K <cit.>. 
Based on tensor-network calculations on this highly frustrated spin model, 
it has been proposed theoretically <cit.> and observed experimentally the long-sought spin supersolidity NBCP, through the
magnetocaloric and neutron diffraction measurements <cit.>. 
Nevertheless, dynamical features of the spin supersolid phase remain 
to be unexplored, both in theory and experiment.

Here we conduct INS measurements on single-crystal samples of NBCP 
down to 55 mK. Complementing these investigations, tensor-network 
computations are performed using the TLAF model, yielding results that are 
in excellent agreement with the INS data. This synergy between experimental 
and theoretical modeling provides a comprehensive understanding of the 
magnetic properties and dynamics of NBCP at low temperatures. In particular,
a substantial downward renormalization is revealed in the magnon dispersion 
near the M point of Brillouin zone (BZ), constituting the magnetic analog of 
roton mode. Meanwhile, a pseudo-Goldstone gap is obtained by computing 
the spin-resolved spectral function. Overall, the peculiar double magnon-roton 
dispersions, one from in-plane spin superfluidity and the other from out-of-plane 
spin solidity, render strong spin fluctuations and large magnetic entropies till low 
temperature and well explain the giant magnetocaloric effect (MCE) observed 
in NBCP.


Samples and neutron scattering measurements.—
Single-crystal samples of  were grown using the flux method as reported 
in Ref. Xiang2024Nature. The INS experiments were performed on 
the cold-neutron time-of-flight spectrometer PELICAN <cit.>, at Australian 
Nuclear Science and Technology Organisation (ANSTO). A total number of 28 
pieces of NBCP single crystals with a total mass of about 3 g were mounted 
on the sample holder made of 6061 aluminum alloy using CYTOP M-type glue. 
The base temperature of 55 mK was achieved using a dilution insert inside a 
7 T vertical cryomagnet. The crystals were co-aligned with their [1, -1, 0] 
direction lying vertically, so that the [H, H, L] scattering plane can be 
mapped out by rotating the sample and an in-plane magnetic field can be 
applied along the [1, -1, 0] direction.

The instrument was configured with an incident neutron wavelength of 5.96 Å, 
providing an incident energy of 2.3 meV with a high energy resolution of about 
0.066 meV at the elastic line. The dataset below 0.1 meV collected at 2.5 T 
was used for background subtraction for the 0 T case, and a standard vanadium 
sample was measured for detector normalization and determination of the energy 
resolution function. The data reductions were performed using the software HORACE
<cit.>.






Spectroscopic evidence for spin supersolid.—

In Fig. <ref>(b) we present the elastic scattering (ω=0±0.05 meV) 
results of the INS measurements, which clearly reveal magnetic Bragg peaks at 
K≡ [1/3, 1/3] points with an incommensurate out-of-plane propagation 
vector of ζ = 0.16, due to sensitivity of the spin supersolid state to weak 
interlayer couplings. Further analysis shows that the magnetic Bragg peaks are 
contributed mainly from the the out-of-plane moments <cit.>, 
consistent with prior neutron diffractions results <cit.>. Besides the magnetic Bragg peaks, the low-energy spin 
fluctuations observed from the INS measurements exhibit a rod-like shape
<cit.>, showing very good two dimensionality of the compound.





In Figs. <ref>(a,b) we present the 
low-energy magnetic excitations observed at 55  mK with zero field and
an in-plane field of B=2.5  T, respectively, and compare them to the model
calculations. At zero field, evidence of gapless Goldstone modes are shown in Fig.
<ref>(a), where the linear spin-wave dispersions emanating from the ordering 
vector K are also plotted. The tensor-network calculations with realistic
model parameters well reproduce the experimental results with similar energy 
resolution. 
On the other hand, with a larger magnetic field B=2.5  T above the in-plane 
critical field <cit.>, a clear magnon dispersion in the nearly 
polarized phase is observed. The theoretical calculations demonstrate an excellent 
match with the experimental results for both zero and 2.5 T field cases, confirming 
once again the validity and accuracy of the effective easy-axis TLAF model for NBCP.

Figure <ref>(a) reveals the presence of an extra intensity that overlays the 
standard linear spin-wave dispersion, notably concentrated around the M≡ [1/2, 1/2] high symmetry point of the Brillouin zone. A downward renormalization of magnon dispersion can be discerned in experimental data in 
Fig. <ref>(a) and simulated results in Fig. <ref>(c), which become more 
evident by improving the energy resolution (see in Fig. <ref>). Below, we 
show such soft mode represents magnetic analog of roton excitation in superfluid 
helium <cit.>. 
 

Dynamical calculations of easy-axis TLAF model.—

Here we consider the easy-axis TLAF model with realistic parameters of
the spin supersolid compound  and compute the dynamical spin structure 
factor S(k, ω) and spectral function A(k, ω). These two quantities
can be obtained from real-time correlation function g^αβ(k,t)= 
1/N e^ i E_0 t∑_i,j e^- i k·(r_i-r_j)⟨ψ_0| 
S^α_i e^- i H t S^β_j |ψ_0⟩, where |ψ_0⟩ (E_0) 
is the ground-state wavefunction (energy), and N is the total site number. The 
ground state |ψ_0⟩ can be obtained with density matrix renormalization 
group <cit.> and the real-time evolution e^- i H t S^β_j |ψ_0⟩ is computed
with time-dependent variational principle approach <cit.>. 

Given the correlation function g^αβ(k,t), the spin-resolved dynamical
structure factor can be computed as S^αβ(k,ω) = ∫_0^t_ max Re [g^αβ(k,t) e^iω t] W(t/t_ max) dt and spectral 
function A^αβ(k,ω)  =  -∫_0^t_ max Im 
[g^αβ(k,t)] sin(ω t) W(t/t_ max) dt, where the energy 
resolution ε is controlled by the maximal evolution time t_ max, 
as ε≃ 8/t_ max with the Parzen window function W(t/t_ max)
<cit.>. In practice, the dynamical calcualtions are performed on 
a YC6 × 15 lattice with the simulated time up to t_max = 100/J_z. 
For zero-field case, the retained bond dimension is D=2000 for all the contour 
plot with momentum scan, and D=3000 for the lines with a fixed momentum [e.g., 
Fig. <ref>(d,e)]. For the B=2.5 T case, we perform real-time evolutions on 
a YC6 × 30 lattice with bond dimension D=600.

In Figs. <ref>(a,b) we show the results of spin-resolved spectral functions 
A^zz(k,ω) and A^+-(k,ω), which exhibit distinct behaviors. 
The spectral functions, rather than the dynamic spin structure factors, are shown 
in Fig. <ref>, which do not include the elastic-scattering peaks at K 
and allow us to concentrate on the low-energy fluctuations. From Fig. <ref>(a) 
we find that the spectral intensities of A^zz(k,ω) that mainly reflect 
the in-plane excitations are significant only for ω≲ 0.15 meV, 
i.e., below about 2 J_xy. On the other hand, the out-of-plane excitations 
reflected in A^+-(k,ω) can extend to higher energies of about 
0.25  meV∼ 2J_z. In addition, in Figs. <ref>(a,b) we find clear 
single-particle excitations for 0 < ω < J_xy, z that we dub as 
magnon-roton dispersions (see discussions below), while for higher energies 
J_xy, z < ω < 2J_xy, z both spectra show excitation continuum
<cit.>.

Goldstone and pseudo-Goldstone magnons.—
To examine the low-energy excitations, we gradually improve the energy resolution 
to about 0.08J_z (t_ max = 100/J_z) in the calculation, and show the 
results of A^zz, +-(k, ω) in Figs. <ref>(a,b). 
where the spectral functions become more coherent as ε decreases. 
The spectral function A^zz(k, ω) reflects the in-plane excitations. As 
the spectral function is parity odd in ω, its peak value Δ_ G 
has been artificially shifted to higher frequencies upon convolution with window 
functions. As shown in Fig. <ref>(d), Δ_ G becomes lowered 
as the resolution ε improves, which extrapolates eventually to
approximately zero energy in the ε=0 limit, indicating the existence 
of gapless Goldstone modes [c.f., Fig. <ref>(c)]. 

On the other hand, in Fig. <ref>(b) we show the low-energy out-of-plane 
excitations modes by computing the spectral function A^+-( K, ω). 
As shown in Fig. <ref>(d), a small but nonzero gap Δ_ PG exists, 
which converge to Δ_ PG≃ 0.063J_z. This is the 
pseudo-Goldstone mode that originates from the modified “Mexican-hat” 
energy landscape with six-fold anisotropy as illustrated in Fig. <ref>(d). 
A complex order parameter Ψ =e^ i θ 
|Ψ| ≡1/N (∑_i∈ A⟨ S^z_i⟩+ ∑_j∈ B⟨ S^z_j⟩ e^ iπ 2/3+ ∑_k ∈ C⟨ S^z_k⟩ 
e^ iπ 4/3) can be introduced, where A,B,C label the three sublattices 
and the U(1) phase θ reveals the “hidden” XY degree of freedom. 
In Fig. <ref>(d), the spin configurations like ↑↑↓ 
and ↑↓↑, etc., correspond to the six-fold degenerate 
ground state with θ = n π/3 (0 ≤ n ≤ 5). The spin excitations in the 
vicinities of 6-fold minima cost a finite amount of energy, and this small pseudo-Goldstone gap is generated by quantum fluctuations via the order-by-quantum-disorder 
mechanism <cit.>.

Double magnon-roton excitations in the spin supersolid phase.—
Besides the conventional phonon dispersion, in superfluid helium-4 there 
exist an anomalous dip — the roton mode — in the excitation spectrum of 
superfluid helium-4. The distinctive phonon-roton dispersion curve was first 
hypothesized by Landau through his seminal work <cit.>, and subsequently substantiated and refined by Feynman
<cit.> by developing a microscopic theory to 
elucidate this feature. The Landau elementary phonon-roton excitations 
<cit.> play an essential role in forming the thermodynamic 
and hydrodynamic characteristics of this quantum fluid <cit.>. 

In isotropic TLAF systems, there have also been theoretical investigations 
on the roton-like minima in spin excitations <cit.>, 
There are theoretical work shedding light on the roton excitations in isotropic 
Heisenberg TLAF <cit.>, and experimental evidence of such magnetic rotons also reported 
<cit.>. This constitutes a 
reminiscent of the phonon-roton spectral characteristics observed in superfluid 
helium-4 <cit.>. The relationship between 
roton excitations and superfluidity remains a subject of intense research, 
with much still to be understood. It is a compelling question to investigate whether 
a magnetic counterpart to roton excitations exists within the spin supersolid state. 

In Figs. <ref>(a,b), we find in both cases there are magnon-roton 
dispersions consisting of linear dispersion and soft quadratic exciations 
near M≡ [1/2, 1/2] points. As summarized in Fig. <ref>(c), 
there are two branches of excitations, where the lower magnon-roton 
dispersion can be associated with the in-plane spin superfluidity, as a 
magnetic analog of phonon-roton dispersion in superfluid helium. Remarkably, 
there is a second magnon-roton dispersion that can be ascribed to the 
out-of-plane fluctuations of spin solidity. Despite a finite pseudo-Goldstone 
gap, the six-fold anisotropy in the pseudo-spin Ψ [see Fig. <ref>(d)] 
becomes irrelevant at elevated temperature and there is an emergent U(1) 
symmetry in the system <cit.>. 
This can give rise to the Berezinskii-Kosterlitz-Thouless (BKT) transition
<cit.> similarly as in triangular lattice quantum
Ising antiferromagnets <cit.>






Thermodynamics of magnon-roton excitations.—
As firstly noticed by Landau <cit.>, rotons are activated at a 
temperature much lower than the roton gap, and thus contribute significantly
to the low-temperature thermodynamics of superfluid helium <cit.>,
due to the very large density of states of roton excitations. In Fig. <ref>(f), 
we show the A(ω) obtained by integrating over the momentum k, 
i.e., A(ω) = 1/𝒩∑_k  A^zz(k,ω) + 
A^+-(k,ω), with 𝒩 a normalization factor such that 
∫ A(ω) dω=1. A prominent peak of A(ω) can be 
observed near the roton mode ω∼Δ_ R1, R2, which 
may strongly influence the low-temperature properties.

In Fig. <ref>, we show the calculated entropy S_ tot and specific 
heat C_ m from thermal tensor-network calculations. By taking 
A(ω) in the particular energy window as the effective density of states 
of the roton excitations <cit.>, we also estimate the roton 
contributions S̃_ R and C̃_ R at low temperature. 
From the results in Fig. <ref>, we find that due to the prominent spectral peak, 
the roton modes have significant contributions even below T_l ≃ 0.2 K, 
despite the considerable roton gap of Δ_ R1, R2≈ 0.78 K. 
Therefore, the double magnon-roton excitations with roton dips significantly 
influence the low-temperature thermodynamics, explaining naturally the giant 
MCE observed in  <cit.>.


Discussion and outlook.—
In this work, we find a coexistence of three-sublattice spin solidity and 
gapless Goldstone magnons, and provide spectroscopic evidence for spin 
supersolidity in NBCP. Predictions have been formulated concerning the 
emergence of double magnon-roton excitations and the pseudo-Goldstone 
gap within the spin supersolid phase. To further explore these phenomena, 
INS measurements with enhanced energy resolution are required.

Beyond NBCP, the main conclusions here on the dynamical properties 
may also apply to other similarly structured triangular-lattice compounds. 
Very recently, emerging evidence suggests the presence of a spin 
supersolid phase within the triangular-lattice cobaltate compound 
K_2Co(SeO_3)_2, despite of a different extent of easy-axis 
anisotropy <cit.>. Owing to the 
substantial spin exchange interactions present in this triangular-lattice 
cobaltate, we anticipate that observing the predicted dynamical signatures, 
including roton modes and pseudo-Goldstone gap, etc, may require 
less stringent experimental conditions compared to .

Note added.— In the finalization of the present work, we get aware 
of two recent studies also on the spin dynamics of spin supersolid phase in  <cit.>, with main conclusions
consistent with our present study.


Acknowledgments.—
W.L., W.J., and Y.G. express their gratitude to Tao Shi, Changle Liu, 
Shang Gao, Bing Li, and Hai-Jun Liao for stimulating discussions. 
This work was supported by the National Natural Science Foundation 
of China (Grant Nos. 12222412, 12047503, and 12074023), and the 
CAS Project for Young Scientists in Basic Research (Grant No. YSBR-057). 
We thank the HPC-ITP for the technical support and generous allocation 
of CPU time. The Australian Center for Neutron Scattering are gratefully 
acknowledged for providing neutron beam time through proposal No. P17086.














Supplementary Material for

 Spin Supersolid Phase and Double Magnon-Roton Excitations in a Cobalt-based Triangular Lattice

 

Yuan Gao, et al.


















§ DATA ANALYSIS OF THE NEUTRON SCATTERING MEASUREMENTS

The low energy range with ω∈  [-0.05, 0.05] meV of the INS data in zero 
field was integrated and treated as the elastic scattering results. As shown in Fig. <ref>(a), the coexistence of bright magnetic Bragg 
peaks with the propagation vector of k = (1/3, 1/3, 0.16) and diffusive rod-like 
scatterings are observed. The latter is along the out-of-plane direction and 
suggest a quasi two-dimensional nature of NBCP.

The integrated intensities of the magnetic Bragg peaks on top of the 
diffuse scatterings were extracted for further analysis. As shown in Fig. <ref>(b), the intensities of the four non-equivalent 
reflections agree with an UUD configuration of the Co^2+ moments
along the c-axis described by the irreducible representation Γ_1. 
The moment sizes on the z = 0 layer are estimated to be 0.606(27), 
-0.303(13) and -0.303(13) μ_B, for the `down', `up', and `up' 
spins in three sublattices, respectively. The results are well consistent with 
our previous neutron diffraction results on NBCP also under zero field
<cit.> and supports the presence of out-of-plane spin 
solidity in the compound.





The low-energy part of the 2.5 T data, with ω∈  [-0.1, 0.1] meV, 
is used for background subtraction for the zero-field INS data. As shown in Fig. <ref>(c), the spin excitations observed with an in-plane field of B=2.5  T are clearly gapped. Therefore, we utilize the low-energy 
part of the 2.5 T data as the intrinsic background to be subtracted for the 
0 T case shown in Fig. <ref>(b), and obtain the result shown in 
Fig. <ref>(a). By integrating the scatterings along the out-of-plane 
[0, 0, ζ] direction, it is found that the spin excitations emanating and 
away from the ordering vector look quite similar, as shown in 
Fig. <ref>(a) for ζ∈  [-0.2, 0.2] and 
Fig. <ref>(b) for ζ∈  [0.3, 0.7], 
also indicating clearly a good two dimensionality of NBCP. 













§ GROUND STATE DYNAMICAL CALCULATIONS



 §.§ Derivation of dynamical spin structure factor and spectral function

In this section, we detail the derivation of ground-state dynamical spin structure 
factor S(k,ω) and spectral function A(k,ω). We start form the 
real-time correlation function 

    g^αβ(k,t)≡1/N e^ i E_0 t∑_i,j e^- i k·(r_i-r_j)⟨ψ_0|S^α_i e^- i H t 
    S^β_j |ψ_0⟩,

where |ψ_0⟩ is the ground state, E_0 is the ground-state energy, 
N is the total site number, and (S_i^α)^† = S_i^β. 
Then we take the complex conjugate and arrive at

    g^αβ(k,t)^* = 1/N e^- i E_0 t∑_i,j e^- i k·(r_j-r_i)⟨ψ_0|S^α_j e^ i H t 
    S^β_i |ψ_0⟩ = g^αβ(k,-t).

With the real-time correlation function, we obtain the dynamical spin structure 
factor S(k,ω)

    S^αβ(k,ω)    ≡∫_-∞^∞ g^αβ(k,t) e^ iω t dt 
    =  ∫_0^∞ g^αβ(k,t) e^ iω t dt + 
    ∫_0^∞ g^αβ(k,-t) e^- iω t dt 
       = 2 ∫_0^∞ Re[g^αβ(k,t) e^ iω t] dt.

Similar, we can calculate the spectral function A(k, ω) ≡ 
-1/π Im [G^ R(k,ω)] where G^ R(k,ω) 
is the retarded Green's function 

    G^ R, αβ(k,ω) ≡    - i/N∑_i,j e^- i k·(r_i-r_j)∫_0^∞ e^ iω t⟨ψ_0|[S^α_i(t), S^β_j(0)]|ψ_0⟩ dt
    
    =    - i/N∑_i,j e^- i k·(r_i-r_j)∫_0^∞ e^ iω t (e^ iE_0 t⟨ψ_0|S^α_i 
    e^- iHt S^β_j|ψ_0⟩ - e^- iE_0 t⟨ψ_0|S^β_j 
    e^ iHt S^α_i|ψ_0⟩) dt
    
    =    - i∫_0^∞ e^ iω t [g^αβ(k,t) - g^βα(-k, -t)] dt.

Note that the Hamiltonian and the ground state is invariant under the space 
reversal transformation r_i → -r_i, we have g^αβ(k,t) 
= g^αβ(-k,t). Besides, with zero magnetic field, the Hamiltonian 
is invariant under the spin flip transformation S_i^z → -S_i^z, S_i^y 
→ -S_i^y, thus we have g^+-(k,t) = g^-+(k,t). Note that we only
consider A^+-(k,ω), A^zz (k,ω) under zero field, we can obtain 
the spectral function as 

    A^αβ(k, ω) ≡    -1/π Im 
    [G^ R,αβ(k,ω)]
    = -2/π Im  [∫_0^∞ [cos(ω t) +  isin(ω t)]   Im[g^αβ(k,t)]]
    
    =    -2/π∫_0^∞ Im[g^αβ(k,t)] sin(ω t) dt.


In the calculations, we first obtain |ψ_0⟩ using density matrix 
renormalization group <cit.> and then perform real-time evolution 
to simulate |ψ'(t)⟩≡ e^- i H t S^β_j |ψ_0⟩ with time 
dependent variational principle (TDVP) <cit.>. 
Having acquired the real-time correlation function, we proceed to calculate the 
dynamic structure factor and the spectral functions, which are convolved with 
an appropriate window function to account for experimental broadening effects,
i.e.,

    S^αβ(k,ω)    = ∫_0^t_ max Re 
    [g^αβ(k,t) e^iω t]   W(t/t_ max) dt,
    
    A^αβ(k,ω)    = -∫_0^t_ max Im 
    [g^αβ(k,t)] sin(ω t)  W(t/t_ max) dt,

where W(t) is the Parzen window function, and t_ max is the 
maximal TDVP evolution time (in natural unit). The energy resolution 
ε≃ 8/t_ max as determined by W(t/t_ max),
the full width at half maximum (FWHM) of its Fourier transform. 
In practical calculations, we perform real-time evolution with retained 
bond dimension up to D=3000 on a YC6× 15 lattice for zero 
field and bond dimension D=600 on a YC6× 30 lattice under 
in-plane field of B=2.5 T, with the involved Y-type cylindrical lattice 
shown in Fig. <ref> below.







 §.§ The estimation of energy resolution
 
The energy resolution ε is determined by W(t/t_ max), 
as the FWHM of its Fourier transform. In practice, we choose 
the Parzen window function following as 

    W(t)=
    
    -2 (-1+t)^3    1/2<t≤ 1 
    
     2 (1+t)^3     -1≤ t<-1/2
    
     1-6 t^2-6 t^3     -1/2≤ t<0 
    
     1-6 t^2+6 t^3     0≤ t≤1/2
    
     0     otherwise
    ,

As the Fourier transformation reads

    F[W(t)](ω) = 6  √(2/π)e^- iω (-1 + e^ iω/2)^4/ω^4,

the energy resolution FWHM can be obtained as 
ε≃ 8/t_ max.



 §.§ Estimation of density of states

Here we show the details of density of states (DOS) estimation. 
We consider the local spectral density A(ω) ≡1/N∑_k∈ BZ A(k,ω), 
where BZ is the first Brillouin zone and N is the system size.
According to Eq. (<ref>), we have

    A^αβ(ω) ≡∑_k∈ BZ 
    A^αβ(k,ω) = -2/π Im  
    [∫_0^∞ e^ iω t Im[1/N∑_k g^αβ(k,t)]].

Note that 

    1/N∑_k g^αβ(k,t) =   1/N
    e^ i E_0 t∑_k⟨ψ_0|S^α_-k e^- i H t 
    S^β_k|ψ_0⟩
    
    =   1/N e^ i (E_0-E_m) t∑_k, m  ||⟨m|S^β_k |ψ_0⟩||^2,

where S^β_k = 1/√(N)∑_i e^ i k r_i S_i^β, 
{|m⟩} are the eigenstates of H with energy E_m 
and we assume (S_-k^α)^†=S_k^β. Substitue 
Eq. (<ref>) into Eq. (<ref>), we have

    A^αβ(ω) =    -2/Nπ∑_k,m  
    ||⟨m|S^β_k |ψ_0⟩||^2
    ∫_0^∞ Im [e^ iω t Im[e^ i(E_0-E_m)t]] dt
    
    =   1/Nπ∑_k,m  
    ||⟨m|S^β_k |ψ_0⟩||^2
    ∫_0^∞ Re[e^ i (ω + E_0 - E_m) t] - 
     Re[e^ i (ω - E_0 + E_m)] dt 
    
    =   1/N∑_k,m ||⟨m|S^β_k |ψ_0⟩||^2 
    (δ(ω + E_0 - E_m) - δ(ω - E_0 + E_m)).

Here we introduce the excitation energy ε_m ≡ E_m-E_0 
and only consider the positive energy part, and we can arrive at 

    A^αβ(ω > 0) = 1/N∑_k,m ||⟨m|S^β_k 
    |ψ_0⟩||^2  δ(ω - ε_m).

Regarding the low-energy excitation states |m⟩ as free gas of 
Bogoliubov quasi-particle with energy ε_k, the Hamiltonian 
can be represented as H ≃∑_k ε_k γ_k^†γ_k. The excitation states can be represented as |m⟩≃Π_k'1/√(n_k'!) (γ_k'^†)^n_k'|ψ_0⟩,
where γ_k is the Bogoliubov quasi-particle operator. Therefore, 
Eq. (<ref>) can be rewritten as 

    A^αβ (ω > 0) ≃1/N∑_k 
    ||⟨ψ_0|γ_k S^β_k |ψ_0⟩||^2  δ(ω - ε_k),

where we assume multi-magnon excitation states have vanishing contributions,
and ||⟨ψ_0|γ_k S^β_k |ψ_0⟩||^2 ∼𝒪(1).
With this (crude) approximation, we arrive at 

    A^αβ (ω > 0) ≈∑_k δ(ω - ε_k),

estimates the density of state of magnon excitations with energy ε_k. 
At low temperature, we treat the magnons and rotons as free boson gas, 
and the entropy and specific heat according to 

    S̃(T) =   ∫_0^∞ A(ω) [ω/T/e^ω/T-1-
     ln(1-e^-ω/T)] dω, 
    C̃(T) =    T∂S̃(T)/∂ T,

where A(ω) ≡ A^+-(ω) + A^zz(ω) is normalized such 
that ∫_0^∞ A(ω) dω= 1. 







 §.§ Data convergence

Below we show in Fig. <ref> the calculated spectral function 
results with different bond dimensions retained. By increasing the bond 
dimension, the magnon-roton excitations become more clear and the 
downward renormalization at M≡ [1/2, 1/2] point becomes
more prominent. Besides, we also show the spectral functions with 
different bond dimensions, momentum points, and energy resolutions. 
As shown in Fig. <ref>, we find the spectral functions are well 
converged with a retained bond dimension D=3000.








 §.§ High-energy excitation continuum

In this section, we discuss about the high-energy continuum observed in both 
INS experiments and tensor-network calculations. In Fig. <ref>(a) we 
show the experimental results, and in Figs. <ref>(b-f) spectral functions 
with different energy resolutions are presented, where the W-shape excitation 
continuum does not sharpen as the energy resolution improves.







 §.§ Fitting of magnon-maxon-roton excitations

Here we show the  magnon-maxon-roton fitting of the calculated spectral 
function. The terminology magnon, maxon, and roton originates from the 
excitations spectrum of superfluid helium-4 <cit.>. Here 
a magnon represents the gapless linear excitation, also known as the 
Goldstone mode arising from U(1) symmetry breaking. The maxon 
corresponds to the quadratic excitation near the peak of the dispersion 
relation, whereas the roton denotes the quadratic excitation at the bottom 
of the dispersion curve, characterized by a finite minimum energy. 
We adopt the undetermined function with the following form

    f(ξ)=
    
    -k_0(ξ-1/3) + Δ_0,    ξ≤ 1/3 
    
    k_0(ξ-1/3) + Δ_0,     1/3 < ξ≤ξ_0 
    
    -a_ M (ξ-ξ_ M)^2 + Δ_ M,    ξ_0 < ξ≤ξ_1
    
    a_ R (ξ-1/2)^2 + Δ_ R,    ξ_1 < ξ

where {k_0, Δ_0, ξ_0, a_ M, ξ_ M, ξ_1, a_ R, 
Δ_ R} are the fitting parameters. The first two lines are the 
(pseudo-)Goldstone part, the third line is the maxon part, and the last 
line is the roton dispersion. Note that f(ξ) should be a continuous 
function, thus the fitting parameters can be reduced to {k_0, Δ_0, 
ξ_0,  ξ_ M, ξ_1, Δ_ R}. With numerical results of 
the dispersion ε_k obtained from the calculated spectral 
functions, these parameters can be fitted and the corresponding results 
are shown in Fig. <ref>.  In examining Goldstone mode excitations 
determined from A^zz(k, ω), we assume Δ_0 to be zero, 
consistent with the anticipated symmetry breaking scenario. Conversely, 
when analyzing the A^+-(k, ω) excitation channel, we allow 
Δ_0 to vary, serving as a free parameter to be precisely fitted to 
the data.








§ LINEAR SPIN WAVE CALCULATIONS
 

Details of the linear spin-wave calculations of easy-axis TLAF model
H = ∑_⟨ i,j⟩ J_xy (S_i^xS_j^x+S_i^yS_j^y) + 
J_z S_i^zS_j^z are shown below, where Δ = J_z/J_xy. 
The ground state of the easy-axis TLAF model with Δ > 1 is a 
Y-shaped state in the x-z plane. There are three sublattices, namely A, 
B, and C, and thus three kinds of Holstein-Primakoff bosons, 
a, b, and c are introduced. For the sublattice A, we have

    S^z = (S - a^† a),
        S^x = √(2S)/2 (a + a^†),
        S^y = √(2S)/2 i (a - a^†).
 
On the other two sublattices, there are angles ±θ between spins 
on A and B(C) sublattices. For sublattice B, the transformation reads

    S^z=cosθ  (S-b^† b)- sinθ √(2S)/2(b+b^†),
       S^x=sinθ  (S-b^† b)+cosθ √(2S)/2(b+b^†),
       S^y=√(2S)/2 i(b-b^†),

and for sublattice C

    S^z=cosθ  (S-c^† c)+ sinθ √(2S)/2(c+c^†), 
       S^x=-sinθ  (S-c^† c)+cosθ √(2S)/2(c+c^†),
       S^y=√(2S)/2 i(c-c^†).


Through the Holstein-Primakoff and Fourier transformations, 
we arrive at a quadratic form of the Hamiltionian 
    H_ HP 
    = ∑_k α_k^† H_0(k) α_k,
 where 
α_k^†  = (a_k^† b_k^† c_k^† a_-k b_-k c_-k) denote the magnon creation operators. 
The quadratic Hamiltonian is H_0(i,j)= S   A_0(i,j)   Z(i,j), 
where the k-independent symmetric part is A_0=A + A^ T, 
i,j label the matrix index, and S=1/2 for the present case. 
A is a 6×6 upper triangular matrix

    A=
    [ M N; 0 M ],

with 

    M=
    [                     -3Δcosθ                    cosθ+1/2                    cosθ+1/2;                           0 3/2(sin^2θ-Δcos^2θ - Δcosθ)         cos^2θ +1-Δsin^2θ/2;                           0                           0 3/2(sin^2θ-Δcos^2θ - Δcosθ);                             ],

and

    N=
    [                   0            cosθ-1/2            cosθ-1/2;            cosθ-1/2                   0 cos^2θ -1-Δsin^2θ/2;            cosθ-1/2 cos^2θ -1-Δsin^2θ/2                   0 ],

Δ = J_z/J_xy is the anisotropic parameter,
and Z(k) is a 6× 6 matrix 

    Z(k) = 
    [   1   z z^*   1   z z^*; z^*   1   z z^*   1   z;   z z^*   1   z z^*   1;   1   z z^*   1   z z^*; z^*   1   z z^*   1   z;   z z^*   1   z z^*   1;     ],

with z=∑_i e^ i k δ_i with δ_1 = (1,0), δ_2 = 
(-1/2, √(3)/2), δ_3 = (-1/2, -√(3)/2).
The angle θ can be obtained by minimizing ground state energy 
E_0 = S^2(2Δcosθ  + Δcos^2θ  - sin^2 θ ),
and for the realistic parameter Δ = 1.68 we have θ = 128.91^∘.
With this, we diagonalize H_HP and obtain the linear spin-wave dispersion
shown in Fig. 2 of the main text.

