Dynamics and Instabilities of the Shastry-Sutherland Model
Zhentao Wang Cristian D. Batista






Department of Physics and Astronomy, The University of Tennessee,
Knoxville, Tennessee 37996, USA



Department of Physics and Astronomy, The University of Tennessee,
Knoxville, Tennessee 37996, USA

Quantum Condensed Matter Division and Shull-Wollan Center, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA



We study the excitation spectrum in the dimer phase of the  Shastry-Sutherland model by using an unbiased variational method
that works in the thermodynamic limit. The method outputs dynamical correlation functions in all possible channels. 


This output is exploited to identify the order parameters with the highest susceptibility (single or multitriplon condensation in a specific channel)
upon approaching a quantum phase transition in the magnetic field versus the J'/J phase diagram.
We find four different instabilities: antiferro spin nematic, plaquette spin nematic, stripe magnetic order,
and plaquette order, two of which have been reported in previous studies.


 

Dynamics and Instabilities of the Shastry-Sutherland Model
    Cristian D. Batista
    December 30, 2023
==========================================================

The Shastry-Sutherland model (SSM) has become a paradigmatic Hamiltonian of frustrated quantum
magnetism <cit.> because it includes an exactly
solvable ground state <cit.>, very heavy low-energy
excitations <cit.>,
exotic phases obtained upon varying the ratio J'/J between two competing exchange constants <cit.>,
and a series of magnetic field induced magnetization plateaux <cit.>.
Its realization in SrCu_2(BO_3)_2 <cit.>
enabled various experimental studies, including magnetization <cit.>,
specific heat <cit.>, inelastic neutron scattering (NS) <cit.>,
far-infrared <cit.>, electron spin resonance (ESR) <cit.>,
Raman scattering <cit.>, and nuclear magnetic resonance
(NMR) <cit.>. These studies revealed that a 
finite Dzyaloshinskii-Moriya (DM) interaction <cit.>
must be added to the SSM in order to account for several properties of SrCu_2(BO_3)_2 <cit.>.

Despite the great theoretical efforts devoted to the SSM, the problem is still far
from being solved. Perturbative approaches are only applicable in  narrow regimes and
conventional numerical methods suffer from severe
size effects. As a consequence,  the nature of the quantum phase diagram  
has been debated for a long time <cit.>.
It is thus desirable to develop and apply alternative approaches. The infinite projected entangled-pair states (iPEPS) is an example of an alternative approach that works in the thermodynamic limit <cit.>. However, it relies heavily on the initial guess of the physical states
and it is difficult to extract dynamical responses.

In this Letter, we introduce an unbiased variational method  to calculate the excitation spectrum
and dynamical responses (susceptibilities) of the SSM in the dimer phase [A numerical method is said to be “variational” when it satisfies the variational principle. Other examples of unbiased variational methods include the density matrix renormalization group <cit.> and matrix product method <cit.>, etc.]. The method 
works in the thermodynamic limit and it complements alternative approaches like iPEPS.  The basic idea  was originally introduced to compute the single-hole 
dispersion of the square lattice t-J model <cit.>. The same
method was  applied to the Shastry-Sutherland lattice t-J model <cit.>.
Here we exploit this idea to compute dynamical correlators and dominant instabilities.  By working in a reduced Hilbert space, which
preserves all model symmetries, we obtain low
energy excitations classified by  quantum numbers. We then predict the character of the neighboring phases by detecting the order parameter with highest susceptibility.
Besides confirming the previously reported plaquette order and antiferro spin-nematic phases,  we  find two new phases; namely,
a plaquette spin-nematic phase and stripe magnetic ordering,  induced by simultaneously  increasing the magnetic field and  J'/J. 
In particular, the plaquette spin-nematic phase explains the nature of the two-triplon states (pinwheels) that crystallize at higher field values <cit.>.


We consider the spin-1/2 SSM under a magnetic field <cit.>:

    ℋ=J∑_⟨ ij⟩S_i·S_j+J^'∑_⟨⟨ ij⟩⟩S_i·S_j-h∑_iS_i^z,

where ⟨ ij⟩ and ⟨⟨ ij⟩⟩
denote intradimer and interdimer neighbors. 

The unit cell has 4 sites (see Fig. <ref>).

 The exact ground state for small enough J^'/J and h 
is a direct product of singlet states on all dimers <cit.>.
The elementary excitations of this “dimer phase” are singlet-triplet excitations known as triplons.
Triplons are dressed by quantum fluctuations with a magnetic correlation length ξ that increases with
J'/J. Methods that can account for the spatial extent of these quantum fluctuations should provide a good description of the low-energy excitation  spectrum of the dimer phase. 










We start the process by creating local excited
states |φ_i⟩ (e.g., single and two triplons). 

We then project these representative
states into subspaces with fixed momentum k,

    |φ_i(k)⟩≡P̂_k |φ_i⟩/√(⟨φ_i|P̂_k|φ_i⟩),

where the projector is defined as P̂_k≡1/N∑_re^ik·rT̂(r).
N→∞ (thermodynamic limit) is the total number of unit cells, and T̂(r)
is the translation operator. Application of ℋ
to  |φ_i(k)⟩ generates
new states that dress the corresponding quasiparticle excitation. 

This procedure can be applied iteratively to systematically improve the variational space. 
After M iterations, we obtain a basis {φ_i(k)}
with good quantum numbers k and S_tot^z. 
The number of iterations determines the spatial range l of the fluctuations that dress the quasiparticle,
so the method is then expected to produce accurate results for l ≳ξ.

The  eigenvalues and eigenvectors of the Hamiltonian restricted to the variational space
are obtained by applying the implicitly restarted Arnoldi method <cit.>. 
The eigenvectors are classified by the Little Group of C_4v
for each momentum k. A continuous phase transition manifests via a  
vanishing gap (condensation)  that  signals a phase transition  into 
a broken symmetry state. The symmetry of the new state is  determined by the irreducible
representation (IREP) of the eigenstate that becomes gapless. 
To keep the method unbiased,   the
initial basis must not break the point group symmetry of ℋ [This means that the basis generates an invariant subspace of the group. In addition, a fully symmetric initial basis leads to better
convergence (as a function of M) in comparison with other choices Ref. <cit.>. ].

For illustration, we first focus on the S_tot^z=0 sector
relevant to h=0.  We include 𝒟=14 S_tot^z=0
initial states to start the iteration <cit.> and then
apply Eq. (<ref>) onto this basis to systematically increase the variational space [The computational cost is exponential in the linear size of the quasiparticle. In the present case, the dimension of the truncated Hilbert space roughly scales as 𝒟∝ 6^M ∝ 6^l. For the S_tot^z=0 sector, the dimension of the variational space is 𝒟=74 for M=1, 𝒟=396 for M=2, …, 𝒟=17  730  330 for M=8.]. After obtaining the lowest energy eigenstates, we use the
eigenfunction to calculate S_tot and its IREP [In principle, one can switch the order of “diagonalization” and “point group symmetry analysis,” to reach larger M and get better convergence. Similar tricks can be found in Refs. <cit.>.].


In contrast to the result obtained with perturbative continuous unitary transformations (CUTs) <cit.>,
we find that the first instability
as a function of J^'/J (for h=0)  takes place in the S_tot=0 channel
with IREP A_2 [The S=1 instability found in Ref. <cit.> is checked by our method both in S_tot^z=0 and S_tot^z=1 sectors. The energy difference between the two sectors is negligible.]. 





Figure <ref> shows the evolution of the gap as a function of M. 
Convergence is reached beyond M=3 for J^'/J≲0.5, but the increase of ξ slows down the convergence for 
larger  J^'/J.  Although Eq. (<ref>) does not conserve
the triplon number, the  state that condenses is adiabatically  connected
with the corresponding  S_tot=0 IREP A_2 pure two-triplon state in the J^'/J → 0 limit (see Fig. <ref>).



We can read out the critical value of J^'/J when this state
condenses. Figure <ref>(a) shows the evolution of the critical
value as we increase M (circles). At M=8, (J_c^'/J)^(M=8)≈0.702.
Previous tensor network based calculations <cit.> showed that the

transition  is actually of first order and the transition point is  at J_c^'/J=0.675 <cit.>. A susceptibility analysis, like the one presented here, is in general inadequate to 
detect first order transitions. However, it can still be used to detect the nature of the order parameter if the system still transitions into
the broken symmetry state with highest susceptibility [A typical example is the case of attractive interaction between the modes (particles) that become soft. Because of the attraction (negative quartic term in a Ginzburg-Landau expansion) the particle density changes discontinuously (first order transition) before the single-particle excitation becomes gapless. ]. Given that the first-order transition takes place when this susceptibility is still finite, J_c^'/J turns out to be smaller than the value at which the susceptibility becomes divergent. This observation explains the difference between the values of J_c^'/J obtained with both approaches. In addition, it  illustrates their complementary nature. The unbiased susceptibility analysis can be used to detect candidates for broken symmetry states. These candidates can then be tested with biased variational techniques, such as iPEPS, which can produce more accurate values of the transition point.





Since the two-triplon bound state has S_tot=0, the
new ground state (condensate) must be nonmagnetic. Furthermore, since the A_2 state is odd (even)
under reflection (rotation) <cit.>, the new ordered state should only break
reflection but not rotation symmetries. These features are consistent with the previously reported plaquette ordering <cit.>.
 Figure <ref>(b) shows a schematic plot of the corresponding bond ordering.
As expected, ⟨S_i·S_j⟩ becomes different 
on different plaquettes and there is no magnetic order. 
In other words, the plaquette order parameter can be defined as ⟨S_i ·S_j  - S_i ·S_j^'⟩, where ij and ij^' are two bonds related by a mirror reflection [see Fig. <ref>(b)].

We consider now  the case of nonzero magnetic field. The energy of excited states with finite 
S_tot^z  decreases linearly in h.
Figure <ref>(b) shows that the dominant instability for  h/J=0.2 corresponds to 
condensation of a state  with S_tot=1, k=(0,0), and IREP E, leading to the stripe magnetic order depicted in 
Fig. <ref>(c). The IREP E is
a two-dimensional representation corresponding to the two possible  directions of the
stripes  (along a_1 or a_2). We note that
the two same-color spins in the same unit cell
are not identical (the corresponding mirror symmetry is broken). 






The stripe state has the highest susceptibility over a  narrow range  0.66≲ J^'/J≲0.70
for M=8 iterations  [see Fig. <ref>(a)].  
Because of the frustrated exchange interactions, the energies of a few other states are not much higher than the stripe state <cit.>.
Among them, the lowest one is a two-triplon state with k=(0,0) and IREP B_1, corresponding to vector chiral order <cit.>.
Although their energies are slightly higher than the stripe magnetic instability within M ≤ 8, the situation may change in the M →∞ limit, or if small perturbations are added to the original model. 



The S_tot=2 excited states take over for higher magnetic field values. 
Figure <ref>(c) shows that for h/J=0.6 the lowest excited state is the S_tot=2 two-triplon bound state with momentum k=(0,0) and IREP A_2.
The fact that this state and the S_tot=0 state that condenses at zero field  belong to the same point group  IREP A_2 indicates that the condensation
of the S_tot=2 A_2 state also leads to “plaquette”
ordering [shown in Fig. <ref>(d)]; the difference being that the S_tot=2 condensate also  breaks the U(1) symmetry group of global spin rotations along the field direction, leading to spin-nematic ordering.
In other words, 
the local bond order parameter is ⟨ S_i^+S_j^+ - S_i^+S_j^'^+⟩ instead of 
⟨S_i·S_j - S_i·S_j'⟩
(ij and ij' denote two bonds connected by a mirror reflection σ_d, see Fig. <ref>).


As indicated in Fig. <ref>(a), the “plaquette spin-nematic” state covers a  wide range 0.40≲ J^'/J≲0.66. It has been shown in Ref. <cit.> that
the 1/8 plateau at slightly higher magnetic field values and
J^'/J=0.63 is induced by crystallization of S_tot^z=2
bound states. A closer scrutiny of the “pinwheel” structure of
these  bound states  shows that they
locally preserve rotational symmetries, while breaking reflection
symmetries <cit.>; i.e.,  they are the same two-triplon  bound states that we are finding in the dilute limit. 

Moving away from the plaquette spin-nematic phase towards the J^'/J ≪ 1 limit, 
it is already known from an early perturbative calculation that S_tot=2 two-triplon bound states with k=(π,π) give the highest susceptibility <cit.>.
This is confirmed by our variational method [see Fig. <ref>(d)].

Since the two-triplon bound state has momentum k=(π,π), the corresponding ordered state also
breaks translational symmetry. As shown in Fig. <ref>(e), ⟨ S_i^+S_j^+⟩
changes sign going from one unit cell to its nearest neighbors. 
Similar to the case of the  stripe ordering (which also comes from condensation of
IREP E states), there are two choices for aligning the bonds. 
We note that the breaking of the C_4 lattice rotational symmetry leads to a modulation of ⟨ S_i^z ⟩ that can be detected with NS experiments.


It is worth mentioning that the two spin-nematic states found in this Letter are different from nematic phases discussed in various other contexts <cit.>.
The so-called “Ising-nematic” ordering corresponds to (discrete)  lattice rotation symmetry breaking. In contrast, “spin-nematic” ordering corresponds to  broken spin rotational symmetry.   The spin-nematic
orderings discussed in this Letter break both the point group symmetry and
spin rotation symmetry <cit.>.
In other words, they are simultaneously Ising nematic and spin nematic.

The frustrated nature of the SSM makes the calculation of dynamical response a difficult task. To date, the only calculation including multitriplon contributions is  the perturbative CUTs, which breaks down for J^'/J ≳ 0.63 <cit.>.
The variational Hilbert space generated by our method thus provides a more reliable access
to dynamical responses via the continued fraction method <cit.>. 

Near the phase boundaries, we expect the
susceptibilities of corresponding order parameters to diverge at ω=0. 
Magnetic orderings, such as the stripe phase, are detected by computing the
dynamic structure factor (DSF) <cit.>:

    S^-+(k,ω)=2π∑_ν|⟨ν|S_k^+|0⟩|^2δ(ω+E_0-E_ν),

which is measured with inelastic NS. As shown Fig. <ref>(b),
the lowest peak of S^-+(k,ω) approaches
ω=0 near the phase boundary indicating condensation
of an S_tot=1 state.

The  divergent susceptibilities of the other phases are revealed by computing two-point dynamical 
correlation functions of the corresponding order parameters. These order parameters are the operators that
create a state that has finite overlap with the one that is condensing. For J^'/J≪1, the lowest energy S_tot=2 eigenstates are known to be a linear
combination of triplons located on nearest (and next-nearest) neighbors <cit.>.
Denoting the order parameter as A_k^S2E, 

the corresponding susceptibility is 

    χ_S2E(k,ω)=2π∑_ν|⟨ν|A_k^S2E|0⟩|^2δ(ω+E_0-E_ν).

Similarly, using the approximate wave functions of the S_tot=2
A_2 state and the S_tot=0 A_2 state <cit.>,
we can also construct the order parameters and compute the corresponding susceptibilities χ_S2A_2(k,ω)
and χ_S0A_2(k,ω).  Figure <ref> shows the nearly divergent susceptibilities in each channel by picking
appropriate Hamiltonian parameters near the phase boundaries.



While the tendency toward stripe ordering can be detected with  inelastic NS,
the experimental detection of the other  phases is nontrivial. Lattice distortions  induced by
the order parameter through magnetostriction  can provide indirect evidence if they are large enough to be 
detected <cit.>.
Experimental knobs, such as pressure, doping, and magnetic field can drive the material into different instabilities <cit.>. 
Thus, the method presented in this Letter provides valuable insight for revealing the nature of the new phases in such experiments.  However, the model relevant to SrCu_2(BO_3)_2 also includes DM interactions that  modify the single triplon dispersion and  can potentially change the phase diagram reported here. In addition, DM interactions reduce the spin rotational symmetry of the model, implying that they can change the nature of the order parameters.


We emphasize that the applicability of the method is not restricted to the SSM considered here.
The same method can be used to detect the instabilities of other quantum paramagnets <cit.>. 
Especially,  it is very difficult to enumerate all the possible instabilities  for highly frustrated systems. 
The low energy spectrum produced by our method provides a valuable educated guess for biased numerical approaches.
Given that the method works in the thermodynamic limit, it can also detect incommensurate instabilities, which cannot be handled by most numerical methods.


In summary, we have used an unbiased variational method  to study the excitation spectrum
of the SSM in the dimer phase. Several instabilities are found next
to the dimer phase corresponding to condensations of single-triplon or two-triplon bound states. Two
of the instabilities (antiferro spin nematic and plaquette) are known
from previous studies and the others (plaquette spin nematic and stripe)
are newly discovered in this Letter. The same method can  be used to compute relevant 
dynamical correlation functions. 

We thank S. Haravifard, B. Shastry, G. Ortiz, S. Zhang, and H. Suwa for helpful
discussions. Z.W. and C.D.B. are supported by funding from the Lincoln
Chair of Excellence in Physics. This work used the Extreme Science
and Engineering Discovery Environment (XSEDE) <cit.> through
allocation TG-DMR170029, which is supported by NSF Grant No. ACI-1548562.
This research used resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

apsrev4-1


 
 
 




   
—Supplemental Material— 





§ CHARACTER TABLE OF C_4V GROUP


The point group of Shastry-Sutherland lattice is C_4v, with 5 classes of symmetry operations:
identity E, two rotations {C_4,C_2} and two reflections
{σ_v,σ_d} (see Table. <ref>). 





§ CHOICE OF INITIAL BASIS


The initial basis is always chosen such that it does not break the
point group symmetries explicitly.
The choices for different values of S_tot^z are listed in Figs. <ref>,
<ref> and <ref>.










§ CONSTRUCTION OF ORDER PARAMETERS


In the J^'/J≪1 limit, the S_tot=2, k=(π,π),
IREP E states are known from Ref. <cit.>.

Besides the S_tot^z=2 initial
states shown in Fig. <ref>, there is another set of relevant states
(see Fig. <ref>) generated at M=2. Here we write down the diagonalized wavefunctions
for Q≡ (π,π):

    |Ψ_1(Q)⟩_S2E
       ≈α[|φ_7^(2)(Q)⟩+|φ_8^(2)(Q)⟩]+β[|φ_2^(2)(Q)⟩-|φ_1^(2)(Q)⟩],
         |Ψ_2(Q)⟩_S2E
       ≈α[|φ_6^(2)(Q)⟩-|φ_5^(2)(Q)⟩]+β[|φ_4^(2)(Q)⟩-|φ_3^(2)(Q)⟩],


where 

    |φ_i^(S_z)(Q) ⟩ = P̂_Q|φ_i^(S_z)⟩/√(⟨φ_i^(S_z) | P̂_Q |φ_i^(S_z)⟩).




The above wavefunctions have
S_tot=2, as long as the normalization condition 2|α|^2+2|β|^2=1
is satisfied. The values of α,β depend on 
J^'/J <cit.>. To construct the nematic 
order parameter, we only need a bound state wavefunction 
that has significant overlap with the exact one.
We then choose α=0.7, β=0.1 to
construct the order parameter 

    A^S2E_Q≡|Ψ_1(Q)⟩_S2E ⟨0|,.

The momentum dependence of the corresponding susceptibility is obtained by computing the two-point correlator after replacing  Q  with k.









We have demonstrated in the main text
that the lowest bound state wave function has S_tot=2,
k=(0,0), IREP A_2 for larger values of J^'/J. Similarly, with the eight states in
Figs. <ref>-<ref>, we can write down the approximate
K≡ (0,0) wavefunction :

    |Ψ(K)⟩_S2A_2
        =α[|φ_5^(2)(K)⟩-|φ_6^(2)(K)⟩-|φ_7^(2)(K)⟩+|φ_8^(2)(K)⟩]
         +β[-|φ_1^(2)(K)⟩-|φ_2^(2)(K)⟩+|φ_3^(2)(K)⟩+|φ_4^(2)(K)⟩].


Once again, S_tot=2 is guaranteed by the normalization condition
4|α|^2+4|β|^2=1. To construct the corresponding  order parameter 
A^S2A_2_k, we use α=0.47, β=0.17 in the main text.

We note that the exact bound state in this region  (larger J^'/J), 
has a larger size than the one in the J^'/J ≪ 1 limit. Consequently, A^S2A_2_k
produces more spectral weight at high energies in comparison to A^S2E_k [see Fig. 5(c)]. 

For even larger J^'/J with h=0, the lowest states have S_tot=0,
k=(0,0), IREP A_2. 
we can similarly construct the approximate wavefunction at K≡ (0,0):

    |Ψ(K)⟩_S0A_2
        =α[|φ_1^(0)(K)⟩+|φ_2^(0)(K)⟩+|φ_3^(0)(K)⟩+|φ_4^(0)(K)⟩
         -|φ_5^(0)(K)⟩-|φ_6^(0)(K)⟩-|φ_7^(0)(K)⟩-|φ_8^(0)(K)⟩]
         -β[|φ_9^(0)(K)⟩+|φ_10^(0)(K)⟩-|φ_11^(0)(K)⟩-|φ_12^(0)(K)⟩].


In this case, the normalization condition 8|α|^2+4|β|^2=1
is different from the S_tot=0 condition, which requires
α=β. Combining both conditions we have α=β=1/2√(3).

Once again, since this approximate wavefunction has a reduced overlap with the exact bound state, we see a significant amount of spectral weight at high energies in Fig. 5(a) in the main text.















§ DOMINANT INSTABILITY AT ZERO FIELD

In the S_tot=0 sector at h=0, we can roughly estimate J_c^'/J ≈ 0.69 in the M→∞ limit (see Fig. 3 in the main text). At this point J^'/J=0.69, we show In Fig. <ref> that the gaps of the other S_tot sectors remain finite in the M→∞ limit, even with a linear extrapolation. Note that the calculation converges for large enough M, implying that the slopes of the curves in Fig. <ref> should finally become zero. In other words, the true value of the gaps should be larger than that obtained from linear extrapolation.
Since all other gaps remain finite while the S_tot=0 gap is closing, we can safely conclude that the order parameter with highest susceptibility has S_tot=0 for h=0. 







§ CONVERGENCE OF DYNAMIC STRUCTURE FACTORS


Generally, the convergence speed of the DSFs is determined by ξ.

Larger ξ  values require larger number of iterations (M) to obtain qualitatively correct results. For instance, 
χ_S2E(k,ω) is almost fully converged at M=3 and  J^'/J=0.2 [see Fig. 5(d) in the main text] 
because ξ is  a few lattice spaces.







To have a more quantative sense of the convergence as a function of
M, we first benchmark the DSF S^-+(k,ω) with J^'/J=0.63,
h=0, which is relevant for the compound SrCu_2(BO_3)_2 <cit.>.
In Fig. <ref>, we show the evolution of the calculated S^-+(k,ω)
as a function of M. The lowest single and two-triplon
bands are already converged at M=5. We note that the method produces not only the bound states,
 but also the multitriplon continuum, which is qualitatively caputured at M=5 and about to
converge at  M≳ 8 for the parameters used in Fig. <ref>.

To measure the relative change between different iterations, we can define

    Δ_M (k)= √(1/Λ∫_0^Λ d ω( S^-+(k,ω)_M - S^-+(k,ω)_M-1/S^-+(k,ω)_M+ S^-+(k,ω)_M-1)^2),

where Λ is the energy cutoff. For the fully converged result, it is expected that Δ_M (k) → 0.



In Fig. <ref>, we show the linecuts of Fig. <ref> at two momenta: k-(4π,0) =  (0,0)  and  (π,π). Correspondingly, the relative errors Δ_M (k) are shown in the bottom panels of the same figure. The relative errors decrease rapidly as we increase M, reaching around 15% at M=8. We note that the value of Δ_M (k) depends on the energy cutoff Λ: the lowest relative errors in Fig. <ref> are found when Λ/J=1, consistent with the fact that low-energy excitations converge faster than the high-energy ones. Note also that Δ_M (k) depends heavily on the choice of Lorentzian broadening factor used in the DSF calculation: larger broadening factors lead to smaller relative errors.




It is also worth benchmarking with Ref. <cit.>, which
carried out the DSF calculation for SSM using the perturbative CUTs. Figure <ref> shows a
good qualitative agreement. However, there are two small differences: the band near 5meV produced
by our method is slightly higher than that obtained from CUT, which
could arise from incomplete convergence for M≲8 in our
method or missing higher order corrections in CUT. We also note that our method captures
more accurately the states in the continuum.






§ COMPETING CHIRAL STATE

The main text describes the S_tot=1 instability that  leads to the stripe magnetic order. 
As shown in Fig. 5(b),  the lowest instability from the DSF arises from 
condensation of the k=(0,0) mode. Figure <ref>
examines this point in more detail. Besides reproducing Fig. 5(b)
in higher resolution, we also denote the symmetries of the k=(4π,0)
peaks [k=(0,0) ≡ (4π,0)].



As we mentioned in the main text, the state that produces the lowest peak at k=(4π,0)
belongs to the IREP E. In Fig. <ref>, we see that the next peak at k=(4π,0)
corresponds to the IREP B_1 [with negligible weight at k=(4π,0)].
Up to M=8, we find
that this B_1 peak is always higher in energy (lower susceptibility) than the one corresponding to the  IREP E.
It is worth mentioning that condensation of the  B_1 state leads to the  vector chiral ordering depicted in Fig. <ref>.







