Field-induced bound-state condensation and spin-nematic phase in SrCu2(BO3)2 revealed by neutron scattering up to 25.9 T

In quantum magnetic materials, ordered phases induced by an applied magnetic field can be described as the Bose-Einstein condensation (BEC) of magnon excitations. In the strongly frustrated system SrCu2(BO3)2, no clear magnon BEC could be observed, pointing to an alternative mechanism, but the high fields required to probe this physics have remained a barrier to detailed investigation. Here we exploit the first purpose-built high-field neutron scattering facility to measure the spin excitations of SrCu2(BO3)2 up to 25.9 T and use cylinder matrix-product-states (MPS) calculations to reproduce the experimental spectra with high accuracy. Multiple unconventional features point to a condensation of S = 2 bound states into a spin-nematic phase, including the gradients of the one-magnon branches and the persistence of a one-magnon spin gap. This gap reflects a direct analogy with superconductivity, suggesting that the spin-nematic phase in SrCu2(BO3)2 is best understood as a condensate of bosonic Cooper pairs.

Bose-Einstein condensation (BEC) underpins exotic forms of order ranging from superconductivity to superfluid 4 He.In quantum magnetic materials, ordered phases induced by an applied magnetic field can be described as the BEC of magnon excitations.With sufficiently strong magnetic frustration, exemplified by the system SrCu 2 (BO 3 ) 2 , no clear magnon BEC is observed and the complex spectrum of multimagnon bound states may allow a different type of condensation, but the high fields required to probe this physics have remained a barrier to detailed investigation.Here we exploit the first purpose-built high-field neutron scattering facility to measure the spin excitations of SrCu 2 (BO 3 ) 2 up to 25.9 T and use cylinder matrix-product-states (MPS) calculations to reproduce the experimental spectra with high accuracy.Multiple unconventional features point to a condensation of S = 2 bound states into a spin-nematic phase, including the gradients of the one-magnon branches, the presence of many novel composite two-and three-triplon excitations and the persistence of a one-magnon spin gap.This gap reflects a direct analogy with superconductivity, suggesting that the spin-nematic phase in SrCu 2 (BO 3 ) 2 is best understood as a condensate of bosonic Cooper pairs.Our results underline the wealth of unconventional states yet to be found in frustrated quantum magnetic materials under extreme conditions.
Condensation of a macroscopic number of particles into a single state is a purely quantum mechanical phe-nomenon exhibited by a wide spectrum of bosonic systems [1].Early and spectacular examples of this Bose-Einstein condensation (BEC) include superfluidity in liquid 4 He and superconductivity in metals, where attractive interactions allow the electrons, which are fermions, to form Cooper pairs, which are bosons.Field-induced gap closure and magnetic order in quantum disordered magnetic materials is conventionally described as a BEC of magnons, which are spin-1 excitations [2,3]; examples of this phenomenon include the dimerized quantum antiferromagnets TlCuCl 3 [4] and BaCuSi 2 O 6 [5,6], and, under quite different experimental conditions, yttrium iron garnet [7,8].
By contrast, the field-induced phase diagram of the highly frustrated quantum magnet SrCu 2 (BO 3 ) 2 (Fig. 1a) [12] shows neither clear gap closure nor induced magnetic order, but instead a spectacular series of magnetization plateaux [12][13][14][15][16]. SrCu 2 (BO 3 ) 2 presents a remarkably faithful realization of the Shastry-Sutherland model (SSM), which is a paradigm for ideal frustration in a two-dimensional (2D) spin system [17].This system consists of S = 1/2 dimer units arranged orthogonally on a square lattice (Fig. 1b), and its properties are governed by the ratio between the intra-dimer interaction, J, and inter-dimer one, J ′ .For J ′ /J < 0.675, the ground state is an exact dimer-product singlet state [18], and the ideally frustrated geometry leads to many unconventional effects, including near-dispersionless one-magnon excitations (to which we refer henceforth as "triplons") and very strongly bound multi-triplon states [19][20][21][22].At ambient pressure and zero field, SrCu 2 (BO 3 ) 2 is very well described by a SSM with J ′ /J ≃ 0.63 [23] and ).c-d Energy levels of a simple J-J ′ plaquette system (inset in panel c) shown as a function of magnetic field.For J ′ /J > 0 (c), the t+ one-triplon branch condenses first, whereas when J ′ /J < 0 (d), the two-triplon bound state condenses first.The two scenarios result in different gradients after condensation.e Excitation spectra measured by inelastic neutron scattering (INS) over a wide range of fields.Data at low fields (4, 8, 9.5 and 12 T) were measured at Q = (1.5, 0.5, 0) as described in the Methods section.Data from our HFM/EXED measurements (0, 15, 18, 20, 23 and 25.9 T) are presented with full momentum integration, as described in Sec.S1 of the SI [9].Colour contours represent the neutron intensity and symbols show the fitted positions of individual excitations (Figs.2g-l).Solid grey lines show the projected locations of the one-triplon branches for a scenario of one-triplon condensation and no DM interactions (panel c), dashed black lines their locations for two-triplon condensation (panel d).The dashed green line shows the Sz = 0 branch of the S = 1 multiplet within the two-triplon bound state.The dotted red line up to 15 T shows the Sz = 2 branch of the S = 2 multiplet as measured by electron spin resonance (ESR) [10], i.e. at Q = (0, 0), where the dispersion of this branch reaches its maximum [11].The red shading represents the estimated bandwidth of this branch, whose minimum determines the binding energy, and our linear extrapolation beyond 15 T indicates its effect on gap closure.f Total momentumintegrated spectral functions obtained from our cylinder MPS calculations on the SSM with additional DM interactions over the same wide field range.
shows anomalous thermodynamic behaviour [24,25] due to the high spectral density of highly localized spin excitations.Under an applied hydrostatic pressure [25][26][27], its behaviour mirrors the phases of the SSM as J ′ /J is changed.
Most remarkable of all is the appearance of fieldinduced magnetization plateaux in SrCu 2 (BO 3 ) 2 : ex-tremely strong magnetic fields applied along the c-axis (Fig. 1a) induce plateaux at 1/8, 2/15, 1/6, 1/4, 1/3, 2/5 and 1/2 of the saturation magnetization, which is reached around 139 T [16].In direct analogy with charge-densitywave (CDW) order in electronic systems, these plateaux can be seen as bosonic CDW phases of magnetic entities.The predicted spin superstructure on the lowest (1/8) plateau is a crystal of S = 2 two-triplon bound states, as opposed to a crystal of unbound S = 1 triplons [28]; this prediction awaits experimental verification.Although the transition to the 1/8 plateau takes place with a jump at µ 0 H c = 27 T, the magnetization is clearly finite above approximately 16 T, and above 21 T it increases linearly until the jump [14].In this crossover region, the gap to the lowest triplon branch does not close as the field is increased, but shows an avoided crossing previously ascribed to the presence of weak Dzyaloshinskii-Moriya (DM) interactions [10,29,30] It was suggested very early that the field-induced gap closure in SrCu 2 (BO 3 ) 2 could be anomalous [31].The zero-field gap of the nearly immobile triplons is ∆ ≃ 3.0 meV, but triplon pairs gain kinetic energy from correlated hopping, such that two-triplon states with S = 0 and 1 lie far below 2∆ and hence are strongly bound (with respective energies 4.0 [32] and 4.8 meV [19]).For the S = 2 multiplet, perturbative (in J ′ /J) calculations on the SSM indicate that its minimum, reached at momentum (π, π), also lies below 2∆ [11].In this situation, as the minimal model of Figs.1c-d makes clear, increasing the field should cause the S = 2 bound state to meet the singlet ground state before the lowest one-triplon branch does.The result would be a BEC of two-triplon bound states, which is a spin-nematic phase [33,34], instead of field-induced magnetic order.However, a conclusive demonstration has not to date been possible: standard momentum-resolved experimental probes do not observe S = 2 excitations, light-scattering methods probe them only at momentum Q = (0, 0) (depicted in Figs 1e) and the DM interactions obscure the situation.
The definitive experiment would be to determine the dynamical structure factor at all Q by inelastic neutron scattering (INS), and to compare this with model calculations.Although INS has historically been limited to fields below 15 T, for a limited period neutron scattering was possible up to 25.9 T (DC) on a specialized facility at the Helmholtz-Zentrum Berlin (HZB) [35][36][37], thus enabling the investigation of SrCu 2 (BO 3 ) 2 over most of the formerly "dark" field regime below H c .Model calculations for frustrated 2D quantum systems have also been a historical impossibility, but recent progress in numerical methods based on matrix-product states (MPS) has made it possible to calculate hitherto unavailable spectral functions with quantitative accuracy.These parallel developments render high-field INS the ideal tool to decode the spin dynamics in SrCu 2 (BO 3 ) 2 and in a range of other systems [6].
Here we harness this progress to demonstrate bound-state condensation forming a spin nematic in SrCu 2 (BO 3 ) 2 .We measure the excitation spectrum using INS at 200 mK with magnetic field strengths up to 25.9 T (Fig. 1e), and we perform MPS calculations of the dynamical structure factor of the SSM with weak DM interactions (Fig. 1f).By tracking the Zeeman splitting of the one-triplon mode, we find field-induced changes of the gradients of the triplon branches and a persisting, DM-independent gap that are fully consistent with bound-state condensation.Our INS and MPS spectra also reveal a rapid, field-induced transfer of weight to multiple novel excitations, which we show are two-and three-triplon composites that can only exist if the condensate contains bound states.

Results
Inelastic neutron scattering.We performed INS experiments at the High-Field Magnet (HFM) facility that operated recently at the HZB.The time-of-flight EXED instrument offered the unique capability of performing neutron diffraction and spectroscopy at fields up to 25.9 T. A further key was a dilution refrigerator enabling a sample temperature of 200 mK, well below the smallest energy scales in SrCu 2 (BO 3 ) 2 .Unless otherwise stated, all the INS data we show were collected with incoming energy E i = 8 meV, yielding an energy resolution of 0.64(1) meV.Further information concerning our experiments is summarized in the Methods section and in Sec.S1 of the Supplementary Information (SI) [9] we provide a detailed description of the HFM/EXED experimental geometry, our data preprocessing and background subtraction.
Figure 1e shows the fully Q-integrated neutron intensity as a function of applied field and energy transfer.In Figs.2a-f we show the neutron intensity measured as a function of Q h , with partial integration around Q k = 0.All the excitations we observe are rather flat over the full reciprocal space, and this justifies integrating the intensity data over a large Q volume, which improves the counting statistics significantly.However, a weak dispersion is visible in most modes (Figs.2a-f), which on integration appears as an additional contribution to the linewidth.The intensity in every branch remains concentrated at the largest values of Q h for all fields, while both the intensities and the relative intensities between branches show weak changes with Q h , which we analyse in Sec.1D of the SI [9].
Integrating over the entire available Q range leads to the results shown in Figs.2g-l.At zero field (Fig. 2g), the one-triplon excitation is found around 3.0 meV, with a number of higher-lying excitations at 4.8, 5.5 and 6.5 meV, all in good agreement with previous INS measurements [19][20][21][22].In an applied field, the one-triplon excitation shows Zeeman splitting into three branches, which are clearly visible at 15 T (Fig. 2h).We label the bottom, middle and top branches respectively as t + , t 0 and t − .Their splitting increases linearly with the field until approximately 18 T (Fig. 2i), as in ESR measurements [10].Above this field, the t + branch curves towards a rather constant value around 0.5 meV, suggesting an avoided crossing, with no indication that the one-triplon gap closes.In addition, we observe the de- A single background was subtracted as detailed in Sec.S1B of the SI [9].The one-triplon excitation (solid blue lines with blue shading) exhibits Zeeman-type splitting in the applied field and its branches were fitted by related Gaussians (Sec.S1C of the SI [9]).The Sz = 0 branch of the S = 1 two-triplon bound state (dark turquoise lines and shading) appears at a constant energy.Additional intensity develops below the t0 one-triplon branch at µ0H ≥ 18 T (solid red lines with red shading) and was fitted with a single Gaussian.Further additional intensity visible below the t− one-triplon branch at µ0H ≥ 15 T (dotted black lines with grey shading) was fitted with multiple Gaussians.
velopment of a shoulder on the low-energy side of the t 0 triplon, which becomes progressively stronger as the field is raised (Figs.2i-l).Around the t − triplon we observe further intensity contributions over a broad energy range, which become stronger at high fields and can be fitted reasonably by three separate Gaussians.Within our experimental resolution, the energies of these new excitations remain largely constant as the magnetic field is increased (Figs. 1e and 2i-l).As Figs. 1c-d showed, the field-induced evolution of the lowest (t + ) triplon mode is important for identifying the nature of condensation.The fact that its gap remains finite throughout the crossover regime improves the prospects for resolving its position close to the elastic line.We performed measurements with an incoming energy of E i = 4 meV, at which the resolution improves to 0.24(1) meV.In Fig. 3a we show the neutron spectra measured at 20 and 25.9 T, where a coherent low-energy feature remains discernible, and we extract the respective positions 0.54 and 0.35 meV (also marked in Fig. 1e), re- inforcing our result that the gap does not close at any field.
To analyse our intensity data in more detail, in Fig. 3b we show the spectral weights accumulated at energy transfers above a 0.7 meV cut-off for each measurement field.At zero field, this weight is zero until the 3.0 meV triplet branch is encountered, after which it is constant again until the S = 1 branch of the two-triplon bound state at 4.8 meV, and beyond this it increases linearly as multiple higher-lying states are encountered.At 15 T, the one-triplon sector is split into three clearly defined branches, each of which gives a significant jump in the accumulated spectral weight.At higher fields, these jumps become increasingly broad as additional states appear at multiple different energies, particularly those directly below the t 0 and t − branches (Figs.1e and Figs. 2i-l).We defer an interpretation of these results (Fig. 3c) to our numerical analysis below.To investigate the physics of the different spectral-weight contributions identified in Figs.2g-l, in Fig. 3d we inspect their evolution as a function of the field.It is clear that the intensity of the additional states (i.e.those below t 0 and below t − ) scales linearly beyond the crossover region, meaning that it follows the magnetization, m, while the one-triplon weight is proportional to 1 − m.Numerical modelling.To model the experimental results, we have performed cylinder MPS calculations of the dynamical spectral function for the SSM supplemented by the weak (3%) DM interactions deduced from experiment [30].To interpret the broadest features of the spectrum we appeal to the plaquette model of Figs.1c-d and for a detailed identification of the excitation types we perform exact diagonalization (ED) calculations for the SSM on clusters of 4×4 spins.A summary of our MPS calculations is provided in the Methods section, and in Sec.S2A of the SI [9] we present additional details including the effects of the MPS cylinder size.MPS results corresponding to the INS spectrum are shown in Fig. 1f, which makes clear that the calculations contain all the features found in the measurements, most notably the additional excitations appearing below the t 0 and t − one-triplon branches beyond 18 T, at energy transfers in quantitative agreement with experiment.
The first striking feature of our MPS results is the minimum in the t + branch at 21 T, with a finite gap of 0.6 meV, beyond which the mode energy increases again with the field.These low energies are difficult to access by INS because of the elastic peak (Fig. 3a), but our MPS calculations confirm the absence of low-energy spectral weight here.The persistence of a gap is not a generic feature of field-induced magnon BEC [2,3], where the gap closes, the lowest magnon condenses and spectral weight appears down to zero energy at all higher fields.This spectral weight should extend up to the magnon bandwidth, which in SrCu 2 (BO 3 ) 2 is only 0.3 meV [38], and thus a minimum around 0.6 meV suggests different physics.
It was assumed previously that this avoided closure is a consequence of the DM interactions in SrCu 2 (BO 3 ) 2 , and indeed DM terms are known to preserve a gap in systems undergoing magnon BEC [39].However, the minimum  should then appear at the field where one-magnon condensation occurs, which by extrapolating both our INS and MPS results (Figs. 1e-f) is clearly 24 T rather than the value µ 0 H * = 21 T that we observe.To confirm that this gap is not caused by the DM interactions, we repeated our MPS calculations with different values of D and D ′ , as we show in Sec.S2B of the SI [9].A robust gap exceeding 0.45 meV remains if the DM interactions are removed, and hence is an intrinsic property of the pure (Heisenberg) SSM.The natural explanation is the condensation of the S = 2 bound state.The minimal two-dimer model (inset, Fig. 1c) is sufficient to illustrate some very clear effects in the one-triplon branches.With antiferromagnetic intradimer interactions (J > 0), the field-induced condensation process is readily controlled using the inter-dimer bonds: if J ′ > 0, the t + branch condenses (Fig. 1c) in a conventional magnon BEC, whereas for J ′ < 0, the lowest branch of the two-triplon multiplet condenses (Fig. 1d).In the latter case, the one-triplon branches show several distinctive features: (i) the minimum energy of the t + branch is half of the binding energy, and the 0.6 meV observed by MPS is consistent with half the binding energy predicted by perturbation theory and iPEPS calculations [28]; (ii) this minimum occurs at a field, H * , lower than the extrapolated t + condensation field; (iii) beyond H * , the t + energy increases linearly with the field; (iv) the t 0 branch has a kink at H * and increases with a gradient that is twice as large; (v) the t − branch also has a kink at H * and increases with a gradient three times as large.These observations are completely generic for the condensation of S = 2 bound states and are fully consistent with both the experimental and the numerical results (Figs. 1e-f).
Turning now to the novel field-induced excitations, we consider first their spectral weight.In Fig. 3c we show the accumulated spectral weight computed by MPS in the same format as the experimental results of Fig. 3b, also extending the MPS results to higher energies (inset).The agreement is excellent, with both datasets reflecting the field-induced smoothing of the steps as increasing amounts of weight are shifted to new excitations that appear at intermediate energies.By benchmarking the total weight at high energies, it is clear that the spectral weight lost from the INS energy window (0.7-7.0 meV) is pushed lower (Fig. 3a).To quantify these weight-shifts, in Fig. 3d we show the MPS intensities integrated over energy ranges spanning the t 0 branch, below this branch and below the t − branch.Beyond the obvious agreement with experiment, the MPS results confirm the linear growth of weight in multi-triplon excitations with fields beyond approximately 16 T and the remarkable (approximately 2/3) loss of one-triplon weight over the field range from there to the 26 T edge of the critical region around H c .
To go further in identifying these multi-triplon states, it is useful to separate the spectral function into its longitudinal [S zz (ω, H)] and transverse [S xx (ω, H)] components, which we show in Figs.4a-b.The t 0 branch appears in the S zz (ω, H) channel and the t ± modes in the S xx (ω, H) channel, as we observe at low fields.However, above 21 T there is a rapid growth of low-energy spectral weight in S zz (ω, H), specifically over a broad range centred near 0.7 meV, and of weight around the t 0 energy in S xx (ω, H).In Figs.4c-d we show the energies and intensities of selected high-weight states obtained around H * in our cluster ED calculations for the SSM (meaning with no DM interactions), which as detailed in Sec.S2D of the SI [9] allow us to make an unambiguous identification of the nature (S and S z quantum numbers) of these states.It is clear that two-triplon composite states appear at low energies in both sectors and that many three-triplon composites are present at energies extending down to the t 0 branch.
The low-energy spectral weight in S zz (ω, H) (Figs. 4a  and 4c) is another strong statement for BEC of bound states: the operator S z can excite either a singlet to a t 0 triplet, at relatively high energy, or a transition inside the band of condensed excitations.As noted above, onetriplon condensation in SrCu 2 (BO 3 ) 2 would be accompanied by weight at and below the triplon bandwidth, 0.3 meV, which is too small to explain the MPS result.By contrast, the internal structure of the S = 2 bound state has four bands, of which the lowest can be considered as the "pinwheel" represented in Fig. 4e [28], and S z can induce "internal excitations" within this structure, a schematic example being shown in Fig. 4f.Using perturbation theory [40] with the parameters of SrCu 2 (BO 3 ) 2 , the energies of the lowest internal excitations should cover the range from 0 to approximately 1.8 meV, in agreement with the energy range and |t + t + ⟩ character we obtain (Figs.4a and 4c).Unlike the t + mode, the band centre of these excitations does not increase linearly with the field beyond H * , and the position and increasing intensity of this additional scattering contribution suggest an origin for the low-energy intensity observed at 25.9 T (Figs. 1e and 3a).
One of the most striking consequences of two-triplon condensation is the presence of low-lying composite three-triplon excitations.The linear increase of spectral weight away from the t 0 and t − branches is not specific to two-triplon condensation, as several composite two-triplon excitations would be observed if single triplons condense [41].A leading example is the |t + t 0 ⟩ state represented schematically in Fig. 4g, but in a scenario of one-triplon (t + ) condensation this would appear  5. a Schematic representations of the ground states (lower panels) and all associated excitations (upper panels) at low fields in the dimer-product phase (left) and at fields beyond H * (right).b Magnetization calculated by MPS in the field regime below the 1/8-plateau and compared with results obtained from NMR and magnetic torque measurements [14] (left axis); for comparison we show the measured intensity of the multi-triplon excitations below t0 (right axis).c Schematic illustration of the field-induced evolution of the primary spectral features identified by combining our neutron spectroscopy data with our MPS and ED results.∆ = 3.0 meV is the one-triplon energy gap at zero field.µ0H * = 21 T is the condensation field for two-triplon bound states (red dashed line).Green, red and magenta shading represent the field-induced increase of spectral weight in novel multi-triplon excitations of the types represented in Figs.4f-h and panel a.
in S zz (ω, H), whereas our discovery of low-energy |t + t 0 ⟩ weight in S xx (ω, H) (Figs. 4b and 4d) underlines its origin in S = 2 bound-state condensation.Quite simply, this scenario makes the number of excitations accessible by INS much larger, and adding three-triplon composites results in the complex spectra we observe both in INS and MPS. Figure 4h shows a schematic representation of a |t + t + t − ⟩ composite obtained by exciting a t − triplon in close proximity to a |t + t + ⟩ bound state, where in general we find that the low-S, low-S z members of the multiplet benefit from sizeable net binding energies.The bound-state population is reflected in the magnetization, and hence the intensity of the two-and three-triplon excitations is also expected to increase linearly with the applied field, as we observed in Fig. 3d.

Discussion
To summarize, we have combined high-field INS experiments and cylinder MPS calculations to reach the unambiguous conclusion that, in the field range between µ 0 H * = 21 T and µ 0 H c = 27 T, the ground state of SrCu 2 (BO 3 ) 2 is not a condensate of single triplets, but a condensate of two-triplon bound states.The consequences of this result, summarized in Fig. 5, are a distinctive one-triplon spectrum with altered gradients beyond a premature kink at H * , a one-triplon gap that persists at all fields and a high field-induced spectral weight of composite two-and three-triplon excitations both internal and external to the two-triplon bound states.We have shown that the weak DM interactions in SrCu 2 (BO 3 ) 2 play no qualitative role in determining this physics, but note nevertheless that they do affect the appearance of our results by creating a broad crossover regime below H * .Figure 5a compares the magnetization computed from the MPS ground state with different experimental measurements, all of which indicate that the crossover from near-zero to linearly increasing magnetization spans the range from 16 to 21 T, consistent with Figs.2h-k and 3d.
Placing this result in perspective, the bound-state condensate is an example of a spin nematic [33,34].This phase has a spontaneously broken rotational symmetry around the magnetic field, but its order parameter is a two-spin operator, and not simply a transverse magnetization.Spin-nematic phases have been predicted to occur just below the saturation field in systems where mixed ferro-antiferromagnetic interactions lead to magnon bound states [42][43][44], and some evidence for this phenomenon has been obtained both by high-field NMR in LiCuVO 4 [45] and by high-field calorimetry in [46].What is new in SrCu 2 (BO 3 ) 2 is that the spin-nematic phase occurs below the first magnetization jump, i.e. at low fields, with dilute triplons binding in a sea of singlet dimers.One primary hallmark of this phase is a gap to one-triplet excitations that persists at all fields and another is the presence of low-lying S = 2 excitations around H * .While NMR or INS provide experimental access to the one-triplet excitations, the weak DM terms ensure that ESR can couple to the S = 2 excitations (at Q = 0).We note that highfield Raman spectroscopy, which is sensitive primarily to S = 0 states and ∆S = 0 processes, has also provided a complementary view of the two-triplon condensation process through the behaviour of those excitations [47].
The spin-nematic phase has a close analogy with the physics of superconductivity.In a superconductor, Cooper pairs of fermions undergo a BEC, resulting in a state with a single-particle gap and collective excitations.In SrCu 2 (BO 3 ) 2 and the SSM, the single particle is a triplon, a bosonic particle with an infinite on-site repulsion (a "hard-core boson").The spin nematic realized in SrCu 2 (BO 3 ) 2 beyond 21 T can be seen as the BEC of Cooper pairs of bosons.This analogy establishes a direct connection between the superconducting gap of BCS theory and the one-triplon gap in SrCu 2 (BO 3 ) 2 , which is equal to half the binding energy at the condensation field (H * ).Building further on the analogy, the BCS-BEC crossover [48] links superconductors with large Cooper pairs that extend over many lattice sites to those with compact pairs, both developing coherence through the condensate.In SrCu 2 (BO 3 ) 2 , the pairs are the strongly localized pinwheel objects of Fig. 4e and the system represents the extreme BEC limit.Also in analogy with unconventional superconductivity, which generically competes with other instabilities of the interacting Fermi sea such as magnetic and charge order, the spin nematic in SrCu 2 (BO 3 ) 2 clearly competes with the bosonic CDW order that is favoured above H c , where the 1/8 magnetization-plateau phase is established.
Finally, our results represent new physics found by combining new capabilities in high-field neutron scattering with state-of-the-art numerical methods, and thus have direct implications in a number of disciplines.They demonstrate the profound importance to correlated condensed matter of scattering facilities able to operate at high magnetic fields, and on the technical side the HZB facility provided valuable information for future generations of high-field experiments.They illustrate directly the value of modern MPS methods in overcoming the long-standing barrier of computing dynamical properties in frustrated systems.For quantum matter, our results reinforce the message that novel composite states and novel forms of many-body order remain to be found at the confluence of strong frustration and controlled extreme conditions.

Methods
Neutron spectroscopy.The inelastic neutron scattering (INS) experiment was performed at the High-Field Magnet (HFM) facility of the HZB.The Extreme Environment Diffractometer (EXED) was a timeof-flight neutron instrument that could be operated either in diffraction and small-angle-scattering modes or as a direct-geometry spectrometer.It was designed specifically to function optimally in combination with the HFM, which provided horizontal magnetic fields up to 25.9 T. This magnet was a hybrid solenoid with a 30 • conical opening and could be rotated in its entirety with respect to the incoming beam in order to extend the accessible region of reciprocal space.
The sample was a single-crystalline rod with a mass of 2.5 g, which was grown by a floating-zone method with its (100) direction orientated approximately along the growth axis [49,50].The magnetic field was applied along the (001) direction in the horizontal scattering plane, a geometry chosen to give access to part of the (Q h , Q k , 0) plane and to avoid the combination of the applied magnetic field with the DM interactions breaking any further symmetries of the system.The incident neutron energies were E i = 4 and 8 meV and the rotation angle of the magnet relative to the incoming neturon beam was φ = −10 • .A dilution refrigerator purpose-built by HZB in collaboration with the University of Birmingham allowed the measurements to be performed at 200 mK.Data were collected for periods between 7 and 20 hours for each selected field strength with E i = 8 meV and between 10 minutes and 6 hours with E i = 4 meV.
Supporting INS data were collected on the tripleaxis spectrometers TASP, at the Paul Scherrer Institute (PSI), and Panda, at the Heinz Maier-Leibnitz Zentrum.On TASP, the same single crystal was orientated with (Q h , Q k , 0) in the horizontal scattering plane and a vertical magnetic field of 9.5 T was applied.The measurement was performed at 2.3 K, with a constant final neutron momentum k f = 1.5 Å−1 and with a liquid-nitrogen-cooled Be filter to suppress the λ/2 contribution.On Panda [22], measurements were performed in the same geometry, at 1.5 K and with final neutron momentum k f = 1.55 Å−1 , for fields of 4, 8 and 12 T. MPS Calculations.Matrix-Product States (MPS) are a variational Ansatz that can be used to represent the wavefunction of 1D quantum systems [51].The representation is formulated in terms of rank-3 tensors and its accuracy controlled by the tensor bond dimension, χ.A 2D system can be described by wrapping the lattice onto a cylinder and reformulating the model as a 1D Hamiltonian with further-neighbour interactions [52].The ground state in the MPS representation is then obtained using the Density-Matrix Renormalization-Group (DMRG) algorithm [51].
The spin physics in SrCu 2 (BO 3 ) 2 is well described by the SSM Hamiltonian with Heisenberg interactions J = 81.5 K and J ′ = 0.63J [53], supplemented by DM interactions D = 0.034J on the intra-dimer bonds [30] and D ′ = −0.02J on the inter-dimer bonds [29].Finite DM interactions were also required to obtain uniform magnetization distributions in our calculations of the ground state in an applied field.All of the spectral functions shown in the main text were computed with a cylinder of length L = 20 and circumference W = 4 sites, as shown in Sec.S2A of the SI [9], where we benchmark the effect of W by computing the energy per site, magnetization distribution and spectral function at two selected fields for a cylinder with W = 6.
To calculate the dynamical properties of the system, we employed the time-dependent variational principle (TDVP) [54].This method makes use of Lie-Suzuki Trotterization not of the Hamiltonian but of projectors to the tangent space of the MPS, which are constructed from the MPS with an accuracy dependent on χ.We used the two-site variant of the TDVP algorithm, which allows the bond dimension of the MPS Ansatz to be readjusted as the entanglement grows upon time-evolution, with a time-step of 0.16/J and keeping a maximum χ of 600.The excitations of the SSM are localized as a consequence of the strong frustration and thus the time-dependence appears as a slowly expanding time-cone.This property allowed us to continue the real-time evolution to rather long times even on short cylinders.
In the magnetized system, ⟨S z ⟩ is uniform while ⟨S x ⟩ and ⟨S y ⟩ are staggered, and hence the magnetic unit cell contains n s = 4 sites.To compute the dynamical structure factor, we therefore started from each of these sites individually and took the Fourier transform where x a and x b are the relative positions of sites a and b within the unit cell and The HFM/EXED instrument that was built at the HZB is described in a series of detailed publications [S1-S5].The sample orientation for our experiment was established by using EXED in diffraction mode with a wavelength range of λ = [0.6,6.3] Å.Three axes are used to define this orientation, Bu, Bv and Bu × Bv, as represented in Fig. S1a.The sample rotation angles with respect to the three axes, (6.0 • , 0.8 • , −2.7 • ), were determined by indexing a number of nuclear Bragg peaks and following their position on the detector at different magnet rotation angles, φ.The scattering geometry, defined in Fig. S1b, was constrained by the cone-shaped 30 • openings of the magnet.The predicted detector image, shown in Fig. S1c, agreed well with the measured image shown in Fig. S1d; these images were made using φ = −8 • , where five nuclear Bragg peaks were observed.The INS experiment was then carried out with φ = −10 • , where only three Bragg peaks were visible, but the larger rotation angle allowed larger scattering angles, which enlarged the accessible portion of reciprocal space.
The range of Q accessible on EXED was restricted by the high-field magnet, as Fig. S2a shows for an energy transfer of 3.0 meV, where almost one Brillouin zone is covered.This area was larger at lower energy transfers, but smaller at higher ones.The comparison in Fig. S2a makes clear that on EXED we could not access the highintensity point at Q = (−0.5,1.5, 0), or any of its equivalent points, found in a previous experiment at zero field [S7].While the structure factors at finite fields are not known, the earlier zero-field data showed that different excitations have markedly different structure factors, but also that this difference is detectable only at higher absolute values of Q.As a result, the only Q-dependence we could observe on EXED was an overall increase of intensity with increasing |Q|.
The range accessible along (0, 0, Q l ) was restricted by the choice of Q h and Q k as well as by the energy transfer, as shown in Fig. S2b.Thus a given point in the spectrum represented a certain range of Q l values and consequently a specific ratio of the measured dynamical structure factors, because S zz ∝ (1 − Q2 l ) whereas S xx ∝ 1 2 (1 + Q2 l ); these two quantities are illustrated in Figs.S2c-d.For the Q-integrated spectra shown in Figs.2g-l of the main text, higher energies contain stronger contributions from S xx (Q, ω) whereas lower energies are more reflective of S zz (Q, ω).We comment that Q k = 0 was used when calculating Q2 l in Fig. S2c-d, although in practice the integration range along (0, Q k , 0) was rather large, as in Figs.2g-l of the main text.

B. Data treatment and background subtraction
The INS data were corrected for the detector efficiency following a standard procedure where individual detector elements are normalized with respect to a vanadium standard.This correction and all subsequent cuts through the data were performed using the Mantid software package [S8].Figures 2g-l of the main text, and the slices of Fig. 1e at the same fields, show the backgroundsubtracted neutron intensity integrated over the entire available area of reciprocal space, Q h = [−0.25,1.25], Q k = [−0.75,0.75] and Q l = [−0.9,1.3], for the data collected with E i = 8 meV.This background was based on the zero-field data, where it is known that the lowestlying excitation is the one-triplon mode at 3.0 meV [S7, S9-S12], and at the low temperature of our experiment one expects a signal only on the neutron energyloss side of the spectrum.The background was defined by a three-step fitting procedure that we illustrate in Fig. S3.(i) Quasi-elastic scattering was approximated by a Lorentzian.(ii) After subtracting the Lorentzian, the elastic lineshape was fitted to six Gaussians.This number is purely empirical and was optimized to describe the shape of the elastic line.(iii) Once both Lorentzian and Gaussian contributions had been subtracted, the timeindependent part of the background was fitted to the form , with E i = 8 meV fixed and the constant I 0 as the only fitting parameter.The final background is then a sum of the three contributions.
For the intensity data measured with E i = 4 meV and shown in Fig. 3a of the main text, the integration range was Q h = [−0.15,0.85], Q k = [−0.5,0.5] and Q l = [−1.5,2.5].The zero-field results in this case did not have good statistics because the count time was only 10 minutes, and instead we based our background estimate on the 25.9 T data, which were collected for 5.5 h.Because the spectrum is not known at this field, the background estimation procedure was based on data points in the energy interval E = [−4.0,0.1] meV.A spurion was observed around 2.1 meV for data with E i = 4 meV and around 5.1 meV for data with E i = 8 meV.This feature was field-independent and appeared close to Q = 0.A Gaussian was fitted and subtracted to remove the spurion before fitting the spectrum of excitations originating in the sample.

C. Peak fitting and intensities
Here we describe in detail the procedure used to fit the data shown in Figs.2g-l of the main text.Previous theoretical analysis at zero field and our own numerical re-sults at finite fields motivate the presence, positions and widths of a significant number of candidate excitations, as a result of which we based our fit on multiple Gaussian profiles.At zero field, the peak width we obtained for the one-triplon mode around 3.0 meV was used to fix the widths of the excitations at 4.8, 5.5 and 6.5 meV, with the ratio of these widths taken from Ref. [S11].The zerofield parameters of the S = 1 two-triplon bound state were used to fix the properties of the S z = 0 branch of this multiplet at finite fields.We did not include Gaussians for the S z = ±1 branches of this multiplet, because their neutron intensities are expected to be weak.We note also that the S = 1 bound state has vanishing in- tensity at Q = (1.5, 0.5, 0), which is why it cannot be observed in some of the data we show in Fig. 1e of the main text, or in Fig. S4 below.
At finite fields, the one-triplon mode undergoes a Zeeman splitting, with the t − and t + branches positioned at energies ±α µ 0 H from the t 0 branch, where α ≃ 0.131 meV/T [S13].We fixed the widths of the t + , t − and S = 1, S z = 0 bound-state branches to that fitted for the t 0 mode.Only at 15 T, where the t 0 peak is close to the S z = +1 branch of the S = 1 multiplet, did we use the fitted width of the t + mode, which is well defined at this field, to fix the others.At all measurement fields from 18 to 25.9 T, we found that a single additional Gaussian placed at an energy below the t 0 branch, and up to three additional Gaussians placed below the t − branch, provided a good fit to the intensity profiles measured as a function of energy transfer.The positions and widths of these Gaussians were allowed to vary freely, except if the widths became much greater or smaller than that of the t 0 branch, in which case their widths were fixed relative to that of this branch.Despite their two-triplon nature, the excitations in this energy and field range are not expected to display a dispersion significantly larger than that due to the DM term [S11], and the width of the momentum-integrated t 0 excitation   In preparing these slices, the intensity data were integrated over the momentum interval Q l = [−0.9,1.3] and over an energy interval of width ∆E = 0.5 meV centred on the fitted positions of the red Gaussians (to obtain the top row) and the t0 Gaussians (to obtain the second row) in Figs.2g-l of the main text.The zero-field intensity of the t0 branch has been divided by 2 for illustration on the same colour scale.Below 18 T, where no extra intensity was detected below t0, the top two colour panels show simply an energy window well below the t0 window.The black boxes in the second row illustrate the region of integration used to produce the panels in the third and fourth rows.These show respectively the integrated intensities and mode positions as a function of Q h .
represents this bandwidth.We stress that our fit does not imply the presence of just one additional excitation below the t 0 branch, only that a single Gaussian is suitable for capturing all of the additional intensity arising from these excitations, a quantity we used in Fig. 3d of the main text and analyse in detail in Sec.S2D.For completeness we state that the position of the highest-energy mode below the t − branch at 23 T was fixed.In Fig. 1e of the main text, the symbols shown at 0, 15, 18, 20, 23 and 25.9 T were taken from the Gaussians obtained in Figs.2g-l.The data slices shown at fields between 0 and 15 T were taken from separate experiments: data at 4, 8 and 12 T were measured on Panda, as noted in the Methods section, at 1.5 K and at constant Q = (1.5, 0.5, 0) [S12, S14]; the data at 9.5 T were measured on TASP, also at constant Q = (1.5, 0.5, 0).The fitting procedure for these data, which are represented with field-scaled separation in Fig. S4, followed a simplified version of the procedure applied at high fields.For the Panda data, three Gaussians were fitted to the three one-triplon branches, with all parameters allowed to vary.For the TASP data, the three one-triplon branches were fitted with three Gaussians of the same width, with the positions of the t − and t + branches constrained to have the same separation from the t 0 position and with their intensities constrained to have half the intensity of the t 0 excitation.Two additional Gaussians were fitted at higher energies, with their individual positions, widths and intensities taken as free parameters.
Figure 3a of the main text shows the low-energy parts of three INS spectra collected with E i = 4 meV.The momentum-integration ranges and background fitting are described in Sec.S1B.To estimate the peak positions close to the elastic line at 20 and 25.9 T, we fitted one Gaussian whose width was fixed to the value obtained for the t 0 mode at zero field in a separate experiment on HFM/EXED.Figure 3d of the main text shows the integrated spectral weights found in energy windows below the t 0 and t − one-triplon branches.To illustrate the extraction of these weights, in Fig. S5 we show the experimental data at a representative field of 18 T.The one-triplon contributions were fitted by three Gaussians following the above procedure.The excess spectral weight induced below t 0 is then defined as that appearing above the onetriplon contribution between t + and t 0 , highlighted with red shading in Fig. S5 and it is integrated numerically.The excess weight below t − is defined in the same way and is shown by the purple shading.

D. Q-dependence
As Figs. 2a-f of the main text make clear, the spin excitations of SrCu 2 (BO 3 ) 2 are extremely insensitive to Q, and this was the basis on which we integrated our EXED results over the full reciprocal space for further analysis.
Here we nevertheless inspect the differences discernible over the Q range covered in our measurements (Fig. S2a).The top two rows of Fig. S6 (colour panels) represent the intensities measured in the (Q h , Q k , 0) plane in energy windows selected to cover the regions below (top row) and spanning (second row) the t 0 excitation branch for all six measurement fields.While the shape of the intensity distribution is indeed very similar in all panels, with the highest intensities at the largest |Q| values, the net intensity reflects the general shift from the t 0 branch to the novel multi-triplon excitations below it (red Gaussians in Figs.2g-l) as the field is increased, as quantified in Fig. 3d of the main text.
The third row of panels in Fig. S6 quantifies the evolution of the two intensity contributions with Q h by integrating over Q k = [−0.1,0.1] (black boxes in second row).At 0 and 15 T, where no statistically significant extra intensity could be discerned below the t 0 branch (Figs.2g and 2h of the main text), the intensity in the selected integration window can be attributed to the sample environment; because different energy intervals correspond to different intervals in Q l , we do not subtract this background contribution.The intensities measured for the t 0 branch at each field show a largely linear dependence on Q h and a steady decrease as the field is increased above 18 T; the intensities of the multi-triplon contribution below t 0 also increase with |Q h | and increase with field until they match the t 0 weight.The bottom row of panels in Fig. S6 shows how the position of the t 0 branch, and the centre position of the spectral weight below t 0 , vary weakly with Q h .The data were discretized in units of ∆Q h = 0.1 for fields up to 20 T and ∆Q h = 0.2 at 23 and 25.9 T. The very minor dispersion observed in all cases is consistent with the bandwidth identified in zero-field studies [S10, S11].

S2. NUMERICAL CALCULATIONS AND ANALYSIS
A. Dependence on cylinder size An MPS Ansatz describes a one-dimensional (1D) quantum many-body wavefunction.A 2D system can be reformulated as a 1D model on a cylinder with longranged interactions [S15], such that MPS-based DMRG [S16, S17] can be used to find the ground state and compute the associated static observables.The dynamical structure factor (DSF) of a 2D system can also be computed using MPS-based methods, and for this work we applied the TDVP [S18-S20], as described in the Methods section.Here we calibrate how the finite length and circumference of the cylinder affects our calculations.All results reported in the main text were obtained using a cylinder circumference of W = 4 sites and a length of L = 20 sites, as represented in Fig. S7a.As a first step towards justifying that these values of W and L return results fully representative of the 2D SSM, we computed the ground-state energy and field-induced static magnetization with W = 6 (Fig. S7b), finding that both quantities change rather little.
TDVP calculations become extremely time-intensive on wider cylinders for two reasons.First the bond dimension of the MPS should increase exponentially with W in order to describe the state with equivalent accuracy.Second, wider cylinders introduce further-neighbour interactions in the effective 1D model that increase the virtual dimension of the Matrix Product Operator in proportion with W .The resulting O(W 2 exp(3W )) scaling of the calculation time therefore prevented us from obtaining the DSF with W = 6 at all values of the magnetic field required for our study, but we did compute it at 15 and 20 T, and in Fig. S8 we compare the Q-integrated DSFs obtained for W = 4 and 6 at both fields.Because the allowed k y values for a W = 4 cylinder are 0 and π, while for a W = 6 cylinder they are 0, 2π/3 and 4π/3, increasing W results in more spectral weight for each mode and thus we normalize the spectra to the same integrated weight.At first sight there is remarkably little change from W = 4 to 6, and on close inspection one may ascertain that the loss of one-triplon peak height in the W = 6 spectrum reappears at energies between and above the three primary branches.This could be expected on the grounds that the wider cylinder allows more space for the development of multi-triplon dynamics.Overall, however, the effects of W are small and the effects of L are negligible beyond L = 20; while we do not show the latter null result, it is readily understood from the extreme localization of all one-and higher-triplon excitations in the ideally frustrated SSM geometry.

B. Role of Dzyaloshinskii-Moriya interactions
As discussed in the main text, a primary consequence of two-triplon bound-state condensation is the presence of a one-triplon gap at all fields.However, in SrCu 2 (BO 3 ) 2 there are weak DM interactions on both the intra-and inter-dimer bonds, as represented in Fig. 1b of the main text, and in particular the intra-dimer D term leads to a direct mixing of S z sectors in the SSM.It is known from studies of a spin ladder with intra-dimer DM terms that an avoided crossing can occur, such that the one-triplon gap reaches a finite minimum at the expected condensation field, and then increases again at higher fields [S21].To eliminate the possibility that the behaviour we observe is a DM-induced avoided crossing, we have performed MPS calculations at different values of D to identify the evolution of the lowest-lying excitations.We note again that a small but finite D is required to obtain a uniformly magnetized ground state in our calculations, and for all the MPS results shown in the main text this was set to the experimentally deduced values D = 0.034J [S22], with D ′ = −0.02J[S23].MPS results for the D, D ′ = 0 limit in the presence of a field may be obtained by varying both parameters and extrapolating the results.
We computed the transverse (S xx ) and longitudinal (S zz ) components of the dynamical structure factor at µ 0 H = 21 T for five different D values, with D ′ fixed to 0.59D, and extracted the lowest-lying excitations in each sector.The one-triplon gap, read from S xx by fitting a Gaussian to the low-energy spectral weight, is shown in Fig. S9.For comparison we show the generic func- tional form ∆(D) = ∆ 2 0 + c 2 D 2 , which interpolates between a constant at D = 0 and a linear increase at large D. This simple form offers a convincing fit, demonstrating that the lower bound for the one-triplon gap as D → 0 is ∆ 0 = 0.45 meV.In the inset of Fig. S9, we show that this gap is also not a consequence of the field choice, as varying H confirms the linear increase away from µ 0 H * = 21 T with decreasing field, and initially with increasing field, as represented qualitatively in Fig. 1e of the main text.For completeness we show also the peak energy determined in the same way from S zz , which appears as a result of internal excitations of the two-triplon bound state (Figs.4a and 4f of the main text), and moves below the S xx peak at the highest fields we study.If one accepts the result of Sec.S2A, that this constant minimum gap of order 0.5 meV at H * is not a finite-size effect, it must therefore be an intrinsic property of the pure SSM (i.e. of the model with only Heisenberg interactions).We repeat the analogy to superconductivity drawn in the main text, where a single-particle excitation formed by pair-breaking in the condensate has an energy gap equal to half the binding energy.This implies that the two-triplon binding energy in SrCu 2 (BO 3 ) 2 is around 1 meV, a value fully consistent with the results of third-order perturbation theory [S24].

C. Q-dependence of MPS results
A full comparison of our cylinder MPS results with the DSF measured in experiment requires the inclusion of the structure factor of SrCu 2 (BO 3 ) 2 , as noted in the Methods section.When the DSF is integrated over a large volume in reciprocal space, as in Figs.1e, 1f, 2, 4a and 4b of the main text, the effect of the structure factor is largely invisible, and for this reason we do not perform extensive comparisons of the Q-dependence in our INS and MPS results.To illustrate the role of the structure factor, for example in preparing constant-energy slices of the type shown in Fig. S6, we integrate the DSF obtained from MPS over energy windows selected to highlight a particular excitation feature and display the results in the space of (Q h , Q k ).Our MPS calculations allow the further advantage of separating the excitations overlapping in a specific energy window by considering the transverse and longitudinal channels separately.In Fig. S10 we show a number of constant-energy slices prepared from our MPS data at a field equivalent to 25 T. The t 0 branch of the one-triplon excitation is found in the longitudinal (S zz ) channel and its Q-dependence is shown in Fig. S10a by selecting an energy window of 0.5 meV around the field-induced centre position of 4.05 meV.Similarly, the internal excitations of the |t + t + ⟩ pinwheel, shown in Fig. 4a of the main text, are found at low energies in S zz (Fig. S10b).Otherwise, stronger spectral contributions from the novel multi-triplon excitations appear in the transverse (S xx ) channel, and in Figs.S10c-e we compare the intensity distributions for energy windows below the t 0 branch (Fig. S10c), around the t 0 branch but in the opposite channel (Fig. S10d) and below the t − branch (Fig. S10e).In each case we observe distinctive patterns that could be analysed in future high-field INS experiments.

D. Identification of multi-triplet excitations
To identify the nature of the excited states found in our cylinder MPS calculations (Figs.1f and 4a-b of the main text), we exploit the fact that the DM interactions of the system are too weak to introduce confusion (by a mixing of spin sectors) and neglect them completely.For identification purposes it is fully sufficient to use a pure SSM and to perform ED on a cluster of only 4×4 spins, where (as noted above) the available k x and k y momentum sectors are 0 and π.The ED spectrum is computed in symmetry sectors labelled by the quantum numbers of total spin (S tot ), total magnetization (S z ) and the momenta k x and k y , and we make use of the separate S zz and S xx components of the DSF to identify states from the action of the S z and S x spin operators.
As discussed in the main text, the two-triplon bound state, |t + t + ⟩, is the first to condense as the applied magnetic field is increased.This state belongs to the symmetry sector S = 2, S z = 2, and its lowest-energy component appears at (k x , k y ) = (0, 0) on the 4×4 cluster.In the limit J ′ /J → 0, the energies of n-triplon states approach nJ, and thus all the states in the spectrum are readily identified by following the evolution of their en- ergies with J ′ /J, as we show in Fig. S11 for the lowest three S tot sectors.It is clear that the different n-triplon sectors remain well separated in energy for smaller values of J ′ /J = 0, but that their energies show significant overlap as the interaction ratio approaches the SrCu 2 (BO 3 ) 2 value of J ′ /J = 0.63.
For more detailed insight into the character of each n-triplon sector, we consider the action of the operators S z and S x on the two-triplon bound state.From the Lehmann representation of the DSF, the excitations of a given initial state lie in symmetry sectors allowed by the selection rules The excitations allowed by these rules do not contribute    , ω, H)) of the DSF.On this cluster, the lowest-energy component of the two-triplon bound state reaches the singlet energy at a field equivalent to µ0H * = 27.7 T when using the g-factor of SrCu2(BO3)2 .We show all of the higher-weight states appearing at low energies in the momentum sectors (0, 0), (π, 0) and (π, π).Because the spin Hamiltonian is of pure Heisenberg type, the field term causes only a linear evolution whose gradient is determined by S z , with a kink and a ∆S z = 2 change of slope at H * .equally to the spectral weight, because all have different matrix elements.To deduce the contribution of each state in the spectrum, we compute the overlap where α = x, z, |ψ 0 ⟩ is the ground (initial) state and |ψ e ⟩ is the excited (final) state.We take N ≥ 0.01 as the criterion to regard |ψ e ⟩ as a state making a significant contribution to the DSF.We note that, while the initial state is always the lowest-energy state of the (S, S z , k x , k y ) = (2, 2, 0, 0) sector, this is not in fact the ground state at H < H * .The final step is to deduce the energies of these highweight levels in the presence of a magnetic field.Because of the purely Heisenberg nature of the Hamiltonian, every state undergoes a simple Zeeman shift determined by its S z quantum number, giving which allows us to track the origin of all the observed states at low fields, as we show in Fig. S12.The small cluster size we have used for this exercise results in nontrivial finite-size effects on the energy levels: in a field, the t + branch crosses the singlet before the lowest |t + t + ⟩ branch, and the crossing field defined from the latter is µ 0 H = 27.7 T.However, these quantitative details are not relevant for the qualitative purpose of identifying the high-weight excitations appearing in an INS experiment on a system with a two-triplon bound-state condensate (i.e. a spin nematic) as the field-induced ground state.Figure S12 provides the full justification for the statements made in the main text that the many excitations visible in the INS and MPS spectral functions are composite two-and three-triplon excitations of several differ-

FIG. 1 .
FIG. 1. a Crystal structure of SrCu2(BO3)2 showing one layer of Cu 2+ ions (half a unit cell along ĉ).b Representation of the superexchange couplings within (J) and between (J ′ ) the Cu-Cu dimers, which reproduce the Shastry-Sutherland model (SSM), and of the corresponding Dzyaloshinskii-Moriya (DM) interactions ( ⃗ D and ⃗ D ′). c-d Energy levels of a simple J-J ′ plaquette system (inset in panel c) shown as a function of magnetic field.For J ′ /J > 0 (c), the t+ one-triplon branch condenses first, whereas when J ′ /J < 0 (d), the two-triplon bound state condenses first.The two scenarios result in different gradients after condensation.e Excitation spectra measured by inelastic neutron scattering (INS) over a wide range of fields.Data at low fields (4, 8, 9.5 and 12 T) were measured at Q = (1.5, 0.5, 0) as described in the Methods section.Data from our HFM/EXED measurements (0, 15, 18, 20, 23 and 25.9 T) are presented with full momentum integration, as described in Sec.S1 of the SI[9].Colour contours represent the neutron intensity and symbols show the fitted positions of individual excitations (Figs.2g-l).Solid grey lines show the projected locations of the one-triplon branches for a scenario of one-triplon condensation and no DM interactions (panel c), dashed black lines their locations for two-triplon condensation (panel d).The dashed green line shows the Sz = 0 branch of the S = 1 multiplet within the two-triplon bound state.The dotted red line up to 15 T shows the Sz = 2 branch of the S = 2 multiplet as measured by electron spin resonance (ESR)[10], i.e. at Q = (0, 0), where the dispersion of this branch reaches its maximum[11].The red shading represents the estimated bandwidth of this branch, whose minimum determines the binding energy, and our linear extrapolation beyond 15 T indicates its effect on gap closure.f Total momentumintegrated spectral functions obtained from our cylinder MPS calculations on the SSM with additional DM interactions over the same wide field range.

FIG. 2 .
FIG.2.a-f Measured neutron intensity (colour contours) shown as a function of energy transfer and of momentum Q h for six different field strengths.These data were integrated over the intervals Q k = [−0.75,0.75] and Q l = [−1.5,2.5] and no background was subtracted.It is clear that all excitations have very little dispersion, justifying the integration of our data over a large region of reciprocal space.g-l Neutron intensity obtained by further integration over Q h = [−0.25,1.25] (black symbols), shown as a function of energy transfer at the same six field strengths.A single background was subtracted as detailed in Sec.S1B of the SI[9].The one-triplon excitation (solid blue lines with blue shading) exhibits Zeeman-type splitting in the applied field and its branches were fitted by related Gaussians (Sec.S1C of the SI[9]).The Sz = 0 branch of the S = 1 two-triplon bound state (dark turquoise lines and shading) appears at a constant energy.Additional intensity develops below the t0 one-triplon branch at µ0H ≥ 18 T (solid red lines with red shading) and was fitted with a single Gaussian.Further additional intensity visible below the t− one-triplon branch at µ0H ≥ 15 T (dotted black lines with grey shading) was fitted with multiple Gaussians.
FIG. 3. a Data collected at zero field (black), 20 T (blue) and 25.9 T (red) with incident energy Ei = 4 meV to resolve the low-energy excitations (solid blue and red lines).The data integration range and determination of the background (solid black line) are detailed in Sec.S1B of the SI [9].The mode positions were determined by fitting to a single Gaussian (dotted blue and red lines and shading).b-c Accumulated spectral weight shown as a function of energy transfer at each field for our INS (b) and MPS data (c).The inset shows an extended energy range calculated by MPS.d Integrated intensity of different low-energy spectral contributions, shown as a function of field for both the INS and MPS data.
FIG. 4. a-b Spectral functions S zz (ω, H) and S xx (ω, H) obtained by cylinder MPS calculations and compared with the one-triplon branches obtained from the one-and two-triplon-condensation scenarios represented in Fig. 1. c-d Energies and ground-state overlaps of states obtained in illustrative ED calculations performed around H * on a 4×4 cluster.e-h Schematic representations of the "pinwheel" structure of the two-triplon bound state (e), an internal excitation of the pinwheel (f), which retains the t+t+ character, a two-triplon excitation consisting of a t+ and a t0 (g) and a t+t+t− three-triplon excitation (h).In each case illustrated, the coefficients of the basis states shown have the same amplitude and differ only by their phases.
S α b,0 (0)⟩dt, with N the number of unit cells, R the spatial position of each cell and spin components α = α ∈ {z, x}.To obtain the results shown in Figs.1f and 4a-b, we integrated these functions over the full reciprocal space.E.F., J.-R.S. and A.A.T. and on Panda by M.E.Z.Data analysis was performed by E.F., M.E.Z. and H.M.R. MPS and ED calculations were performed by M.N. and their theoretical interpretation was provided by M.N., B.N. and F.M. The manuscript was written by E.F., H.M.R., M.N., B.N. and F.M. with contributions from all the authors.

(
FIG. S1. a Schematic representation of the instrument geometry, showing the magnet axis (black line), the neutron beam (thick blue line) and the detector positions.The axes Bu and Bv (red arrows) define the sample orientation with respect to the field.The inset shows a photo of the sample mounted in the dilution refrigerator, indicating the orientation.b Scattering triangle showing the wavevectors ki of the incoming neutrons, k f of the outgoing neutrons and the scattering vector, Q = (Q h , Q k , Q l ).The cones represent the opening angles of the magnet, whose rotation with respect to the incoming beam is φ.c Simulated detector view for φ = −8 • , showing the predicted positions of five different nuclear Bragg peaks.d Detector view from the experiment with the same settings as in panel c.Panels a and c were generated using the EXEQ prediction tool [S6].

FIG
FIG. S2. a The area inside the solid circle shows the portion of reciprocal space accessible in our EXED experiment and the area inside the dashed lines its comparison with the range covered in a typical time-of-flight INS experiment with Ei = 8 meV and detector coverage 10-100 • .This slice in (Q h , Q k , 0) was made for energies in the range [2.79, 3.29] meV around the one-triplon excitation.The pattern shown with faded colour is the corresponding structure factor calculated from the zero-field MPS spectrum.b Representation of the connection between the energy transfer and the Q l component of the scattering vector.The data shown were taken at zero field and integrated over the range Q k = [−0.75,0.75].c, d Structure-factor components (1 − Q2 l ) and 1 2 (1 + Q2 l ), shown as functions of energy transfer and of Q h , again with integration over Q k = [−0.75,0.75].

FIG. S3 .
FIG.S3.Three-step process for background estimation based on the zero-field data.The black data points are those used for each step of the fitting procedure and the grey ones are those which were disregarded.The red lines show the fit at each step and the final background is their sum.

FIG. S4 .
FIG. S4.Neutron intensities measured as functions of energytransfer at constant Q = (1.5, 0.5, 0) on Panda and TASP at selected field strengths.The data points are offset along the vertical axis by a distance proportional to the corresponding applied field value, as noted at right.The blue shading shows fits to the three one-triplon branches.The grey shading shows additional excitations at higher energies.The data from TASP were scaled to those from Panda by matching the integrated intensities of the measured t0 branches.

FIG. S5 .
FIG. S5.Illustration of INS data fitting, shown for an applied field of 18 T.The one-triplon branches are fitted by three Gaussians of fixed width (black lines and grey shading).Red shading indicates the extra intensity found below the t0 branch and purple shading the extra intensity appearing below the t− branch.

FIG. S6 .
FIG.S6.Data slices showing the Q-dependence of the neutron intensity below the t0 excitation branch (top row) and around the t0 branch (second row).In preparing these slices, the intensity data were integrated over the momentum interval Q l = [−0.9,1.3] and over an energy interval of width ∆E = 0.5 meV centred on the fitted positions of the red Gaussians (to obtain the top row) and the t0 Gaussians (to obtain the second row) in Figs.2g-lof the main text.The zero-field intensity of the t0 branch has been divided by 2 for illustration on the same colour scale.Below 18 T, where no extra intensity was detected below t0, the top two colour panels show simply an energy window well below the t0 window.The black boxes in the second row illustrate the region of integration used to produce the panels in the third and fourth rows.These show respectively the integrated intensities and mode positions as a function of Q h .

6 FIG. S8 .
FIG. S7.Illustration of the cylinder geometries used for our MPS calculations.a Cylinder width W = 4. b W = 6.The red (positive) and blue (negative) circles indicate the magnetization per site obtained at H = 20 T, and the degree to which it is uniform across the central regions of the cylinder.The dashed rectangle indicates the region used for the spatial Fourier transform.

FIG. S9 .
FIG. S9.Excitation energies obtained by MPS calculations on cylinders of width W = 4 to illustrate the effect of the DM interactions at a field equivalent to µ0H = 21 T. On the x-axis, D and D ′ are scaled to their experimental values.On the y-axis, the position of the t+ branch was obtained from the lowest peak in the S xx component of the dynamical structure factor and the vertical bar indicates the width of this peak.(Inset) Variation of the t+ gap with magnetic field at fixed D = 0.034J and D ′ = −0.02J.For comparison we show the peak energy of the low-energy spectral weight in the S zz component, which arises from internal excitations within the |t+t+⟩ bands (Fig. 4a of the main text).

8 FIG= 1 bS tot = 2 cFIG. S11 .
FIG. S10.Constant-energy slices at zero field displayed in(Q h , Q k ).a Slice including the one-triplon t 0 branch extracted from S zz .a Slice including the internal excitations of the two-triplon bound state, extracted from S zz .c-e Slices including different multi-triplon bound states extracted from S xx for energy intervals 2.6-3.1 c, 3.8-4.3d and 5.5-6.0 meV e.We comment that the weak breaking of fourfold symmetry visible in these panels is a consequence of the cylinder geometry (20×4, Fig.S7a).

S
FIG.S12.Identification of composite multi-triplon states in the spectrum of the 4×4 SSM spin cluster in a magnetic field strong enough to condense the lowest |t+t+⟩ state.a Excitations accessible in the longitudinal component (S zz (k, ω, H)) of the DSF.b Excitations accessible in the transverse component (S xx (k, ω, H)) of the DSF.On this cluster, the lowest-energy component of the two-triplon bound state reaches the singlet energy at a field equivalent to µ0H * = 27.7 T when using the g-factor of SrCu2(BO3)2 .We show all of the higher-weight states appearing at low energies in the momentum sectors (0, 0), (π, 0) and (π, π).Because the spin Hamiltonian is of pure Heisenberg type, the field term causes only a linear evolution whose gradient is determined by S z , with a kink and a ∆S z = 2 change of slope at H * .