compat=1.1.0 /tikzfeynman/momentum/arrow shorten = 0.3 /tikzfeynman/warn
luatex = false

# Elastic signatures of a spin-nematic

Monica Bapna<EMAIL_ADDRESS>International Centre for Theoretical
Sciences, Tata Institute of Fundamental Research, Bangalore 560089, India
Subhro Bhattacharjee<EMAIL_ADDRESS>International Centre for Theoretical
Sciences, Tata Institute of Fundamental Research, Bangalore 560089, India

###### Abstract

We study the elastic signatures– renormalisation of sound velocity and
magnetostriction – of the spin-nematic phase of a spin-$1$ magnet on a
triangular lattice described by the bilinear-biquadratic spin Hamiltonian. We
show that at low temperatures, the scattering of the acoustic phonons from the
Goldstone modes of the nematic phase lead to a powerlaw renormalisation of the
fractional change in the sound velocity, $v$, as a function of temperature,
$T$, i.e. $\Delta v/v\propto T^{3}$ as opposed to the same in the high
temperature paramagnet where $\Delta v/v\propto T^{-1}$. At the generically
discontinuous- nematic transition, there is a jump in magnetostriction as well
as $\Delta v/v$ along with enhanced $h^{4}$ dependence on the magnetic field,
$h$, near the nematic transition. These signatures can help positively
characterise the spin-nematic in general and in particular the realisation of
such a phase in the candidate material NiGa2S4.

##  I Introduction

Interplay of symmetries and competing interactions can stabilise a plethora of
magnetic phases in spin systems with inconclusive experimental signatures for
conventional probes. This not only include issues of finding smoking-gun
signatures of fractionalised quasi-particles in quantum spin liquids, but a
much broader context pertaining to many other unconventional phases such as
higher (than dipole) multi-pole orders in a variety of candidate systems with
spin moments $S>1/2$.Santini and Amoretti (2000); Chandra  _et al._ (2002);
Sato  _et al._ (2012); Patri  _et al._ (2019); Harter  _et al._ (2017); Voleti
_et al._ (2020) A possible resolution to the above conundrum is to invoke a
combination of probes to gather complementary insights. An important class of
such experimental probes like vibrational Raman and infrared scatterings as
well as ultrasonic spectroscopy aim to exploit the ubiquitous magnetoelastic
coupling to reveal the properties of the unconventional magnetic ground states
and low energy excitations via the phonons.Chapon  _et al._ (2004);
Bhattacharjee  _et al._ (2011, 2016); Seth  _et al._ (2022)

In this paper, we shall develop the theory of magnetoelastic coupling for the
rather elusive spin-nematic phases Blume and Hsieh (1969); Chen and Levy
(1971); Matveev (1973); Chandra and Coleman (1991) in spin-1 magnets and apply
it to predict its elastic signatures. In a spin-nematic, the ground state has
no magnetic dipole moment but has a finite expectation value for the the
quadrupole moment which is a bilinear of spins. Here, we shall consider only
on-site spin-nematic characterised by the symmetric traceless
operatorBhattacharjee  _et al._ (2006); Penc and Läuchli (2011) :

$Q_{i}^{\mu\nu}=\frac{(S_{i}^{\mu}S_{i}^{\nu}+S_{i}^{\nu}S_{i}^{\mu})}{2}-\frac{S(S+1)\delta^{\mu\nu}}{3},~{}$
(1)

where $S^{\alpha}_{i}$ ($\alpha=x,y,z$) are spin-1 operators at lattice site
$i$. Such order has been proposed for the triangular lattice magnet NiGa2S4
Nakatsuji  _et al._ (2005); Tsunetsugu and Arikawa (2006); Bhattacharjee  _et
al._ (2006); Läuchli  _et al._ (2006); Stoudenmire  _et al._ (2009) where it
is stabilized by a sizeable biquadratic term $\sim({\bf S}_{i}\cdot{\bf
S}_{j})^{2}$ that can arise from spin-lattice coupling.Kittel (1960); Barma
(1975); Bhattacharjee  _et al._ (2006)

${\text{NiGa}_{2}\text{S}_{4}}$ is a layered material where $\text{Ni}^{2+}$
forms an isotropic triangular lattice with $S=1$ at each site. The system
fails to show any conventional long range magnetic order in neutron scattering
experiments to the lowest temperature measured (Curie Weiss temperature,
$\theta_{W}=-80K$).Nakatsuji  _et al._ (2005, 2007) However the state below
$T\sim 10$ K shows $\sim T^{2}$ magnetic specific heat and constant magnetic
susceptibility indicating the presence of low energy linearly dispersing
excitations. It was subsequently proposed that this spin-1 triangular lattice
magnet could possibly have spin ferronematic Bhattacharjee  _et al._ (2006);
Läuchli  _et al._ (2006); Stoudenmire  _et al._ (2009) or three sublattice
nematicTsunetsugu and Arikawa (2006) ordering (see Fig. 1). This spin-nematic
state forms the right starting point Stoudenmire  _et al._ (2009) to
understand the relevance of the third nearest neighbour Heisenberg exchange
Pradines  _et al._ (2018); Stock  _et al._ (2010); Mazin (2007); Takubo  _et
al._ (2007) as well as spin freezing below $10$ K Nambu  _et al._ (2015)– both
relevant for the material. However, the most pertinent experiments for our
work are the recent Raman measurements probing the phonons in NiGa2S4
Valentine  _et al._ (2020) which shows substantial spin-lattice coupling. We
build on the above experimental indication of sizeable magnetoelastic coupling
to show that the possible spin-nematic can have a strong signature in the
elastic sector – namely the strain and the sound velocity – which can form
further experimental probes to such unconventional spin-quadrupolar order.
These lattice signatures can be measured very accurately and have recently
proved very useful in obtaining information about unconventional magnetic
phases and phase transitions.Bhattacharjee  _et al._ (2011); Tang  _et al._
(2022)

Being bilinear in spins, the spin-nematic is time-reversal (TR) invariant,
but, breaks the spin-rotation symmetry and has been dubbed as moment-free
magnetism.Chandra and Coleman (1991) Unlike dipolar ordering, such quadrupolar
ordering is hard to detect via neutron scattering.Barzykin and Gor’kov (1993)
However, the same TR even order parameter is expected to couple strongly to
the lattice vibrations and hence provides a possible way to probe it. Further,
even in absence of static dipole moment, the breaking of spin-rotation
symmetry lead to gapless Goldstone modes – spin-nematic waves. The coupling of
such gapless modes to the acoustic phonons further gives a way to detect the
former in ultrasound experiments.

Here, we report the effect of magnetoelastic coupling in a uniaxial spin-
nematic by deriving the microscopics of such coupling and in particular the
coupling of the acoustic phonons to the nematic Goldstone modes to obtain the
renormalisation of sound speed at low temperatures deep inside the spin-
nematic phase. We complement the microscopics approach with a phenomenological
Landau-Ginzburg theory for the long wavelength dipolar, quadrupolar and strain
modes to capture the effect of the thermal transition out of the spin-nematic
phase on the magnetostriction and sound speed renormalisation. While we focus
on NiGa2S4 Nakatsuji  _et al._ (2005, 2007), our results are easily
generalised to other cases of spin-nematic order.

The rest of this paper is organised as follows. In section II, we introduce
the microscopic setting of magnetoelastic coupling (Eq. 7) in spin-1 magnets
on triangular lattice with bilinear-biquadratic exchanges (Eq. 2) and
summarize its different phases– motivated by the low energy physics of
NiGa2S4– and discuss the effect of such coupling in fractional change of sound
speed (Eq. 13). The above formalism is used to calculate the coupling between
the nematic Goldstone modes and the acoustic phonons and hence the temperature
dependence of fractional change in sound speed deep inside the both the
ferronematic and the three sublattice nematic phase in Section III. Typically
the temperature dependence of the sound speed is given by a powerlaw, $\Delta
v/v\sim T^{a}$ where $a=3$ (Eq. 24) for a ferronematic and can vary between
$1-3$ (Eq. 32) for the three sublattice nematic depending on the temperature
range and the ratio of bilinear and biquadratic coupling $J/K$ (see Eq. 2).
This is unlike the case of the thermal paramagnet where $\Delta v/v\propto
1/T$ (Eq. 36). Complementary to the microscopics, we study the
phenomenological Landau-Ginzburg theory for the long wavelength dipolar,
quadrupolar and strain modes in Section IV. We find the symmetry allowed
irreducible representations to obtain the mean field free energy which we use
in Section V to examine the elastic signatures – fractional change in length
(magnetostriction), $\Delta L/L$, and $\Delta v/v$ – of the dipolar and
quadrupolar ordering via symmetry allowed magnetoelastic couplings. We find
that magnetic(nematic) ordering would lead to continuous change(jumps) in
$\Delta L/L$ and $\Delta v/v$ along with an enhanced $h^{4}$ dependence on the
magnetic field, $h$, near the nematic transitions. Observing jumps in
fractional change in length and fractional change in sound speed in the
absence of any magnetization would strongly favour the case for a spin-
nematic. Various details are summarised in appendices.

##  II Spin-lattice coupling in Spin-1 magnet

The starting point is the nearest neighbour minimal spin-1 bilinear-
biquadratic model on the triangular lattice Bhattacharjee  _et al._ (2006);
Penc and Läuchli (2011) given by

$H=J\sum_{\langle ij\rangle}(\mathbf{S}_{i}\cdot\mathbf{S}_{j})+K\sum_{\langle
ij\rangle}(\mathbf{S}_{i}\cdot\mathbf{S}_{j})^{2}$ (2)

where in addition to the usual (first) Heisenberg term, we also have the
(second) biquadratic spin exchanges. These higher order exchanges can be
obtained from an underlying higher energy multi-orbital Hubbard model.Fazekas
(1999) Notably, the biquadratic term can be further renormalised by
magnetoelastic coupling by integrating out the phonons Bhattacharjee  _et al._
(2006); Tchernyshyov  _et al._ (2002) whose importance is evident in the
recent Raman scattering experiments.Valentine  _et al._ (2020) Both virtual
hopping and phonon effects naturally give rise to $K<0$ while $J>0$ from the
former.

It is useful to re-write the above Hamiltonian, up to a constant, using the
spin-$1$ operator identity Penc and Läuchli (2011) :
$(\mathbf{S}_{i}\cdot\mathbf{S}_{j})^{2}=\frac{1}{3}S^{2}(S+1)^{2}-\frac{1}{2}(\mathbf{S}_{i}\cdot\mathbf{S}_{j})+\sum_{\mu,\nu}Q_{i}^{\mu\nu}Q_{j}^{\mu\nu}$
as $H=\left(J-\frac{K}{2}\right)\sum_{\langle
ij\rangle}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+K\sum_{\langle
ij\rangle;\mu,\nu}Q^{\mu\nu}_{i}Q^{\mu\nu}_{j}$. Eq. 2 is clearly invariant
under global spin rotations as well as various symmetries of the triangular
lattice and TR. The ground states, however, spontaneously break different
symmetries depending on the ratio of $K/J$.

Most important to us is the large $K/J$ regime where an uniaxial uniform
(ferro) spin-nematic is stabilised Bhattacharjee  _et al._ (2006) for $K<0$.
For spin-1, an on-site uniaxial nematic state is stabilized (this is clear
from the discussion on the spin-1 wave functions in Refs. Läuchli  _et al._ ,
2006; Penc and Läuchli, 2011; Tóth, 2011 as summarized in Appendix A) where
the director of the nematic is uniform at all the lattice sites as shown in
Fig. 1. The order parameter that characterise such a ferronematic order is

$\langle{Q}^{\mu\nu}\rangle=\mathcal{Q}_{FN}\bigg{(}n^{\mu}n^{\nu}-\frac{1}{3}\delta^{\mu\nu}\bigg{)}~{}.$
(3)

where $\mathcal{Q}_{FN}$ and $\hat{\mathbf{n}}$ are respectively the magnitude
and director of the ferronematic. Indeed, for $J>0$ and $K<0$, the
ferronematic (spiral) phase is stable for $|K|/J>2(<2)$ within mean-field
analysis.Stoudenmire  _et al._ (2009)

Figure 1: Schematic : The uniaxial ferro (Eq. 3) and the three sublattice (Eq.
97) spin-nematic orders on the triangular lattice. The directors are shown in
deep blue.

While the microscopic mechanisms as discussed above, favours $K<0$, it is
useful to consider the case of $K>0$, whence for large $K/J$ a three
sublattice nematic is stabilised in the triangular lattice where the directors
of the spin-nematic are orthogonal to each other on the three sublattices
Tsunetsugu and Arikawa (2006) as shown in Fig. 1. The relevant mean-field
wave-functions for the three sublattice nematic orders are briefly summarised
in Appendix A.

Turning the the limit of $J/|K|\gg 1$, on a triangular lattice, for $J>0$, a
natural competing (with the above two nematic) phase is the 120∘ coplanar
spiral with a non-zero spin expectation given by

$\displaystyle\langle\mathbf{S}_{i}\rangle$
$\displaystyle=m\hat{\mathbf{s}}_{i}~{}~{}{\rm
with}~{}~{}\hat{\mathbf{s}}_{i}=\text{cos}(\mathbf{q}\cdot\mathbf{r}_{i})\hat{\mathbf{e}}_{x}+\text{sin}(\mathbf{q}\cdot\mathbf{r}_{i})\hat{\mathbf{e}}_{y}$
(4)

where $m$ is the magnitude of magnetization,
$\mathbf{q}=\frac{2\pi}{3}\hat{\mathbf{x}}+\frac{2\pi}{\sqrt{3}}\hat{\mathbf{y}}$
is the spiral wave vector and $\hat{\mathbf{e}}_{x}$ and
$\hat{\mathbf{e}}_{y}$ are orthogonal unit vectors in the plane of spin
ordering. Note that the spin ordering induces a parasitic quadrupole moment,
$\mathcal{Q}_{par}\sim m^{2}$, which should be distinguished from the pure
spin-nematic ordering discussed above.

###  II.1 Spin-phonon coupling

The Raman scattering experiments Valentine  _et al._ (2020) indicate the
presence of substantial magnetoelastic coupling in NiGa2S4, possibly arising
from the modulation of the spin exchange coupling constants ($J$ and $K$) in
the Hamiltonian in Eq. 2 by the phonons. This is obtained via Taylor expansion
of the exchange constants in lattice displacements about their equilibrium
positions as Bhattacharjee  _et al._ (2011)

$\displaystyle J_{ij}$ $\displaystyle=$ $\displaystyle J_{0}+\frac{\partial
J_{ij}}{\partial\mathbf{R}_{ij}}\cdot\mathbf{R}_{ij}+\frac{1}{2}\mathbf{R}_{ij}\cdot\frac{\partial^{2}J_{ij}}{\partial\mathbf{R}_{ij}^{2}}\cdot\mathbf{R}_{ij}$
$\displaystyle K_{ij}$ $\displaystyle=$ $\displaystyle K_{0}+\frac{\partial
K_{ij}}{\partial\mathbf{R}_{ij}}\cdot\mathbf{R}_{ij}+\frac{1}{2}\mathbf{R}_{ij}\cdot\frac{\partial^{2}K_{ij}}{\partial\mathbf{R}_{ij}^{2}}\cdot\mathbf{R}_{ij}~{},$
(5)

where $\mathbf{R}_{ij}=\mathbf{R}_{i}-\mathbf{R}_{j}$ and $\mathbf{R}_{i}$ is
the displacement of site $i$ from its equilibrium position
$\mathbf{R}_{i}^{0}$. Eq. 2 then becomes

$\displaystyle H$ $\displaystyle=H_{sp}+H_{sp-ph}$ (6)

where $H_{sp}$ is the spin Hamiltonian in Eq. 2 with the exchanges now being
given by the equilibrium values (first term in Eq. 5) and $H_{sp-ph}$ being
the spin-phonon coupling Hamiltonian of the form,

$\displaystyle H_{sp-ph}$ $\displaystyle=$ $\displaystyle H_{1}+H_{2},$ (7)

where

$\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{k}}U_{\mathbf{k}}^{(1)}A_{\mathbf{k}}+\sum_{\mathbf{k}}\widetilde{U}_{\mathbf{k}}^{(1)}A_{\mathbf{k}}$
(8) $\displaystyle H_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\mathbf{k},\mathbf{k}^{\prime}}U_{\mathbf{k},\mathbf{k}^{\prime}}^{(2)}A_{\mathbf{k}}A_{-\mathbf{k}^{\prime}}+\frac{1}{2}\sum_{\mathbf{k},\mathbf{k}^{\prime}}\widetilde{U}_{\mathbf{k},\mathbf{k}^{\prime}}^{(2)}A_{\mathbf{k}}A_{-\mathbf{k}^{\prime}}$
(9)

are respectively linear and quadratic in lattice displacement operator given
by $A_{\mathbf{k}}=a_{\mathbf{k}}+a_{-\mathbf{k}}^{\dagger}$ with
$a^{\dagger}_{\mathbf{k}}$ being the bosonic phonon creation operator.
$U_{\mathbf{k}}^{(1)}$, $\widetilde{U}_{\mathbf{k}}^{(1)}$,
$U_{\mathbf{k},\mathbf{k}^{\prime}}^{(2)}$ and
$\widetilde{U}_{\mathbf{k},\mathbf{k}^{\prime}}^{(2)}$ depend on the spin
operators and are given by

$\displaystyle U_{\mathbf{k}}^{(1)}=\sum_{\langle
ij\rangle}(\mathbf{S}_{i}\cdot\mathbf{S}_{j})\frac{(e^{i\mathbf{k}\cdot{\mathbf{R}_{i}^{0}}}-e^{i\mathbf{k}\cdot{\mathbf{R}_{j}^{0}}})}{\sqrt{2MN\omega_{0,\mathbf{k}}^{ph}}}\left(\frac{\partial
J_{ij}}{\partial\mathbf{R}_{ij}}\cdot\mathbf{e}_{\mathbf{k}}\right)$
$\displaystyle\widetilde{U}_{\mathbf{k}}^{(1)}=\sum_{\langle
ij\rangle}(\mathbf{S}_{i}\cdot\mathbf{S}_{j})^{2}\frac{(e^{i\mathbf{k}\cdot{\mathbf{R}_{i}^{0}}}-e^{i\mathbf{k}\cdot{\mathbf{R}_{j}^{0}}})}{\sqrt{2MN\omega_{0,\mathbf{k}}^{ph}}}\left(\frac{\partial
K_{ij}}{\partial\mathbf{R}_{ij}}\cdot\mathbf{e}_{\mathbf{k}}\right)$ (10)

and

$\displaystyle U_{\mathbf{k},\mathbf{k}^{\prime}}^{(2)}=\sum_{\langle
ij\rangle}\frac{(\mathbf{S}_{i}\cdot\mathbf{S}_{j})}{2MN\sqrt{\omega_{0,\mathbf{k}}^{ph}\omega_{0,-\mathbf{k}^{\prime}}^{ph}}}(e^{i\mathbf{k}\cdot{\mathbf{R}_{i}^{0}}}-e^{i\mathbf{k}\cdot{\mathbf{R}_{j}^{0}}})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\left(\mathbf{e}_{-\mathbf{k}^{\prime}}\cdot\frac{\partial^{2}J_{ij}}{\partial{\mathbf{R}_{ij}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}\right)(e^{-i\mathbf{k}^{\prime}\cdot{\mathbf{R}_{i}^{0}}}-e^{-i\mathbf{k}^{\prime}\cdot{\mathbf{R}_{j}^{0}}})$
$\displaystyle\widetilde{U}_{\mathbf{k},\mathbf{k}^{\prime}}^{(2)}=\sum_{\langle
ij\rangle}\frac{(\mathbf{S}_{i}\cdot\mathbf{S}_{j})^{2}}{2MN\sqrt{\omega_{0,\mathbf{k}}^{ph}\omega_{0,-\mathbf{k}^{\prime}}^{ph}}}(e^{i\mathbf{k}\cdot{\mathbf{R}_{i}^{0}}}-e^{i\mathbf{k}\cdot{\mathbf{R}_{j}^{0}}})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\left(\mathbf{e}_{-\mathbf{k}^{\prime}}\cdot\frac{\partial^{2}K_{ij}}{\partial{\mathbf{R}_{ij}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}\right)(e^{-i\mathbf{k}^{\prime}\cdot{\mathbf{R}_{i}^{0}}}-e^{-i\mathbf{k}^{\prime}\cdot{\mathbf{R}_{j}^{0}}})$
(11)

where $M$ and $N$ are the mass and number of a nickel ions respectively,
$\mathbf{e}_{\mathbf{k}}$ is the phonon polarization vector,
$\omega_{0,\mathbf{k}}^{ph}$ is the bare phonon frequency corresponding to the
Harmonic phonon Hamiltonian,

$\displaystyle H_{ph}=$
$\displaystyle\sum_{\mathbf{k}}\omega_{0,\mathbf{k}}^{ph}a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}$
(12)

For acoustic phonons at long wavelengths (${\bf k}\rightarrow 0$) we have
$\omega_{0,\mathbf{k}}^{ph}\approx v_{0}|\mathbf{k}|$ where $v_{0}$ is the
bare sound speed. Above and in the rest of this work, we have set $\hbar=1$.

###  II.2 Fractional Change in the sound speed

The effect of the spin-phonon coupling on the fractional change in sound speed
is obtained from the real part of the phonon self-energy,
$\Sigma(\mathbf{q},\omega)$, as,Bhattacharjee  _et al._ (2011)

$\frac{\Delta
v}{v}=\lim_{\mathbf{q}\xrightarrow{}0}\frac{\Delta\omega_{\mathbf{q}}^{ph}}{\omega_{0,\mathbf{q}}^{ph}}\approx\lim_{\mathbf{q}\xrightarrow{}0}\frac{Re\Sigma(\mathbf{q},\omega_{0,\mathbf{q}}^{ph})}{\omega_{0,\mathbf{q}}^{ph}}.$
(13)

The phonon self-energy due to the interactions with the spins, in turn, can be
obtained from the phonon propagator
$G^{ph}(\mathbf{q},\tau-\tau^{\prime})=-\langle
T_{\tau}(A_{\mathbf{q}}(\tau)A_{-\mathbf{q}}(\tau^{\prime}))\rangle$ in
Matsubara space and is given by the Dyson equation

$\displaystyle
G^{ph}(\mathbf{q},i\Omega_{n})=\frac{1}{\left(G^{ph}_{0}(\mathbf{q},i\Omega_{n})\right)^{-1}-\Sigma(\mathbf{q},i\Omega_{n})}$
(14)

where $\Omega_{n}$ are the bosonic Matsubara frequencies and the bare phonon
propagator given by

$\displaystyle
G^{ph}_{0}(\mathbf{q},i\Omega_{n})=\frac{2\omega_{0,\mathbf{q}}^{ph}}{(i\Omega_{n})^{2}-(\omega_{0,\mathbf{q}}^{ph})^{2}}.$
(15)

The phonon self energy due to the magnetoelastic coupling would depend on
correlations of the spins and hence fractional change in sound speed can be
used to probe the spin physics. In the following sections we use this
formalism to calculate the temperature dependence of $\Delta v/v$ in the
ferronematic, three sublattice nematic and the thermal paramagnet.

##  III Fractional change in sound speed in spin-nematic phase

Deep inside the spin-nematic phase – both ferro and three sublattice, the
renormalisation of the sound speed is brought about by the interaction between
the acoustic phonons and the spin-nematic waves Barma (1974) via the coupling
given by Eq. 7. To this end we use the spin-nematic wave theory for the
ferronematic (Matveev et. al.Matveev (1973)) and three sublattice nematic
(Tsunetsugu et. al.Tsunetsugu and Arikawa (2006)) to write the spins in terms
of the low energy Goldstone bosons of the respective nematic orders –
summarised in Appendix C for completeness.

###  III.1 The Ferronematic phase

To get the temperature dependence of $\Delta v/v$ (Eq. 13) in the
ferronematic, we calculate phonon self energy starting with the following
Hamiltonian for the phonon and linear spin-nematic waves

$\displaystyle\mathcal{H}_{FN}$ $\displaystyle=H_{ph}+H_{sp,FN}+H_{sp-
ph,FN}~{},$ (16)

obtained from Eqs. 6 and 12. Here $H_{ph}$ is the harmonic phonon Hamiltonian
(Eq. 12), $H_{sp,FN}$ and $H_{sp-ph,FN}$ are respectively the linear spin-
nematic wave Hamiltonian for the ferronematic phase and the phonon-spin-
nematic wave coupling Hamiltonian respectively whose form we now discuss.

$H_{sp,FN}$ can be obtained from Eq. 2 by expressing the spins in terms of
boson operators that create spin-nematic wave.Matveev (1973) This is done for
the ferronematic state by introducing two bosons at every site $i$, with
creation operators given by $b_{i1}^{\dagger}$ and $b_{i2}^{\dagger}$ that
capture deviations from the mean field ferronematic ground state. For a mean
field state with the director along the $\hat{\bf z}$-direction, the wave
function is $\otimes_{i}\ket{0}_{i}$ (Appendix A) and
$b_{i1}^{\dagger}\ket{0}_{i}=\ket{1}_{i},b_{i1}\ket{1}_{i}=\ket{0}_{i},b_{i2}^{\dagger}\ket{0}_{i}=\ket{\bar{1}}_{i}$
and $b_{i2}\ket{\bar{1}}_{i}=\ket{0}_{i}$. The relation between the spin
operators and the above bosons are given in Appendix C.1. Expressing the spin
operators in terms of the bosons in Eq. 2 and approximating to the harmonic
order we get

$\displaystyle H_{sp,FN}$
$\displaystyle=\sum_{\mathbf{k}}\omega^{s}_{\mathbf{k}}\psi_{\mathbf{k}}^{\dagger}\psi_{\mathbf{k}}~{},$
(17)

with $\psi_{\mathbf{k}}=(d_{\mathbf{k},1},d^{\dagger}_{-\mathbf{k},2})$ where
$d_{\mathbf{k},1}$ and $d_{\mathbf{k},2}$ are related to
$(b_{\mathbf{k},1},b_{\mathbf{k},2})$ via Bogoliubov transformation (Eq. 101).
The dispersion is given by

$\displaystyle\omega^{s}_{\mathbf{k}}$
$\displaystyle=6K_{0}\sqrt{(1-\gamma_{\mathbf{k}})\left(1+\gamma_{\mathbf{k}}-\frac{2J_{0}}{K_{0}}\gamma_{\mathbf{k}}\right)}~{},$
(18)

where
$\gamma_{\mathbf{k}}=\frac{1}{6}\sum_{\bm{\delta}}e^{i\mathbf{k}\cdot\bm{\delta}}$
and $\bm{\delta}$ is the distance to the six nearest neighbours. In the long
wavelength limit, $\omega^{s}_{\mathbf{k}}\approx\bar{c}_{s}|\mathbf{k}|$ with
$\bar{c}_{s}$ being the ferro spin-nematic wave speed.

The coupling between the phonons and the ferro spin-nematic waves is obtained
by using the same bosonic representation of the spin operators in Eq. 7. The
resultant Hamiltonian is given by

$\displaystyle H_{sp-ph,FN}=H_{1,FN}+H_{2,FN}$ (19)

where

$\displaystyle H_{1,FN}=$
$\displaystyle\sum_{\mathbf{k},\mathbf{q_{2}}}\psi^{\dagger}_{\mathbf{k}+\mathbf{q_{2}}}\mathcal{M}_{\mathbf{k},\mathbf{q_{2}}}^{(1)}\psi_{\mathbf{q_{2}}}A_{\mathbf{k}}$
$\displaystyle H_{2,FN}=$
$\displaystyle\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}\psi^{\dagger}_{\mathbf{k}-\mathbf{k}^{\prime}+\mathbf{q_{2}}}\mathcal{M}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}^{(2)}\psi_{\mathbf{q_{2}}}A_{\mathbf{k}}A_{-\mathbf{k}^{\prime}}$
$\displaystyle+\sum_{\mathbf{k}}\mathcal{M}_{\mathbf{k}}^{(2),0}A_{\mathbf{k}}A_{-\mathbf{k}}$
(20)

represent respectively the linear (Eq. 8) and quadratic (Eq. 9) coupling with
the phonons as shown in Fig. 2. The detailed expression of the scattering
vertices $\mathcal{M}_{\mathbf{k},\mathbf{q_{2}}}^{(1)}$,
$\mathcal{M}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}^{(2)}$ and
$\mathcal{M}_{\mathbf{k}}^{(2),0}$ are given in Appendix D.1. Notably the
second term in $H_{2,FN}$, given by $\mathcal{M}^{(2),0}_{\bf k}$ only depends
on the displacement operators via a quadratic form similar to that of the
Harmonic potential. This corresponds to the renormalisation of the bare phonon
frequency and hence to sound speed due to ferronematic ordering.

The contributions from the spin-nematic waves of the ferronematic to $\Delta
v/v$ can be captured to the leading order in magnetoelastic coupling by
calculating the phonon free energy due to
$\mathcal{M}_{\mathbf{k},\mathbf{q_{2}}}^{(1)}$ and
$\mathcal{M}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}^{(2)}$ terms
given by Feynman diagrams in Fig. 3.

$\mathbf{k}$$\mathbf{{k+q_{2}}}$$\mathbf{{q_{2}}}$

$\mathbf{k}$$\mathbf{-k^{\prime}}$$\mathbf{{k-k^{\prime}+q_{2}}}$$\mathbf{{q_{2}}}$

Figure 2: Interactions between the phonon and ferronematic spin waves (Eq.
20). The left diagram comes from $H_{1,FN}$ and right diagram from first term
in $H_{2,FN}$. Solid lines represent phonons and dashed lines represent
ferronematic spin waves ($\psi$ in Eq. 17)

.

$\mathbf{q}$$\mathbf{{q_{2}}}$$\mathbf{{q+q_{2}}}$$\mathbf{q}$

$\mathbf{q}$$\mathbf{q}$$\mathbf{{q_{2}}}$

Figure 3: Diagrams used to calculate contribution to the phonon self energy
due to nematic fluctuations (Fig. 2). The left diagram comes from vertices due
to $H_{1,FN}$ and right diagram from $H_{2,FN}$ (Eq. 20).

.

The phonon self energy (in Eq. 14) due to the above two contributions is then
given by

$\displaystyle\Sigma(\mathbf{q},i\Omega_{n})=\Sigma_{1}(\mathbf{q},i\Omega_{n})+\Sigma_{2}(\mathbf{q},i\Omega_{n})$
(21)

where

$\displaystyle\Sigma_{1}(\mathbf{q},i\Omega_{n})$ $\displaystyle\approx$
$\displaystyle-\frac{1}{\beta}\sum_{\begin{subarray}{c}\mathbf{q}_{2},\lambda,\\\
\rho,\omega_{1}\end{subarray}}\mathcal{M}_{-\mathbf{q},\mathbf{q}+\mathbf{q}_{2}}^{(1),\lambda\rho}\mathcal{M}_{\mathbf{q},\mathbf{q}_{2}}^{(1),\rho\lambda}G^{m,\lambda\lambda}_{0}(\mathbf{q}_{2},i\omega_{1})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times
G^{m,\rho\rho}_{0}(\mathbf{q}+\mathbf{q}_{2},i\Omega_{n}+i\omega_{1})$
$\displaystyle\Sigma_{2}(\mathbf{q},i\Omega_{n})$ $\displaystyle\approx$
$\displaystyle
2\mathcal{M}_{\mathbf{q}}^{(2),0}+\sum_{\mathbf{q}_{2},\lambda}2\mathcal{M}_{\mathbf{q},\mathbf{q},\mathbf{q_{2}}}^{(2),\lambda\lambda}\langle\psi^{\dagger\lambda}_{\mathbf{q}_{2}}\psi^{\lambda}_{\mathbf{q}_{2}}\rangle~{},$
(22)

with $\lambda,\rho=1,2$ and $G^{m,\lambda\rho}_{0}$ denotes the bare Matrix
ferronematic spin-wave Green’s function defined as
$G^{m,\lambda\rho}_{0}(\mathbf{q},\tau-\tau^{\prime})=-\langle
T_{\tau}(\psi^{\lambda}_{\mathbf{q}}(\tau)\psi^{\dagger\rho}_{\mathbf{q}}(\tau^{\prime}))\rangle$.
Due to the diagonal form of the spin wave Hamiltonian (Eq. 17), the matrix is
given by

$G^{m}_{0}(\mathbf{q},i\Omega_{n})=\left(\begin{matrix}\frac{1}{i\Omega_{n}-\omega^{s}_{\mathbf{q}}}&0\\\
0&\frac{-1}{i\Omega_{n}+\omega^{s}_{-\mathbf{q}}}\end{matrix}\right)~{}.$ (23)

The sum over the bosonic Matsubara frequencies in Eq. 22 is evaluated using
the standard techniques Bruus and Flensberg (2004) and using Eq. 13, the
temperature dependence of the fractional change in sound speed is obtained as

$\frac{\Delta v}{v}=\bar{c}_{1}+\bar{c}_{2}T^{3}~{},$ (24)

The factors $\bar{c}_{1}$ and $\bar{c}_{2}$ (see expressions in Appendix D.1)
have contributions from first and second derivatives of both the bilinear and
the biquadratic couplings and depend on the details of the lattice via the
phonon spectrum and polarisation. However, the above temperature dependence is
generically valid for other two dimensional lattices. The sound attenuation,
$\propto Im\Sigma/v_{0}$, can similarly be calculated from the imaginary part
of the phonon self energy. At this order only $\Sigma_{1}$ contributes and
such sound attenuation generically proportional, at low frequencies, the to
the phonon frequency.Cottam (1974); Gen  _et al._ (2019)

###  III.2 The three sublattice nematic phase

The $\Delta v/v$ for the three sublattice nematic can be obtained in a similar
way. In analogy with the ferronematic case (Eq. 16), the relevant Hamiltonian
is given by

$\displaystyle\mathcal{H}_{3SN}$ $\displaystyle=H_{ph}+H_{sp,3SN}+H_{sp-
ph,3SN}~{},$ (25)

where $H_{ph}$ is the harmonic phonon Hamiltonian (Eq. 12), $H_{sp,3SN}$ and
$H_{sp-ph,3SN}$ are respectively the linear spin-nematic wave Hamiltonian for
the three sublattice nematic phase and the phonon-three sublattice spin-
nematic wave coupling Hamiltonian.

Due to the non-uniform structure of this nematic phase, the calculations are
somewhat more tedious and the relevant parts are relegated to Appendix D.2.
The difference, however, in this case stems from the three sublattice spin-
nematic wave spectrum. The spin-nematic wave theory about a three sublattice
nematic state, Tsunetsugu and Arikawa (2006) such as
$\prod_{\mathbf{\widetilde{R}}}\ket{x}_{\mathbf{\widetilde{R}},1}\ket{y}_{\mathbf{\widetilde{R}},2}\ket{z}_{\mathbf{\widetilde{R}},3}$
is obtained as follows. Here $\mathbf{\widetilde{R}}$ stands for the three
sublattice unit cell and 1, 2 and 3 denote the sublattices of each unit cell
with their nematic directors being along the $\hat{\mathbf{x}}$,
$\hat{\mathbf{y}}$ and $\hat{\mathbf{z}}$ directions respectively as shown in
Fig. 1. We introduce two bosons $\widetilde{\alpha}^{\dagger}$ and
$\widetilde{\beta}^{\dagger}$ at each sublattice to capture the deviations
from the ground state, e.g., for sublattice-3,
$\ket{z}_{\mathbf{\widetilde{R}},3}\equiv\ket{\text{vac}}_{\mathbf{\widetilde{R}},3}$
and $\ket{S_{z}=\pm
1}_{\mathbf{\widetilde{R}},3}=\frac{1}{\sqrt{2}}(\widetilde{\alpha}^{\dagger}_{\mathbf{\widetilde{R}},3}\pm
i\widetilde{\beta}^{\dagger}_{\mathbf{\widetilde{R}},3})\ket{\text{vac}}_{\mathbf{\widetilde{R}},3}$.
The details are summarised in Appendix C.2.

The resultant diagonalised harmonic Hamiltonian (similar to Eq. 17 for the
ferronematic) is given by

$H_{sp,3SN}=\sum_{\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}\bigg{[}\widetilde{\omega}^{s}_{+,\mathbf{k}}\text{\
}\widetilde{d}_{+,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}^{~{}\dagger}\widetilde{d}_{+,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}+\widetilde{\omega}^{s}_{-,\mathbf{k}}\text{\
}\widetilde{d}_{-,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}^{~{}\dagger}\widetilde{d}_{-,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}\bigg{]}$
(26)

where $\widetilde{\lambda}\widetilde{\rho}=\\{12,23,31\\}$,
$\widetilde{d}_{\pm,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}$ are the
diagonalised bosonic annihilation operators with dispersions

$\displaystyle\widetilde{\omega}^{s}_{\pm,\mathbf{k}}=3K_{0}\sqrt{(1\pm|\widetilde{\gamma}_{\mathbf{k}}|)\left(1\pm\bigg{(}1-\frac{2J_{0}}{K_{0}}\bigg{)}|\widetilde{\gamma}_{\mathbf{k}}|\right)}~{},$
(27)

where
$\widetilde{\gamma}_{\mathbf{k}}=\frac{1}{3}\sum_{\bm{\widetilde{\delta}}}e^{i\mathbf{k}.\bm{\widetilde{\delta}}}$
and
$\bm{\widetilde{\delta}}=\\{-\hat{\mathbf{x}}/2-\sqrt{3}\hat{\mathbf{y}}/2,\hat{\mathbf{x}},-\hat{\mathbf{x}}/2+\sqrt{3}\hat{\mathbf{y}}/2\\}$
(lattice constant set to unity).

For $J_{0}\neq K_{0}$, the $\widetilde{\omega}^{s}_{+,\mathbf{k}}$ branch is
gapped. For $J_{0}/K_{0}\ll 1$ the $\widetilde{\omega}^{s}_{+,\mathbf{k}}$
branch lies entirely above the $\widetilde{\omega}^{s}_{-,\mathbf{k}}$ branch
(touching $\widetilde{\omega}^{s}_{-,\mathbf{k}}$ at the Brillouin zone
corners.Tsunetsugu and Arikawa (2006)) Hence its effect can be neglected at
low temperatures which is dominated by the long wavelength behaviour of the
gapless mode $\widetilde{\omega}^{s}_{-,\mathbf{k}}$. The latter shows, in
turn, a crossover depending on the ratio of $J_{0}/K_{0}$. Eq. 27, for generic
$J_{0}/K_{0}$, at long wavelengths,
$\widetilde{\omega}^{s}_{-,\mathbf{k}}\approx c_{s}|{\bf k}|$ (with
$c_{s}=3\sqrt{\frac{J_{0}K_{0}}{2}}$). However for $J_{0}=0$, the long
wavelength scaling of the dispersion changes to
$\widetilde{\omega}^{s}_{-,\mathbf{k}}\approx\widetilde{c}_{s}|{\bf k}|^{2}$
(with $\widetilde{c}_{s}=\frac{3K_{0}}{4}$ )Tsunetsugu and Arikawa (2006) such
that for $J_{0}/K_{0}\ll 1$, there is a crossover scale
$k^{*}=\sqrt{\frac{8J_{0}}{K_{0}}}$ (in inverse units of lattice length scale)
below which the dispersion is approximately linear and above which it is
approximate quadratic. This change of the dispersion as well as the crossover
affects the temperature dependence of the $\Delta v/v$ in the three sublattice
nematic depending on the ratio of $J_{0}/K_{0}$ (see below).

Finally, the three sublattice nematic spin wave-phonon Hamiltonian, $H_{sp-
ph,3SN}$ (in eq. 25) is obtained, similar to the ferronematic case, by writing
the spin operators (Eq. 103) in terms of the
$\widetilde{d}_{\pm,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}$ bosons of
Eq. 26 and substituting in Eq. 8, 9, 10 and 11 (analogous to Eq. 19). This
leads to Eq. 112. The fractional change of sound speed due to the three
sublattice ordering is then obtained similar to the the ferronematic (details
in Appendix D.2) and is given by

$\displaystyle\frac{\Delta v}{v}$ $\displaystyle\approx$
$\displaystyle\lim_{\mathbf{q}\xrightarrow{}0}\Bigg{[}\int_{0}^{k^{*}}dk\int
d\theta_{\mathbf{\hat{k}}}\frac{Re\widetilde{\Sigma}_{\mathbf{k}}\left(\widetilde{\omega}^{s}_{-,\mathbf{k}}\propto|\mathbf{k}|\right)}{\omega_{0,\mathbf{q}}^{ph}}$
(28) $\displaystyle+$ $\displaystyle\int_{k^{*}}^{\infty}dk\int
d\theta_{\mathbf{\hat{k}}}\frac{Re\widetilde{\Sigma}_{\mathbf{k}}\left(\widetilde{\omega}^{s}_{-,\mathbf{k}}\propto|\mathbf{k}|^{2}\right)}{\omega_{0,\mathbf{q}}^{ph}}\Bigg{]}~{},$

accounting for the crossover between the linear and quadratic dispersions at
$k^{*}$. This approximately leads to, in the three sublattice nematic,

$\displaystyle\frac{\Delta
v}{v}=\left\\{\begin{array}[]{l}\widetilde{c}_{1}+\widetilde{c}_{2}T^{3}~{}~{}~{}~{}{\rm
for}~{}\beta J_{0}\gg 1\\\
\widetilde{c}_{3}+\widetilde{c}_{4}T^{2}~{}~{}~{}~{}{\rm for}~{}\beta
J_{0}\sim 1\\\ \widetilde{c}_{5}+\widetilde{c}_{6}T~{}~{}~{}~{}~{}{\rm
for}~{}\beta J_{0}<1\\\ \end{array}\right.$ (32)

where $\beta=1/(k_{B}T)$ and the details of the pre-factors are given in
Appendix D.2. No such crossover behaviour is observed for the ferronematic
case since the spin wave dispersion remains linear in $k$ even on setting
$J_{0}=0$.

###  III.3  The thermal paramagnet

In contrast to the low temperature nematic phases discussed above, for the
high temperature paramagnetic phase, the spin dynamics is faster than the
acoustic phonons and integrating out spins, gives the following effective
interaction Hamiltonian for phonons,Bhattacharjee  _et al._ (2011)

$H_{\text{eff}}=\sum_{\mathbf{k}}\omega_{0,\mathbf{k}}^{ph}a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}+\frac{1}{2}\sum_{\mathbf{k},\mathbf{k}^{\prime}}V_{\mathbf{k}\mathbf{k}^{\prime}}A_{\mathbf{k}}A_{-\mathbf{k}^{\prime}}$
(33)

where

$V_{\mathbf{k}\mathbf{k}^{\prime}}=\langle\bar{\bar{U}}_{\mathbf{k}\mathbf{k}^{\prime}}^{(2)}\rangle-\beta\langle\langle\bar{\bar{U}}_{\mathbf{k}}^{(1)}\bar{\bar{U}}_{-\mathbf{k}^{\prime}}^{(1)}\rangle\rangle$
(34)

Here
$\bar{\bar{U}}_{\mathbf{k}}^{(1)}=U_{\mathbf{k}}^{(1)}+\widetilde{U}_{\mathbf{k}}^{(1)}$
and
$\bar{\bar{U}}_{\mathbf{k}\mathbf{k}^{\prime}}^{(2)}=U_{\mathbf{k}\mathbf{k}^{\prime}}^{(2)}+\widetilde{U}_{\mathbf{k}\mathbf{k}^{\prime}}^{(2)}$
(which have been defined in Eq. 10 and Eq. 11) and
$\langle\langle....\rangle\rangle$ denotes the connected correlators for the
spins, averaged over a thermal ensemble. The self energy is given by :

$\Sigma(\mathbf{q},i\Omega_{n})=V_{\mathbf{q}\mathbf{q}}$ (35)

So, calculating the leading order temperature dependence of $\Delta v/v$ would
involve various spin correlators which are then calculated using high
temperature series expansion. To the leading order in $1/T$, this gives

$\frac{\Delta v}{v}=a_{1}+\widetilde{a}_{2}\beta~{},$ (36)

where, notably, a constant contribution arises from the biquadratic term such
that

$a_{1}=\sum_{\bm{\delta}}\frac{1}{3Mv_{0}^{2}}(\hat{\mathbf{q}}.\bm{\delta})^{2}\left(\mathbf{e}_{-\mathbf{q}}.\frac{\partial^{2}K}{\partial{\bm{\delta}^{2}}}.\mathbf{e}_{\mathbf{q}}\right)~{},$
(37)

where $\bm{\delta}$ runs over the six nearest neighbours. It is evident that
the constant arises only due to the biquadratic coupling and would be absent
for the case of pure Heisenberg model.

##  IV The Landau-Ginzburg Theory

The above microscopic approach works deep inside the nematic at low
temperatures or in the thermal paramagnet at high temperatures. For general
elastic responses for the entire phase diagram, a more phenomenological Landau
-Ginzburg theory– that accounts for the 120∘ spiral and the spin-nematic (both
ferro and three sub-lattice) ordering as well as the elastic degrees of
freedom– is useful to account for the spin-lattice physics of the system.Patri
_et al._ (2019) Below, we construct this theory to derive the elastic
signatures of a spin-nematic in context of the NiGa2S4. Our calculations are
easily generalised to other situations, in particular the spin-orbit coupled
multi-polar orders.Voleti  _et al._ (2022); Patri  _et al._ (2019)

We start with the symmetry analysis for the three fields – the magnetic, the
spin-nematic and elastic – for the point group symmetry of
${\text{NiGa}_{2}\text{S}_{4}}$ which is $\text{D}_{3d}$. For our present
calculations, we choose the largest unit cell for all the type of orderings
discussed above– a single triangle– and systematically isolate the relevant
symmetry allowed terms for the three fields including their interactions. We
choose an up triangle that consists of three sites of $\text{Ni}^{2+}$ (the
triangular unit and site labels are shown in Fig.1 and their positions being
$\mathbf{r}_{1}=0$, $\mathbf{r}_{2}=\hat{\mathbf{x}}$ and
$\mathbf{r}_{3}=\frac{1}{2}\hat{\mathbf{x}}+\frac{\sqrt{3}}{2}\hat{\mathbf{y}}$)
and impose inversion symmetry to obtain the normal modes. The non-trivial
transformations for the up triangle are $C_{3}$ and $\sigma_{h}C_{2}^{\prime}$
(these transformations keep the centre of the triangle fixed, $\sigma_{h}$ is
required with the $2$-fold rotation to bring the crystal field environment
back to itself, details in Appendix E). Combining this with inversion
generates all the $6$ conjugacy classes of $\text{D}_{3d}$. We can then
decompose the dipole, quadrupole fields and the elastic modes into the
irreducible representations (irrep) to construct the Landau-Ginzburg free
energy.

###  IV.1 The magnetic (dipolar) and nematic (quadrupolar) modes

Eq. 2 the system has full SU(2) spin rotation symmetry. Hence the spin
operators remain un-rotated (in spin space) under various lattice
transformations while the site indices transform. Using these latter
transformations, the non-trivial irreducible representations are constructed.

Irrep | Expression
---|---
$m_{a}^{\alpha}$ | $(S_{1}^{\alpha}+S_{2}^{\alpha}+S_{3}^{\alpha})/\sqrt{3}$
$m_{e1}^{\alpha}$ | $(S_{1}^{\alpha}-S_{2}^{\alpha})/\sqrt{2}$
$m_{e2}^{\alpha}$ | $(S_{1}^{\alpha}+S_{2}^{\alpha}-2S_{3}^{\alpha})/\sqrt{6}$
$Q_{a}^{\alpha\beta}$ | $(Q_{1}^{\alpha\beta}+Q_{2}^{\alpha\beta}+Q_{3}^{\alpha\beta})/\sqrt{3}$
$Q_{e1}^{\alpha\beta}$ | $(Q_{1}^{\alpha\beta}-Q_{2}^{\alpha\beta})/\sqrt{2}$
$Q_{e2}^{\alpha\beta}$ | $(Q_{1}^{\alpha\beta}+Q_{2}^{\alpha\beta}-2Q_{3}^{\alpha\beta})/\sqrt{6}$
Table 1: Irreducible representations for spins and quadrupoles
($\alpha,\beta=x,y,z$).

The irreducible representations for dipoles and quadrupole on the up triangle
with site labels as shown in Fig. 1 are given in Table 1. Note that for both
the fields the irreducible representations consist of a singlet ($a$) and a
doublet ($e=e1,e2$). However, it is useful to note that while the dipolar
field is odd under time reversal, the quadrupolar field is even. The details
of the symmetry transformations of the irreps are listed in Appendix E. All
the relevant orders can be represented as different combinations of the
irreducible representations.

####  IV.1.1 Magnetic orders

The relevant magnetic orders are :

##### Ferromagnetic order :

This is given by $\langle\mathbf{m}_{a}\rangle\neq 0$,
$\langle\mathbf{m}_{e1}\rangle=0$ and $\langle\mathbf{m}_{e2}\rangle=0$ (where
bold font is used to suppress the spin index). Inverting the relations for
spins in Table 1, we get :

$\displaystyle\langle\mathbf{S}_{1}\rangle$
$\displaystyle=\frac{1}{6}\left(2\sqrt{3}\langle\mathbf{m}_{a}\rangle+3\sqrt{2}\langle\mathbf{m}_{e1}\rangle+\sqrt{6}\langle\mathbf{m}_{e2}\rangle\right)$
$\displaystyle\langle\mathbf{S}_{2}\rangle$
$\displaystyle=\frac{1}{6}\left(2\sqrt{3}\langle\mathbf{m}_{a}\rangle-3\sqrt{2}\langle\mathbf{m}_{e1}\rangle+\sqrt{6}\langle\mathbf{m}_{e2}\rangle\right)$
$\displaystyle\langle\mathbf{S}_{3}\rangle$
$\displaystyle=\frac{1}{3}\left(\sqrt{3}\langle\mathbf{m}_{a}\rangle-\sqrt{6}\langle\mathbf{m}_{e2}\rangle\right)$
(38)

which is evidently a ferromagnet for only $\langle\mathbf{m}_{a}\rangle\neq
0$. The ferromagnetic order also has a non-vanishing parasitic quadrupolar
moment (see Appendix A).

#####  $120^{\circ}$ spiral order :

The spiral order is given by the three simultaneous conditions

$\displaystyle\langle\mathbf{m}_{a}\rangle=0,~{}~{}~{}~{}~{}\langle\mathbf{m}_{e1}\rangle.\langle\mathbf{m}_{e2}\rangle=0$
$\displaystyle\langle\mathbf{m}_{e1}\rangle.\langle\mathbf{m}_{e1}\rangle=\langle\mathbf{m}_{e2}\rangle.\langle\mathbf{m}_{e2}\rangle$
(39)

The first condition (using Eq. 38) translates into zero magnetisation
condition per triangle , while the second and third conditions results in the
equality of magnitude of the moments at different sites, i.e.,
$\langle\mathbf{S}_{1}\rangle.\langle\mathbf{S}_{1}\rangle=\langle\mathbf{S}_{2}\rangle.\langle\mathbf{S}_{2}\rangle=\langle\mathbf{S}_{3}\rangle.\langle\mathbf{S}_{3}\rangle$
and equal angle between the ordered moments, i.e.,
$\langle\mathbf{S}_{1}\rangle.\langle\mathbf{S}_{2}\rangle=\langle\mathbf{S}_{2}\rangle.\langle\mathbf{S}_{3}\rangle=\langle\mathbf{S}_{3}\rangle.\langle\mathbf{S}_{1}\rangle$.
From this is is fairly easy to show that the angle between any two nearest
neighbour moments is $120^{\circ}$ such that

$\displaystyle\langle\mathbf{S}_{i}\rangle=$
$\displaystyle\frac{\sqrt{6}}{3}\bigg{(}\langle\mathbf{m}_{e1}\rangle\text{cos}\bigg{(}\mathbf{q}.\mathbf{r}_{i}+\frac{\pi}{6}\bigg{)}+\langle\mathbf{m}_{e2}\rangle\text{sin}\bigg{(}\mathbf{q}.\mathbf{r}_{i}+\frac{\pi}{6}\bigg{)}\bigg{)}$
(40)

with $i=1,2,3$ and
$\mathbf{q}=\frac{2\pi}{3}\hat{\mathbf{x}}+\frac{2\pi}{\sqrt{3}}\hat{\mathbf{y}}$.
The phase factor of $\pi/6$ is due our choice of doublet modes of spins in
Table. 1. Similar to the ferromagnet, the spiral order also has non-vanishing
parasitic quadrupolar moments.

####  IV.1.2 Nematic orders

For pure spin nematic orders, the expectation value of the magnetic moments
vanishes, i.e. $\langle{\bf S}_{i}\rangle=0$. The on-site quadrupolar tensor
is symmetric and can be diagonalized and using Eq. 1, one can obtain the
expectation values $\langle(\mathbf{S}.\hat{\mathbf{e}}_{1})^{2}\rangle$,
$\langle(\mathbf{S}.\hat{\mathbf{e}}_{2})^{2}\rangle$ and
$\langle(\mathbf{S}.\hat{\mathbf{e}}_{3})^{2}\rangle$ along the principal axes
$\hat{\mathbf{e}}_{1}$, $\hat{\mathbf{e}}_{2}$ and $\hat{\mathbf{e}}_{3}$ of
the quadrupole ellipsoid.Kosmachev  _et al._ (2015) In terms of the quadrupole
ellipsoid, there are three possibilities : (1) No nematic order :
$\langle(\mathbf{S}.\hat{\mathbf{e}}_{1})^{2}\rangle=\langle(\mathbf{S}.\hat{\mathbf{e}}_{2})^{2}\rangle=\langle(\mathbf{S}.\hat{\mathbf{e}}_{3})^{2}\rangle$,
(2) Biaxial nematic :
$\langle(\mathbf{S}.\hat{\mathbf{e}}_{1})^{2}\rangle\neq\langle(\mathbf{S}.\hat{\mathbf{e}}_{2})^{2}\rangle\neq\langle(\mathbf{S}.\hat{\mathbf{e}}_{3})^{2}\rangle$,
and, (3) uniaxial nematic : This is obtained when
$\langle(\mathbf{S}.\hat{\mathbf{e}}_{i})^{2}\rangle=\langle(\mathbf{S}.\hat{\mathbf{e}}_{j})^{2}\rangle\neq\langle(\mathbf{S}.\hat{\mathbf{e}}_{k})^{2}\rangle$
for $i\neq j\neq k$. For spin-1, the first case implies a paramagnet, while
the biaxial nematic is not relevant for us. Hence we focus on the uniaxial
nematic where, if the unequal expectation value is greater (lesser) than the
equal ones, we have a rod (disc)-like uniaxial nematic. For spin-1, only a
uniaxial nematic of the disc type is allowed.Penc and Läuchli (2011); Tóth
(2011) Therefore for the on-site quadrupolar order relevant to our
calculations, the order parameter is characterised by

$\displaystyle\langle
Q^{\alpha\beta}_{i}\rangle=\mathcal{Q}_{i,N}\big{(}n^{\alpha}_{i}n^{\beta}_{i}-\delta^{\alpha\beta}/3\big{)}$
(41)

with $\mathcal{Q}_{i,N}<0$ and $\mathbf{n}_{i}$ being the director of the
nematic.

All possible three sublattice disk-like uniaxial nematic orders can be
obtained from the quadrupolar irreps in Table 1. The on-site expectation
values obtained by inverting the relations are given by :

$\displaystyle\langle Q_{1}^{\alpha\beta}\rangle$
$\displaystyle=\frac{1}{6}\left(2\sqrt{3}\langle
Q_{a}^{\alpha\beta}\rangle+3\sqrt{2}\langle
Q_{e1}^{\alpha\beta}\rangle+\sqrt{6}\langle
Q_{e2}^{\alpha\beta}\rangle\right)$ $\displaystyle\langle
Q_{2}^{\alpha\beta}\rangle$ $\displaystyle=\frac{1}{6}\left(2\sqrt{3}\langle
Q_{a}^{\alpha\beta}\rangle-3\sqrt{2}\langle
Q_{e1}^{\alpha\beta}\rangle+\sqrt{6}\langle
Q_{e2}^{\alpha\beta}\rangle\right)$ $\displaystyle\langle
Q_{3}^{\alpha\beta}\rangle$ $\displaystyle=\frac{1}{3}\left(\sqrt{3}\langle
Q_{a}^{\alpha\beta}\rangle-\sqrt{6}\langle Q_{e2}^{\alpha\beta}\rangle\right)$
(42)

#####  Ferronematic order:

This is given by

$\displaystyle\langle
Q_{a}^{\alpha\beta}\rangle=\sqrt{3}\mathcal{Q}_{FN}\big{(}n^{\alpha}n^{\beta}-\frac{1}{3}\delta^{\alpha\beta}\big{)},~{}~{}~{}~{}\langle{\bf
Q}_{\bf e}\rangle=0$ (43)

The on-site expectation values can be obtained from Eq. 42 and are given by
$\langle Q_{1}^{\alpha\beta}\rangle=\langle Q_{2}^{\alpha\beta}\rangle=\langle
Q_{3}^{\alpha\beta}\rangle=\frac{1}{\sqrt{3}}~{}\langle
Q_{a}^{\alpha\beta}\rangle$, which is clearly a ferronematic (Eq. 3) that is
energetically favoured Bhattacharjee  _et al._ (2006) when $K<0$ in Eq. 2. For
$J>0$ and $K<0$, the ferronematic phase is stable for $|K|/J>2$ within mean-
field analysis.Läuchli  _et al._ (2006)

##### Three sublattice nematic order:

This is given by:

$\displaystyle\langle Q_{a}^{\alpha\beta}\rangle$ $\displaystyle=0$ (44)
$\displaystyle\langle Q_{e1}^{\alpha\beta}\rangle$
$\displaystyle=\frac{\mathcal{Q}_{3SN}}{\sqrt{2}}(n_{1}^{\alpha}n_{1}^{\beta}-n_{2}^{\alpha}n_{2}^{\beta})$
(45) $\displaystyle\langle Q_{e2}^{\alpha\beta}\rangle$
$\displaystyle=\frac{\mathcal{Q}_{3SN}}{\sqrt{6}}(n_{1}^{\alpha}n_{1}^{\beta}+n_{2}^{\alpha}n_{2}^{\beta}-2n_{3}^{\alpha}n_{3}^{\beta})~{},$
(46)

where $\mathbf{n}_{1},\mathbf{n}_{2},\mathbf{n}_{3}$ are mutually orthogonal
unit vectors and $\mathcal{Q}_{3SN}<0$. Again, using Eq. 42, we can see that
the above conditions correspond to the three sublattice nematic order with the
directors on sites $1$, $2$ and $3$ being $\mathbf{n}_{1}$, $\mathbf{n}_{2}$
and $\mathbf{n}_{3}$ respectively. Such three sublattice ordering is expected
to be stabilised when the sign of the biquadratic term in the spin Hamiltonian
(Eq. 2), $K>0$. For $J>0$ and $K>0$, the three sublattice nematic is obtained
for $K/J>1$.Läuchli  _et al._ (2006)

####  IV.1.3 Coexistence of nematic order and collinear sinusoidal dipolar
order

In passing, we point out the possibility of an interesting phase of coexisting
nematic and dipole order with a three site unit-cell and hence captured within
the above formulation. This is given by
$\langle\mathbf{m}_{a}\rangle=0~{},~{}~{}~{}~{}\langle\mathbf{m}_{e1}\rangle\parallel\langle\mathbf{m}_{e2}\rangle$.
The resultant dipolar order is collinear and sinusoidal Kawamura (1998) with
the spin configuration given by
$\langle\mathbf{S}_{i}\rangle=\frac{\sqrt{6}}{3}\langle\mathbf{m}_{e1}\rangle\bigg{(}\sqrt{1+\widetilde{\lambda}^{2}}\bigg{)}\text{sin}\bigg{(}\mathbf{q}.\mathbf{r}_{i}+\frac{\pi}{6}+\phi\bigg{)}$,
where
$\mathbf{q}=\frac{2\pi}{3}\hat{\mathbf{x}}+\frac{2\pi}{\sqrt{3}}\hat{\mathbf{y}}$.
The above form can be obtained from Eq. 40 by using
$\langle\mathbf{m}_{e2}\rangle=\widetilde{\lambda}\langle\mathbf{m}_{e1}\rangle$.
Thus the state corresponds to collinear sinusoidal magnetic order. Since the
spins are not completely polarised it allows for (non-parasitic) nematic
ordering. Using Eq. 90 we find that the sinusoidal dipole ordering is
accompanied by bi-axial nematic ordering with the two orthogonal principal
directions of the nematic directors being along
$\hat{n}=\langle\mathbf{m}_{e1}\rangle/|\langle\mathbf{m}_{e1}\rangle|$ and
perpendicular to it respectively. The quadrupole moment and the magnetic
moment do not have the same symmetry implying the coexistence of nematic and
collinear sinusoidal dipolar order. In the present Hamiltonian we do not
expect this order to be stabilised. However, they may be relevant for more
generic models such as the one studied in Ref. Seifert and Savary, 2022.

###  IV.2 The Elastic modes

Having discussed the dipole and the quadrupole modes of interest, we now turn
to the elastic modes that they can couple to. The normal modes of a single
triangle are given in Appendix E.2. This consists of

$\displaystyle\text{a
singlet}:\epsilon_{a}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{a
doublet}:(\epsilon_{e1},\epsilon_{e2})$ (47)

They are linearly related with the Cartesian strain tensors,
$\varepsilon_{ij}$ (where $\varepsilon_{ij}=\frac{1}{2}\left(\frac{\partial
u_{i}}{\partial R_{j}}+\frac{\partial u_{j}}{\partial R_{i}}\right)$ ; $u_{i}$
is the $i^{th}$ component of displacement from equilibrium position, $R_{i}$)
as (see Appendix E.2)

$\displaystyle\varepsilon_{a}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\varepsilon_{xx}+\varepsilon_{yy});$
$\displaystyle\varepsilon_{e1}$ $\displaystyle=$
$\displaystyle\varepsilon_{xy},~{}~{}~{}~{}~{}~{}~{}~{}\varepsilon_{e2}=\frac{1}{2}(\varepsilon_{xx}-\varepsilon_{yy})$
(48)

We now study their coupling with the dipole and the quadrupole modes to
understand the nature of magnetoelastic response of the system.

###  IV.3 The Landau free energy

To write the Landau-Ginzburg free energy, we consider the long wavelength
symmetry allowed terms for the dipolar, quadrupolar and the elastic modes. At
the mean field level we only consider uniform terms and drop all spatially
fluctuating ones.

####  IV.3.1  Dipole-quadrupole free energy

The dipole-quadrupole free energy is

$\mathcal{F}_{mQ}=\mathcal{F}_{m}+\mathcal{F}_{Q}+\mathcal{F}_{mQ}^{\text{int}}~{},$
(49)

where the three terms denote pure dipolar and pure quadrupolar free energies;
and the interaction between the dipoles and the quadrupoles.

##### Dipolar terms:

The dipolar modes are odd under time reversal so they appear in even powers in
the free energy. The dipoles have $a$ and $e$ irreps so the dipolar free
energy can be written as the sum of free energy for individual modes and the
interaction terms between them

$\mathcal{F}_{m}=\mathcal{F}_{m,a}+\mathcal{F}_{m,e}+\mathcal{F}_{m,ae}^{int}~{},$
(50)

where upto quartic orders

$\mathcal{F}_{m,a}=r_{m_{a}}(\mathbf{m}_{a}.\mathbf{m}_{a})+u_{m_{a}}(\mathbf{m}_{a}.\mathbf{m}_{a})^{2}~{},$
(51) $\displaystyle\mathcal{F}_{m,e}$ $\displaystyle=$ $\displaystyle
r_{m_{e}}(\mathbf{m}_{\bf e}.\mathbf{m}_{\bf e})+u_{m_{e}}(\mathbf{m}_{\bf
e}.\mathbf{m}_{\bf e})^{2}$ (52) $\displaystyle+$ $\displaystyle
v_{m_{e}}[(\mathbf{m}_{e1}.\mathbf{m}_{e2})^{2}-(\mathbf{m}_{e1}.\mathbf{m}_{e1})(\mathbf{m}_{e2}.\mathbf{m}_{e2})],$

are the free energies of singlet and doublet modes while

$\displaystyle\mathcal{F}_{m,ae}^{int}=$ $\displaystyle
v_{m}(\mathbf{m}_{a}.\mathbf{m}_{a})(\mathbf{m}_{\bf e}.\mathbf{m}_{\bf
e})+\tilde{v}_{m}(m_{a}^{\alpha}m_{a}^{\beta})(m_{e}^{\alpha}.m_{e}^{\beta})~{},$
(53)

represents the interaction between them.

As discussed in section IV.1.1, ferromagnetic order corresponds to
$\langle\mathbf{m}_{a}\rangle\neq 0$ with $\langle\mathbf{m}_{\bf
e}\rangle=0$. The mean field theory leads to a continuous transition between
paramagnetic (for $r_{m_{a}}>0$) and ferromagnetic (for $r_{m_{a}}<0$) ordered
state which takes place at $r_{m_{a}}=0$. However, for the microscopic model
in Eq. 2, we expect that $r_{m_{a}}>0$.

For the doublet mode, the continuous transition between the thermal paramagnet
(for $r_{m_{e}}>0$) and the dipole (for $r_{m_{e}}<0$) ordered phase occurs at
$r_{m_{e}}=0$ for Eq. 52. The details of the dipole ordered phase is
controlled by $v_{m_{e}}$.Kawamura (1998) For $v_{m_{e}}>0(<0)$, the 120∘
spiral (collinear sinusoidal) state is stabilised. Note that the dipolar free
energy for the doublet mode, Eq. 52, is invariant under :

$\displaystyle\mathbf{m}_{e1}$ $\displaystyle\xrightarrow{}$
$\displaystyle\mathbf{m}_{e1}^{\prime}=\cos{\theta}\mathbf{m}_{e1}-\sin{\theta}\mathbf{m}_{e2}$
$\displaystyle\mathbf{m}_{e2}$ $\displaystyle\xrightarrow{}$
$\displaystyle\mathbf{m}_{e2}^{\prime}=\pm(\sin{\theta}\mathbf{m}_{e1}+\cos{\theta}\mathbf{m}_{e2})~{},$
(54)

for $\theta\in(0,2\pi]$. The origin of this enhanced symmetry is due to the
enlargement of the translation symmetry from $\mathbb{Z}_{3}$ to $U(1)$ as can
be quickly checked by writing the spin configurations in Eq. 40. This enhanced
symmetry, absent in the microscopic model is broken down at the sixth order by
the term
$(\mathbf{m}_{e1}.\mathbf{m}_{e1}-\mathbf{m}_{e2}.\mathbf{m}_{e2})[(\mathbf{m}_{e1}.\mathbf{m}_{e1}-\mathbf{m}_{e2}.\mathbf{m}_{e2})^{2}-12(\mathbf{m}_{e1}.\mathbf{m}_{e2})^{2}]$.Kawamura
(1990)

Finally the interaction term in Eq. 53 leads to effective attraction (for
$v_{m},\tilde{v}_{m}<0$) or repulsion (for $v_{m},\tilde{v}_{m}>0$) between
the ferromagnetic and the spiral or the sinusoidal orders. In particular,
attraction between the ferromagnetic and collinear sinusoidal order can open
up a regime of ferrimagnetic order where both these order parameters are non-
zero.

##### Quadrupolar terms:

The quadrupolar modes are even under time reversal and can be compactly
written in an SU(2) spin rotation invariant form in terms of traces over spin
indices giving rise to

$\mathcal{F}_{Q}=\mathcal{F}_{Q,a}+\mathcal{F}_{Q,e}+\mathcal{F}_{Q,ae}^{int}~{}.$
(55)

where up to quartic orders,

$\mathcal{F}_{Q,a}=r_{Q_{a}}\text{Tr}\mathbf{Q}_{a}^{2}+w_{Q_{a}}\text{Tr}\mathbf{Q}_{a}^{3}+u_{Q_{a}}\text{Tr}\mathbf{Q}_{a}^{4}$
(56)
$\displaystyle\mathcal{F}_{Q,e}=r_{Q_{e}}\text{Tr}(\mathbf{Q}_{e1}^{2}+\mathbf{Q}_{e2}^{2})+s_{Q_{e}}\text{Tr}(\mathbf{Q}_{e2}^{3}-3\mathbf{Q}_{e1}^{2}\mathbf{Q}_{e2})$
$\displaystyle+u_{Q_{e},1}\text{Tr}(\mathbf{Q}_{e1}^{2}+\mathbf{Q}_{e2}^{2})^{2}+u_{Q_{e},2}[\text{Tr}(\mathbf{Q}_{e1}^{2}+\mathbf{Q}^{2}_{e2})]^{2}$
$\displaystyle+u_{Q_{e},3}[\text{Tr}(\mathbf{Q}_{e1}\mathbf{Q}_{e2})\text{Tr}(\mathbf{Q}_{e1}\mathbf{Q}_{e2})-\text{Tr}(\mathbf{Q}_{e1}^{2})\text{Tr}(\mathbf{Q}_{e2}^{2})]~{}.$

for the singlet and doublet modes and

$\displaystyle\mathcal{F}_{Q,ae}^{int}=$ $\displaystyle
s_{Q_{a}Q_{e}}Q_{a}^{\alpha\beta}\bigg{[}\bigg{(}Q_{e1}^{\beta\gamma}Q_{e1}^{\gamma\alpha}+Q_{e2}^{\beta\gamma}Q_{e2}^{\gamma\alpha}\bigg{)}$
$\displaystyle-\frac{1}{3}\delta^{\alpha\beta}\text{Tr}(\mathbf{Q}_{e1}^{2}+\mathbf{Q}_{e2}^{2})\bigg{]}~{},$
(58)

represents the leading order symmetry allowed interaction between them. On
general grounds all associated transitions out of the nematic orders described
by Eq. 55, are first order (without fine-tuning) due to the presence of the
third order terms. Understanding these different nematic orders is tedious due
to the matrix order parameter $\mathbf{Q}$ and here we restrict ourselves to
the simple nematic orders relevant to Eq. 2.

As discussed in section IV.1.2, the ferronematic order corresponds to Eq. 43.
This is obtained in the regime $r_{Q_{a}}<0$ and $r_{Q_{e}}>0$ in Eq. 55 where
the above solution is gotten by considering a general singlet quadrupolar
tensor in its eigen-basis
($\mathbf{Q}_{a}=\text{diag}(l_{1},l_{2},-l_{1}-l_{2})$) and extremising with
respect to the eigenvalues.

For the three sublattice nematic which is an ordering in the doublet mode, we
substitute the ansatz of Eq. (45) and (46) in the doublet free energy
Eq.(LABEL:fqe) to get,

$\displaystyle\mathcal{F}_{Q,e}=$ $\displaystyle
2r_{Q_{e}}\mathcal{Q}_{3SN}^{2}-\frac{4}{\sqrt{6}}s_{Q_{e}}\mathcal{Q}_{3SN}^{3}$
$\displaystyle+(4u_{Q_{e},1}/3+4u_{Q_{e},2}-u_{Q_{e},3})\mathcal{Q}_{3SN}^{4}~{}.$
(59)

For positive quartic term, we get a three sub-lattice nematic via a first
order phase transition.

The interaction term between the singlet and doublet quadrupolar modes (Eq.
58) indicates that ordering of the doublet mode would give rise to a parasitic
singlet quadrupole expectation value of

$Q_{a}^{\alpha\beta}\approx-\frac{s_{Q_{a}Q_{e}}}{2r_{Q_{a}}}[(Q_{e1}^{\beta\gamma}Q_{e1}^{\gamma\alpha}+Q_{e2}^{\beta\gamma}Q_{e2}^{\gamma\alpha})-\frac{Tr(\mathbf{Q}_{\bf
e}^{2})\delta^{\alpha\beta}}{3}].$

However, for the particular three sublattice nematic (Eqs. 45 and 46) such
parasitic ferronematic order vanishes.

##### Dipole-quadrupole interaction :

The general form of interaction terms between dipoles and quadrupoles, due to
time reversal symmetry, consists of even powers of dipoles and any powers for
quadrupoles consistent with other symmetries. Therefore the lowest order term
has the form $\sim m^{2}Q$. This implies that pure quadrupolar ordering
renormalises the mass for dipoles while dipolar ordering gives rise to
parasitic quadrupole moment. The nature of the parasitic quadrupole moments is
obtained by extremising the leading order dipole-quadrupole interaction :

$\displaystyle\mathcal{F}_{mQ}^{\text{int}}=$ $\displaystyle
s_{m_{a}Q_{a}}Q_{a}^{\alpha\beta}\bigg{(}m_{a}^{\alpha}m_{a}^{\beta}-\frac{1}{3}(\mathbf{m}_{a}.\mathbf{m}_{a})\delta^{\alpha\beta}\bigg{)}$
$\displaystyle+s_{m_{e}Q_{a}}Q_{a}^{\alpha\beta}\bigg{(}m_{e1}^{\alpha}m_{e1}^{\beta}+m_{e2}^{\alpha}m_{e2}^{\beta}$
$\displaystyle-\frac{1}{3}(\mathbf{m}_{\bf e}.\mathbf{m}_{\bf
e})\delta^{\alpha\beta}\bigg{)}+s_{m_{e}Q_{e}}\bigg{[}-2Q_{e1}^{\alpha\beta}\bigg{(}\frac{m_{e1}^{\alpha}m_{e2}^{\beta}}{2}$
$\displaystyle+\frac{m_{e2}^{\alpha}m_{e1}^{\beta}}{2}-\frac{1}{3}(\mathbf{m}_{e1}.\mathbf{m}_{e2})\delta^{\alpha\beta}\bigg{)}+Q_{e2}^{\alpha\beta}\bigg{(}m_{e2}^{\alpha}m_{e2}^{\beta}$
$\displaystyle-
m_{e1}^{\alpha}m_{e1}^{\beta}-\frac{1}{3}(\mathbf{m}_{e2}.\mathbf{m}_{e2}-\mathbf{m}_{e1}.\mathbf{m}_{e1})\delta^{\alpha\beta}\bigg{)}\bigg{]}$
$\displaystyle+u_{m_{a}Q_{e}}m_{a}^{\alpha}m_{a}^{\beta}(Q_{e1}^{\alpha\gamma}Q_{e1}^{\gamma\beta}+Q_{e2}^{\alpha\gamma}Q_{e2}^{\gamma\beta})$
$\displaystyle+s_{m_{a}m_{e}Q_{e}}\bigg{[}Q_{e1}^{\alpha\beta}\bigg{(}\frac{m_{a}^{\alpha}m_{e1}^{\beta}+m_{e1}^{\alpha}m_{a}^{\beta}}{2}$
$\displaystyle-\frac{1}{3}(\mathbf{m}_{a}.\mathbf{m}_{e1})\delta^{\alpha\beta}\bigg{)}+Q_{e2}^{\alpha\beta}\bigg{(}\frac{m_{a}^{\alpha}m_{e2}^{\beta}+m_{e2}^{\alpha}m_{a}^{\beta}}{2}$
$\displaystyle-\frac{1}{3}(\mathbf{m}_{a}.\mathbf{m}_{e2})\delta^{\alpha\beta}\bigg{)}\bigg{]}~{}.$
(60)

For the ferromagnetic ordering, this leads to

$Q_{a}^{\alpha\beta}\approx-\frac{s_{m_{a}Q_{a}}}{2r_{Q_{a}}}\bigg{(}m_{a}^{\alpha}m_{a}^{\beta}-\frac{1}{3}(\mathbf{m}_{a}.\mathbf{m}_{a})\delta^{\alpha\beta}\bigg{)}~{},$
(61)

while for the $120^{\circ}$ spiral order we get

$\displaystyle
Q_{a}^{\alpha\beta}\approx-\frac{s_{m_{e}Q_{a}}}{2r_{Q_{a}}}\bigg{(}m_{e1}^{\alpha}m_{e1}^{\beta}+m_{e2}^{\alpha}m_{e2}^{\beta}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{1}{3}(\mathbf{m}_{e1}.\mathbf{m}_{e1}+\mathbf{m}_{e2}.\mathbf{m}_{e2})\delta^{\alpha\beta}\bigg{)}$
$\displaystyle
Q_{e1}^{\alpha\beta}\approx\frac{s_{m_{e}Q_{e}}}{r_{Q_{e}}}\bigg{(}\frac{m_{e1}^{\alpha}m_{e2}^{\beta}+m_{e2}^{\alpha}m_{e1}^{\beta}}{2}-\frac{1}{3}(\mathbf{m}_{e1}.\mathbf{m}_{e2})\delta^{\alpha\beta}\bigg{)}$
$\displaystyle
Q_{e2}^{\alpha\beta}\approx-\frac{s_{m_{e}Q_{e}}}{2r_{Q_{e}}}\bigg{(}m_{e2}^{\alpha}m_{e2}^{\beta}-m_{e1}^{\alpha}m_{e1}^{\beta}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{1}{3}(\mathbf{m}_{e2}.\mathbf{m}_{e2}-\mathbf{m}_{e1}.\mathbf{m}_{e1})\delta^{\alpha\beta}\bigg{)}$
(62)

The above parasitic moments for the ferromagnetic and $120^{\circ}$ spiral
order are in accordance with what we expect from the wave functions listed in
appendix A.

####  IV.3.2 Coupling with the magnetic field

Magnetic field ($\mathbf{h}=h\hat{\mathbf{h}}$) couples to the spins via the
usual Zeeman term

$\mathcal{F}_{hm}=g_{m}\mathbf{h}.\mathbf{m}_{a}$ (63)

such that only the ${\bf m}_{a}$ mode couples linearly to the uniform magnetic
field. Such linear coupling evidently favours the polarised phase along the
magnetic field.

As quadrupoles are even under time reversal, the coupling takes the form:

$\mathcal{F}_{hQ}=g_{Q}Q_{a}^{\alpha\beta}\bigg{(}h^{\alpha}h^{\beta}-\frac{1}{3}\delta^{\alpha\beta}(\mathbf{h}.\mathbf{h})\bigg{)}$
(64)

Note again, only the $Q_{a}$ modes couple to bilinears of the uniform magnetic
field. This can be seen from the microscopic term of the form
$h^{\alpha}h^{\beta}\sum_{i}Q_{i}^{\alpha\beta}$ (where $i$ is the site index)
and writing it in terms of the irreducible representations in Table 1. The
effect (to linear order) of small uniform magnetic field on the ferronematic
and the three sublattice nematic order is summarised below.

##### Ferronematic order:

For the ferronematic order (Eq. 43) due to $\mathcal{F}_{hm}$ (Eq. 63), a
uniform magnetic moment proportional to the magnetic field develops
$m_{a}^{\alpha}\approx-\frac{g_{m}h^{\alpha}}{2r_{m_{a}}}$. Due to this
uniform magnetic moment, the coupling between the ferronematic order parameter
and the singlet dipolar mode in Eq. 60 becomes,

$\mathcal{F}_{mQ,FN}^{\text{int}}=\frac{\sqrt{3}s_{m_{a}Q_{a}}\mathcal{Q}_{FN}g_{m}^{2}h^{2}}{4r_{m_{a}}^{2}}\bigg{(}(\hat{\mathbf{h}}.\mathbf{n})^{2}-\frac{1}{3}\bigg{)}.$
(65)

When $g_{Q}\approx 0$ in Eq. 64, the above term alone decides whether the
ferronematic directors turn perpendicular(parallel) to the magnetic field for
$s_{m_{a}Q_{a}}<0$ ($s_{m_{a}Q_{a}}>0$) (since $\mathcal{Q}_{FN}<0$ for disk
like ferronematic order.Tóth (2011); Penc and Läuchli (2011)) With $g_{Q}\neq
0$, $\mathcal{F}_{hQ}$ becomes,

$\mathcal{F}_{hQ,FN}=\sqrt{3}\mathcal{Q}_{FN}g_{Q}h^{2}\bigg{(}(\hat{\mathbf{h}}.\mathbf{n})^{2}-\frac{1}{3}\bigg{)}~{},$
(66)

which competes with Eq. 65 in deciding whether director is parallel or
perpendicular to the magnetic field. Extremising the Gaussian free energy (Eq.
56) along with these terms, we get

$\displaystyle\mathcal{Q}_{FN}=-\frac{\sqrt{3}(g_{Q}+s_{m_{a}Q_{a}}g_{m}^{2}/4r_{m_{a}}^{2})}{6r_{Q_{a}}}h^{2}$
(67)

while is the magnetic field induced ferronematic ordering where we have
assumed for concreteness $s_{m_{a}Q_{a}},g_{Q}>0$.

##### Three sublattice nematic:

For the three sublattice nematic too the uniform magnetic moment proportional
to and in the direction of the magnetic field develops with
$m_{a}^{\alpha}\approx-\frac{g_{m}h^{\alpha}}{2r_{m_{a}}}$. In contrast to the
ferronematic case, however, here the application of magnetic field also
results in doublet dipolar modes via ${\bf m}_{a}$. This can be seen by
extremising the $s_{m_{a}m_{e}Q_{e}}$ (Eq. 60) :

$\displaystyle m_{\bf e}^{\alpha}$
$\displaystyle\approx-\frac{s_{m_{a}m_{e}Q_{e}}Q_{\bf
e}^{\alpha\beta}m_{a}^{\beta}}{2r_{m_{e}}}=g_{m}\frac{s_{m_{a}m_{e}Q_{e}}Q_{\bf
e}^{\alpha\beta}h^{\beta}}{4r_{m_{e}}r_{m_{a}}}$

for ${\bf e}=(e1,e2)$. It is evident that the doublet modes depend on the
orientation of the magnetic field relative to the directors. When the magnetic
field is along any one of the orthogonal directors of the three sublattice
nematic, we see from Eq. 45, Eq. 46 and Eq. LABEL:me1me2duetoh that
$\mathbf{m}_{e1}$ and $\mathbf{m}_{e2}$ are also in the direction of the
magnetic field. This gives rise to a collinear sinusoidal magnetization in
addition to the uniform magnetization (due to $\mathbf{m}_{a}$) along the
direction of the magnetic field. Considering the general case where the
directors are $\mathbf{n}_{1}=\hat{\mathbf{x}}$,
$\mathbf{n}_{2}=\hat{\mathbf{y}}$ and $\mathbf{n}_{3}=\hat{\mathbf{z}}$ and
the magnetic field is at a general inclination to the directors,
$\mathbf{h}=h(\sin{\theta}\cos{\phi}\hat{\mathbf{x}}+\sin{\theta}\sin{\phi}\hat{\mathbf{y}}+\cos{\theta}\hat{\mathbf{z}})$,
using Eq. 38 we see that the resulting doublet and singlet magnetisations
correspond to,
$\langle\mathbf{S}_{1}\rangle=(a_{m}\sin{\theta}\cos{\phi}\hat{\mathbf{x}}+b_{m}\sin{\theta}\sin{\phi}\hat{\mathbf{y}}+b_{m}\cos{\theta}\hat{\mathbf{z}})$,
$\langle\mathbf{S}_{2}\rangle=(b_{m}\sin{\theta}\cos{\phi}\hat{\mathbf{x}}+a_{m}\sin{\theta}\sin{\phi}\hat{\mathbf{y}}+b_{m}\cos{\theta}\hat{\mathbf{z}})$
and
$\langle\mathbf{S}_{3}\rangle=(b_{m}\sin{\theta}\cos{\phi}\hat{\mathbf{x}}+b_{m}\sin{\theta}\sin{\phi}\hat{\mathbf{y}}+a_{m}\cos{\theta}\hat{\mathbf{z}})$
where
$a_{m}=\frac{g_{m}h}{3}\left(-\frac{\sqrt{3}}{2r_{m_{a}}}+\frac{2s_{m_{a}m_{e}Q_{e}}Q_{3SN}}{4r_{m_{e}}r_{m_{a}}}\right)$
and
$b_{m}=\frac{g_{m}h}{3}\left(-\frac{\sqrt{3}}{2r_{m_{a}}}-\frac{s_{m_{a}m_{e}Q_{e}}Q_{3SN}}{4r_{m_{e}}r_{m_{a}}}\right)$.
The resultant phase is a combination of uniform polarization along the field
as well as sub-lattice dependent polarisation along the three orthogonal
direction.

####  IV.3.3 Elastic term and coupling of strains to dipoles and quadrupoles

Finally, we turn to the contributions to the free energy due to elastic fields
as well as spin-phonon coupling.

##### Elastic energy:

The harmonic elastic energy for the singlet and the doublet strain modes is

$\mathcal{F}_{\text{elastic}}=\frac{1}{2}(c_{a}\varepsilon_{a}^{2}+c_{e}\varepsilon_{e}^{2})~{},$
(69)

where $c_{a}$ and $c_{e}$ are two independent elastic constants.

##### Coupling between dipoles, quadrupoles and strains :

The coupling term between dipoles and strain fields is given by

$\mathcal{F}_{\varepsilon m}=\mathcal{F}_{\varepsilon
ma}+\mathcal{F}_{\varepsilon me}~{},$ (70)

where

$\displaystyle\mathcal{F}_{\varepsilon ma}$ $\displaystyle=$
$\displaystyle[\mu_{a1}\varepsilon_{a}+\frac{\mu_{a}}{2}\varepsilon_{a}^{2}+\frac{\mu_{e}}{2}\varepsilon_{e}^{2}]p_{1}(\mathbf{m}_{a}.\mathbf{m}_{a})$
(71) $\displaystyle\mathcal{F}_{\varepsilon me}$ $\displaystyle=$
$\displaystyle[\mu_{a1}\varepsilon_{a}+\frac{\mu_{a}}{2}\varepsilon_{a}^{2}+\frac{\mu_{e}}{2}\varepsilon_{e}^{2}]p_{2}(\mathbf{m}_{\bf
e}.\mathbf{m}_{\bf e})$ (72)

denotes the coupling between the singlet and doublet dipolar modes.

Similarly, the coupling term between quadrupoles and strains is given by

$\mathcal{F}_{\varepsilon Q}=\mathcal{F}_{\varepsilon
Qa}+\mathcal{F}_{\varepsilon Qe}~{},$ (73)

where

$\displaystyle\mathcal{F}_{\varepsilon Qa}$ $\displaystyle=$
$\displaystyle[\tilde{\mu}_{a1}\varepsilon_{a}+\frac{\tilde{\mu}_{a}}{2}\varepsilon_{a}^{2}+\frac{\tilde{\mu}_{e}}{2}\varepsilon_{e}^{2}]p_{3}\text{Tr}\mathbf{Q}_{a}^{2}$
(74) $\displaystyle\mathcal{F}_{\varepsilon Qe}$ $\displaystyle=$
$\displaystyle[\tilde{\mu}_{a1}\varepsilon_{a}+\frac{\tilde{\mu}_{a}}{2}\varepsilon_{a}^{2}+\frac{\tilde{\mu}_{e}}{2}\varepsilon_{e}^{2}]p_{4}\text{Tr}(\mathbf{Q}_{\bf
e}.\mathbf{Q}_{\bf e})$ (75)

In both the cases of dipolas and quadrupoles, the free energy allows for
linear coupling with the singlet strain field with bilinear of the
dipole/quadrupole fields. This would lead to lattice distortions upon dipolar
or quadrupolar ordering while the quadratic terms (in elastic fields) lead to
renormalisation of elastic constants and hence the sound speed.

##  V Fractional change in sound speed and fractional change in length

###  V.1 Fractional change in length

Fractional change in length along direction
$\hat{r}\equiv(\cos\theta,\sin\theta)$ (where $\theta$ is the angle with
respect to the Cartesian $x$-axis in Fig. 1) can be obtained from the
Cartesian strain fields $\varepsilon_{ij}$ :

$\displaystyle\bigg{(}\frac{\Delta L}{L}\bigg{)}_{\hat{r}}$
$\displaystyle=\sum_{ij}\varepsilon_{ij}\hat{r}_{i}\hat{r}_{j}=\varepsilon_{a}+\varepsilon_{e2}\cos(2\theta)+\varepsilon_{e1}\sin(2\theta).$
(76)

Due to the linear coupling term between the uniform strain field,
$\varepsilon_{a}$ and the dipole/nematic bilinears (Eqs. 71, 72, 74 and 75),
this leads to isotropic magnetostriction of the triangular lattice given by

$\displaystyle\bigg{(}\frac{\Delta L}{L}\bigg{)}_{\hat{r}}=\varepsilon_{a}$
$\displaystyle\approx$
$\displaystyle-\frac{1}{c_{a}}\bigg{[}\mu_{a1}\bigg{(}p_{1}(\mathbf{m}_{a}.\mathbf{m}_{a})+p_{2}(\mathbf{m}_{\bf
e}.\mathbf{m}_{\bf e})\bigg{)}$ (77) $\displaystyle+$
$\displaystyle\tilde{\mu}_{a1}\bigg{(}p_{3}\text{Tr}\mathbf{Q}_{a}^{2}+p_{4}\text{Tr}(\mathbf{Q}_{\bf
e}.\mathbf{Q}_{\bf e})\bigg{)}\bigg{]}$

Focussing on the $120^{\circ}$ spiral order, Eq. 76 reduces to

$\displaystyle\bigg{(}\frac{\Delta
L}{L}\bigg{)}_{\hat{r}}=-\frac{\mu_{a1}p_{2}}{c_{a}}(\mathbf{m}_{\bf
e}.\mathbf{m}_{\bf e})$ (78)

such that below the critical point ($r_{m_{e}}<0$ in Eq. 52) when ${\bf
m}_{e}\sim\sqrt{-r_{m_{e}}/u_{m_{e}}}$,

$\displaystyle\bigg{(}\frac{\Delta L}{L}\bigg{)}_{\hat{r}}\propto|r_{m_{e}}|$
(79)

as expected from the lowest symmetry allowed coupling. Thus for a thermal
phase transition where $r_{m_{e}}\propto(T-T_{c})$, one expects a linear
turning on of the magnetiostrictive distortion with measurable consequence for
thermal expansion experiments.

A more startling effect occurs across a spin nematic transition. In particular
for a ferronematic, the above expression reduces to

$\displaystyle\bigg{(}\frac{\Delta
L}{L}\bigg{)}_{\hat{r}}=-\frac{\mu_{a1}p_{3}}{c_{a}}Tr{\bf
Q}_{a}^{2}=-\frac{2\mu_{a1}p_{3}}{c_{a}}\mathcal{Q}_{FN}^{2}$ (80)

such that across the discontinuous nematic transition, there is a simultaneous
jump in the lattice volume.

On turning on the magnetic field inside either the spiral or the ferronematic,
there are added contributions due to the Zeeman term (Eq. 63) resulting in
non-zero $m_{a}^{\alpha}=-\frac{g_{m}h^{\alpha}}{2r_{m_{a}}}$ and
$\mathcal{Q}_{FN}$ (Eq. 67) such that we have

$\displaystyle\bigg{(}\frac{\Delta L}{L}\bigg{)}_{\hat{r}}=$
$\displaystyle-\frac{\mu_{a1}p_{2}}{c_{a}}(\mathbf{m}_{\bf e}.\mathbf{m}_{\bf
e})-\frac{\mu_{a1}p_{1}g_{m}^{2}}{4c_{a}r_{m_{a}}^{2}}h^{2}$
$\displaystyle-\frac{\tilde{\mu}_{a1}p_{3}(g_{Q}+s_{m_{a}Q_{a}}g_{m}^{2}/4r_{m_{a}}^{2})^{2}}{6r_{Q_{a}}^{2}c_{a}}~{}~{}h^{4}$
(81)

for the spiral and

$\displaystyle\bigg{(}\frac{\Delta L}{L}\bigg{)}_{\hat{r}}=$
$\displaystyle-\frac{2\mu_{a1}p_{3}}{c_{a}}\mathcal{Q}_{FN}^{2}-\frac{\mu_{a1}p_{1}g_{m}^{2}}{4c_{a}r_{m_{a}}^{2}}h^{2}$
$\displaystyle-\frac{\tilde{\mu}_{a1}p_{3}(g_{Q}+s_{m_{a}Q_{a}}g_{m}^{2}/4r_{m_{a}}^{2})^{2}}{6r_{Q_{a}}^{2}c_{a}}~{}~{}h^{4}$
(82)

for the ferronematic.

While both the forms predict a similar mixture of $h^{2}$ and $h^{4}$
dependence, one expects that near the ferronematic phase transition
$r_{Q_{a}}\rightarrow 0$ such that the $h^{4}$ dependence becomes more
pronounced on approaching the ferronematic phase transition. Similar results
hold for the three sublattice nematic.

###  V.2 Fractional change in sound speed

The fractional change in sound speed, $v$, is related to the change in elastic
constants, $c$, (Eq. 69) as :

$\frac{\Delta v}{v}\propto\frac{\Delta c}{c}~{},$ (83)

The component of the elastic tensor depends on the propagation direction and
polarization of the sound and hence is generically a linear combination of
$c_{a}$ and $c_{e}$.

The renormalisation of these two elastic constants are due to coupling with
the dipole and nematic orders are readily obtained from Eqs. 70 and 73 as :

$\displaystyle\frac{\Delta c_{a}}{c_{a}}$ $\displaystyle=$
$\displaystyle\frac{1}{c_{a}}\bigg{[}\mu_{a}\bigg{(}p_{1}(\mathbf{m}_{a}.\mathbf{m}_{a})+p_{2}(\mathbf{m}_{\bf
e}.\mathbf{m}_{\bf e})\bigg{)}$ $\displaystyle+$
$\displaystyle\tilde{\mu}_{a}\bigg{(}p_{3}\text{Tr}\mathbf{Q}_{a}^{2}+p_{4}\text{Tr}(\mathbf{Q}_{\bf
e}.\mathbf{Q}_{\bf e})\bigg{)}\bigg{]}$ $\displaystyle\frac{\Delta
c_{e}}{c_{e}}$ $\displaystyle=$
$\displaystyle\frac{1}{c_{e}}\bigg{[}\mu_{e}\bigg{(}p_{1}(\mathbf{m}_{a}.\mathbf{m}_{a})+p_{2}(\mathbf{m}_{\bf
e}.\mathbf{m}_{\bf e})\bigg{)}$ (84) $\displaystyle+$
$\displaystyle\tilde{\mu}_{e}\bigg{(}p_{3}\text{Tr}\mathbf{Q}_{a}^{2}+p_{4}\text{Tr}(\mathbf{Q}_{\bf
e}.\mathbf{Q}_{\bf e})\bigg{)}\bigg{]}$

Therefore similar to the case of fractional change in length, the fraction
change in sound speed would be continuous across the $120^{\circ}$ spiral
ordering transition and discontinuous across the nematic transitions. The
magnetic field dependence below the ordering transitions will have a similar
dependence.

Observing jumps in fractional change in length and fractional change in sound
speed in the absence of any magnetization, along with the enhanced $h^{4}$
magnetic field dependence would strongly favour the case for spin-nematic
orders.

##  VI Summary and Outlook

In this work, we have studied the possible elastic signatures of a spin-
nematic, possibly realised in the triangular lattice compound NiGa2S4.
Inspired by recent Raman scattering experiments Valentine  _et al._ (2020)
which indicate substantial spin-phonon coupling in the material, and the
general expectation that the time reversal symmetric spin-nematic order
parameter can linearly couple to the elastic degrees of freedom, we show that–
(1) inside such a spin-nematic phase, the interaction between the nematic
Goldstone boson in a ferronematic phase and the phonons lead to a powerlaw
dependence of the fractional change in sound speed, i.e., $\Delta v/v\propto
T^{3}$, (2) across the spin-nematic phase transition, there is a finite jump
in both magnetostriction as well as $\Delta v/v$ stemming from the
discontinuous nature of the nematic transition, and (3) possible enhanced
$h^{4}$ dependence of the magnetostriction and $\Delta v/v$ on the magnetic
field just below the nematic transition. These signatures, along with
thermodynamic measurements such as low temperature powerlaw specific heat
Nakatsuji  _et al._ (2005); Tsunetsugu and Arikawa (2006); Bhattacharjee  _et
al._ (2006); Läuchli  _et al._ (2006); Stoudenmire  _et al._ (2009) and
absence of long range spin correlations provides important steps in
characterising the spin-nematic phase in particular and paves a way for
concrete experimental signatures for higher-moment magnetic orderings that
have proved quite challenging in spite of several candidate materials. It is
useful to note that light scattering can also act as a complementary
spectroscopic probe for spin-nematic as shown in Ref. Michaud  _et al._ ,
2011.

In context of NiGa2S4, we have used the minimal nearest neighbour bilinear-
biquadratic spin-1 Hamiltonian (Eq. 2) to understand the the low temperature
elastic signatures. In the actual material, in addition to nearest neighbour
bilinear and biquadratic terms, a third neighbour antiferromagnetic bilinear
spin interaction appears to be relevant to understand the detailed physics.
Here, however, we have neglected such terms as we are interested in the
elastic response of the spin-nematic where the effect of such terms are
expected to be secondary and as far as the temperature dependence is
concerned. We further note that in addition to spin-rotation invariant terms,
the microscopics of NiGa2S4 can also admit small single-ion anisotropies as
well as DM interactions. While the former can pin the nematic
director,Bhattacharjee  _et al._ (2006) the latter can pin the chirality of
the spiral state for the DM vector pointing perpendicular to the plane (see
Appendix B). Therefore the former can gap out the nematic Goldstone modes and
change the powerlaw dependence to an exponential damping for $\Delta v/v$.
However, since such effects have not been measured in the low temperature
specific heat,Nakatsuji  _et al._ (2005); Tsunetsugu and Arikawa (2006);
Bhattacharjee  _et al._ (2006); Läuchli  _et al._ (2006); Stoudenmire  _et
al._ (2009) we neglect such single-ion anisotropies.

Several studies on ultrasound properties on candidates for multipolar order in
SOC coupled systems Patri  _et al._ (2019) or in-field spin-nematic Gen  _et
al._ (2019) exists. We hope our results will generate interest in probing
ultrasound renormalisation in NiGa2S4 which is a rare spin-nematic candidate
in a system described by a spin rotation symmetric spin Hamiltonian to a very
good approximation.

###### Acknowledgements.

The authors thank R. Moessner, S. Nakatsuji, N. Drichko and S. Zherlitsyn for
discussion. SB acknowledges adjunct fellow program at SNBNCBS, Kolkata for
hospitality. The authors acknowledge funding from Max Planck Partner group
Grant at ICTS, Swarna Jayanti fellowship grant of SERB-DST (India) Grant No.
SB/SJF/2021-22/12 and the Department of Atomic Energy, Government of India,
under Project No. RTI4001.

##  Appendix A  Wave function for various spin-1 orders

The mean-field product wave functions for spin-1 magnets for the dipole and
nematic orders Läuchli  _et al._ (2006); Penc and Läuchli (2011); Tóth (2011)
relevant to this work is best understood in the basis :

$\displaystyle\ket{x}=i\frac{\ket{1}-\ket{\bar{1}}}{\sqrt{2}},~{}\ket{y}=\frac{\ket{1}+\ket{\bar{1}}}{\sqrt{2}},~{}\ket{z}=-i\ket{0}.$
(85)

The on-site spin and the quadrupole operators are

$\displaystyle
S^{\bar{\mu}}=-i\sum_{\bar{\nu}\bar{\lambda}}\epsilon_{\bar{\mu}\bar{\nu}\bar{\lambda}}\ket{\bar{\nu}}\bra{\bar{\lambda}};~{}Q^{\bar{\mu}\bar{\nu}}=\frac{1}{3}\delta^{\bar{\mu}\bar{\nu}}-\frac{\ket{\bar{\nu}}\bra{\bar{\mu}}}{2}-\frac{\ket{\bar{\mu}}\bra{\bar{\nu}}}{2}$
(86)

where $\epsilon_{\bar{\mu}\bar{\nu}\bar{\lambda}}$ denotes the Levi Civita
tensor and $\bar{\mu}$, $\bar{\nu}$, $\bar{\lambda}=x,y,z$. The most general
single spin-1 state is

$\ket{\mathbf{w}}=\sum_{\bar{\mu}}w_{\bar{\mu}}\ket{\bar{\mu}}~{},$ (87)

where $w_{\bar{\mu}}$ are the three components of a complex vector

$\mathbf{w}=\mathbf{u}+i\mathbf{v}$ (88)

with ${\bf u},{\bf v}\in\mathcal{R}^{3}$ and
$\mathbf{u}\cdot\mathbf{u}+\mathbf{v}\cdot\mathbf{v}=1$ with the overall phase
fixed by $\mathbf{u}\cdot\mathbf{v}=0$. The expectation of spin and quadrupole
operators are :

$\displaystyle\bra{\mathbf{w}}\mathbf{S}\ket{\mathbf{w}}$ $\displaystyle=$
$\displaystyle 2\mathbf{u}\times\mathbf{v}$ (89)
$\displaystyle\bra{\mathbf{w}}Q^{\bar{\mu}\bar{\nu}}\ket{\mathbf{w}}$
$\displaystyle=$
$\displaystyle\frac{1}{3}\delta^{\bar{\mu}\bar{\nu}}-u^{\bar{\mu}}u^{\bar{\nu}}-v^{\bar{\mu}}v^{\bar{\nu}}~{}.$
(90)

The mean field wave function for different magnetic and nematic ordered states
are direct products of the onsite wave functions,
$\ket{\psi}_{MF}=\otimes_{i}\ket{\mathbf{w}_{i}}$.

##### Ferromagnetic order:

The ferromagnetic order in the $\hat{\mathbf{z}}$ direction is given by

$\displaystyle\mathbf{u}_{i}$
$\displaystyle=\frac{1}{\sqrt{2}}(\text{sin}\theta\text{
}\hat{\mathbf{x}}+\text{cos}\theta\text{
}\hat{\mathbf{y}}),~{}~{}\mathbf{v}_{i}=\frac{-1}{\sqrt{2}}(\text{cos}\theta\text{
}\hat{\mathbf{x}}-\text{sin}\theta\text{ }\hat{\mathbf{y}}),$ (91)

such that the spin expectation value (Eq. 89) is given by

$\bra{\mathbf{w}}\mathbf{S}\ket{\mathbf{w}}=\text{ }\hat{\mathbf{z}}$ (92)

which characterises the fully polarised ferromagnetic state. The on-site
parasitic quadrupole order is

$\bra{\mathbf{w}}Q^{\bar{\mu}\bar{\nu}}\ket{\mathbf{w}}=\begin{pmatrix}-\frac{1}{6}&0&0\\\
0&-\frac{1}{6}&0\\\ 0&0&\frac{1}{3}\\\ \end{pmatrix}~{},$ (93)

in the $\\{\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}}\\}$ basis and
corresponds to uniaxial parasitic ferronematic order with director along
$\hat{\mathbf{z}}$.

#####  $120^{\circ}$ spiral order:

For a $120^{\circ}$ spiral order :

$\displaystyle\mathbf{u}_{1}$
$\displaystyle=\frac{1}{\sqrt{2}}\hat{\mathbf{z}}\text{ ,
}\mathbf{v}_{1}=-\frac{1}{\sqrt{2}}\hat{\mathbf{y}}$
$\displaystyle\mathbf{u}_{2}$
$\displaystyle=\frac{1}{\sqrt{2}}\hat{\mathbf{z}}\text{ ,
}\mathbf{v}_{2}=-\frac{1}{\sqrt{2}}\bigg{(}\frac{\sqrt{3}}{2}\hat{\mathbf{x}}-\frac{1}{2}\hat{\mathbf{y}}\bigg{)}$
$\displaystyle\mathbf{u}_{3}$
$\displaystyle=\frac{1}{\sqrt{2}}\hat{\mathbf{z}}\text{ ,
}\mathbf{v}_{3}=-\frac{1}{\sqrt{2}}\bigg{(}-\frac{\sqrt{3}}{2}\hat{\mathbf{x}}-\frac{1}{2}\hat{\mathbf{y}}\bigg{)}$
(94)

such that the spin expectation value is given by Eq. 4 and corresponds to Fig.
4. Again, each on-site parasitic quadrupolar moment would be uniaxial with the
director being in the direction of the respective magnetic moment.

Pure nematic orders require
$\bra{\mathbf{w}}\mathbf{S}\ket{\mathbf{w}}=2\mathbf{u}\times\mathbf{v}=0$ in
addition to $\mathbf{u}\cdot\mathbf{v}=0$ which correspond to either
$\mathbf{v}=0$ or $\mathbf{u}=0$ i.e., $\mathbf{w}$ is real or imaginary.
Considering $\mathbf{v}=0$, the expectation value of the quadrupole operator
(Eq. 90) is,

$\bra{\mathbf{w}}Q^{\bar{\mu}\bar{\nu}}\ket{\mathbf{w}}=\frac{1}{3}\delta^{\bar{\mu}\bar{\nu}}-u^{\bar{\mu}}u^{\bar{\nu}}~{},$
(95)

which is a uniaxial nematic with the director along $\mathbf{u}$.

##### Ferronematic order:

All sites have

$\displaystyle\mathbf{u}$
$\displaystyle=\hat{\mathbf{n}},~{}~{}~{}~{}\mathbf{v}=0$ (96)

with $\hat{\mathbf{n}}$ being the director (Fig.1).

##### Three sublattice nematic order:

The three sites of the triangular unit (see Fig. 1)) have :

$\displaystyle\mathbf{u}_{1}$ $\displaystyle=\hat{\mathbf{n}}_{1}\text{ ,
}\mathbf{v}_{1}=0$ $\displaystyle\mathbf{u}_{2}$
$\displaystyle=\hat{\mathbf{n}}_{2}\text{ , }\mathbf{v}_{2}=0$
$\displaystyle\mathbf{u}_{3}$ $\displaystyle=\hat{\mathbf{n}}_{3}\text{ ,
}\mathbf{v}_{3}=0~{}$ (97)

where $\hat{\mathbf{n}}_{1}$ , $\hat{\mathbf{n}}_{2}$ and
$\hat{\mathbf{n}}_{3}$ are mutually orthogonal.

##### Co-existing nematic order and collinear sinusoidal dipolar order:

The spins for a three sublattice collinear sinusoidal dipolar order (with co-
existing nematic order discussed in Section IV.1.2) can be expressed as
$\langle\mathbf{S}_{i}\rangle=\text{sin}\eta_{i}\text{ }\hat{m}$ (with
$\eta_{i}=\mathbf{q}.\mathbf{r}_{i}+\pi/6+\phi$). Using Eq. 91, the wave
functions with the above order are :

$\displaystyle\mathbf{u}_{i}$
$\displaystyle=\text{cos}(\eta_{i}/2)(\text{sin}\theta_{i}\text{
}\hat{\mathbf{m}}_{\perp_{1}}+\text{cos}\theta_{i}\text{
}\hat{\mathbf{m}}_{\perp_{2}})$ $\displaystyle\mathbf{v}_{i}$
$\displaystyle=-\text{sin}(\eta_{i}/2)(\text{cos}\theta_{i}\text{
}\hat{\mathbf{m}}_{\perp_{1}}-\text{sin}\theta_{i}\text{
}\hat{\mathbf{m}}_{\perp_{2}})$ (98)

where
$\hat{\mathbf{m}}_{\perp_{1}}\times\hat{\mathbf{m}}_{\perp_{2}}=\hat{\mathbf{m}}$,
$i=1,2,3$ and $\eta_{2}=\eta_{1}+2\pi/3$ , $\eta_{3}=\eta_{1}+4\pi/3$. Such an
order would feature sites having partially polarized spins. Unlike quadrupole
moments for fully polarized spins (Eq. 93), for the case of partially
polarized spins one can see that the expectation values of the on-site
quadrupolar moment (using Eq. 90) would depend on the choice of $\theta_{i}$
for the site and the quadrupole moment would be bi-axial - with the
eigenvalues being
$\\{1/3,1/3-|\mathbf{u}_{i}|^{2},1/3-|\mathbf{v}_{i}|^{2}\\}$ for the
eigenvectors being
$\\{\hat{\mathbf{m}},\hat{\mathbf{u}}_{i},\hat{\mathbf{v}}_{i}\\}$
respectively.

##  Appendix B  Mean field T=0 phase diagram and DM term

The mean field phase diagram of Eq. 2 was obtained in Ref. Bhattacharjee  _et
al._ , 2006. On adding to it, a DM term given by
$H_{DM}=\mathbf{D}_{ij}\cdot(\mathbf{S}_{i}\times\mathbf{S}_{j})$, where
$\mathbf{D}_{ij}=-\mathbf{D}_{ji}$ and $\mathbf{D}_{ij}=D\hat{\mathbf{z}}$
vector pointing out of the triangular lattice plane, the mean field phase
diagram is given by Fig. 4 for $J>0$ and $K<0$. The main role of the DM term
is to lift the degeneracy of the clockwise and anti-clockwise spirals. The
different order parameters are plotted in Fig. 5.

Figure 4: Mean field $T=0$ phase for bi-linear biquadratic spin-1 ($J>0$ and
$K<0$) Hamiltonian with the DM interaction. Dashed line denotes first order
transition and continuous lines are for continuous transition. Orientation 1
has spiral
$\mathbf{q}=-\frac{2\pi}{3}\hat{\mathbf{x}}-\frac{2\pi}{\sqrt{3}}\hat{\mathbf{y}}$,
Orientation 2 has spiral
$\mathbf{q}=\frac{2\pi}{3}\hat{\mathbf{x}}+\frac{2\pi}{\sqrt{3}}\hat{\mathbf{y}}$
(lattice constant taken to be unity). Figure 5: Results of the self
consistency calculations (for $J>0$ and $K<0$) : a),b): Energy and order
parameters for $D=0$, c),d): Energy and order parameters for $D/J=0.42$

##  Appendix C Linear Spin-nematic wave theory

###  C.1 The ferronematic

As discussed in Sec. III.1, the Goldstone modes out of a ferronematic with
director along the $\hat{\bf z}$ direction can be captured by defining two
bosons at each site $i$, $b_{i1}^{\dagger}$ and $b_{i2}^{\dagger}$. For a
faithful description of the spin-1 Hilbert space, the bosons are subject to
the constraints –
$b_{i1}^{\dagger}b_{i1}^{\dagger}=b_{i1}b_{i1}=b_{i2}^{\dagger}b_{i2}^{\dagger}=b_{i2}b_{i2}=b_{i1}^{\dagger}b_{i2}^{\dagger}=b_{i1}b_{i2}=b_{i1}b_{i2}^{\dagger}=b_{i2}b_{i1}^{\dagger}=0$
and
$b_{i2}b_{i2}^{\dagger}=b_{i1}b_{i1}^{\dagger}=1-b_{i1}^{\dagger}b_{i1}-b_{i2}^{\dagger}b_{i2}$.
The spin operators can be written in terms of these bosons as

$\displaystyle
S_{i}^{+}=\sqrt{2}[b_{i1}^{\dagger}+b_{i2}],~{}S_{i}^{-}=\sqrt{2}[b_{i1}+b_{i2}^{\dagger}],~{}S_{i}^{z}=n^{\prime}_{i1}-n^{\prime}_{i2},$
(99)

where, $n^{\prime}_{i1(2)}=b_{i1(2)}^{\dagger}b_{i1(2)}$. It is useful to note
that the above two bosons are related to the three bosons $a_{x},a_{y},a_{z}$
of the SU(3) flavour wave theory Läuchli  _et al._ (2006) as :
$b_{1}=(-a_{x}+ia_{y})/\sqrt{2}$ and $b_{2}=(a_{x}+ia_{y})/\sqrt{2}$ while
$a_{z}$ is condensed. Using Eq. 99 in Eq. 2 and using the constraints, we get

$\displaystyle H_{sp,FN}$ $\displaystyle=$ $\displaystyle
6K_{0}N-6K_{0}\sum_{i,\mu^{\prime}}n^{\prime}_{i\mu^{\prime}}+2J_{0}\sum_{\langle
ij\rangle}(b_{i1}^{\dagger}b_{j1}+b_{i2}^{\dagger}b_{j2})$ (100)
$\displaystyle+2(J_{0}-K_{0})\sum_{\langle
ij\rangle}(b_{i1}^{\dagger}b_{j2}^{\dagger}+b_{i1}b_{j2})$

where $\mu^{\prime}=1,2$. The quadratic spin-nematic wave Hamiltonian above
can be diagonalized via Fourier transform followed by Bogoliubov
transformation,

$\displaystyle\left(\begin{matrix}b_{\mathbf{k},1}\\\
b_{-\mathbf{k},2}^{\dagger}\end{matrix}\right)$ $\displaystyle=$
$\displaystyle\mathcal{R}_{\mathbf{k}}\psi_{\mathbf{k}}~{},$ (101)

with

$\displaystyle\mathcal{R}_{\mathbf{k}}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}u_{\mathbf{k}}&v_{\mathbf{k}}\\\
v_{\mathbf{k}}&u_{\mathbf{k}}\end{matrix}\right)\text{ and
}\psi_{\mathbf{k}}=\left(\begin{matrix}d_{\mathbf{k},1}\\\
d_{-\mathbf{k},2}^{\dagger}\end{matrix}\right)$ (102)

such that $u_{\mathbf{k}}=\text{cosh}\bar{\theta}_{\mathbf{k}}$,
$v_{\mathbf{k}}=\text{sinh}\bar{\theta}_{\mathbf{k}}$ and
$\text{tanh}2\bar{\theta}_{\mathbf{k}}=\frac{(K_{0}-J_{0})\gamma_{\mathbf{k}}}{J_{0}\gamma_{\mathbf{k}}-K_{0}}$.
The diagonalized Hamiltonian and dispersion are given in Eq. 17 and Eq. 18
respectively.

###  C.2 The three sublattice nematic

Spin-nematic wave theory for the three sublattice nematic,Tsunetsugu and
Arikawa (2006) as mentioned in the main text, requires two bosons at every
sublattice, $\tilde{\alpha}_{i}$ and $\tilde{\beta}_{i}$, that account for the
deviations about the respective nematic director. In terms of these bosons,
the spin operators at sublattice $3$ (say) (see Fig. 1) are given by

$\displaystyle S^{x}_{\mathbf{\widetilde{R}},3}$ $\displaystyle=$
$\displaystyle\widetilde{\alpha}^{\dagger}_{\mathbf{\widetilde{R}},3}+\widetilde{\alpha}_{\mathbf{\widetilde{R}},3},~{}~{}~{}S^{y}_{\mathbf{\widetilde{R}},3}=\widetilde{\beta}^{\dagger}_{\mathbf{\widetilde{R}},3}+\widetilde{\beta}_{\mathbf{\widetilde{R}},3}$
$\displaystyle S^{z}_{\mathbf{\widetilde{R}},3}$ $\displaystyle=$
$\displaystyle-i(\widetilde{\alpha}^{\dagger}_{\mathbf{\widetilde{R}},3}\widetilde{\beta}_{\mathbf{\widetilde{R}},3}-\widetilde{\beta}^{\dagger}_{\mathbf{\widetilde{R}},3}\widetilde{\alpha}_{\mathbf{\widetilde{R}},3})~{},$
(103)

with constraints
$(\widetilde{\alpha}_{\mathbf{\widetilde{R}},3})^{2}=(\widetilde{\beta}_{\mathbf{\widetilde{R}},3})^{2}=\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}\widetilde{\beta}_{\mathbf{\widetilde{R}},3}=\widetilde{\beta}_{\mathbf{\widetilde{R}},3}\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}=\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}\widetilde{\beta}_{\mathbf{\widetilde{R}},3}^{\dagger}=\widetilde{\beta}_{\mathbf{\widetilde{R}},3}\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}^{\dagger}=0$
and
$\widetilde{\beta}_{\mathbf{\widetilde{R}},3}\widetilde{\beta}_{\mathbf{\widetilde{R}},3}^{\dagger}=\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}^{\dagger}=1-\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}^{\dagger}\widetilde{\alpha}_{\mathbf{\widetilde{R}},3}-\widetilde{\beta}_{\mathbf{\widetilde{R}},3}^{\dagger}\widetilde{\beta}_{\mathbf{\widetilde{R}},3}$
for faithful representation of the Hilbert space. Likewise, the spin operators
can be written for the other two sublattices.Tsunetsugu and Arikawa (2006) Eq.
2, up to harmonic orders in the bosons is,

$H_{sp}=3K_{0}N+H(\widetilde{\beta}_{1},\widetilde{\alpha}_{2})+H(\widetilde{\beta}_{2},\widetilde{\alpha}_{3})+H(\widetilde{\beta}_{3},\widetilde{\alpha}_{1})~{},$
(104)

with the following form of the Hamiltonian for each pair featuring above:

$\displaystyle
H(\widetilde{\beta}_{\widetilde{\lambda}},\widetilde{\alpha}_{\widetilde{\rho}})$
$\displaystyle=$ $\displaystyle
3\sum_{\mathbf{k}}\big{[}\widetilde{\gamma}_{\mathbf{k}}(J_{0}\widetilde{\beta}^{\dagger}_{\mathbf{k},\widetilde{\lambda}}\widetilde{\alpha}^{\dagger}_{-\mathbf{k},\widetilde{\rho}}+(J_{0}-K_{0})\widetilde{\beta}^{\dagger}_{\mathbf{k},\widetilde{\lambda}}\widetilde{\alpha}_{\mathbf{k},\widetilde{\rho}})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\text{h.c}\big{]}+3K_{0}\sum_{\mathbf{k}}(\bar{n}_{\widetilde{\beta},\mathbf{k},\widetilde{\lambda}}+\bar{n}_{\widetilde{\alpha},\mathbf{k},\widetilde{\rho}})$

where $\widetilde{\gamma}_{\mathbf{k}}$ is as defined below Eq. 27 and
$\bar{n}_{\widetilde{\beta},\mathbf{k},\widetilde{\lambda}}=\widetilde{\beta}^{\dagger}_{\mathbf{k},\widetilde{\lambda}}\widetilde{\beta}_{\mathbf{k},\widetilde{\lambda}}$,
$\bar{n}_{\widetilde{\alpha},\mathbf{k},\widetilde{\rho}}=\widetilde{\alpha}^{\dagger}_{\mathbf{k},\widetilde{\rho}}\widetilde{\alpha}_{\mathbf{k},\widetilde{\rho}}$.
The above Hamiltonian can be diagonalized using Bogoliubov
transformation,Kondo  _et al._ (2020)

$\left(\begin{matrix}\widetilde{\alpha}_{\mathbf{k},\widetilde{\rho}}\\\
\widetilde{\beta}_{\mathbf{k},\widetilde{\lambda}}\\\
\widetilde{\alpha}_{-\mathbf{k},\widetilde{\rho}}^{\dagger}\\\
\widetilde{\beta}_{-\mathbf{k},\widetilde{\lambda}}^{\dagger}\end{matrix}\right)=T_{\mathbf{k}}\left(\begin{matrix}\widetilde{d}_{+,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}\\\
\widetilde{d}_{-,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}\\\
\widetilde{d}_{+,-\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}^{~{}\dagger}\\\
\widetilde{d}_{-,-\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}^{~{}\dagger}\end{matrix}\right)\text{
with
}T_{\mathbf{k}}=\left(\begin{matrix}\widetilde{T}_{1,{\mathbf{k}}}&\widetilde{T}_{2,{\mathbf{k}}}\\\
\widetilde{T}_{2,{\mathbf{k}}}&\widetilde{T}_{1,{\mathbf{k}}}\end{matrix}\right)$
(106)

where

$\displaystyle\widetilde{T}_{1,{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\frac{\widetilde{t}_{+,\mathbf{k}}}{\widetilde{\gamma}_{\mathbf{k}}\widetilde{N}_{+,\mathbf{k}}}&\frac{\widetilde{t}_{-,\mathbf{k}}}{\widetilde{\gamma}_{\mathbf{k}}\widetilde{N}_{-,\mathbf{k}}}\\\
-\frac{\widetilde{t}_{+,\mathbf{k}}}{|\widetilde{\gamma}_{\mathbf{k}}|\widetilde{N}_{+,\mathbf{k}}}&\frac{\widetilde{t}_{-,\mathbf{k}}}{|\widetilde{\gamma}_{\mathbf{k}}|\widetilde{N}_{-,\mathbf{k}}}\end{matrix}\right)$
$\displaystyle\widetilde{T}_{2,{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}-\frac{\widetilde{\gamma}^{*}_{\mathbf{k}}J_{0}}{|\widetilde{\gamma}_{\mathbf{k}}|\widetilde{N}_{+,\mathbf{k}}}&\frac{\widetilde{\gamma}^{*}_{\mathbf{k}}J_{0}}{|\widetilde{\gamma}_{\mathbf{k}}|\widetilde{N}_{-,\mathbf{k}}}\\\
\frac{J_{0}}{\widetilde{N}_{+,\mathbf{k}}}&\frac{J_{0}}{\widetilde{N}_{-,\mathbf{k}}}\end{matrix}\right)$
(107)

with
$\widetilde{t}_{\pm,\mathbf{k}}=\pm(J_{0}-K_{0})|\widetilde{\gamma}_{\mathbf{k}}|-\left(K_{0}+\frac{\widetilde{\omega}^{s}_{\pm,\mathbf{k}}}{3}\right)$,
$\widetilde{N}_{\pm,\mathbf{k}}=\frac{2K_{0}}{|\widetilde{\gamma}_{\mathbf{k}}|}\sqrt{1+{\omega}^{\prime}_{\pm,\mathbf{k}}\pm(1-J^{\prime})|\widetilde{\gamma}_{\mathbf{k}}|(2+{\omega}^{\prime}_{\pm,\mathbf{k}})+(1-2J^{\prime})|\widetilde{\gamma}_{\mathbf{k}}|^{2}}$
where
${\omega}^{\prime}_{\pm,\mathbf{k}}=\frac{\widetilde{\omega}^{s}_{\pm,\mathbf{k}}}{3K_{0}}$,
$J^{\prime}=\frac{J_{0}}{K_{0}}$ and the Hamiltonian and dispersions are given
in Eq. 26 and Eq. 27.

##  Appendix D  Temperature dependence of fractional change in sound speed in
the spin nematic phases

###  D.1  Ferronematic phase

The coupling Hamiltonian for the ferro spin-nematic waves and the phonons are
given by Eq. 20 where the scattering vertices,
$\mathcal{M}_{\mathbf{k},\mathbf{q_{2}}}^{(1)}$,
$\mathcal{M}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}^{(2)}$ and
$\mathcal{M}_{\mathbf{k}}^{(2),0}$ are given by

$\displaystyle\mathcal{M}_{\mathbf{k},\mathbf{q_{2}}}^{(1)}=\frac{1}{2}\sum_{\bm{\delta}}\frac{(1-e^{i\mathbf{k}\cdot{\bm{\delta}}})}{\sqrt{2MN\omega_{0,\mathbf{k}}^{ph}}}~{}\overline{\mathcal{M}}_{\mathbf{k},\mathbf{q_{2}}}^{(1)}(\bm{\delta})$
(108)

with

$\displaystyle\overline{\mathcal{M}}_{\mathbf{k},\mathbf{q_{2}}}^{(1)}(\bm{\delta})=\mathcal{R}_{\mathbf{k}+\mathbf{q_{2}}}\left[(e^{i\mathbf{q}_{2}\cdot{\bm{\delta}}}+e^{-i(\mathbf{k}+\mathbf{q}_{2})\cdot{\bm{\delta}}})\left[-\frac{\partial
J_{\bm{\delta}}}{\partial\bm{\delta}}\cdot\mathbf{e}_{\mathbf{k}}\left[\begin{matrix}1&1\\\
1&1\end{matrix}\right]+\frac{\partial
K_{\bm{\delta}}}{\partial\bm{\delta}}\cdot\mathbf{e}_{\mathbf{k}}\left[\begin{matrix}0&1\\\
1&0\end{matrix}\right]\right]+(1+e^{-i\mathbf{k}\cdot{\bm{\delta}}})\frac{\partial
K_{\bm{\delta}}}{\partial\bm{\delta}}\cdot\mathbf{e}_{\mathbf{k}}\left[\begin{matrix}1&0\\\
0&1\end{matrix}\right]\right]\mathcal{R}_{\mathbf{q_{2}}}$ (109)

and

$\displaystyle\mathcal{M}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}^{(2)}$
$\displaystyle=$
$\displaystyle\sum_{\bm{\delta}}\frac{(1-e^{i\mathbf{k}\cdot{\bm{\delta}}})(1-e^{-i\mathbf{k}^{\prime}\cdot{\bm{\delta}}})}{8MN\sqrt{\omega_{0,\mathbf{k}}^{ph}\omega_{0,-\mathbf{k}^{\prime}}^{ph}}}\mathcal{R}_{\mathbf{k}-\mathbf{k}^{\prime}+\mathbf{q_{2}}}\bigg{[}\left(\mathbf{e}_{-\mathbf{k}^{\prime}}\cdot\frac{\partial^{2}J_{\bm{\delta}}}{\partial{\bm{\delta}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}\right)(e^{i\mathbf{q_{2}}\cdot\bm{\delta}}+e^{-i(\mathbf{k}-\mathbf{k}^{\prime}+\mathbf{q_{2}})\cdot\bm{\delta}})\left(\begin{matrix}1&1\\\
1&1\end{matrix}\right)$ (110)
$\displaystyle~{}~{}~{}~{}~{}-\left(\mathbf{e}_{-\mathbf{k}^{\prime}}\cdot\frac{\partial^{2}K_{\bm{\delta}}}{\partial{\bm{\delta}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}\right)\bigg{(}(e^{i\mathbf{q}_{2}\cdot{\bm{\delta}}}+e^{-i(\mathbf{k}-\mathbf{k}^{\prime}+\mathbf{q}_{2})\cdot{\bm{\delta}}})\left(\begin{matrix}0&1\\\
1&0\end{matrix}\right)+(1+e^{-i(\mathbf{k}-\mathbf{k}^{\prime})\cdot{\bm{\delta}}})\left(\begin{matrix}1&0\\\
0&1\end{matrix}\right)\bigg{)}\bigg{]}\mathcal{R}_{\mathbf{q}_{2}}$
$\displaystyle\mathcal{M}_{\mathbf{k}}^{(2),0}$ $\displaystyle=$
$\displaystyle\sum_{\bm{\delta}}\frac{1}{2M\omega_{0,\mathbf{k}}^{ph}}(1-e^{i\mathbf{k}\cdot{\bm{\delta}}})(1-e^{-i\mathbf{k}\cdot{\bm{\delta}}})\left(\mathbf{e}_{-\mathbf{k}}\cdot\frac{\partial^{2}K_{\bm{\delta}}}{\partial{\bm{\delta}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}\right).$
(111)

The matrix $\mathcal{R}_{\mathbf{k}}$ is defined in Eq. 102. $J_{\bm{\delta}}$
and $K_{\bm{\delta}}$ denote the coupling constants depend on the relative
separation ($\bm{\delta}$) between the two nearest neighbour spins. The
resultant temperature dependence of the fractional change in sound speed is
given by Eq. 24 with

$\displaystyle\bar{c}_{1}=$
$\displaystyle\sum_{\bm{\delta}}\frac{1}{Mv_{0}^{2}}(\hat{\mathbf{q}}\cdot\bm{\delta})^{2}\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}K_{\bm{\delta}}}{\partial{\bm{\delta}^{2}}}\cdot\mathbf{e}_{\mathbf{q}}\right)-\frac{V}{32\pi^{2}MNv_{0}^{2}}\int_{\text{BZ}}dq_{2}d\theta_{\hat{\mathbf{q}}_{2}}q_{2}\sum_{\bm{\delta},\bm{\delta}^{\prime}}\frac{(\hat{\mathbf{q}}\cdot\bm{\delta})(\hat{\mathbf{q}}\cdot\bm{\delta}^{\prime})}{\omega^{s}_{\mathbf{q}_{2}}}\overline{\mathcal{M}}_{0,\mathbf{q}_{2}}^{(1),12}(\bm{\delta})\overline{\mathcal{M}}_{0,\mathbf{q}_{2}}^{(1),21}(\bm{\delta}^{\prime})$
$\displaystyle+\frac{V}{8\pi^{2}MNv_{0}^{2}}\int_{\text{BZ}}dq_{2}d\theta_{\hat{\mathbf{q}}_{2}}q_{2}\sum_{\bm{\delta}}(\hat{\mathbf{q}}\cdot\bm{\delta})^{2}\bigg{[}\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}J_{\bm{\delta}}}{\partial\bm{\delta}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)\text{cos}(\mathbf{q}_{2}\cdot\bm{\delta})(\text{cosh}2\bar{\theta}_{\mathbf{q}_{2}}+\text{sinh}2\bar{\theta}_{\mathbf{q}_{2}})$
$\displaystyle-\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}K_{\bm{\delta}}}{\partial\bm{\delta}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)(\text{cosh}2\bar{\theta}_{\mathbf{q}_{2}}+\text{cos}(\mathbf{q}_{2}\cdot\bm{\delta})\text{sinh}2\bar{\theta}_{\mathbf{q}_{2}})\bigg{]}$
$\displaystyle\bar{c}_{2}=$
$\displaystyle\frac{Vk_{B}^{3}}{32\pi^{2}MNv_{0}^{2}\bar{c}_{s}^{4}}\int_{0}^{\infty}dx\int_{0}^{2\pi}d\theta_{\hat{\mathbf{q}}_{2}}\bigg{[}\frac{2x^{3}e^{x}\bar{\mathcal{A}}_{+}(\hat{\mathbf{q}},\hat{\mathbf{q}_{2}})\bar{c}_{s}^{2}\cos^{2}{\theta_{\mathbf{q}\mathbf{q}_{2}}}}{(e^{x}-1)^{2}(v_{0}^{2}-\bar{c}_{s}^{2}\cos^{2}{\theta_{\mathbf{q}\mathbf{q}_{2}}})}+\frac{2x^{2}}{(e^{x}-1)}\bigg{(}-\bar{\mathcal{A}}_{-}(\hat{\mathbf{q}},\hat{\mathbf{q}}_{2})+2\bar{c}_{s}\bar{\mathcal{C}}(\hat{\mathbf{q}},\hat{\mathbf{q}}_{2})\bigg{)}\bigg{]}$

where $q_{2}=|\mathbf{q}_{2}|$, $V$ is the total volume of the lattice,
$\theta_{\mathbf{q}\mathbf{q}_{2}}$ is the angle between $\hat{\mathbf{q}}$
and $\hat{\mathbf{q}}_{2}$, the integration limits in $\bar{c}_{2}$ can be
extended upto infinity due to exponential damping factor and

$\displaystyle\bar{\mathcal{A}}_{\pm}(\hat{\mathbf{q}},\hat{\mathbf{q}_{2}})=$
$\displaystyle\sum_{\bm{\delta},\bm{\delta}^{\prime}}(\hat{\mathbf{q}}\cdot\bm{\delta})(\hat{\mathbf{q}}\cdot\bm{\delta}^{\prime})\bigg{[}\left(\frac{\partial
J_{\bm{\delta}}}{\partial\bm{\delta}}\cdot\mathbf{e}_{-\mathbf{q}}\right)\frac{3K_{0}}{\bar{c}_{s}}-\left(\frac{\partial
K_{\bm{\delta}}}{\partial\bm{\delta}}\cdot\mathbf{e}_{-\mathbf{q}}\right)\left(\frac{3K_{0}}{\bar{c}_{s}}\pm\frac{(\hat{\mathbf{q}}_{2}\cdot\bm{\delta})^{2}}{\bar{c}_{s}}6(K_{0}-J_{0})\right)\bigg{]}$
$\displaystyle~{}~{}~{}~{}\times\bigg{[}\left(\frac{\partial
J_{\bm{\delta}^{\prime}}}{\partial\bm{\delta}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}\right)\frac{3K_{0}}{\bar{c}_{s}}-\left(\frac{\partial
K_{\bm{\delta}^{\prime}}}{\partial\bm{\delta}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}\right)\left(\frac{3K_{0}}{\bar{c}_{s}}\pm\frac{(\hat{\mathbf{q}}_{2}\cdot\bm{\delta^{\prime}})^{2}}{\bar{c}_{s}}6(K_{0}-J_{0})\right)\bigg{]}$
$\displaystyle\bar{\mathcal{C}}(\hat{\mathbf{q}},\hat{\mathbf{q}_{2}})=$
$\displaystyle\sum_{\bm{\delta}}(\hat{\mathbf{q}}\cdot\bm{\delta})^{2}\bigg{[}-\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}J_{\bm{\delta}}}{\partial\bm{\delta}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)\frac{3K_{0}}{\bar{c}_{s}}+\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}K_{\bm{\delta}}}{\partial\bm{\delta}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)\left(\frac{3K_{0}}{\bar{c}_{s}}+\frac{(\hat{\mathbf{q}}_{2}\cdot\bm{\delta})^{2}}{\bar{c}_{s}}6(K_{0}-J_{0})\right)\bigg{]}$

To obtain the sound attenuation $\propto\Sigma(\mathbf{q},\omega)/v_{0}$, we
note that at the leading order only $\Sigma_{1}$ has an imaginary part. From
Eq. 108, it is fairly straight forward to show that it is proportional to
$\omega_{0,\mathbf{q}}^{ph}$ and hence disappears for the acoustic phonons in
the limit of ${\bf q}\rightarrow 0$.Cottam (1974)

###  D.2  Three sublattice nematic phase

For the three sublattice nematic, the spin-phonon coupling Hamiltonian in Eq.
25 is given by

$\displaystyle H_{sp-ph,3SN}=H_{1,3SN}+H_{2,3SN}$ (112)

where the two terms are similar to Eqs. 20 are given by

$\displaystyle H_{1,3SN}=$
$\displaystyle\sum_{\mathbf{k},\mathbf{q_{2}},\widetilde{\lambda}\widetilde{\rho}}\widetilde{\psi}^{\dagger}_{\mathbf{k}+\mathbf{q_{2}},\widetilde{\lambda}\widetilde{\rho}}\widetilde{\mathcal{M}}_{\mathbf{k}+\mathbf{q_{2}},\mathbf{q_{2}}}^{(1)}\widetilde{\psi}_{\mathbf{q_{2}},\widetilde{\lambda}\widetilde{\rho}}A_{\mathbf{k}}$
(113)

$\displaystyle H_{2,3SN}=$
$\displaystyle\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}},\widetilde{\lambda}\widetilde{\rho}}\widetilde{\psi}^{\dagger}_{\mathbf{k}-\mathbf{k}^{\prime}+\mathbf{q_{2}},\widetilde{\lambda}\widetilde{\rho}}\widetilde{\mathcal{M}}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}^{(2)}\widetilde{\psi}_{\mathbf{q_{2}},\widetilde{\lambda}\widetilde{\rho}}A_{\mathbf{k}}A_{-\mathbf{k}^{\prime}}$
$\displaystyle+\sum_{\mathbf{k}}\widetilde{\mathcal{M}}_{\mathbf{k}}^{(2),0}A_{\mathbf{k}}A_{-\mathbf{k}}~{},$
(114)

where
$\widetilde{\psi}_{\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}=\bigg{(}\widetilde{d}_{+,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}},\text{
}\widetilde{d}_{-,\mathbf{k},\widetilde{\lambda}\widetilde{\rho}},\text{
}\widetilde{d}^{~{}\dagger}_{+,-\mathbf{k},\widetilde{\lambda}\widetilde{\rho}},\text{
}\widetilde{d}^{~{}\dagger}_{-,-\mathbf{k},\widetilde{\lambda}\widetilde{\rho}}\bigg{)}^{T}$
and

$\displaystyle\widetilde{\mathcal{M}}_{\mathbf{k}+\mathbf{q_{2}},\mathbf{q_{2}}}^{(1)}=$
$\displaystyle\sum_{\bm{\widetilde{\delta}}}\sqrt{\frac{1}{2MN\omega_{0,\mathbf{k}}^{ph}}}(1-e^{i\mathbf{k}\cdot{\bm{\widetilde{\delta}}}})\overline{\widetilde{\mathcal{M}}}_{\mathbf{k}+\mathbf{q_{2}},\mathbf{q_{2}}}^{(1)}(\bm{\widetilde{\delta}})$
(115)

where

$\displaystyle\overline{\widetilde{\mathcal{M}}}_{\mathbf{k}+\mathbf{q_{2}},\mathbf{q_{2}}}^{(1)}(\bm{\widetilde{\delta}})=$
$\displaystyle
T^{\dagger}_{\mathbf{k}+\mathbf{q_{2}}}\Bigg{[}\left(\frac{\partial
J_{\bm{\widetilde{\delta}}}}{\partial(-\bm{\widetilde{\delta}})}\cdot\mathbf{e}_{\mathbf{k}}\right)e^{i\mathbf{q}_{2}\cdot{\bm{\widetilde{\delta}}}}\left(\begin{matrix}0&0&0&0\\\
1&0&1&0\\\ 0&0&0&0\\\ 1&0&1&0\end{matrix}\right)+\left(\frac{\partial
K_{\bm{\widetilde{\delta}}}}{\partial(-\bm{\widetilde{\delta}})}\cdot\mathbf{e}_{\mathbf{k}}\right)\left(\begin{matrix}e^{-i\mathbf{k}\cdot{\bm{\widetilde{\delta}}}}&-e^{-i(\mathbf{k}+\mathbf{q}_{2})\cdot{\bm{\widetilde{\delta}}}}&0&0\\\
-e^{i\mathbf{q}_{2}\cdot{\bm{\widetilde{\delta}}}}&1&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{matrix}\right)\Bigg{]}T_{\mathbf{q}_{2}}$ (116)

and

$\displaystyle\widetilde{\mathcal{M}}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}^{(2)}=$
$\displaystyle\sum_{\bm{\widetilde{\delta}}}\frac{1}{4MN\sqrt{\omega_{0,\mathbf{k}}^{ph}\omega_{0,-\mathbf{k}^{\prime}}^{ph}}}(1-e^{i\mathbf{k}\cdot{\bm{\widetilde{\delta}}}})(1-e^{-i\mathbf{k}^{\prime}\cdot{\bm{\widetilde{\delta}}}})\overline{\widetilde{\mathcal{M}}}^{(2)}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}(\bm{\widetilde{\delta}})$
(117)

where

$\displaystyle\overline{\widetilde{\mathcal{M}}}^{(2)}_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q_{2}}}(\bm{\widetilde{\delta}})=$
$\displaystyle
T^{\dagger}_{\mathbf{k}-\mathbf{k}^{\prime}+\mathbf{q}_{2}}\Bigg{[}\mathbf{e}_{-\mathbf{k}^{\prime}}\cdot\frac{\partial^{2}J_{\bm{\widetilde{\delta}}}}{\partial{\bm{\widetilde{\delta}}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}e^{i\mathbf{q}_{2}\cdot{\bm{\widetilde{\delta}}}}\left[\begin{matrix}0&0&0&0\\\
1&0&1&0\\\ 0&0&0&0\\\
1&0&1&0\end{matrix}\right]+\mathbf{e}_{-\mathbf{k}^{\prime}}\cdot\frac{\partial^{2}K_{\bm{\widetilde{\delta}}}}{\partial{\bm{\widetilde{\delta}}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}\left[\begin{matrix}e^{-i(\mathbf{k}-\mathbf{k}^{\prime})\cdot{\bm{\widetilde{\delta}}}}&-e^{-i(\mathbf{k}-\mathbf{k}^{\prime}+\mathbf{q}_{2})\cdot{\bm{\widetilde{\delta}}}}&0&0\\\
-e^{i\mathbf{q}_{2}\cdot{\bm{\widetilde{\delta}}}}&1&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{matrix}\right]\Bigg{]}T_{\mathbf{q}_{2}}$ (118)

$\displaystyle\widetilde{\mathcal{M}}_{\mathbf{k}}^{(2),0}=$
$\displaystyle\sum_{\bm{\widetilde{\delta}}}\frac{1}{4M\omega_{0,\mathbf{k}}^{ph}}(1-e^{i\mathbf{k}\cdot{\bm{\widetilde{\delta}}}})(1-e^{-i\mathbf{k}\cdot{\bm{\widetilde{\delta}}}})\left(\mathbf{e}_{-\mathbf{k}}\cdot\frac{\partial^{2}K_{\bm{\widetilde{\delta}}}}{\partial{\bm{\widetilde{\delta}}^{2}}}\cdot\mathbf{e}_{\mathbf{k}}\right)$
(119)

Here, the matrix $T_{\mathbf{k}}$ is as defined in Eq. 106 and
$\bm{\widetilde{\delta}}$ are the displacements mentioned below Eq. 27.

The phonon self energy is then calculated perturbatively (in spin-phonon
interaction) and the leading contributions are given by:

$\displaystyle\Sigma(\mathbf{q},i\Omega_{n})$
$\displaystyle=\sum_{\mathbf{k}}\widetilde{\Sigma}_{\mathbf{k}}(\mathbf{q},i\Omega_{n})$
$\displaystyle=\sum_{\mathbf{k}}[\widetilde{\Sigma}_{1\mathbf{k}}(\mathbf{q},i\Omega_{n})+\widetilde{\Sigma}_{2\mathbf{k}}(\mathbf{q},i\Omega_{n})]~{},$
(120)

where

$\displaystyle\widetilde{\Sigma}_{1\mathbf{k}}(\mathbf{q},i\Omega_{n})\approx\frac{3}{2}\sum_{\widetilde{\mu},\widetilde{\nu}\in\\{+,-\\}}$
$\displaystyle\bigg{[}-\frac{1}{\beta}\mathcal{N}_{2,\mathbf{q},\mathbf{k}}^{\widetilde{\mu}\widetilde{\nu}}\sum_{\omega_{1}}\widetilde{G}^{m,\widetilde{\mu}}_{0}(\mathbf{q}+\mathbf{k},i\Omega_{n}+i\omega_{1})\widetilde{G}^{m,\widetilde{\nu}}_{0}(\mathbf{k},i\omega_{1})$
$\displaystyle-\frac{1}{\beta}\mathcal{N}_{3,\mathbf{q},\mathbf{k}}^{\widetilde{\mu}\widetilde{\nu}}\sum_{\omega_{1}}\widetilde{G}^{m,\widetilde{\mu}}_{0}(-\mathbf{q}+\mathbf{k},-i\Omega_{n}-i\omega_{1})\widetilde{G}^{m,\widetilde{\nu}}_{0}(-\mathbf{k},i\omega_{1})$
$\displaystyle-\frac{1}{\beta}\mathcal{N}_{1,\mathbf{q},\mathbf{k}}^{\widetilde{\mu}\widetilde{\nu}}\sum_{\omega_{1}}\widetilde{G}^{m,\widetilde{\mu}}_{0}(\mathbf{q}+\mathbf{k},i\Omega_{n}-i\omega_{1})\widetilde{G}^{m,\widetilde{\nu}}_{0}(-\mathbf{k},i\omega_{1})\bigg{]}$
(121)
$\displaystyle\widetilde{\Sigma}_{2\mathbf{k}}(\mathbf{q},i\Omega_{n})\approx
2\widetilde{\mathcal{M}}_{\mathbf{q}}^{(2),0}+6\sum_{\tilde{i}=\\{1\text{ to
}4\\}}\widetilde{\mathcal{M}}_{\mathbf{q},\mathbf{q},\mathbf{k}}^{(2),\tilde{i}\tilde{i}}\langle\widetilde{\psi}^{\dagger\tilde{i}}_{\mathbf{k}}\widetilde{\psi}^{~{}\tilde{i}}_{\mathbf{k}}\rangle~{},$
(122)

where the $\widetilde{\lambda}\widetilde{\rho}$ indices are suppressed since
each pair gives the same contribution and has been accounted by appropriate
multiplicative factors in the expressions above,
$\widetilde{G}^{m,+}_{0}(\mathbf{q},i\Omega_{n})$ and
$\widetilde{G}^{m,-}_{0}(\mathbf{q},i\Omega_{n})$ are the bare propagators for
the $\widetilde{d}_{+,\mathbf{q}}$ and $\widetilde{d}_{-,\mathbf{q}}$ bosons:
$\widetilde{G}^{m,+}_{0}(\mathbf{q},i\Omega_{n})=\frac{1}{i\Omega_{n}-\widetilde{\omega}^{s}_{+,\mathbf{q}}}$
and
$\widetilde{G}^{m,-}_{0}(\mathbf{q},i\Omega_{n})=\frac{1}{i\Omega_{n}-\widetilde{\omega}^{s}_{-,\mathbf{q}}}$
and $\mathcal{N}$ in $\Sigma_{1}$ stands for the appropriate matrix elements
of
$\mathcal{O}((\widetilde{\mathcal{M}}_{\mathbf{k},\mathbf{q_{2}}}^{(1)})^{2})$.
The sum over Matsubara frequencies can be performed using the techniques in
Ref. Bruus and Flensberg, 2004. Converting the momentum sums into integrals
one obtains Eq. 28 in the main text.

For $J_{0}/K_{0}\ll 1$ and temperatures much less than the gap scale of
$\widetilde{\omega}^{s}_{+,\mathbf{q}}$, the temperature dependent
contribution to Eq. 121 due to both the internal lines in Fig. 3 belonging to
the $\widetilde{\omega}^{s}_{+,\mathbf{q}}$ branch, the contribution is
exponentially damped due to the finite gap. When one branch is
$\widetilde{\omega}^{s}_{+,\mathbf{q}}$ and the other is
$\widetilde{\omega}^{s}_{-,\mathbf{q}}$, the temperature dependence obtained
is similar to that obtained for both the branches being
$\widetilde{\omega}^{s}_{-,\mathbf{q}}$ albeit with different prefactors.
Likewise, the contribution to the second term of Eq. 122 due to the
$\widetilde{\omega}^{s}_{+,\mathbf{q}}$ would be exponentially damped.

We now list the detailed expressions obtained by considering only the
$\widetilde{\omega}^{s}_{-,\mathbf{q}}$ branch for obtaining the self energy.
Now, to scale out the temperature dependence due to the Bose distribution
function, we perform the standard change of variables
$x=\beta\widetilde{\omega}^{s}_{-,\mathbf{k}}$ (where $\beta=1/(k_{B}T)$) and
obtain the temperature dependence. This changes the upper limit of the first
term in Eq. 28 and the lower limit of the second term in Eq. 28 to $\beta
c_{s}k^{*}=\beta\widetilde{c_{s}}(k^{*})^{2}=6\beta J_{0}$ for integration
over the variable $x$. In the limit where $\beta J_{0}\gg 1$, only the
contribution from the first integral in Eq. 28 with the linear dispersion
would be important, giving the following expressions for the pre-factors in
Eq. 32,

$\displaystyle\widetilde{c}_{1}=$
$\displaystyle\sum_{\bm{\widetilde{\delta}}}\frac{1}{2Mv_{0}^{2}}(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}})^{2}\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}K_{\bm{\widetilde{\delta}}}}{\partial{\bm{\widetilde{\delta}}^{2}}}.\mathbf{e}_{\mathbf{q}}\right)-\frac{3V}{16\pi^{2}}\lim_{\mathbf{q}\xrightarrow{}0}\int_{\text{BZ}}dkd\theta_{\hat{\mathbf{k}}}k\frac{(\mathcal{N}_{1,\mathbf{q},\mathbf{k}}^{--}+\mathcal{N}_{3,\mathbf{q},\mathbf{k}}^{--})}{\widetilde{\omega}^{s}_{-,\mathbf{k}}\omega_{0,\mathbf{q}}^{ph}}$
$\displaystyle+\frac{3V}{8\pi^{2}MNv_{0}^{2}}\int_{\text{BZ}}dkd\theta_{\hat{\mathbf{k}}}k\sum_{\bm{\widetilde{\delta}}}(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}})^{2}\overline{\widetilde{\mathcal{M}}}^{(2),44}_{0,0,\mathbf{k}}(\bm{\widetilde{\delta}})$
(123)

$\displaystyle\widetilde{c}_{2}=$
$\displaystyle\frac{3Vk_{B}^{3}}{4\pi^{2}MNv_{0}^{2}c_{s}^{4}}\int_{0}^{\infty}dx\int_{0}^{2\pi}d\theta_{\hat{\mathbf{k}}}\bigg{[}\frac{x^{3}e^{x}c_{s}\cos{\theta_{\mathbf{q}\mathbf{k}}}\mathcal{A}_{+}(\hat{\mathbf{q}},\hat{\mathbf{k}})K_{0}}{(e^{x}-1)^{2}(v_{0}-c_{s}\cos{\theta_{\mathbf{q}\mathbf{k}}})J_{0}}+\frac{x^{2}}{(e^{x}-1)}\big{(}-\mathcal{A}_{-}(\hat{\mathbf{q}},\hat{\mathbf{k}})\frac{K_{0}}{J_{0}}+c_{s}\sqrt{\frac{K_{0}}{2J_{0}}}\mathcal{C}(\hat{\mathbf{q}},\hat{\mathbf{k}})\big{)}\bigg{]}$
(124)

where $k=|\mathbf{k}|$, $\theta_{\mathbf{q}\mathbf{k}}$ is the angle between
$\hat{\mathbf{q}}$ and $\hat{\mathbf{k}}$, the integration limits in
$\widetilde{c}_{2}$ have been extended to infinity in the $\beta J_{0}\gg 1$
limit and

$\displaystyle\mathcal{A}_{\pm}(\hat{\mathbf{q}},\hat{\mathbf{k}})=$
$\displaystyle\sum_{\bm{\widetilde{\delta}},\bm{\widetilde{\delta}}^{\prime}}(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}})(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}}^{\prime})\bigg{[}\frac{1}{4}\frac{\partial
J_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}}\cdot\mathbf{e}_{-\mathbf{q}}\pm\frac{\partial
K_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}}\cdot\mathbf{e}_{-\mathbf{q}}\left(\frac{J_{0}}{2K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}})^{2}\bigg{]}\bigg{[}\frac{1}{4}\frac{\partial
J_{\bm{\widetilde{\delta}}^{\prime}}}{\partial\bm{\widetilde{\delta}}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}\pm\frac{\partial
K_{\bm{\widetilde{\delta}}^{\prime}}}{\partial\bm{\widetilde{\delta}}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}\left(\frac{J_{0}}{2K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}}^{\prime})^{2}\bigg{]}$
(125)

$\displaystyle\mathcal{C}(\hat{\mathbf{q}},\hat{\mathbf{k}})=$
$\displaystyle\sum_{\bm{\widetilde{\delta}}}(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}})^{2}\bigg{[}\frac{1}{4}\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}J_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)+\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}K_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)\left(\frac{J_{0}}{2K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}})^{2}\bigg{]}$
(126)

For $\beta J_{0}\lesssim 1$, the leading temperature dependence is obtained
from the second integral in Eq. 28 with the quadratic dispersion,
$\frac{\Delta
v}{v}=\widetilde{c}_{1}+\mathcal{I}_{1}+\mathcal{I}_{2}+\mathcal{I}_{3}$ where
$\widetilde{c}_{1}$ is given in Eq. 123 and

$\displaystyle\mathcal{I}_{1}=$
$\displaystyle-\frac{3V}{8\pi^{2}MNv_{0}^{2}\widetilde{c}_{s}}\int_{6\beta
J_{0}}^{\infty}\frac{dx}{x(e^{x}-1)}\int
d\theta_{\hat{\mathbf{k}}}\widetilde{\mathcal{A}}_{-}(\hat{\mathbf{q}},\hat{\mathbf{k}})$
$\displaystyle\mathcal{I}_{2}=$
$\displaystyle\frac{27V}{4\pi^{2}MNv_{0}^{2}}\int_{6\beta
J_{0}}^{\infty}\frac{dx\sqrt{x}e^{x}}{\sqrt{\beta}(e^{x}-1)^{2}}\int
d\theta_{\hat{\mathbf{k}}}\frac{\widetilde{\mathcal{A}}_{+}(\hat{\mathbf{q}},\hat{\mathbf{k}})\cos{\theta_{\mathbf{q}\mathbf{k}}}}{v_{0}\sqrt{\widetilde{c}_{s}}}$
$\displaystyle\mathcal{I}_{3}=$
$\displaystyle\frac{9V}{8\sqrt{2}\pi^{2}MNv_{0}^{2}\widetilde{c}_{s}}\int_{6\beta
J_{0}}^{\infty}\frac{dx}{\beta(e^{x}-1)}\int
d\theta_{\hat{\mathbf{k}}}\widetilde{\mathcal{C}}(\hat{\mathbf{q}},\hat{\mathbf{k}})$
(127)

where the upper limit can be taken to infinity due to exponential damping and

$\displaystyle\widetilde{\mathcal{A}}_{+}(\hat{\mathbf{q}},\hat{\mathbf{k}})=$
$\displaystyle\sum_{\bm{\widetilde{\delta}},\bm{\widetilde{\delta}}^{\prime}}(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}})(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}}^{\prime})\bigg{[}\frac{1}{6}\frac{\partial
J_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}}\cdot\mathbf{e}_{-\mathbf{q}}+\frac{\partial
K_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}}\cdot\mathbf{e}_{-\mathbf{q}}\left(\frac{J_{0}}{K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}})^{2}\bigg{]}\bigg{[}\frac{1}{6}\frac{\partial
J_{\bm{\widetilde{\delta}}^{\prime}}}{\partial\bm{\widetilde{\delta}}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}+\frac{\partial
K_{\bm{\widetilde{\delta}}^{\prime}}}{\partial\bm{\widetilde{\delta}}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}\left(\frac{J_{0}}{K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}}^{\prime})^{2}\bigg{]}$
(128)
$\displaystyle\widetilde{\mathcal{A}}_{-}(\hat{\mathbf{q}},\hat{\mathbf{k}})=$
$\displaystyle\sum_{\bm{\widetilde{\delta}},\bm{\widetilde{\delta}}^{\prime}}(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}})(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}}^{\prime})\bigg{[}\frac{1}{2}\frac{\partial
J_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}}\cdot\mathbf{e}_{-\mathbf{q}}-\frac{\partial
K_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}}\cdot\mathbf{e}_{-\mathbf{q}}\left(\frac{J_{0}}{K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}})^{2}\bigg{]}\bigg{[}\frac{1}{2}\frac{\partial
J_{\bm{\widetilde{\delta}}^{\prime}}}{\partial\bm{\widetilde{\delta}}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}-\frac{\partial
K_{\bm{\widetilde{\delta}}^{\prime}}}{\partial\bm{\widetilde{\delta}}^{\prime}}\cdot\mathbf{e}_{\mathbf{q}}\left(\frac{J_{0}}{K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}}^{\prime})^{2}\bigg{]}$
(129)

$\displaystyle\widetilde{\mathcal{C}}(\hat{\mathbf{q}},\hat{\mathbf{k}})=$
$\displaystyle\sum_{\bm{\widetilde{\delta}}}(\hat{\mathbf{q}}\cdot\bm{\widetilde{\delta}})^{2}\bigg{[}\frac{1}{6}\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}J_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)+\left(\mathbf{e}_{-\mathbf{q}}\cdot\frac{\partial^{2}K_{\bm{\widetilde{\delta}}}}{\partial\bm{\widetilde{\delta}}^{2}}\cdot\mathbf{e}_{\mathbf{q}}\right)\left(\frac{J_{0}}{K_{0}}\right)(\hat{\mathbf{k}}\cdot\bm{\widetilde{\delta}})^{2}\bigg{]}$
(130)

At exactly $J_{0}=0$, the lower limit of the integrals in Eq. 127 is zero and
the integrals diverge suggesting an instability at the $J_{0}=0$ point. At
finite but small $J_{0}$, on performing numerical integration, we see the
approximate temperature dependence as mentioned in Eq. 32.

##  Appendix E Details of the Landau-Ginzburg Theory

The point group symmetry of ${\text{NiGa}_{2}\text{S}_{4}}$ is
$\text{D}_{3d}$. This, along with translations by lattice vectors
($\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ in Fig. 6 in the plane of the
triangular lattice of magnetic atoms) generate the set of relevant symmetry
operations. Notably, the point group $\text{D}_{3d}$ is $\text{D}_{3}\otimes
C_{i}$ (where $C_{i}$ consists of identity($E$) and inversion($I$)
operations).Hamermesh (1989) $\text{D}_{3}$ consists of a $3$-fold rotation
axis and a $2$-fold rotation axis perpendicular to the $3$-fold axis (The
$2$-fold rotation axes $a$,$b$ and $c$ are shown in Fig.6). Combining with
inversion this gives $12$ elements divided into $6$ classes for
$\text{D}_{3d}$ : ($E$), ($C_{3},C_{3}^{2}$),
($C_{2,a}^{\prime},C_{2,b}^{\prime},C_{2,c}^{\prime}$), ($I$),
($IC_{3},IC_{3}^{2}$) and
($IC_{2,a}^{\prime},IC_{2,b}^{\prime},IC_{2,c}^{\prime}$).

For the 3-sublattice dipolar and quadrupolar orders that we are interested in,
we consider a single triangle as the unit cell. The non-trivial
transformations for the triangle formed by the points $1$,$2$,$3$ in Fig. 6
which keep its centre (point x) fixed are $C_{3}$ and
$\sigma_{h}C_{2}^{\prime}$ where $\sigma_{h}$ denotes reflection about the
horizontal plane (the plane containing the triangular lattice). The axes
$a^{\prime}$,$b^{\prime}$ and $c^{\prime}$ for the $C_{2}^{\prime}$ rotation
are shown in Fig.6. $\sigma_{h}$ is required with the $2$-fold rotation to
bring the crystal field environment back to itself. Combining these with
inversion about a lattice point generates the $12$ elements of $D_{3d}$ by up
to translations by lattice vectors.

Figure 6: Triangular lattice of magnetic $\text{Ni}^{2+}$. Axes for the
different symmetry transformations of $D_{3d}$ are shown

Denoting the transformations with an argument indicating the coordinate point
that is left unchanged by the transformations ( eg. $C_{3}(1)$ indicates the
operation of $C_{3}$ rotation keeping point $1$ fixed, $I(1)C_{3}$(x)
indicates the operation of $C_{3}$ rotation keeping point x fixed followed by
inversion about the site now at location $1$), for the elements apart from
identity and inversion, we have:

$\displaystyle T_{-\mathbf{a}_{1}}C_{3}(\text{x})=$ $\displaystyle C_{3}(1)$
$\displaystyle T_{-\mathbf{a}_{2}}C_{3}^{2}(\text{x})=$ $\displaystyle
C_{3}^{2}(1)$ $\displaystyle
T_{-\mathbf{a}_{1}}\sigma_{h}C_{2,a^{\prime}}^{\prime}(\text{x})=$
$\displaystyle I(1)C_{2,a}^{\prime}(1)$
$\displaystyle\sigma_{h}C_{2,b^{\prime}}^{\prime}(\text{x})=$ $\displaystyle
I(1)C_{2,b}^{\prime}(1)$ $\displaystyle
T_{-\mathbf{a}_{2}}\sigma_{h}C_{2,c^{\prime}}^{\prime}(\text{x})=$
$\displaystyle I(1)C_{2,c}^{\prime}(1)$ $\displaystyle
T_{\mathbf{a}_{1}}I(1)C_{3}(\text{x})=$ $\displaystyle I(1)C_{3}(1)$
$\displaystyle T_{\mathbf{a}_{2}}I(1)C_{3}^{2}(\text{x})=$ $\displaystyle
I(1)C_{3}^{2}(1)$ $\displaystyle
T_{\mathbf{a}_{1}}I(1)\sigma_{h}C_{2,a^{\prime}}^{\prime}(\text{x})=$
$\displaystyle C_{2,a}^{\prime}(1)$ $\displaystyle
I(1)\sigma_{h}C_{2,b^{\prime}}^{\prime}(\text{x})=$ $\displaystyle
C_{2,b}^{\prime}(1)$ $\displaystyle
T_{\mathbf{a}_{2}}I(1)\sigma_{h}C_{2,c^{\prime}}^{\prime}(\text{x})=$
$\displaystyle C_{2,c}^{\prime}(1)~{}.$

###  E.1 Transformation of dipoles and quadrupoles under point group
symmetries

The symmetry transformations for irreps of the dipoles on the triangle of Fig.
1 listed in Table 1:

$\displaystyle C_{3}:\quad m_{a}^{\alpha}\xrightarrow{}$ $\displaystyle
m_{a}^{\alpha}$ $\displaystyle m_{e1}^{\alpha}\xrightarrow{}$
$\displaystyle-\frac{1}{2}m_{e1}^{\alpha}+\frac{\sqrt{3}}{2}m_{e2}^{\alpha}$
$\displaystyle m_{e2}^{\alpha}\xrightarrow{}$
$\displaystyle-\frac{\sqrt{3}}{2}m_{e1}^{\alpha}-\frac{1}{2}m_{e2}^{\alpha}$
$\displaystyle\sigma_{h}C_{2,a^{\prime}}^{\prime}:\quad
m_{a}^{\alpha}\xrightarrow{}$ $\displaystyle m_{a}^{\alpha}$ $\displaystyle
m_{e1}^{\alpha}\xrightarrow{}$ $\displaystyle-m_{e1}^{\alpha}$ $\displaystyle
m_{e2}^{\alpha}\xrightarrow{}$ $\displaystyle m_{e2}^{\alpha}$

Further, inversion about a lattice point of a sublattice interchanges the
other two sublattices. This action of interchange of two sublattices is the
same as that of $\sigma_{h}C_{2}^{\prime}$ about the appropriate two fold
axis. Under translation $T_{\mathbf{a}_{1}}$, sublattice
$1\xrightarrow{}2\xrightarrow{}3\xrightarrow{}1$ which is also the action of
$C_{3}$. Under translation $T_{\mathbf{a}_{2}}$, sublattice
$1\xrightarrow{}3\xrightarrow{}2\xrightarrow{}1$ which is also the action of
$C_{3}^{2}$. Under global SU(2) spin rotation, such that we have
$m_{\widetilde{I}}^{\alpha}\xrightarrow{}R^{\alpha\beta}m_{\widetilde{I}}^{\beta}$
for $\widetilde{I}=a,e1,e2$, where $R^{\alpha\beta}$ are the SU(2) spin
rotation matrices for spin-1. The dipoles are odd under time reversal, so
$m_{\widetilde{I}}^{\alpha}\xrightarrow{}-m_{\widetilde{I}}^{\alpha}$.

The symmetry transformations for irreps of the quadrupoles on the up triangle
of Fig. 1 listed in Table 1:

$\displaystyle C_{3}:\quad Q_{a}^{\alpha\beta}$ $\displaystyle\xrightarrow{}$
$\displaystyle Q_{a}^{\alpha\beta}$ (131) $\displaystyle Q_{e1}^{\alpha\beta}$
$\displaystyle\xrightarrow{}$
$\displaystyle-\frac{1}{2}Q_{e1}^{\alpha\beta}+\frac{\sqrt{3}}{2}Q_{e2}^{\alpha\beta}$
$\displaystyle Q_{e2}^{\alpha\beta}$ $\displaystyle\xrightarrow{}$
$\displaystyle-\frac{\sqrt{3}}{2}Q_{e1}^{\alpha\beta}-\frac{1}{2}Q_{e2}^{\alpha\beta}$
$\displaystyle\sigma_{h}C_{2,a^{\prime}}^{\prime}:\quad Q_{a}^{\alpha\beta}$
$\displaystyle\xrightarrow{}$ $\displaystyle Q_{a}^{\alpha\beta}$
$\displaystyle Q_{e1}^{\alpha\beta}$ $\displaystyle\xrightarrow{}$
$\displaystyle-Q_{e1}^{\alpha\beta}$ $\displaystyle Q_{e2}^{\alpha\beta}$
$\displaystyle\xrightarrow{}$ $\displaystyle Q_{e2}^{\alpha\beta}$

Like the case for dipolar modes, inversion about a lattice point of a
sublattice interchanges the other two sublattices. This action of interchange
of two sublattices is the same as that of $\sigma_{h}C_{2}^{\prime}$ about the
appropriate two fold axis. Under translation $T_{\mathbf{a}_{1}}$, sublattice
$1\xrightarrow{}2\xrightarrow{}3\xrightarrow{}1$ which is also the action of
$C_{3}$. Under translation $T_{\mathbf{a}_{2}}$, sublattice
$1\xrightarrow{}3\xrightarrow{}2\xrightarrow{}1$ which is also the action of
$C_{3}^{2}$. Similarly, under spin rotations and time reversal, we
respectively have
$Q_{\widetilde{I}}^{\alpha\beta}\xrightarrow{}R^{\alpha\rho}R^{\beta\lambda}Q_{\widetilde{I}}^{\rho\lambda}$
and
$Q_{\widetilde{I}}^{\alpha\beta}\xrightarrow{}Q_{\widetilde{I}}^{\alpha\beta}$
for $\widetilde{I}=a,e1,e2$.

###  E.2 Normal elastic modes and their transformations

We consider the up triangle of Fig.1 with $\mathbf{u}_{1}$,$\mathbf{u}_{2}$
and $\mathbf{u}_{3}$ being the displacements from mean position. The three
non-zero energy normal modes (in the basis {
$u_{1x}$,$u_{1y}$,$u_{2x}$,$u_{2y}$,$u_{3x}$,$u_{3y}$}) arranged as
irreducible representations are given in Table 2.

Irrep | Eigenvector
---|---
$\epsilon_{a}$ | {$\frac{-1}{2}$,$\frac{-1}{2\sqrt{3}}$,$\frac{1}{2}$,$\frac{-1}{2\sqrt{3}}$,$0$,$\frac{1}{\sqrt{3}}$}
$\epsilon_{e1}$ | {$\frac{-1}{2\sqrt{3}}$,$\frac{-1}{2}$,$\frac{-1}{2\sqrt{3}}$,$\frac{1}{2}$,$\frac{1}{\sqrt{3}}$,$0$}
$\epsilon_{e2}$ | {$\frac{-1}{2}$,$\frac{1}{2\sqrt{3}}$,$\frac{1}{2}$,$\frac{1}{2\sqrt{3}}$,$0$,$\frac{-1}{\sqrt{3}}$}
Table 2: Normal modes of an up triangle

Under the lattice symmetries, the normal modes transform as:

$\displaystyle C_{3}:\quad\epsilon_{a}$ $\displaystyle\xrightarrow{}$
$\displaystyle\epsilon_{a}$ (132) $\displaystyle\epsilon_{e1}$
$\displaystyle\xrightarrow{}$
$\displaystyle-\frac{1}{2}\epsilon_{e1}+\frac{\sqrt{3}}{2}\epsilon_{e2}$
$\displaystyle\epsilon_{e2}$ $\displaystyle\xrightarrow{}$
$\displaystyle-\frac{\sqrt{3}}{2}\epsilon_{e1}-\frac{1}{2}\epsilon_{e2}$
$\displaystyle\sigma_{h}C_{2,a^{\prime}}^{\prime}:\quad\epsilon_{a}$
$\displaystyle\xrightarrow{}$ $\displaystyle\epsilon_{a}$
$\displaystyle\epsilon_{e1}$ $\displaystyle\xrightarrow{}$
$\displaystyle-\epsilon_{e1}$ $\displaystyle\epsilon_{e2}$
$\displaystyle\xrightarrow{}$ $\displaystyle\epsilon_{e2}$

The long wavelength strain modes,
$\varepsilon_{a},\varepsilon_{e1},\varepsilon_{e2}$ (where
$\varepsilon_{a}=\epsilon_{a}/a$, $\varepsilon_{e1}=\epsilon_{e1}/a$ and
$\varepsilon_{e2}=\epsilon_{e2}/a$ ; $a=1$ is the lattice spacing), can be
expressed in terms of Cartesian strains as given by Eq. 48. There are two
independent elastic constants associated with the two modes as given in Eq.
69. Using the standard relations for the elastic constants for $D_{3d}$
symmetry Lüthi (2007) and comparing the harmonic elastic energy in written in
terms of Cartesian strains with the Eq. 69 we obtain,

$\displaystyle c_{a}$ $\displaystyle=2(c_{11}+c_{12})$ $\displaystyle c_{e}$
$\displaystyle=2(c_{11}-c_{12})~{},$ (133)

where the elastic constants on the right hand side are written in Voigt
notation.

## References

  * Santini and Amoretti (2000) P. Santini and G. Amoretti, Phys. Rev. Lett. 85, 2188 (2000).
  * Chandra  _et al._ (2002) P. Chandra, P. Coleman, J. Mydosh, and V. Tripathi, Nature 417, 831 (2002).
  * Sato  _et al._ (2012) T. J. Sato, S. Ibuka, Y. Nambu, T. Yamazaki, T. Hong, A. Sakai, and S. Nakatsuji, Phys. Rev. B 86, 184419 (2012).
  * Patri  _et al._ (2019) A. S. Patri, A. Sakai, S. Lee, A. Paramekanti, S. Nakatsuji, and Y. B. Kim, Nature Communications 10, 4092 (2019).
  * Harter  _et al._ (2017) J. Harter, Z. Zhao, J.-Q. Yan, D. Mandrus, and D. Hsieh, Science 356, 295 (2017).
  * Voleti  _et al._ (2020) S. Voleti, D. D. Maharaj, B. D. Gaulin, G. Luke, and A. Paramekanti, Phys. Rev. B 101, 155118 (2020).
  * Chapon  _et al._ (2004) L. C. Chapon, G. R. Blake, M. J. Gutmann, S. Park, N. Hur, P. G. Radaelli, and S.-W. Cheong, Phys. Rev. Lett. 93, 177402 (2004).
  * Bhattacharjee  _et al._ (2011) S. Bhattacharjee, S. Zherlitsyn, O. Chiatti, A. Sytcheva, J. Wosnitza, R. Moessner, M. E. Zhitomirsky, P. Lemmens, V. Tsurkan, and A. Loidl, Phys. Rev. B 83, 184421 (2011).
  * Bhattacharjee  _et al._ (2016) S. Bhattacharjee, S. Erfanifam, E. L. Green, M. Naumann, Z. Wang, S. Granovsky, M. Doerr, J. Wosnitza, A. A. Zvyagin, R. Moessner, A. Maljuk, S. Wurmehl, B. Büchner, and S. Zherlitsyn, Phys. Rev. B 93, 144412 (2016).
  * Seth  _et al._ (2022) A. Seth, S. Bhattacharjee, and R. Moessner, Phys. Rev. B 106, 054507 (2022).
  * Blume and Hsieh (1969) M. Blume and Y. Hsieh, Journal of Applied Physics 40, 1249 (1969).
  * Chen and Levy (1971) H. H. Chen and P. M. Levy, Phys. Rev. Lett. 27, 1383 (1971).
  * Matveev (1973) V. Matveev, Zh. Eksp. Teor. Fiz 65, 1626 (1973).
  * Chandra and Coleman (1991) P. Chandra and P. Coleman, Phys. Rev. Lett. 66, 100 (1991).
  * Bhattacharjee  _et al._ (2006) S. Bhattacharjee, V. B. Shenoy, and T. Senthil, Phys. Rev. B 74, 092406 (2006).
  * Penc and Läuchli (2011) K. Penc and A. M. Läuchli, in  _Introduction to Frustrated Magnetism_ (Springer, 2011) pp. 331–362.
  * Nakatsuji  _et al._ (2005) S. Nakatsuji, Y. Nambu, H. Tonomura, O. Sakai, S. Jonas, C. Broholm, H. Tsunetsugu, Y. Qiu, and Y. Maeno, Science 309, 1697 (2005).
  * Tsunetsugu and Arikawa (2006) H. Tsunetsugu and M. Arikawa, Journal of the Physical Society of Japan 75, 083701 (2006).
  * Läuchli  _et al._ (2006) A. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. 97, 087205 (2006).
  * Stoudenmire  _et al._ (2009) E. M. Stoudenmire, S. Trebst, and L. Balents, Phys. Rev. B 79, 214436 (2009).
  * Kittel (1960) C. Kittel, Phys. Rev. 120, 335 (1960).
  * Barma (1975) M. Barma, Phys. Rev. B 12, 2710 (1975).
  * Nakatsuji  _et al._ (2007) S. Nakatsuji, Y. Nambu, K. Onuma, S. Jonas, C. Broholm, and Y. Maeno, Journal of Physics: Condensed Matter 19, 145232 (2007).