Complete field-induced spectral response of the spin-1/2 triangular-lattice antiferromagnet CsYbSe2

Fifty years after Anderson’s resonating valence-bond proposal, the spin-1/2 triangular-lattice Heisenberg antiferromagnet (TLHAF) remains the ultimate platform to explore highly entangled quantum spin states in proximity to magnetic order. Yb-based delafossites are ideal candidate TLHAF materials, which allow experimental access to the full range of applied in-plane magnetic fields. We perform a systematic neutron scattering study of CsYbSe2, first proving the Heisenberg character of the interactions and quantifying the second-neighbor coupling. We then measure the complex evolution of the excitation spectrum, finding extensive continuum features near the 120°-ordered state, throughout the 1/3-magnetization plateau and beyond this up to saturation. We perform cylinder matrix-product-state (MPS) calculations to obtain an unbiased numerical benchmark for the TLHAF and spectacular agreement with the experimental spectra. The measured and calculated longitudinal spectral functions reflect the role of multi-magnon bound and scattering states. These results provide valuable insight into unconventional field-induced spin excitations in frustrated quantum materials.


INTRODUCTION
Frustrated quantum magnets provide an intriguing playground for investigating novel many-body phenomena in condensed matter 1,2 .The triangular-lattice Heisenberg antiferromagnet (TLHAF) is a prototypical example of geometrical frustration, and its ground state in the quantum limit of S = 1/2 spins has the 120 • AF order of the classical (large-S) case 3 , albeit with an ordered moment strongly suppressed by quantum fluctuation effects ( 4 and references therein).Proposals to capture these effects include the resonating valencebond (RVB) paradigm 5 , and the addition of a weak nextneighbour HAF interaction (0.06 ≲ J 2 /J 1 ≲ 0.15) does drive the S = 1/2 TLHAF into a quantum spin-liquid (QSL) phase [6][7][8][9] of some type [10][11][12] .Over a finite range of applied magnetic fields, AF quantum fluctuations favour a collinear up-up-down (UUD) ordered phase and thus stabilize a magnetization plateau with M = M Sat /3 (where M Sat is the saturation magnetization) 13,14 .
In this context, the Yb-based delafossite family AYbQ 2 , with A an alkali metal and Q a chalcogenide, has attracted widespread attention 28,29 .The Yb ions form perfect and well separated TLs (Fig. 1a) [30][31][32] without the structural disorder intrinsic to YbMgGaO 4 33,34 .The combination of strong spin-orbit coupling and the crystalline electric field (CEF) creates a ground-state doublet that gives an effective S = 1/2 pseudospin at low temperatures 30,31,35,36 .Although the J = 7/2 CEF level structure is manifest in a strong spatial anisotropy of the response to applied magnetic fields 37,38 , initial scattering studies provided no evidence for a strongly non-Heisenberg pseudospin Hamiltonian 39,40 .Early specific heat, magnetization, muon spin-rotation spectroscopy and neutron diffraction studies of multiple AYbQ 2 materials found no magnetic order at zero field down to their base temperatures 31,37,[39][40][41] , but recent studies, including our own (Fig. 1b), indicate its presence in some materials at the 0.1 K scale.The 1/3-magnetization plateau is found at in-plane fields in the 3-5 T range 37,39,40,42 (Fig. 1c), with robust UUD order up to 1 K. Clearly the delafossite family offers an excellent platform to study the field-controlled magnetic states of the S = 1/2 TLAF.The solid and open circles represent respectively the temperatures of sharp peaks and broad humps in the corresponding specific-heat curves, as described in Supplementary Note 1D.The arrows indicate schematically the spin order of the five phases of the classical TLHAF, which are consistent with the magnetic peaks we observe.e INS spectrum measured under a magnetic field of 5 T, showing an absence of well defined ∆S = 1 excitations away from the Γ and K points but extensive continuum features (Fig. 3).The grey horizontal bar at low energy masks the elastic-line contribution.f Comparison with a matrix-product-state (MPS) calculation of the spectrum at the same field.
single-crystalline Yb delafossites at zero field 40,43 suggested a gapless excitation continuum, which was interpreted as originating from a QSL ground state, but appears to persist even in the presence of weak magnetic order 44 .Early INS studies of the spin dynamics in the field-induced phases were limited by their polycrystalline samples 39,45 , but the 1/3-magnetization plateau has recently been analyzed in some detail 46 .As in Ba 3 CoSb 2 O 9 , where the plateau has been reached despite the higher energy scales in this family of materials 47 , the magnetic excitations were captured largely by semiclassical nonlinear spin-wave theory (SWT).The lower in-plane energy scale and vanishing inter-plane coupling in delafossites present more experimental challenges, but also the key advantage of reaching saturation within laboratory-available magnetic fields.
On the theoretical side, the challenge of computing the dynamical spectral functions of frustrated models lies in the absence of analytical methods that capture all the physics of non-collinear magnetic states with a field-controlled ratio of weak order to strong quantum fluctuations.The application of unbiased numerical methods, meaning those whose truncation methodology can be extended systematically to convergence, has in the past been impossible, but continuous progress in dynamical quantum Monte Carlo techniques and matrix-product-state (MPS) representations is placing this goal within reach.For the TLHAF, the zero-field spectral function has been obtained by a number of biased methods, by which we mean those based on initial assumptions that have to be assesed a posteriori; these include series expan-sions 48 , interacting spin waves 49,50 , Schwinger bosons 51,52 , bond operators 53 and variational Monte Carlo 54 .Despite recent progress with MPS calculations in a cylinder geometry 11,12 , the full field-induced dynamics has remained an unsolved problem.All of these methods produce scattering continua whose origin may lie in fractional excitations, multimagnon states or possibly neither.
In this work we perform high-resolution neutron spectroscopy on CsYbSe 2 using two different spectrometers to span the full range of applied in-plane fields, meaning from zero to beyond saturation.Our measurements reveal the pronounced changes in the magnetic excitation spectrum as it evolves with the magnetic field, and we associate these with the field-driven phase transitions of the ground state (Fig. 1d).In parallel we perform large-cylinder MPS calculations of the full TLHAF spectral function at all fields to obtain a hitherto unavailable benchmark for the model, semi-quantitative agreement with experiment (Figs.1e-f) and a robust foundation for any effective quasiparticle descriptions of the spin dynamics.

Ground state at zero field
The growth and structural characterization of our single crystals are summarized in the Methods section and detailed in Supplementary Note 1A.Neutron diffraction at zero magnetic field (Supplementary Note 2) reveals a series of weak magnetic intensity peaks at Q = (1/3, 1/3, L) for odd-integer L (Fig. 1b).The (1/3, 1/3, 1) peak develops a finite intensity below T ≃ 0.4 K, which increases on cooling.The propagation wavevector matches the 120 • state of the TL-HAF with AF out-of-plane correlations (represented by the arrows in Fig. 1a) and the low-temperature ordered moment is m Yb ≃ 0.1 µ B .In Supplementary Note 2B we extract the in-and out-of-plane correlation lengths, ξ ab = 60 (7)  Å and ξ c = 23(5) Å, which are not resolution-limited, meaning that CsYbSe 2 does not exhibit true, long-ranged AF order at zero field down to T = 0.02 K.However, the presence of the magnetic peak clearly excludes a QSL, as in KYbSe 2 44 but in contrast to NaYbSe 2 43 .Given that CsYbSe 2 (space group P6 3 /mmc) has AA layer stacking (Fig. 1a), which should favour an unfrustrated collinear c-axis order, whereas the Na, K and Rb materials (space group R3m) have an ABC stacking that should produce interlayer frustration, we suggest in Supplementary Note 1B that the origin of this behaviour may instead lie in the next-nearest neighbour coupling, J 2 (below).

Magnetic phase diagram of CsYbSe 2
We performed low-temperature magnetization, specific-heat and neutron diffraction measurements over a wide field range, as described in Supplementary Notes 1 and 2. Figure 1c shows isothermal magnetization data, with evidence of a plateau at M Sat /3 corresponding to the UUD phase 13,14 .To interpret these data despite their finite-temperature rounding, we estimate B Sat = 9.6(2) T from the TLHAF model parameters obtained by high-field INS (Fig. 2) and perform a grand canonical DMRG calculation 55 of the magnetization (Supplementary Note 4) that allows us to deduce the boundaries of the 1/3-magnetization plateau.We compare these data with the integrated intensity of the (1/3, 1/3, 1) magnetic peak, which increases strongly from 0 T to the UUD state, then remains maximal and almost constant over the field range of the plateau (3 T ≤ B ≤ 4.5 T), before decreasing strongly towards the fully polarized (FP) state.Our thermodynamic and neutron diffraction results yield the field-temperature phase diagram shown in Fig. 1d, where we indicate the spin alignments of the classical TLHAF.

Parameters of the magnetic Hamiltonian
We made INS measurements up to 5 T on the time-of-flight (ToF) spectrometer CNCS at ORNL and up to 11 T on the multiplexing spectrometer CAMEA at PSI, as detailed in the Methods section.To quantify the parameters of the spin Hamiltonian, we exploit our ability to perform INS at fields B > B Sat , where the magnetic excitations of the FP phase can be described by linear SWT. Figure 2b shows that the spectrum measured along [H H 3] at 11 T consists of a single, sharp magnon mode with a cosinusoidal dispersion above a field-induced gap at the K point.In Supplementary Note 2D we show CNCS data indicating a complete lack of outof-plane dispersion over the whole field range, and hence that a two-dimensional TL model is appropriate.We fit this dispersion using the SPINW package 56 by considering a Hamiltonian with an anisotropic XXZ-type J 1 term and a Heisenberg J 2 term (Supplementary Note 3A).With the in-plane gfactor fixed (below), the optimal fit (solid red line in Fig. 2b) yields two essential pieces of information.First, the nearestneighbour interaction has no XXZ anisotropy within the precision of the measurement, i.e. despite the strongly anisotropic field response 29,38 , the spin dynamics are of Heisenberg type; in Ref. 35 it was shown how these contrasting forms of behaviour can appear simultaneously in edge-sharing octahedral Yb 3+ systems.Second, the next-neighbour interaction is sufficiently weak, J 2 /J 1 ≃ 0.03, that CsYbSe 2 remains on the ordered side of the phase boundary separating the 120 • and QSL states in the S = 1/2 J 1 -J 2 TLHAF [6][7][8][9] .We therefore conclude that a J 1 -J 2 Heisenberg Hamiltonian with J 1 = 0.395 meV, provides a complete description of the low-energy magnetic behaviour in CsYbSe 2 .
Turning to fields below saturation, we begin in Fig. 2d by considering constant-Q cuts at the Γ-point for each field.These show a clear, single-peak feature for B ≥ 2 T, whose energy obeys the field-linear form ℏω(B) = ℏω 0 + g ab µ B BS with g ab = 3.20 (6) and ℏω 0 = 0.00(2) meV.This fieldinduced behaviour at Γ is generic for the Heisenberg model 57 , reinforcing our conclusion concerning the absence of XXZ anisotropy, and in the TLHAF is also present at K (Fig. 2d).To characterize the anisotropic field response, we have performed electron spin resonance (ESR) measurements that determine the g-tensor parameters shown in Fig. 2c.The narrow and well-defined ESR spectrum reflects the high quality of our crystal.The best fit in Fig. 2c, g ab = 3.25(0) and g c = 0.3(0), completes our determination of the Hamiltonian  parameters for CsYbSe 2 in any applied field and also shows the consistency of our INS result for g ab .Our measurements also demonstrate that g ab is isotropic in the ab plane to within the experimental accuracy (data not shown).The very small g c is a consequence of strong hybridization with the first excited CEF doublet 38 , and we comment below on its role in our scattering study.

MPS calculations of spectral functions
To interpret the measured spin dynamics, we have performed MPS calculations on a finite cylinder to obtain the dynamical spectral function of the TLHAF with J 2 = 0.03J 1 , where the energy unit is fixed to J 1 = 0.395 meV.The cylinder size, matrix bond dimensions and time-evolution procedures are summarized in the Methods section and their convergence to the properties of the TLHAF is benchmarked in Supplementary Note 5. We compute the spin correlation functions both transverse and longitudinal to the applied field, which in the notation of Fig. 2a are respectively S xx (Q, ω) and S zz (Q, ω).In experiment, the strong g-tensor anisotropy (g ab ≫ g c ) means that the component fluctuating parallel to the c axis (S yy ) is hidden [Fig.2a defines the (HKL) and (xyz) coordinate frames].The measured intensities then represent the sum of S xx and S zz weighted by the polarization factor, which ensures that spectra taken along [H H 0] have no contribution from S xx , i.e.only from the component longitudinal to the applied field.In order to sample both components, in Fig. 3 we integrate our INS and MPS spectra over a wide range of the out-of-plane momentum, L.

Field-induced evolution of the spectrum
As expected from the phase diagram (Fig. 1d), both the observed and computed spectra in Fig. 3 are readily classified by their field-induced evolution into four regimes, to which we refer as Y, UUD, V and FP (the last analysed in Fig. 2).Starting with Y, we have shown (Fig. 1b) that the zero-field ground state of CsYbSe 2 is consistent with threesublattice 120 • order, and thus the spectrum should contain three excitation branches.However, both the INS and MPS spectra exhibit only a broad, V-shaped continuum around K (Figs. 3a-b), similar to the spectra observed in NaYbSe 2

43
and KYbSe 2 44 .Because this clear signature of strong quantum fluctuations on top of weak magnetic order has received extensive theoretical 49,51,58 and numerical analysis 11,12,44,48,54 , which our results confirm but do not extend, we focus rather on adding to the understanding of the finite-field spectra.At 2 T (Figs. 3c-d), most of the spectral intensity shifts upwards, forming the broad feature, with a gap around 0.4 meV at Γ and K, seen in Fig. 2d.Away from Γ and K, however, it has no well defined magnonic form and the spectrum appears as a weak and highly dispersed continuum.Nevertheless, we mark the maximum intensity of this feature by the circles in Figs.3c-d and refer to it as mode I, observing its bandwidth falling rapidly across the Y regime.Another mode is present at lower energies, whose gapless nature is clearly visible in the MPS spectrum.Linear SWT captures rather well the maximum of mode I, and also finds two gapless branches, but cannot reproduce the extreme broadening of the measured modes away from Γ and K, their intensity distribution or the contin-uum scattering above mode I.
The data collected at 3 and 4 T represent respectively the lower edge and upper middle of the UUD regime (Fig. 1c).In contrast to the Y phase, a number of rather sharp excitations extend across the full Brillouin zone, and at 3 T we identify four distinct features (Figs.3e-f).Mode I shifts upwards, becomes resolution-limited and has intensity over a large Q range.A weak and very low-lying feature II is visible only around its maximum near 0.4 meV.A broad feature III is concentrated around the K point and disperses upwards to touch mode I.A continuum feature IV disperses from around 0.8 meV at K to 1.3 meV at X and M. At 4 T (Figs. 3g-h), features II and III have almost merged to be-come a strong, spin-wave-like branch.Mode I continues its upward shift while continuum IV remains almost unchanged in position and intensity.Modes I-III have been observed in Ba 3 CoSb 2 O 9 47 , and very recently all four features were measured in KYbSe 2 46 .Both studies used a nonlinear SWT to obtain a good account of modes I-III, and in KYbSe 2 it was suggested that feature IV is a two-magnon continuum.Given that the 1/3 plateau is absent in linear SWT, it is not surprising that the orange lines in Figs.3f and 3h provide at best partial agreement with some observed branches.By contrast our MPS calculations provide a quantitatively excellent description of every feature in the observed UUD spectra, which will allow a deeper analysis in Fig. 4.
Proceeding into the V phase at 5 T causes a qualitative modification of the spectrum (Figs.3i-j).As in the Y phase, no clear magnon branches are visible away from Γ and K.More specifically, mode II softens, broadens and decreases in intensity, while mode III merges fully with it.Mode I decays into a broad continuum with sharp intensity peaks only at Γ and K, and obscures continuum IV.Linear SWT traces only the lower boundary of the mode-I continuum and the dispersion of mode II.The 8 T dataset in Fig. 3k shows a sharp band maximum at Γ (mode I) together with a weak replica at K.Here our MPS results (Fig. 3l) clarify how mode I becomes very broad and mode II becomes very soft; linear SWT provides an acceptable guide to the positions, but absolutely not to the emerging mid-zone continuum nature, of these features.In Supplementary Note 6 we present cuts through the data displayed in Figs.3a-j that confirm the near-quantitative agreement between the INS and MPS spectra at almost all points in Q and ω.

Two-magnon bound and scattering states
To analyse the spectra in Fig. 3 we begin in the UUD phase, where all the ordered moments are orientated (anti)parallel to the field and thus the S zz channel contains purely those spin fluctuations longitudinal to the field.Figure 4a shows the longitudinal excitation spectrum obtained from a data slice in the [H H 0] direction and Fig. 4b the analogous MPS calculation.Both spectra show a weakly dispersive, low-energy branch running from X to M, and above this the entirety of continuum IV.To understand the origin of these longitudinal features, we appeal first to the Ising limit, where in the UUD phase a single spin-flip against the field direction costs no energy, whereas the opposite flip costs 3J 1 .If both processes occur on neighbouring spins, the energy cost is only 2J 1 (Fig. 4c), forming a localized two-magnon bound state.The spectrum close to the Ising limit then contains a nearly flat bound-state mode at an energy of 2J 1 , which is clearly split off from a continuum of states that starts around 3J 1 (Supplementary Note 7).In Fig. 4d we show an MPS calculation at rather strong Ising anisotropy that nevertheless shows many properties of the Heisenberg case (Fig. 4b), and in Supplementary Figure 16 we show a more complete interpolation.These results demonstrate clearly the evolution of the lowest localized modes into a split-off and weakly dispersive longitudinal two-magnon bound state, while the upper localized modes evolve into a scattering resonance that forms the characteristic shape of continuum IV.Thus the Ising picture of these features remains valid even at the Heisenberg point.
Increasing the field into the V phase (Figs.3i-j) causes the longitudinal spectrum to show little change, whereas the transverse magnons disintegrate rapidly.Increasing noncollinearity leads to a mixing of transverse and longitudinal character, such that both sets of excitations merge into narrow continua (on the scale of the bandwidth) with strong intensity concentrated only at the Γ and K points.These continuum features become both sharper and more dispersive with increasing field, regaining their single-magnon character above B Sat (Fig. 2b).By contrast, as the field is decreased into the Y phase (Figs.3c-d), the effects of non-collinearity and dominant quantum fluctuations lead to a rapid loss of one-magnon character (again intensity is concentrated only at the Γ and K points) and the emergence of wide excitation continua in both S zz and S xx .

DISCUSSION
Our INS measurements demonstrate unambiguously that the excitations of the TLHAF at all fields consist of magnon-like features only around the Γ and K points that merge into extensive continua across the rest of the Brillouin zone.Despite the presence of at least short-ranged magnetic order at all fields, only in the UUD (1/3-plateau) phase can single magnons provide an adequate basis for describing the spectrum.Capturing the effects of strong quantum fluctuations on such weak order remains a major challenge, which we address by cylinder MPS calculations of the spectrum.Deploying such an unbiased numerical method allows one to divide the process of obtaining physical understanding into a two-step exercise of 'expression' and 'interpretation,' but brings into focus a dichotomy between the two.The expressibility of the MPS method is excellent, in that it captures all the features of the measured excitations with semi-quantitative accuracy, but as a numerical experiment its interpretability is limited.The primary contribution of our study is at the first step, namely providing an unbiased approach that confirms the true spectral content of a paradigm model.At least for the TLHAF, several different biased methods exist that interpret some of the observed spectral features, but to date have lacked a benchmark.
Here it is the agreement between our INS and MPS results which allows us to assert that we have delivered the required benchmark.
We have in addition provided a modern standard for theoretical methods by employing the applied field as a control parameter to access four different, but continuously connected, physical regimes.Thus the ability to separate the transverse and longitudinal spectral functions in one regime affords some key insight that we use to interpret the longitudinal response in the other regimes.When we do consider one biased approach to interpreting the measured and calculated spectra, we find the hallmarks of bound and resonant states of magnon pairs.In the literature it has been argued that scattering continua can arise either from fractionalization (into bosonic 51,52,58 or fermionic 54 components) or from the formation of two-particle and higher-order bound and scattering states of spin-1 excitations 48,49,59 , a subset of the latter being the magnon-breakdown scenario 60,61 .Although none of our present results necessitate a fractionalization scenario to explain the observed spectra, we certainly cannot exclude that fully quantitative analyses of the low-field limit could yet reveal the presence of deconfining S = 1/2 entities in the TLHAF.
To place our results in perspective, to our knowledge the complete field-induced spectrum of a 2D Heisenberg system has not previously been determined in experiment, and here we provide it for the TLHAF realized in CsYbSe 2 .Methodologically, we have used our spectral data to benchmark cylinder MPS calculations of the dynamical spectral function at all applied fields, demonstrating that these now provide a powerful numerical method delivering near-quantitative accuracy.The next-neighbour TLHAF with J 2 on the cusp of the QSL phase provides a microcosm of all the key questions in quantum magnetism, arising where strong quantum spin fluctuations cause a partial or total suppression of magnetic order, whose extent can be controlled by an applied field.We believe that the combination of the three themes of our study, namely neutron spectroscopy in quantum materials, magnetic field-induced phenomena and MPS methods of accessing the complete spectral response of arbitrary locally interacting spin models, offers an exciting near-term future for quantum magnetism.

Experimental information
High-quality single crystals of CsYbSe 2 were prepared using the flux method 62 .Refinement of single-crystal X-ray diffraction data demonstrated the complete absence of Cs/Yb site mixing, as detailed in Supplementary Note 1A.For the characterization of our crystals, we measured the magnetization up to 60 T in pulsed magnetic fields at the National High Magnetic Field Laboratory (MagLab) and the specific heat in a dilution refrigerator at temperatures down to 0.05 K and magnetic fields up to 9 T (results shown in Supplementary Notes 1C and 1D).Our electron spin resonance measurements were performed using a continuous-wave ESR spectrometer, collecting data at X-band frequencies (ν = 9.4 GHz) and at T = 15 K.The resonance signal was measured from the field-derivative, dP/dB of the power, P , absorbed in a transverse microwave magnetic field and the spectra were fitted to a Lorentzian lineshape.
Approximately 200 single-crystalline pieces totalling around 0.5 g of material were co-aligned on copper plates to obtain a mosaic sample shown in Supplementary Figure 4. Our neutron scattering experiments were performed on the time-of-flight (ToF) Cold Neutron Chopper Spectrometer (CNCS) 63 at the Spallation Neutron Source at Oak Ridge National Laboratory (ORNL), the multiplexing Continuous Angle Multiple Energy Analysis spectrometer (CAMEA) 64 , and the cold-neutron Triple-Axis Spectrometer (TASP), the latter both located at the Swiss Spallation Neutron Source (SINQ) at the Paul Scherrer Institut (PSI).The measurements at CNCS were performed with an incident neutron energy E i = 3.32 meV, providing an energy resolution of 0.11 meV.A cryomagnet equipped with a dilution refrigerator was used to provide a maximum magnetic field of B = 5 T at temperatures down to 0.07 K. Measurements at CAMEA were performed with incident neutron energies E i = 5.2 and 6.2 meV (giving an energy resolution of 0.18 meV) and those at TASP with fixed k i = k f = 1.5 Å−1 , both using an 11 T cryomagnet reaching a base temperature of T ≃ 0.02 K.In all three experiments, the sample was orientated in the (HHL) scattering plane, such that the vertical magnetic field was applied along the [−1 1 0] direction in the ab plane.The software packages MANTID-PLOT 65 and HORACE 66 were employed for the data reduction and analysis at CNCS, while the data collected at CAMEA were analysed with the MJOLNIR software package 67 .

MPS calculations
We applied a cylinder MPS method to compute the dynamical spectral function of the isotropic spin-1/2 TLHAF in a magnetic field, as defined in Eq. ( 1), with J 2 /J 1 = 0.03.The MPS method proceeds by computing the time-dependent spin-spin correlation function where r is the site at which the initial spin operator is applied, x is the vector separation in the two-point correlator, and α, β ∈ {x, y, z}.The cylinder size, the bond dimension of the matrices used in the representation and the timeevolution procedures required to obtain well converged spectral functions at all fields are discussed and benchmarked in Supplementary Note 5.The calculations were implemented in Python using the package TENPY 68 .The dynamical spin spectral function was obtained from the Fourier transform The subscript r is retained because the correlation function is computed on a finite cylinder, with respect to the site r at its centre, and by time-evolving a ground state that breaks the translational symmetry of the Heisenberg model.To restore this symmetry in the spectral function, we average over three distinct time-evolved states, each corresponding to a site in the central unit cell, as explained in Supplementary Note 5.This procedure offers a strong reduction of the computational cost when compared with the MPS calculation of time evolution for a single spatially symmetric state.To account for artifacts in the spectral function caused by the finite cylinder length and time series in the Fourier transform, we convolve C αβ r (x, t) in Eq. ( 3) with a Gaussian envelope, as described in Supplementary Note 5.This results in an effective energy resolution of 0.1J 1 ≡ 0.038 meV and a momentum resolution of 0.032/a ≡ 0.006 r.l.u.
For comparison with the measured INS data, the calculated components of the dynamical structure factor were converted into a cross-section using the relation where F (Q) is the magnetic form factor of the Yb 3+ ion, is the neutron scattering polarization factor and g αβ specifies the components of the g-tensor determined by ESR.A detailed comparison of our INS and MPS results is shown in Supplementary Note 6.
One of the key questions in the analysis of every compound in the Yb-delafossite family is whether the ground state of the S = 1/2 system can be a quantum spin liquid (QSL) at zero applied magnetic field.In our work we have confirmed the presence of 120 • order in CsYbSe 2 , at temperatures below 0.4 K and at least over a spatial range far exceeding the lattice constant [Fig.1c of the main text and Supplementary Figure 6(e) below].The same order is found in KYbSe 2 [3], but notably not in NaYbSe 2 [4].Factors destabilizing this order include a possible frustrated coupling between triangular-lattice (TL) planes and a possible next-neighbour interaction within the planes, where (as noted in the main text) recent numerical studies have achieved partial agreement that the ground state is a QSL in the regime 0.06 ≲ J 2 /J 1 ≲ 0.15 [5-10].
Addressing first the issue of layer stacking, the majority of Yb delafossites have a structure described by the R3m space group, as listed in Supplementary Table 2, which has an ABC stacking that does indeed suggest frustration of antiferromagnetic (AF) interlayer couplings.However, some members of the series with larger alkali-metal ions, including Cs (Supplementary Table 2), have a structure with space group P6 3 /mmc that has AA stacking, and hence no interlayer frustration.Thus the available materials examples tend to suggest that the stacking of TL layers is not a relevant factor.Turning to J 2 , in our work we have deduced a weak next-Supplementary Table 1.Crystallographic data for CsYbSe2 determined by single-crystal X-ray diffraction.

Empirical Formula
CsYbSe2 Formula weight (g/mol) 463.87 neighbour interaction in CsYbSe 2 , J 2 = 0.03J 1 .The authors of Ref. [3] deduced a value J 2 = 0.05J 1 in KYbSe 2 , placing it closer to the QSL regime.While the authors of Ref. [4] were not able to determine a J 2 value, one may certainly suggest that the larger lattice constant inherent to the members with larger alkali-metal ions (Supplementary Table 2) causes a reduction of J 2 and hence an increasing stabilization of the 120 • -ordered ground state at base temperature in the Yb-delafossite materials.

C. Magnetization in Pulsed Magnetic Field
We measured the isothermal magnetization at T = 0.4 K in pulsed magnetic fields up to 60 T, and for calibration in a MPMS-7 with a 3 He insert at fields up to 7 T. We note that the magnetic field in all our experiments is applied in the ab plane, unless otherwise stated, and that we do not distinguish between in-plane directions.In the raw magnetization data shown by the black solid line in Supplementary Figure 2, we observe a clear and continuous increase above a saturation field of approximately 10 T. This contribution results from van Vleck paramagnetism and can be approximated in lowest order by a linear function (the dashed red line), from which the van Vleck susceptibility may be estimated as χ V V ≈ 0.0152 µ B /T. Subtracting this contribution from the raw data [12-14] leaves an intrinsic magnetization [solid blue line in Supplementary Figure 2, shown also in Fig. 1c of the main text] displaying near-perfect saturation above a field B Sat that we estimate most accurately from the excitation spectrum of the fully polarized phase [Fig.2b of the main text].In Supplementary Figure 2 we show also the bulk magnetization obtained from the field-induced enhancement of the (0 0 4) Bragg peak measured on TASP at the Paul Scherrer Institute (PSI), and these results agree quantitatively with the di- rect measurements.We defer the accurate modelling of these data, which we performed to obtain the orange line in Fig. 1c of the main text, to Supplementary Note 4.

D. Specific Heat
We measured the specific heat in a dilution refrigerator at temperatures down to 50 mK and magnetic fields up to 9 T. A sharp peak in C(T ) can be used to establish the presence of long-range order in the system.In Supplementary Figure 3(a) we group the applied fields, namely those below 3 T and above 6 T, in which no clear peak appears in C(T ) at all, only a broad hump at the low fields and a weak shoulder at the fields.In Supplementary Figure we show that relatively sharp, λ-shaped peaks can be found 3, 4 and 5 T, in partial correspondence with the long-range order of the UUD phase found by neutron diffraction.The temperatures of these sharp peaks are indicated as solid circles in the phase diagram shown in Fig. 1d of the main text and the temperatures of the features found at B ≤ 2 T and B ≥ 6 T as open circles.
In more detail, the specific heat at zero field in Supplementary Figure 3(a) shows a broad hump with no trace of a phase transition.This should be contrasted with the neutron diffraction measurements in Fig. 1b of the main text, which display a series of weak magnetic Bragg peaks with a clear onset at T ≃ 0.4 K.However, the analysis of the correlation lengths presented in Supplementary Note 2B leads to the result that these are finite at zero field, signalling a lack of true, longranged AF order even at T = 0.02 K, which is consistent with the absence of a sharp peak in the corresponding C(T ) data.The specific heat does show sharp, λ-shaped peaks in the field range 3-5 T [Supplementary Figure 3(b)] that reflect phase transitions into a long-range-ordered phase on, and apparently near, the 1/3-plateau state.These thermodynamic data are confirmed by the field-dependence of the (1/3, 1/3, 1) magnetic peak presented in Fig. 1c of the main text and analyzed in detail in Supplementary Figure 7 below.

Supplementary Note 2. NEUTRON SCATTERING EXPERIMENTS A. Sample Preparation
We coaligned around 200 single crystallites on copper plates, as shown in Supplementary Figure 4(a), to obtain a mosaic sample with a total mass of approximately 0.5 g.The rocking curve obtained by a neutron diffraction measurement of one selected magnetic peak [Supplementary Figure 4(b)] shows that the mosaicity of this sample is approximately 2 • FWHM.

B. Elastic Neutron Scattering at CNCS
Supplementary Figure 5 shows two-dimensional (2D) slices through the measured intensity dataset at zero energy in the (H H L) plane at fields from 0 to 5 T. At 0 T we observe weak but nonetheless clear intensity peaks at Q = (1/3, 1/3, L) for odd-integer L, which become weaker and almost invisible at 1 T.For B ≥ 2 T we find rods of magnetic intensity elongated along the [1/3 1/3 L] direction, with peaks at the same Q values.The intensities of these peaks signal effectively long-ranged order around 3-5 T, their very broad nature in L reflects the 2D nature of the magnetic system (as opposed to new peaks at intermediate L values) and their (H, K) position confirms the persistence of threefold periodicity in the plane of the triangular lattice at all applied fields.In the context of Supplementary Figure 5(a), we state for completeness that the data shown in the upper panel of Fig. 1b of the main text were integrated over the interval K = [−0.05,0.05], symmetrized according to the crystal symmetry and unfolded for visual clarity, meaning that the intensities at ±L are equivalent.
For a specific magnetic peak, the spin correlation length can be estimated using the formula ξ = 2π/ √ w 2 − R 2 [15], where w is the full width at half maximum (FWHM) height of the magnetic peak and R is the instrumental momentum resolution at this peak.To estimate this resolution, we first (a) To obtain the FWHM of the magnetic peaks for a quantitative determination of the correlation length, we fitted 1D cuts through the magnetic peaks to the function, which is defined as where the * denotes a convolution.In (1), is the Lorentz function, with peak centre x c and FWHM w L , and ) is a Gaussian with peak centre x = 0, unit peak area and FWHM w G , which should be fixed to the instrumental resolution, R(|Q|), for a specific Q.The FWHM of the Voigt function is then w V = 0.5346w L + 0.2166 • w 2 L + w 2 G .By inserting w V for w and R(|Q|) for R into the above expression for ξ, we estimate the correlation lengths at Q = (1/3, 1/3, 1) for each different field.Supplementary Figure 6(e) shows the correlation lengths we deduce from all of our CNCS data, which from our choice of the (0 0 L) series of nuclear peaks correspond to the in-and out-of-plane correlations.At B = 0 we obtain the values ξ ab = 60(7) Å and ξ c = 23(5) Å quoted in the main text.We observe that the correlation lengths dip at 1 T before ξ ab rises strongly, to values in excess of 200 Å, on and just above the 1/3 plateau.To gauge the meaning of this number, in Supplementary Figure 6(f) we show a 1D cut along [H H 1] at B = 4 T: the horizontal blue bar is the |Q|resolution at the magnetic Bragg peak (1/3, 1/3, 1), which is clearly resolution-limited, and thus the intrinsic ξ values diverge as expected for the long-range-ordered UUD state.

C. Elastic Neutron Scattering at TASP
The dependence of the integrated intensity of the magnetic peak (1/3, 1/3, 1) on the applied magnetic field is summarized in Fig. 1c of the main text and its dependence on temperature in Fig. 1b.Supplementary Figure 7(a) shows the raw elastic neutron scattering data measured on TASP in scans along the [H H 1] direction at base temperature for a series of fields.The peak intensity is clearly suppressed by a weak field, with a minimum around 1 T, before becoming sharper and as the field is increased, reaching a maximum around 4 T.This intensity then decreases beyond 5 T to small values higher the V Supplementary Figure 7(b) shows data obtained for the same scan at zero field for multiple temperatures below 1 K, where the robust peak present at T = 0.35 K clearly becomes sharper and stronger towards base temperature.

D. Two-Dimensional Excitation Character
The spin excitations we observe in our measurements retain their highly 2D nature over the full energy range and under all applied fields.In the representative 2D constant-energy slices presented in Supplementary Figure 8, the intensity distribution in all cases takes the form of multiple rods extending largely unchanged along the L direction.This definitive proof of 2D spin excitations is fully consistent with previous reports on CsYbSe 2 [16] and related delafossites [3, 4].As a consequence we may integrate our intensity data over a wide L range to analyse the in-plane spin dynamics, and for the data shown both below and in Fig. 3 of the main this integration range, 1.2 ≤ L ≤ 3.8, is denoted by the L label 2.5.

E. Background Definition and Subtraction for INS Spectra
The INS spectrum is always contaminated by a background contribution that arises from incoherent neutron scattering and from scattering due to the sample environment.Raw spectral data obtained on CNCS for the [H H 2.5] direction (meaning with the broad L integration described above) are shown in the left panels of Supplementary Figure 9(a,c,d) for three different magnetic fields at our base temperature of 70 mK.For an accurate characterization of the background, we make use of the fact that the 5 T spectrum at the Γ point, I(Q = 0, E), consists of a single and well-defined inelastic peak near 1 meV.We assume that an appropriate 1D intensity cut, prepared from the shaded region in Supplementary Figure 9(a) and shown by the black points in Supplementary Figure 9(b), can be described within the energy window 0.5 meV ≤ E ≤ 1.5 meV by a Lorentzian peak and a linear background.Thus we fitted the measured signal to the form where I 0 , W 0 and E 0 characterize respectively the intensity, width and centre of the inelastic peak.We subtract the fitted Lorentzian, shown by the blue line in Supplementary Figure 9(b), from the raw spectrum to obtain a residual intensity that we consider as fully representative of the background [shown by the red points in Supplementary Figure 9(b)].The results of subtracting this background are shown in the right panel of Supplementary Figure 9(a), where clear limits to the extent of each continuum become visible.Because our data provide no evidence that the background varies with Q, we subtracted the same form from our cuts everywhere in reciprocal space.

Supplementary Note 3. LINEAR SPIN-WAVE THEORY
In our analysis we use linear spin-wave theory (SWT) to achieve two separate goals.Working above the saturation field, B Sat , we make use of the fact that the excitations are guaranteed to be well defined spin waves to obtain the most accurate available fit of the magnetic parameters of CsYbSe 2 (Supplementary Note 3A).Working below B Sat , we use linear SWT for a preliminary indication of the locations and energy scales of putative ∆S = 1 excitations in the TL-HAF, and hence of departures from semiclassical magnetism arising due to quantum corrections.In Supplementary Note 3B we outline how the SWT results shown in Fig. 3 of the main text were obtained and use these to obtain an overview of possible two-magnon contributions.

A. Fitting of Magnetic Interactions
As described in the main text, we determine the magnetic interactions by using the 11 T dataset from CAMEA, which shows a single, sharp magnon mode with its maximum at the Γ point [Fig.2b of the main text].We quantified the location of this mode at 17 points in reciprocal space and used SPINW [17] to fit these to the Hamiltonian The optimal fit in the 3D parameter space of J 1 , ∆ and J 2 is J 1 = 0.395 (7) meV, ∆J 1 = 0.39(2) meV and J 2 = 0.011(4) meV.In Supplementary Figure 10 we illustrate the quality of this fit by showing two 2D cross-sections: fixing J 2 to its optimal value [Supplementary Figure 10(a)] shows that J 1 = ∆J 1 , i.e. the nearest-neighbour interaction is isotropic within the precision of our measurements; set-ting ∆ = 1 [Supplementary Figure 10(b)] allows us to optimize J 2 , and hence to conclude that magnetic Hamiltonian of CsYbSe 2 is described by the J 1 -J 2 Heisenberg model with J 1 = 0.395(8) meV and J 2 = 0.011(4) meV (J 2 /J 1 = 0.03).

B. SWT Spectra
Linear SWT ceases to provide a complete description of the spectrum of the TLHAF as soon as the field is lowered below B Sat , when quantum corrections become finite.The TL-HAF in an in-plane magnetic field has long been known [18]  to exhibit a deformed 120 • phase (Y), the 1/3 plateau phase (UUD) and a V (or 2:1) phase, with two of the spins parallel, below saturation, as illustrated schematically in Fig. 1d of the main text.At the mean-field level, the magnetization is perfectly linear in field from 0 to B Sat = 9.6 T (Fig. Supplementary Figure 12), which is obtained exactly due to the classical nature of the problem at B ≥ B Sat .However, M (B) has no plateau around B Sat /3, as this phase is stabilized only at higher order in 1/S.In linear SWT, there is a large family of classically degenerate states at any field B < B Sat , and the precise orientation of the Y, UUD or V state at any given field is set manually in our analysis.The spin excitations around this fixed ground state were then calculated using SPINW and the resulting spectra at all selected fields in our experimental range (B = 0-11 T) are shown as the orange lines in Fig. 3 of the main text.
For completeness, in Supplementary Figure 11 we show the densities of two-magnon states computed using the linear SWT one-magnon branches at fields of 0, 2, 3, 4, 5 and 8 T. The full two-magnon spectral functions are reweighted versions of these densities of states, which for the Y and V phases are more complex to compute and to normalize to the one-magnon branches.The density of states therefore serves as a useful semi-quantitative guide to the location of the two-magnon continuum in wavevector and energy, and to its qualitative shape.The densities of states in Supplementary Figure 11 each reflect the nine different continua arising from the three one-magnon branches, whose overlap results in the edge structures that appear throughout the Brillouin zone.The primary difference between linear SWT and our INS and MPS results is clearly the absence of well defined one-magnon branches extending over most of the zone at all fields outside the UUD phase.While the two-magnon continua of linear SWT occupy the entire zone over a broad energy window centred at 1 meV, they are in general too flat (except at 8 T), too uniform and have too many edges to bear a close resemblance to the continuum features in our measured and calculated spectral functions.plementary Note 3A) and perform a grand canonical DMRG calculation of the full magnetization curves that allows us to deduce the boundaries of the 1/3-magnetization plateau.
For a state-of-the-art numerical determination of the magnetization response of the S = 1/2 TLHAF, we apply the grand canonical DMRG method [19].This technique computes the infinitesimally small magnetization response to a change in the applied field, and thus the physical quantities it provides mimic the thermodynamic limit to approximately 1 part in 10 3 for a 2D system.The method is based on a graded division of a finite-sized cluster into a centre region and edge regions, such that the centre reproduces the continuous bulk response by using the near-zero-energy states of the edge as a type of buffer.In a system of fixed size and shape, the grading is introduced by modulating the energy scale with the externally imposed function which deforms the Hamiltonian smoothly from its standard energy at the centre of the system (r = 0) to zero at the open cluster edges (r = R).After obtaining the lowest eigenwavefunction of the deformed Hamiltonian, the magnetization can be read as M = ⟨S z (r = 0)⟩, because this wavefunction optimizes the expectation value to its thermodynamic limit at any given magnetic field.To mimic the results of a bulk measurement, appropriately weighted for the different possible directions of the magnetic field with respect to the triangular lattice, we performed our calculations using the hexagonal 75-site cluster shown in the inset of Supplementary Figure 12.
The calculations were performed for the J 1 -J 2 TLHAF using the interaction parameters deduced in Supplementary Note 3A and an in-plane g-factor of g ab = 3.2, as deduced from ESR and INS in the main text.The result shown in Supplementary Figure 12 is that reproduced in Fig. 1c of the main text, and was used for our extraction of the lower and upper boundaries (B l and B u ) of the 1/3 plateau.gle symmetry-broken ground state (Supplementary Figure 13) and we take the equally weighted average to obtain the symmetric spectral function This procedure is repeated for all αβ ∈ {zz, +−, −+}, meaning that in total nine time-evolution runs are required at every value of the magnetic field, B.
In addition, symmetries of the Heisenberg Hamiltonian allow us to obtain the negative-time correlation functions at no additional computational cost as the complex conjugate of the positive-time ones, i.e.C αβ r (x, −t) = C αβ r (x, t).Putting these together, the symmetric spectral function of Eq. ( 7) simplifies to the form which results in an effective broadening of the spectral function.For this we used σ t = 0.005J 2 1 and σ x = 0.02/a 2 , where a is the unit-cell size of the TL and e L the unit vector along the cylinder axis, leading to the effective resolution in energy and momentum given in the Methods section of the main text.
In Supplementary Figure 14 we show spectra obtained using this procedure in the UUD (1/3-plateau) phase, at a magnetic field equivalent to B = 3 T, to illustrate the degree to which our calculations have converged to the spectral function of the infinite system.We observe that increasing the cylinder circumference from C = 6 to 8, the length from L = 30 to 60 or the bond dimension of the MPS from χ = 256 to 512 make only minor changes to the quality of the spectra at the values of σ t and σ x chosen to match the experimental resolution.We comment that higher-resolution experimental data would require higher C, L or χ values to achieve convergence, and that our present MPS studies allow a factor-2 reduction in σ t where h = µ B g ab B. At m = 1/3, the ground state in the Ising limit, J z /J xy ≫ 1, is the UUD product state illustrated in Fig. 4c of the main text.An excitation at the same S z consists at least of one pair of local spin-flips, U → D and D → U .If the two flipped spins are spatially well separated, they cost an energy 3J z in the Ising limit, but if they are nearest neighbours then they cost only 2J z .We note that, because the total S z is unchanged by two opposing spin flips, the z-axis magnetic field has no effect on the energies of these two-spin excitations.
As the next step we consider the perturbative limit with a finite, but small, J xy .The nearest-neighbour bound states now show a spatially localized yet partially mobile structure.A schematic representation is shown in Fig. 4c of the main text, where the flipped red spin is confined, but can hop around the flipped blue spin along a finite chain of six sites, i.e. on the hexagon around the blue spin.The discrete energy spectrum of this six-site chain shows a lifting of the bound-state degeneracy to six flat bands with energies {2J z + J xy , 2J z + J xy /2, 2J z − J xy /2, 2J z − J xy }, with the second and third energies each twofold degenerate.The manifold of excitations arising from pairs of flipped spins far away from each other is also split at first order in J xy , and forms a two-particle continuum at an energy around 3J whose support can be derived from the dispersion of the single spin-flips.
The next step is to determine whether this structure of bound states is visible in the longitudinal dynamical structure factor, S zz (Q, ω).For this we performed MPS calculations for the nearest-neighbour model at different values of J z /J xy , and in Supplementary Figure 16(a) we show the spectral function at J z /J xy = 100 along the usual Γ-K-M path in the Brillouin zone.Indeed one observes that the spectral weight becomes Q-dependent, and that the excitation energies are confined to the perturbative branches (horizontal dashed red lines) discussed above.As J z /J xy is lowered to 5 [Supplementary Figure 16(b)], 2 [Supplementary Figure 16(c)] and the Heisenberg point [Supplementary Figure 16(d)], the boundstate branches acquire an increasing dispersion and broaden in energy; in particular, while the lower bound state remains relatively sharp as J z /J xy → 1, the upper scattering states broaden into the bow-tie structure of continuum IV.
To track the origin of this broadening, in Supplementary Figure 16(e-h) we show four spectral functions on a logarithmic intensity scale.In Supplementary Figure 16(e), which corresponds to Supplementary Figure 16(b), the two-magnon scattering continuum begins to become visible, centred at an energy around 3J z (15J xy ).As J z /J xy is lowered, this continuum rises in intensity and overlaps increasingly with the energy window of the bound states, until at the Heisenberg point [Supplementary Figure 16(h)] only a part of the lowest bound-state branch remains as the isolated and well defined mode observed in experiment (Fig. 4a of the main text).

FIG. 1 .
FIG. 1. Structure, properties, phase diagram and finite-field spectra of CsYbSe2.a Crystal structure of CsYbSe2 and representation of the ideal Yb 3+ TL layer.The red arrows represent the ordered spins of the weak 120 • AF order at zero field.b Upper panel: elastic scattering intensity in the (HHL) plane for 0 T, with the data presented as described in Supplementary Note 2B.Red arrows indicate weak magnetic intensity peaks at Q = (1/3, 1/3, ±L) with L = 1, 3, 5. Lower panel: temperature-dependence of the (1/3, 1/3, 1) peak area at zero field.The red solid line is a fit to an order-parameter form.The grey shaded area represents the approximate sensitivity limit of our measurement.c Upper panel: isothermal magnetization (blue symbols) measured at T = 0.4 K as a function of magnetic field applied in the ab plane and with the van Vleck contribution subtracted as described in Supplementary Note 1C.The solid orange line shows a grand canonical density-matrix renormalization-group (DMRG) calculation of the magnetization performed for the TLHAF using parameters deduced in Fig. 2, from which we determined the saturation field, BSat = 9.6(2) T (vertical solid line), and the lower and upper boundaries of the 1/3 plateau as B l = 2.95(6) T and Bu = 4.5(1) T (vertical dashed lines).Lower panel: integrated intensity of the (1/3, 1/3, 1) magnetic Bragg peak measured at T < 0.05 K. d Phase diagram of CsYbSe2.The blue diamond is the phase-transition temperature obtained from neutron diffraction at zero field.The two open squares indicate the region (3 T ≤ B ≤ 4.5 T) where the peak intensity in the lower half of panel c remains almost unchanged.The solid and open circles represent respectively the temperatures of sharp peaks and broad humps in the corresponding specific-heat curves, as described in Supplementary Note 1D.The arrows indicate schematically the spin order of the five phases of the classical TLHAF, which are consistent with the magnetic peaks we observe.e INS spectrum measured under a magnetic field of 5 T, showing an absence of well defined ∆S = 1 excitations away from the Γ and K points but extensive continuum features (Fig.3).The grey horizontal bar at low energy masks the elastic-line contribution.f Comparison with a matrix-product-state (MPS) calculation of the spectrum at the same field.

FIG. 2 .
FIG. 2. INS and ESR determination of the spin Hamiltonian.a Representation of the Brillouin zone in the (H K 0) plane and definition of directions x ∥ [1 1 0], ŷ ∥ [0 0 1] and ẑ ∥ B, which is orthogonal to x and ŷ. b Spin excitations along the [H H 3] direction measured on CAMEA at T = 0.05 K in the field-polarized regime (B = 11 T).The solid line shows the dispersion calculated using linear spin-wave theory (SWT) with J1 = 0.395(8) meV, J2 = 0.011(4) meV and g = 3.2.The grey horizontal bar at low energy masks the elastic-line contribution.c Angular dependence of the g-factor measured by ESR at T = 15 K.The solid line indicates the form g(θ) = g 2 ab sin 2 θ + g 2 c cos 2 θ, which allows the extraction of the strongly anisotropic in-and out-of-plane coefficients g ab = 3.25 and gc = 0.3.Inset: representative ESR spectrum with a Lorentzian fit shown by the dashed red line.d Field-dependence of the INS signal at the Γ point; data for B ≤ 5 T were measured on CNCS and data for B = 8 and 11 T on CAMEA.Black and red points show the respective positions of the magnon mode as extracted from the CNCS and CAMEA datasets, for both the Γ and K points.The solid line shows a linear fit that yields g ab = 3.20(6).

FIG. 3 .
FIG. 3. Complete field-induced spectral response of CsYbSe2.Panels a, c, e, g, i and k show spin excitation spectra measured under different magnetic fields at T = 0.07 K.The open circles, squares, triangles and diamonds indicate respectively excitation features in the categories I, II, III and IV described in the text.All data have been symmetrized according to the crystal symmetry.The orthogonal in-plane integration range along the [−K K 0] direction is K = [−0.05,0.05] and the out-of-plane range is 1.2 ≤ L ≤ 3.8 for our CNCS data (0-5 T) and 2 ≤ L ≤ 4 on CAMEA (8 T).The background subtraction is described in Supplementary Note 2E.The narrow horizontal grey regions mask the elastic line.Panels b, d, f, h, j and l present dynamical spin structure factors obtained at different magnetic fields by cylinder MPS calculations (Supplementary Note 5).Colour bars represent both the measured and calculated intensities in a single set of arbitrary units (i.e. the same units are used at all fields).Orange lines show the mode positions and intensities given by linear SWT with the same interaction parameters.The open points are identical to those shown in panels a, c, e, g, i and k.
FIG. 4. Longitudinal spin excitations.a Longitudinal component of the excitation spectrum extracted for the [H H 0] direction at B = 4 T. The orthogonal in-plane integration range is K = [−0.05,0.05] and the out-of-plane range is L = [−0.5,0.5].The narrow grey region masks the elastic line.b Corresponding MPS calculation of the longitudinal component, Szz(Q, ω), of the dynamical structure factor.c Schematic representation of spin-flip processes in the UUD phase: red and blue circles represent respectively U and D spins, the grey dashes highlight flipped spins (U → D and D → U) and the green circle delineates the hexagon on which the blue flipped spin may propagate at no energy cost in the Ising limit, ensuring its localization.d MPS calculation of Szz(Q, ω) for the strongly Ising-type parameter choice Jz = 5Jxy (Supplementary Note 7); the mode energies are shown in units of Jxy to illustrate the role of this interaction in setting the splitting and dispersion of the bound states centred at 2J1 = 2Jz.

Supplementary Figure 3 .
Temperature-dependence of the specific heat measured in a range of applied in-plane magnetic fields (B ⊥ ĉ).
Supplementary Figure 4. (a) Photographs of the CsYbSe2 crystals coaligned and assembled for neutron scattering experiments.(b) Rocking curve of the coaligned sample measured by neutron diffraction, with a Gaussian fit shown by the solid line.

Supplementary Figure 5 .
Constant-energy slices in the (H H L) plane taken at zero energy transfer under different magnetic fields, with the data symmetrized about the [H H 0] axis.The strong spots at Q = (0, 0, L) for even-integer L are nuclear Bragg peaks.

Supplementary Figure 6 .
(a,b) 1D intensity cuts at zero energy along the [0 0 L] and [H H 2] directions at the (0, 0, 2) nuclear Bragg peak at 0 T. (c) FWHM of nuclear peaks obtained from cuts along the [H H 0] direction (black) and [0 0 L] direction (red) at Q = (0, 0, 2), (0, 0, 4) and (0, 0, 6).Solid lines in both panels show polynomial fits.(d) Perpendicular and parallel |Q|-resolution estimated from the FWHMs shown in panel (c).(e) Perpendicular (ξ ⊥ ) and parallel (ξ ∥ ) correlation lengths calculated from the magnetic Bragg peaks shown in Supplementary Figure5.For this geometry ξ ⊥ ≡ ξ ab , the in-plane correlation length, and ξ ∥ ≡ ξc, the out-of-plane correlation length.(f) An example 1D cut along the [H H 1] direction through the magnetic Bragg peak (1/3, 1/3, 1) at 4 T. The horizontal blue bar represents the instrumental |Q|-resolution, which indicates a resolution-limited peak at 4 T. prepared 1D intensity cuts along the [H H 0] and [0 0 L] directions at the nuclear peaks (0, 0, 2) [Supplementary Figure6(a,b)], (0, 0, 4) and (0, 0, 6), then fitted a Gaussian function to them.The FWHM values of these three nuclear peaks reflect the momentum resolution of the instrument perpendicular (for [H H 0]) and parallel (for [0 0 L]) to Q at each reciprocallattice point.Polynomial fits to the FWHM as a function of L [Supplementary Figure6(c)] yield approximate estimates of the resolution as a function of |Q| for an arbitrarily cho-sen momentum transfer, Q, in the perpendicular and parallel directions [Supplementary Figure6(d)].

Supplementary Figure 7 .
(a) Intensity of the (1/3, 1/3, 1) magnetic Bragg peak measured at T < 0.05 K for a series of magnetic fields, which are offset by 50 units for clarity.Solid lines and shading show Gaussian fits.(b) Intensity of this peak measured at zero field for a series of temperatures below 1 K, which are offset by 20 units for clarity.

Supplementary
Figure 10.(a) Fit quality shown as a function of assumed XXZ-model parameters J1 and ∆ for the dominant nearestneighbour interaction.Here J2 has been fixed to its optimal value of 0.011 meV.(b) Fit quality shown as a function of J1 and J2, assuming ∆ = 1.Red crosses in both panels indicate the position of the global minimum.

11 .
Two-magnon densities of states computed from linear SWT at six different magnetic fields for the J1-J2 TLHAF with J2/J1 = 0.03.L denotes integration over the full range of L, equivalent to an ideally two-dimensional system.

Supplementary Figure 14 .
Comparison of spectral functions computed using the MPS method for different system sizes and bond dimensions at the fixed values of the Gaussian envelope parameters stated in the text.Here we show the quantity Sxx(Q, ω) + Szz(Q, ω) at a magnetic field corresponding to the experimental value B = 3 T.The columns show increasing values of the bond dimension, χ, while the rows show cylinders of different sizes, as labelled in the panels in the first column.

e
−i x•Q cos(ωt) Re C αβ r (x, t)−sin(ωt) Im C αβ r (x, t), where the outer sum runs over the three centre sites.To compensate for the finite cylinder length and time-step series, it is standard before Fourier transformation to convolve the correlation function with a Gaussian envelope, C αβ r (x, t) → e −σtt 2 e −σx(e L •x) 2 C αβ r (x, t),