Dynamics of quantum spin-nematics: Comparisons with canted antiferromagnets

The identification of spin-nematic states is challenging due to the absence of Bragg peaks. However, the study of dynamical physical quantities provides a promising avenue for characterizing these states. In this study, we investigate the dynamical properties of spin-nematic states in three-dimensional quantum spin systems in a magnetic field, using a two-component boson theory that incorporates magnons and bi-magnons. Our particular focus lies on the dynamical spin structure factor at zero temperature and the nuclear magnetic resonance (NMR) relaxation rate at finite temperatures. Our findings reveal that the dynamical structure factor does not exhibit any diverging singularity across momentum and frequency while providing valuable information about the form factor of bi-magnon states and the underlying structure of spin-nematic order. Furthermore, we find a temperature dependence in the NMR relaxation rate proportional to $T^3$ at low temperatures, similar to canted antiferromagnets. A clear distinction arises as there is no critical divergence of the NMR relaxation rate at the spin-nematic transition temperature. Our theoretical framework provides a comprehensive understanding of the excitation spectrum and the dynamical properties of spin-nematic states, covering arbitrary spin values $S$ and encompassing site and bond nematic orders. Additionally, we apply the same methodology to analyze these dynamical quantities in a canted antiferromagnetic state and compare the results with those in the spin-nematic states.


I. INTRODUCTION
Exploring hidden order phases in condensed matter is challenging due to their elusive nature in conventional static measurements.An intriguing hidden order is a spin-nematic order in spin systems, characterized by the absence of a magnetic Bragg peak and the presence of a broken partial spin rotation symmetry resulting from spin quadrupolar order [1].The emergence of spinnematic order arises from the Bose-Einstein condensation of bound magnon pairs [2].Theoretical investigations have predicted the appearance of spin-nematic order in various systems, including spin-S bilinear-biquadratic systems [3][4][5][6], frustrated ferromagnets [2,[7][8][9][10][11][12][13][14][15], and quantum dimer systems [16][17][18][19].Importantly, the symmetry of magnon pairing varies across systems, with Swave symmetric magnon pairs leading to quadrupolar order on each lattice site and magnon pairs with other symmetries, such as d-wave symmetry, leading to bond nematic order with quadrupolar order on bonds.Subsequent to these predictions, there has been a growing interest in experimental investigations aimed at verifying the existence of spin-nematic phases [20][21][22][23][24][25][26][27][28].However, conclusive experimental confirmation requires further insights into the intrinsic characteristics of these phases.
Detecting spin-nematic phases is challenging, but dynamical quantities show promise as valuable tools for identification.
In one-dimensional spin-nematic Tomonaga-Luttinger liquids, the temperature dependence of the nuclear-magnetic-resonance (NMR) relaxation rate 1/T 1 shows a slower decay compared to conventional one-dimensional antiferromagnets [29], providing a means to detect spin-nematic liquids (see, for example, Ref. [23]).However, in three-dimensional spin-nematic ordered phases, different methods have yielded inconsis-tent results regarding the temperature dependence of the NMR relaxation rate at low temperatures [30,31], leading to a lack of comprehensive understanding.In a 1/N expansion, the temperature dependence of 1/T 1 was estimated as T 5 [30], whereas a field-theoretical approach concluded T 7 [31].Both results deviate from the temperature dependence of T 3 for the conventional antiferromagnets derived by Moriya [32][33][34].Another distinct difference from antiferromagnets was argued to be the absence of a critical divergence in 1/T 1 at the transition temperature of the spin-nematic phase [31].
In this paper, we investigated the dynamical properties of spin-nematic phases in three-dimensional systems in a magnetic field to address this confusing issue.To accomplish this, we developed a methodology for describing spin-nematic states using a boson representation, which enables us to use the established standard interacting boson theory.Expanding on our previous study [28], where we employed a two-component boson theory including both magnons and bi-magnons to investigate the thermodynamic properties of the spin-gap phase near a spinnematic phase, the present work extends this methodology.Specifically, we investigated the dynamical quantities in the spin-nematic and paramagnetic phases, considering the influence of interactions and bi-magnon condensation.Our approach incorporates the structure of spin-nematic order, as captured by the form factor of bound magnon pairs, enabling the analysis of various types of spin-nematic order, including site-nematic and bond-nematic orders, within a unified theoretical framework.Consequently, this method enables a comprehensive study of the spin-nematic phases.
In particular, we calculated the dynamical spin structure factor at T = 0 and the NMR relaxation rate at finite temperatures.We find that the dynamical struc-arXiv:2308.12569v2[cond-mat.str-el]17 Feb 2024 ture factor shows no diverging singularity at any momentum and frequency, consistent with previous studies on specific spin models [5,6,30,35].Our results uncover that the dynamical spin structure factor carries information about the form factor of bi-magnon states, providing insights into the underlying structure of spin-nematic order.Furthermore, in contrast to previous studies, we observe that the NMR relaxation rate in the spin-nematic phase exhibits a temperature dependence proportional to T 3 , similar to canted antiferromagnets.A notable distinction from antiferromagnets is the absence of a critical divergence in the NMR relaxation rate at the spinnematic transition temperature, consistent with a prior field-theoretical study [31].Among these results, explicit dependence on pairing symmetry appears only in the intensity of the one-magnon mode in the dynamical structure factor.The remaining results constitute universal features of spin-nematic phases that do not explicitly rely on pairing symmetry.Additionally, using the same methodology, we reanalyzed these dynamical quantities in a canted antiferromagnetic phase, facilitating comparisons between spin-nematic and antiferromagnetic systems.
This paper is organized as follows: Section II details the formulation for describing spin-nematic phases in a magnetic field.Section III presents the results for the dynamical quantities of spin-nematic phases.Section IV describes an application of our results to the spin-nematic phase in the low magnetic field regime of spin-dimer systems.Section V focuses on the analyses of canted antiferromagnets.Finally, we summarize our main findings in Sec.VI, along with discussions.

II. BEC OF BI-MAGNONS
When magnetic excitations above a spin-gapped ground state form stable two-particle bound states with lower energy than two scattering excited states, these bound states close the gap as the magnetic field varies.Subsequently, at low temperatures, a macroscopic number of magnetic bound particles are induced by the field, interacting with each other, which leads to a spin-nematic phase [2].In this section, we introduce an interacting boson theory to describe the spin-nematic phase in a magnetic field.In the Supplemental Materials of our previous paper [28], we discussed the characteristics of the spingapped phase adjacent to the spin-nematic phase.Extending the previous analysis, we investigate the effects of interactions and Bose-Einstein condensation (BEC) of bi-magnons.

A. Interacting boson theory for spin nematics
In the following theory, we treat the spin-nematic phases near saturation in spin-S systems on a threedimensional lattice Λ, where S can be arbitrary.One example of such a system is a layered square-lattice frustrated ferromagnet, where spins on each layer couple with the ferromagnetic 1st-neighbor and antiferromagnetic 2nd-neighbor interactions [2,12].These arguments are equally applicable to the spin-nematic phase observed in low fields in frustrated spin-dimer systems, including the two-dimensional Shastry-Sutherland model [16,17,28,36], as discussed in Sec.IV.

Boson mapping
To describe dilute magnetic excitations on top of the fully polarized state, we set up Fock spaces for both single magnon excitations (S z = −1) and bound pairs (S z = −2) as the Hilbert space.To manage the state overlap between two magnons and a bound pair, we extend the Hilbert space and introduce an effective repulsive interaction in the Hamiltonian.The creation operator for a bound pair is expressed as b Here a † k denotes the bosonic creation operator of a magnon with momentum k, N Λ is the number of the lattice sites, and g q is a form factor satisfying the normalization 2(N Λ ) −1 q |g q | 2 = 1.We use the form factor derived from the exact bound two-magnon eigenstate in the lowest-energy mode above the fully polarized state, available through various methods [10,14,37,38].In spin-nematic phases in frustrated ferromagnets on the zigzag chain [39], the square lattice [2], and the body-centered-cubic lattice [14], the function g q is real, whereas, in a frustrated ferromagnet on the triangular lattice, g q can be complex due to a chiral degeneracy [13].
We investigate an interacting boson system with these two types of magnetic excitations, characterized by the Hamiltonian Here, ε 1,k (ε 2,k ) and v 1,q (v 2,q ) denote the excitation energy and the repulsive potential of one-magnons (bound bi-magnons), respectively, while u q represents the repulsion between magnons and bi-magnons.These interactions encompass both microscopic interactions and the effect of removing the extended excess boson space.A similar approach is adopted in the flavor wave expansion, where repulsive interactions are included by expanding Holstein-Primakoff type square root operators [4,35].The energies ε 1,k and ε 2,k include the Zeeman energies h and 2h, respectively, with explicit expressions: 2,k + 2h.The Hamiltonian with only bi-magnon operators b has been previously discussed in studies of spin-nematic states [8,35,39], while the Hamiltonian involving both types of bosons has been utilized to describe the spin-nematic Tomonaga-Luttinger liquid [40].
Here, we briefly delve into the generality of the Hamiltonian form.The original spin system maintains spin U(1) rotation symmetry around a magnetic field.During rotation, the boson operators a and b undergo a phase change to ae iθ and be 2iθ , respectively.Due to the U(1) symmetry, interactions involving the same kind of bosons must have an equal number of creation and annihilation operators.Interactions between different boson species are allowed, provided they satisfy U(1) invariance.For example, a † ab † b and a † a † b.In our study, we omitted pair creation and annihilation terms such as aab † and a † a † b for simplicity, as they only quantitatively alter the physical quantities in our mean-field treatment.
To simplify the analysis, we disregard the momentum dependences in the potentials, specifically setting v 1,q = v 1 , v 2,q = v 2 , and u q = u.Thus, microscopic interactions are approximated with local on-site effective repulsions.Since we focus on the low-density regime, where bosonic particles are well separated each other, microscopic interaction details are not crucial; instead, they can be represented by on-site effective couplings.The remaining information from the microscopic models are incorporated in the energy spectra ε 1,k and ε 2,k , as well as in the form factor g q of bi-magnons in this framework.Additionally, for simplicity, we assume the space inversion symmetry, implying In the extended boson Hilbert space, we have determined the matrix elements of spin operators in a lowdensity regime.For S = 1/2 [28], we find that while for arbitrary spin S, where g(0) = N Λ −1 q g q denotes the on-site form factor of bi-magnons, and g(0) = 0 for S = 1/2 due to the hard-core on-site repulsion.In the expression for S − k , we have omitted terms with three or more operators.The terms with b † k+q a q represent transitions from a single magnon to a bound pair.An illustration of these transitions is shown in Fig. 1.The mapping between spin and boson operators in S z k was used to construct the bosonic low-energy effective theory for the spin-nematic Tomonaga-Luttinger liquid in one dimension [39], where the field-theoretical prediction showed an excellent agreement with the numerical results obtained from the corresponding S = 1/2 spin model [8].Additionally, for the S ≥ 1 case, approximating the bi-magnon state as consisting solely of two magnons on the same site (i.e., g(k) = 1/ √ 2) makes the mappings (4) and ( 5) equivalent to one form of the flavor-wave expansions [4,35].However, we note that, in the exact two-magnon bound states at the saturation field of quadrupolar phases, the two magnons are usually spatially spread out and distributed at nearby sites [41,42].
Energy spectrum for single-particle excitations with S z = −1 and bound excitation pairs with S z = −2 in a gapped phase.The arrows represent selected matrix elements of the spin lowering operator S − q .
We also find that a spin-pair lowering operator, important in describing the spin-nematic order parameter, operates as in the low-density limit.Establishing a mapping between spin and boson operators enables us to use various tools and existing results in the framework of interacting boson theory to study spin-nematic states.

Bi-magnon BEC
We consider the case in which, as the magnetic field h is decreased, bound pairs first close the excitation gap at a critical field h c at the Γ point k = 0.At this critical field, we have ε (0) 2,0 +2h c = 0, while the single particle excitations still have a positive energy gap with ε (0) 1,k + h c > 0 for any k.Following the closure of the excitation energy gap, the bound pairs show the Bose-Einstein condensation (BEC).In this regime, the bound-pair creation operator is written as where N c denotes the number of condensed bound pairs.By using this relation, the spin operators are expressed, for example, in the case of S = 1/2, as where n c denotes the condensate density To describe the dilute Bose gas composed of bimagnons and magnons, as well as bi-magnon BEC, we use the self-consistent Hartree-Fock-Popov (HFP) approximation [43,44], as used in Ref. [45].Within this approximation, we obtain the mean-field Hamiltonian H HFP in a quadratic form.By applying the Bogoliubov transformation, we diagonalize the Hamiltonian as where Here n 1 and n 2 denote the densities given by the thermal averages n 1 = ⟨a † i a i ⟩ and n 2 = ⟨b † i b i ⟩, which are evaluated with the Hamiltonian H HFP at a finite temperature.The bi-magnon and condensate densities, n 2 and n c , satisfy The function n B (E) denotes the Bose distribution function [exp(E/k B T )−1] −1 .Furthermore, the HFP approximation [44] imposes the condition Field dependence of the lowest energy of one-magnon excitations (E 1,k 0 ) and bi-magnon excitations (E2,0).Inside the spin-nematic phase, the slope of the one-magnon gap changes depending on the interaction ratio u/v2.
Consequently, by solving these relations self-consistently, we determine the values of n 1 , n 2 , and n c at a finite temperature.
The one-magnon excitation energy E 1,k undergoes corrections arising from finite densities and interactions.A similar calculation has been done in a one-dimensional system [40].In the condensed phase (n c > 0), E 2,k becomes gapless at k = 0 due to relation (15), and it exhibits a k-linear dispersion relation of the Nambu-Goldstone mode near k = 0 with the velocity v NG = v 2 n c /M , where M is given by ε 2,k ≃ (2M ) −1 k 2 + 2(h − h c ) in the limit of small wave vectors.In the hightemperature non-condensed phase (n c = 0), E 2,k exhibits a positive energy gap, which vanishes at a transition temperature with decreasing temperature.Hereafter, we focus on situations where the dispersion relation E 1,k possesses a positive energy gap ∆ 1 , i.e., E 1,k ≥ ∆ 1 > 0.
The field dependence of the excitation gaps at T = 0 is shown in Fig. 2. In the spin-gap ground state, the onemagnon excitation energy, Eq. (11), results in E 1,k = ε (0) 1,k + h since both n 1 and n 2 are zero.In the spinnematic phase, the energy includes a correction from the finite density n 2 and interactions; at T = 0, the onemagnon density n 1 is found to vanish since one-magnons have a gap.The densities n 2 and n c , which usually satisfy n 2 > n c due to zero-point oscillation, can be calculated from Eqs. ( 14) and ( 15) at zero temperature.In the limit of low density n 2 → 0+ for bi-magnons, as h → h c −, the condensate density approaches the total density, resulting in n c /n 2 → 1.In the case n c = n 2 , one finds ε 2,0 + n 2 v 2 = 0, leading to the field dependence of the one-magnon energy gap expressed as 1,k0 + h c is the one-magnon gap at the critical field h c , and below h c , the slope as a function of h changes to 1−2u/v 2 .Hence, the interaction effect in the spin-nematic phase impedes the rapid closure of the one-magnon gap and may even widen the gap.
In the case of u < v 2 /2, the one-magnon gap closes as the magnetic field decreases in the spin-nematic phase.This closure leads to an additional phase transition, resulting in magnon BEC and a transition to an antiferromagnetic phase.This phenomenon has been observed in a two-dimensional hard-core boson system with correlated hopping [46,47] and a two-dimensional ferromagnetic J 1 -antiferromagnetic J 2 model [12], indicating that these systems are effectively described by the model with u < v 2 /2.Conversely, in the case of u > v 2 /2, the onemagnon gap begins to increase at the critical field.Numerical calculations of the zigzag chain showed a weak gap increase below saturation with decreasing field [11], implying u > v 2 /2.Experimental observations on triplon excitations in SrCu 2 (BO 3 ) 2 , a potential material for spinnematic order, showed an upward turn with a kink near the magnetic field where the magnetization begins to increase [36,48].This suggests an effective description by the model with u > v 2 /2.(The application of our theory to the spin-nematic phase at low magnetic fields in spindimer systems is discussed in Sec.IV.)These results from theoretical model calculations may vary in the presence of interlayer or interchain coupling.

B. Static quantities
The above two-component boson theory describes a spin-nematic state in a magnetic field.The spin expectation values are expressed as at a finite temperature.Notably, the transverse spin components (in the x and y directions) show no dipolar spin order at any temperature.This result does not depend on the form factor g q .In contrast, in the low-temperature condensed phase (n c > 0), the following quadrupolar order appears on bonds as for i ̸ = j, where g(r) = N Λ −1 q g q e −iq•r , and on the same sites for S ≥ 1 as Note that ⟨S x i S x j − S y i S y j ⟩ and ⟨S x i S y j + S y i S x j ⟩ are spinnematic order parameters, capturing a breakdown of the spin U(1) symmetry about the applied field [2].
In the case the bi-magnons have s-wave symmetry, e.g., g(k) = 1/ √ 2, a dominant quadrupolar order appears on the lattice sites.On the other hand, if the bi-magnons exhibit d-wave symmetry, such as g(k) = (cos k x − cos k y )/ √ 2, a dominant quadrupolar order appears on the bonds, resulting in an orthogonal director structure on two-dimensional planes.See Fig. 3 for an illustration.

III. DYNAMICAL QUANTITIES IN THE SPIN-NEMATIC STATE
Our effective theory can describe characteristic behaviors in dynamical quantities of the spin-nematic state.This section presents the results of the dynamical spin structure factor at zero temperature and the NMR relaxation rate at finite temperatures.These analyses are conducted without assuming any symmetry of the form factor of bi-magnons.
A. Dynamical spin structure factor at T = 0 At zero temperature, the dynamical spin structure factor is written as where µ = x, y, z.Here, |n⟩ and E n denote eigenstates and eigenenergies, with n = 0 signifying the lowest energy level.This quantity has been calculated for the field-induced spin-nematic state in the spin-1 bilinearbiquadratic model in Ref. [35], where the computation is performed under the assumption of a solely on-site swave pairing.For comparison, we show the spin-1/2 case with magnon pairing of arbitrary symmetry as follows: The transverse component S xx (k, ω) comprises a gapful one-magnon excitation mode with energy E 1,k and a continuous spectrum of two-magnon states above it.Due to the absence of transverse spin order, this component does not exhibit diverging singularity.Importantly, the intensity of the one-magnon mode depends on the form factor g k of two-magnon bound states and the order parameter of bi-magnon BEC |⟨b † i ⟩| = √ n c .Thus, the intensity contains valuable information about the form factor.In the case of an s-wave symmetric form factor, the intensity of the one-magnon mode appears isotropic in momentum space around k = 0. Otherwise, the intensity displays oscillations in momentum space, reflecting the symmetry of the bound magnon pairs.Hence, this momentum-dependent signal can be a valuable tool for detecting the spin-nematic order and identifying the symmetry of the bound magnon pairs.
The longitudinal component S zz (k, ω) contains, in the lowest energy mode, a gapless k-linear Nambu-Goldstone mode with energy E 2,k , which is gapless only at k = 0. Near the gapless point k = 0, this component behaves as S zz (k, ω) ≃ 2π n c /M v 2 ∥k∥ δ(ω − v NG ∥k∥), and the intensity of this mode vanishes linearly.The loss of intensity behavior at the gapless point comes from the U(1) spin symmetry of the system, and hence, it is common in systems showing magnon BEC (see Sec. V B) and bi-magnon BEC [35].The second term in S zz (k, ω) shows a continuum above the mode, originating from two Nambu-Goldstone bosons created by the operators β † k+q β † −q .We note that this continuum part does not satisfy the exact frequency sum rule lim |k|→0 ∞ 0 ωS zz (k, ω)dω = 0, presumably due to the Hartree-Fock approximation.
Overall, the dynamical spin structure factors do not show any diverging singularity, and the intensity decreases linearly in the vicinity of gapless point k = 0.This behavior has also been observed in quadrupolar phases in the spin-1 bilinear-biquadratic models in zero [5,41,49] and applied magnetic fields [35].

Formula
Dynamical spin correlations at finite temperatures give the NMR relaxation rate 1/T 1 as where A ⊥ k and A ∥ k denote the form factors describing the coupling between nuclear and electronic spins, ω 0 the resonance frequency of nuclear spins, and δS z k = S z k − ⟨S z k ⟩.We neglect the momentum dependence in the form factors A ⊥k and A ∥k , setting |A ⊥k | 2 = c ⊥ and |A ∥k | 2 = c ∥ , and take the limit ω 0 → 0.

Transverse component
The transverse component 1/T 1,⊥ comes from the dynamics of excitations created by the operator S − k in Eq. ( 8).Since the one-magnon excitation created by a † k has an energy gap, it can not receive the small energy ℏω 0 from nuclear spins and does not affect the relaxation.If there is an energy overlap between one-magnon excitations and gapless collective modes, their scattering process, such as those involving integral ∞ −∞ dt⟨β † k+q (t)a q (t)a † q β k+q ⟩, contributes to the relaxation rate [28].Because of the one-magnon excitation gap ∆ 1 , this process shows an exponential decay form exp(−∆ 1 /k B T ) for low temperatures k B T ≪ ∆ 1 .Three magnon terms with the form a † a † a or a † aa in S x k , omitted in Eqs. ( 2) and ( 4), also result in the exponential form in the relaxation rate.

Longitudinal component
The longitudinal component 1/T 1,∥ originates from excitations created by the operator S z k in Eq. ( 5).This operator includes a linear term of bi-magnon operators, contributing to a one-particle propagator in the formula for 1/T 1,∥ .However, since the linear term in terms of the operators β and β † has a very weak coefficient near the contribution from the one-particle propagator in 1/T 1,∥ vanishes in the limit of ω 0 → 0.
Among the various terms in the longitudinal component, the dominant contribution arises from the Raman process induced by the operator β † k+q β q in S z k , similar to the case of antiferromagnets [32].This process leads to a relaxation rate given by where D 2 (E) represents the density of states per site for bosons created by β † .We note that this relaxation rate equation has the same form as Eq. ( 41) derived from the longitudinal component for the canted antiferromagnetic phase.In the long-wavelength approximation, the leading term results in a T 3 temperature dependence, in the low-temperature range ℏω 0 ≪ k B T ≪ k B T SN , reflecting the k-linear behavior of the collective mode.Here T SN denotes the spin-nematic phase transition temperature.a result, the NMR relaxation rate 1/T 1 in the spin-nematic phase exhibits a T 3 dependence at low temperatures.Notably, this finding differs from the previous studies [30,31].This most dominant behavior arises from a term that was previously omitted.
Based on the preceding discussion, the dominant contribution to the NMR relaxation rate in the spin-nematic phase comes from the longitudinal component presented in Eq. ( 23).To further investigate the temperature dependence, we numerically computed Eq. ( 23) on a threedimensional lattice, solving self-consistent equations at finite temperatures while employing a long-wavelength approximation for the density of states.The results of these computations are presented in Fig. 4.
Throughout the low-temperature spin-nematic phase, the evaluated values clearly demonstrate a distinct T 3 temperature dependence.At the transition temperature T SN , the energy gap in Eq. ( 12) closes, indicating the achievement of the critical point as the temperature is decreased from above T SN .However, when the temperature is increased from below T SN , a first-order phase transition occurs at the transition point.Importantly, the evaluated value of 1/T 1 does not exhibit a diverging singularity at the transition temperature, even as the temperature is decreased.Instead, it shows a cusp-like behavior consistent with earlier field-theoretical research [31].This absence of a diverging singularity represents a characteristic feature of the spin-nematic phase.

IV. APPLICATION TO SPIN-DIMER SYSTEMS IN A LOW MAGNETIC FIELD
In this section, we discuss an application of our results to the spin-nematic phase in spin-dimer systems, especially in the low magnetic field regime adjacent to the zero-field spin-gap phase.For this investigation, we use the two-dimensional orthogonal dimer model, known as the Shastry-Sutherland model.In the absence of a  23).The system continuously changes to the critical point (T = TSN) with decreasing the temperature and, just below the critical temperature, the system shows a first-order phase transition to a spin-nematic ordered state.The value does not diverge at T = TSN.
magnetic field, the ground state is the spin-singlet dimer state [50].The excitations in this state comprise triplons, characterized by low mobility [51], and bound triplon pairs, which show high mobility [52,53].Upon applying a magnetic field, triplons with S z = 1 and bound triplon pairs with S z = 2 become prominent low-energy degrees of freedom, and the bound pairs close the energy gap [16,17] at a magnetic field.Hereafter, we focus on these two types of degrees of freedom.Due to the inherent structure of the Shastry-Sutherland lattice, a refinement in the definition of operators a and b is necessary.This lattice contains two types of orthogonal dimer bonds, and each unit cell hosts two distinct dimer bonds.The labeling of spin operators involves the unit cell position index j, the dimer-type index η = A, B in each unit cell, and the site index m = 1, 2 in each dimer, denoted as S α j,η,m (α = x, y, z).A single triplon occupies a dimer bond, and the creation operator for triplons with S z = 1 on an η-type dimer is written as a † k,η .Despite the existence of two branches of triplon energy modes, we focus on the lowest energy mode, represented by a creation operator in the form √ 2 with σ = 1 or −1 depending on the model parameter.Here we assume the low-est mode is non-degenerate.We omit the higher energy mode created by (a † k,A − σa † k,B )/ √ 2. Additionally, we consider the lowest-energy model of the bound triplon pairs with S z = 2, and its creation operator is represented as b with a form factor g q,η,η ′ .In the extended boson Hilbert space, spanned by a † k and b † k , the low-energy effective Hamiltonian retains the same form as in Eq. ( 1), but with the opposite magnetic field dependence: Consequently, the observed results of the excitation spectrum, as depicted in Fig. 2, can applied to the Shastry-Sutherland model by reversing the sign of the magnetic field dependence.In SrCu 2 (BO 3 ) 2 , a material with a Shastry-Sutherland lattice structure sharing many magnetic properties with the Shastry-Sutherland model, the one-triplon excitation energy exhibits an upward trend with a kink near the magnetic field where the magnetization starts to increase [36,48].This suggests that the system enters into a spin-nematic phase at this magnetic field, and the effective couplings for SrCu To evaluate dynamical spin quantities, a refinement in the mapping between spin and boson operators is also necessary.In the dilute regime, the matrix elements of spin operators can be expressed as follows: in the boson Hilbert space.Here, 1 represents the bond operator creating a triplon with S z = 1 on the η-type dimer in the j-th unit cell.These expressions enable the extraction of information about dynamical spin quantities in the Shastry-Sutherland model using interacting boson theory.

V. COMPARISONS WITH CANTED ANTIFERROMAGNETS
To compare with the dynamical quantities in the spinnematic phases, we investigate the dynamical structure factors and the NMR relaxation rate in canted antiferromagnets in a magnetic field within the present theoretical framework.

A. Interacting boson theory of magnons
We consider a spin-S antiferromagnetic system near saturation.The magnons, which are bosonic particles with quantum number S z = −1, arise as excitations above the fully polarized state.To study this system, we employ the interacting boson theory to describe a magnon Bose gas and magnon BEC [45,54,55].We follow the approach used in Ref. [45].This method succeeded in explaining various thermodynamic properties, including magnetization [45], specific heat [56], and thermal Hall conductivity [57].
The Hamiltonian for interacting bosons can be expressed as where ε k denotes the excitation energy, including the Zeeman energy k + h, and a † k (a k ) denotes the bosonic creation (annihilation) operator for a magnon with momentum k.The interaction v q denotes the repulsive potential between two magnons.We omit the momentum dependence of repulsive coupling, setting v q = v.
The spin operators are written in terms of the boson operators as in the dilute limit.In Eq. ( 29), we have expanded the Holstein-Primakoff transformation [58].
Let us assume that the energy spectrum ε k has the lowest energy at a specific momentum k 0 .We consider the case where bosons condense solely at a single momentum k 0 , i.e., 2k 0 lies within the reciprocal lattice space.Above the saturation field h c , the ground state has a positive energy gap ε k0 > 0. This gap closes at h = h c , satisfying the condition ε (0) k0 + h c = 0.In the condensed phase (h < h c ), the interactions lead to BEC of magnons, and the boson operators can be written as where N c is the number of condensed bosons.By treating interactions with the HF approximation and applying the Bogoliubov transformation, we obtain the effective Hamiltonian with where n denotes the particle density per site, n c represents the condensate density n c = N c /N Λ , and θ k satisfies cosh 2θ k = (ε k + 2vn)/E k and sinh 2θ k = vn c /E k .
In the presence of finite condensates (n c > 0), the particle and condensate densities satisfy the relation in the HFP approximation [44], leading to a gapless energy spectrum with E k0 = 0.In the absence of a condensate (n c = 0), there is no similar constraint on n.The values of n and n c also satisfy the self-consistent equation where ⟨• • • ⟩ denotes the thermal average with the Hamiltonian H caf-HF .The particle and condensate densities can be calculated by solving these self-consistent equations at a finite temperature, as in Refs.[45,56].
We set θ = 0 without loss of generality.By using the operators ãk , we can express the spin operators: For S = 1/2, we have where we have omitted three-body terms in S − k .Both spin operators consist of linear and quadratic terms.The transverse spin components align in the x direction, and their expectation values are given by ⟨S 2 )e −ik0•ri and ⟨S y i ⟩ = 0, where the zero-point oscillation from the quadratic terms in Eq. ( 36) has been disregarded.

B. Dynamical structure factors at T = 0
We obtain the dynamical structure factors in the S = 1/2 canted antiferromagnet at zero temperature using this framework.Here, we focus on the one-magnon excitation mode, which behaves as In the low-energy regime, the transverse component S yy (k, ω) shows a diverging singularity at (k, ω) = (k 0 , 0) as S yy (k, ω) and v caf = vn c /m represents the velocity of magnon excitations.
In contrast, the intensities of the low-energy regions in S xx (k, ω) near (k, ω) = (k 0 , 0) and S zz (k, ω) near (k, ω) = (0, 0) decay linearly as S xx (k, ω) ≃ , respectively.The linear decay of intensity in the longitudinal component S zz (k, ω) near the gapless point k = 0 is due to the conservation of S z k in the limit of k = 0 and hence is a common feature observed in systems with spin U(1) symmetry.
Lastly, we note that there have been discussions regarding the decay of magnons in canted antiferromagnets [59].This decay mechanism contributes to the broadening of the one-magnon excitation peak.However, the treatment of this issue needs further study beyond the HF approximation and hence lies beyond the scope of the present study.

C. NMR relaxation rate
Moriya argued that the two-magnon Raman process contributes most dominantly to the NMR relaxation rate in zero-field antiferromagnets [32,33].We derive the same result for canted antiferromagnets using the interacting boson theory described above.Our conclusion is summarized in Table I.TABLE I. Temperature dependence of the NMR relaxation rate in the field-induced spin-nematic (SN) phase and the canted antiferromagnetic (CAF) phase.The NMR relaxation rate consists of transverse and longitudinal components corresponding to the first and second terms in Eq. ( 22), respectively.Here, ∆1 denotes the gap of one-magnon excitations in the spin-nematic state.The temperature is within the range of ℏω0 ≪ kBT ≪ kBTc, where Tc denotes the transition temperature.For the CAF phase, the behavior depends on the spin anisotropy: (a) when the spin anisotropy gap ∆aniso exceeds the NMR frequency ω0, ω0 < ∆aniso, and (b) either in the presence of perfect U(1) spin symmetry or when ∆aniso < ω0.The operator S z k shown in Eq. ( 37) includes a linear term with operators a and a † .The contribution from this term to 1/T 1,∥ vanishes in the limit of ω 0 → 0, similar to the extinction of the effect of the bi-magnon propagator in the longitudinal component in the spin-nematic phase.Among the various operators present in S z k , the quadratic term with the operator α † k+q α q contributes to the Raman process.After the time integration, only the Raman process remains non-zero in the calculation of Here D(E) represents the density of states per site.This equation has the same form as the well-known result obtained by Moriya using the spin-wave expansion [33], except that n c is self-consistently determined in our approach.Hence, the longitudinal part 1/T 1,∥ shows a T 3 temperature dependence at low temperatures.

Transverse component
The transverse component 1/T 1,⊥ in Eq. ( 22) comes from propagations of excitations created by the operators in S − k [Eq.(36)].It has been argued that the one-magnon propagator gives 1/T 1,∥ a T -linear temperature dependence at low temperatures [55].This result can be seen in our calculations; since the coupling of the magnon creation and annihilation operators in S y k has a diverging singularity at k = k 0 as the contribution from the one-magnon propagator to 1/T 1,⊥ survives and results in a T -linear temperature dependence even in the limit of ω 0 → 0. However, in real materials, this effect is typically eliminated by a spin anisotropy gap [32].
To demonstrate the spin anisotropy dependence of the T -linear term in 1/T 1,⊥ , we maintain a small but finite value for ω 0 instead of taking the limit ω 0 → 0. We introduce a weak uniaxial anisotropy along the spin x axis, which is added to the Hamiltonian H caf as the term η k ε k (a † k a † −k + h.c.), with η being a small positive real constant.In the ordered phase, the energy is minimized with θ = 0.In the HFP approximation, the relation between n and n c in Eq. ( 35) is replaced by (1 + 2η)ε k0 + (2n − n c )v = 0. Consequently, the energy spectrum of magnons is given by E aniso (k) = (ε k + 2vn) 2 − (vn c + 2ηε k ) 2 , and it has an energy gap ∆ aniso = 2η(h c − h)vn c , which closes at the critical temperature due to n c = 0. Using these results, we find that the contribution from the one-magnon propagator to 1/T 1,⊥ behaves as for ℏω 0 ≪ k B T , where D aniso (ω) denotes the density of states with the energy spectrum E aniso (k) and ).Though n and n c implicitly depend on the temperature, they can be well approximated at low temperatures with the values at zero temperature.Hence, the dominant temperature dependence in Eq. ( 42) is T -linear at low temperatures.
Note that ω 0 is typically much smaller than the energy scales of spin interactions.Only in cases where ω 0 > ∆ aniso , does Eq.( 42) exhibit a T -linear behavior.In many magnetic materials, where ω 0 < ∆ aniso , the term in Eq. ( 42) disappears due to D aniso (ω 0 ) = 0. Consequently, no T -linear behavior is observed in nearly the entire temperature range of the canted antiferromagnetic phase.In such cases, the operator α † k+q α k0+q in S − k , which gives the Raman process, results in the most dominant T 3 temperature dependence in 1/T 1,⊥ .
Even in the presence of anisotropy, the critical divergence of 1/T 1 arises from the one-particle propagator.The anisotropy gap closes near the critical region, including the quantum critical point at T = 0 and the phase transition line at finite temperatures.Consequently, the influence of the one-magnon propagator becomes pronounced in the critical region, contributing to the diverging singularity at T = T c shown in Ref. [60].

VI. SUMMARY AND DISCUSSIONS
We have successfully formulated an interacting boson theory that incorporates both magnon and bi-magnon degrees of freedom, establishing a comprehensive framework for describing spin-nematic states.Our theory accommodates a wide range of spin-nematic states, encompassing site to bond spin-nematic states.The characteristics of the spin-nematic order structure are embedded into the form factor of the bound bi-magnons.By utilizing well-established methods in the field of interacting boson theory, our versatile method enables the calculation of various physical quantities.In this paper, we applied the self-consistent HFP approximation and successfully obtained the dynamical quantities in the spinnematic phases.
The dynamical quantities in the spin-nematic phases exhibit distinct behaviors compared to canted antiferromagnets.The dynamical spin structure factors in the spin-nematic phases do not show diverging singularities at any combination of momentum and frequency.This behavior, which has also been previously observed in some specific models [5,30,35,41,49], is a general characteristic of spin-nematic phases.Furthermore, the transverse component of the dynamical structure factor displays a gapful structure, indicating the presence of a finite binding energy for bi-magnons.Notably, the intensity of one-magnon excitations in this transverse component is influenced by the form factor governing the structure of bound magnon pairs.As a result, the intensity can exhibit oscillatory patterns around k = 0, providing information about the symmetry characteristics, such as d-wave, of the spin-nematic order parameter.
The NMR relaxation rate 1/T 1 shows another characteristic in the spin-nematic phases.Specifically, the relaxation rate does not show any diverging singularity; instead it shows a cusp at the critical temperature of these phases.This observation aligns with a previous fieldtheoretical argument [31].Moreover, through our reanalysis of the low-temperature behavior of the relaxation rate, we revealed a previously unreported T 3 temperature dependence to be the most dominant one in 1/T 1 .This behavior is similar to that observed in canted antiferromagnets with weak spin anisotropy, which is common in typical compounds.Notably, in systems with perfect U(1) spin symmetry, canted antiferromagnets show a T -linear temperature dependence at low temperatures, while the spin-nematic phases maintain the T 3 behavior.
Among the obtained results concerning the dynamical quantities, only the intensity of one-magnon excitation mode in the dynamical structure factor explicitly depends on the form factor of the bound magnon pairs.Conversely, the remaining results exhibit universal features of the spin-nematic states, indicating their independence from the symmetry of the magnon pairs.It would be worth future work to explore additional quantities beyond those analyzed in this study that hold the potential to reveal and identify magnon-pair symmetry.
In this paper, we focused on the case where the wavenumber vector of spin-nematic order is given as k SN = 0. Now, we discuss the case where the wavevector is non-zero.If 2k SN points on one of the vertices of the reciprocal lattice, bi-magnons only condense at one wavenumber, making our discussion mostly applicable.In this case, the excitation energies take the same form as in Eq. ( 10), but the HFP condition, Eq. ( 15), changes to ε 2,kSN + (2n 2 − n c )v 2 + n 1 u = 0, and the bi-magnon excitation energy E 2,k has a gapless point at wavenumber k SN .However, the temperature dependence of the NMR relaxation rate remains unchanged, indicating that the results shown in Table I exhibit universal behaviors in spin-nematic ordered states.If the wavenumber k SN is incommensurate, as observed in the J 1 -J 2 zigzag chain [39], the discussion goes beyond the analysis in this paper since bi-magnons can condense at two wavenumbers ±k SN .
Lastly, it is crucial to consider spin anisotropy when comparing with experiments on spin-nematic states.In systems with spin anisotropy allowed by crystal structures, such as Dzyaloshinskii-Moriya interactions, the symmetry breaking of spin-nematic order is either discrete or absent, as discussed in Refs.[61,62].This significantly influences the critical properties of the spinnematic phase transition.Excitations have an anisotropy gap, and the temperature dependence of thermodynamic quantities varies at very low temperatures.When no symmetry is broken, the system undergoes a crossover from a paramagnetic state to a spin-nematic state, and the anisotropy gap remains open, even near the crossover temperature.Whether symmetry breaking occurs or not needs to be considered for each individual crystalline structure of the material.

FIG. 3 .
FIG. 3. Illustration of quadrupolar orders on two-dimensional planes.The directors represent the arrangements of quadrupolar moments.(a) On-site ferro-quadrupolar moments are formed by s-wave magnon pairing.(b) Antiferroquadrupolar moments are formed on bonds through d-wave paring.

1 FIG. 4 .
FIG. 4. Temperature dependence of the NMR relaxation rate 1/T1 near the transition temperature TSN of a spin-nematic phase in (a) logarithmic scales and (b) linear scales.The plot is evaluated from Eq. (23).The system continuously changes to the critical point (T = TSN) with decreasing the temperature and, just below the critical temperature, the system shows a first-order phase transition to a spin-nematic ordered state.The value does not diverge at T = TSN.