Nearly order from quantum disorder phenomena and its observation in a bosonic quantum anomalous Hall system
Fadi Sun Jinwu Ye



Substantially revised version of arXiv:1903.11134 



Department of Physics and Astronomy,
Mississippi State University, Mississippi State, Mississippi 39762, USA

Department of Physics, Capital Normal University,

Beijing 100048, China







We report a new many body phenomena called " Nearly order from quantum disorder phenomena" (NOFQD).

We demonstrate the NOFQD  in the experimentally realized weakly interacting Quantum Anomalous Hall system
of spinor bosons in an optical lattice.

At a zero Zeeman field h=0, there is a classically infinite degeneracy
due to a spurious spin SU(2) symmetry.
By developing a systematic order from quantum disorder (OFQD) analysis,
we find the quantum ground state to be a N=2 XY-AFM superfluid (SF) state and also evaluate its excitation spectrum
due to the effective potential generated by the OFQD at h=0.


At a high h, the system is in a Z-FM SF state.
At a small h>0, even there is no spurious symmetry anymore,
there is still a delicate competition between the Zeeman field effect at h >0 and
the effective potential



leading to a quantum phase transition between the two states.

We name this novel many body phenomena as NOFQD which describes how a OFQD state response non-trivially to a small deformation.


Starting from symmetry principle, we construct a Ginsburg-Landau (GL) type of effective action to
investigate the nature of the transition.



We establish intrinsic  connections between the phenomenological GL theory and the microscopic calculations on

the effective potential.


Connections with the bilayer quantum Hall system with a total filling factor ν_T=1 are made.
Some insightful analogy with NAdS_2/NCFT_1 ( where N also means nearly ) correspondence in the context of Sachdev-Ye-Kitaev models are hinted.


Two types of OFQDs are classified, one response trivially, another non-trivially to a small deformation to the Hamiltonian leading to NOFQD.
The NOFQD can be detected in the current cold atom  bosonic quantum anomalous Hall experiments and may also appear
in many other frustrated systems.





 Nearly order from quantum disorder phenomena and
its observation in the bosonic quantum anomalous Hall system 
     Fadi Sun and Jinwu Ye  
    December 30, 2023
===============================================================================================================





1. Introduction 

The investigation and control of spin-orbit coupling (SOC)
have become subjects of intensive research
in both condensed matter and cold atom systems
after the discovery of the topological insulators <cit.>.
In materials side, the SOC plays crucial roles
in various 2d or layered insulators, semi-conductor systems,
metals and superconductors without inversion symmetry.
The quantum anomalous Hall (QAH) effect was
experimentally realized in Cr doped Bi(Sb)_2Te_3 thin films
<cit.>.
In the cold atom side, using Raman schemes,
several experimental groups <cit.>
generated 2d SOC for neutral cold atoms
in both continuum and in optical lattices.
Especially, the bosonic analog of the QAH
for spinor bosons ^87Rb was realized in <cit.>,
and the lifetime of SOC ^87Rb Bose-Einstein condensation (BEC) have been improved from 300 ms to 900 ms recently.


Motivated by the recent experiment on bosonic QAH <cit.>,
we investigate possible quantum many body phenomena in the bosonic quantum anomalous Hall (BQAH) model in a square lattice.
In contrast to its fermonic counterpart, the bosonic QAH system is necessarily interacting.
Most importantly, we discover a new many body phenomena we name " Nearly order from quantum disorder phenomena" (NOFQD).
Our main results are presented in Fig.1. We achieve Fig.1 by the combination of specific microscopic calculations and
the symmetry based phenomenological classification approaches.

We first perform microscopic calculations at  weak interaction limit and also identify its limitations.
At a high Zeeman field h, we find the system is in a spin-polarized  SF state along the z direction ( named as Z-FM SF state ). We evaluate the SF Goldstone mode at k⃗=(0,0) and the roton mode at Q⃗=(π,π).
We delineate the behaviour of the roton mode as h decreases: the roton mode continues to drop
and eventually becomes a gapless quadratic mode ∼ v^2k^2 as h → 0^+.
At zero Zeeman field h=0, we identify a spurious spin SU(2) symmetry
which leads to  a classically infinitely degenerate ground state manifold.
We develop a systematic order from quantum disorder (OFQD) analysis to determine
an effective potential generated by the the OFQD effects. It selects the true quantum ground state from such
a manifold to be a N=2 XY-AFM SF state where N=2 stands for the number of BEC momenta.
It also opens a gap to the spurious quadratic roton mode.
However, at any small h >0, the spurious spin SU(2) symmetry disappears, a simple classical analysis leads to a unique classical ground state which is the Z-FM SF. Although  at a high h, the system is indeed in this Z-FM SF state, the classical picture may break down
at a small h>0. At a small h >0, even there is not any spurious symmetries,
there is still a delicate competition between the classical  Z-FM SF state at a small h >0 and
the N=2 XY-AFM superfluid (SF) state selected by the OFQD mechanism at h=0
which may lead to a quantum phase transition between the two states.
We name this new quantum many body phenomena as " Nearly order from quantum disorder phenomena “ (NOFQD) which
happens slightly away from the conventional OFDQ in the parameter space of a Hamiltonian.

At the mean field level, we find that at a finite small h, the competition between the Zeeman energy and the effective potential
first leads to a canted N=2 XY-AFM SF state (CAFM), then drives a transition to the Z-FM state.
Unfortunately, it becomes difficult to study the nature of the quantum phase transition
from the microscopic weak coupling  approach.

Then we change our strategy to symmetry based classification approach. Just from the constraints of the symmetries of the Hamiltonian,
we construct a phenomenological Ginsburg-Landau (GL) type of effective action to investigate
the NOFQD phenomena and the associated quantum phase transition from the a weak to strong Zeeman field.


Using the GL action, we first find the excitations in the CAFM and Z-FM at low and high Zeeman field respectively which include
the SF mode and the roton mode. It is the roton gap closing which drives the quantum phase transition from the CAFM to the Z-FM.
Then by performing a renormalization group (RG) analysis, we find the transition is in the same universality class as
the zero density Mott to superfluid transition,
subject to a marginally and dangerously irrelevant U(1) symmetry breaking  ( cubic ) anisotropic term.
So it has the critical exponents z=2, ν=1/2, η=0 subject to new logarithmic corrections.
This U(1) symmetry breaking cubic anisotropic term leads to experimentally observable logarithmic corrections to
various scaling functions in the quantum critical regime.

By comparing our symmetry based GL action with the microscopic calculations at h=0 and a high h, we establish the intrinsic and
profound connections between the phenomenological GL theory and the microscopic calculations, especially the
effective potential generated by the OFQD at h=0, therefore endow the broad impacts of OFQD from a new and profound  perspective.
The microscopic calculations and the phenomenological GL theory are two different, but complementary approaches.
Each has its own advantages and limitations. The combination of the two approaches lead to rather deep and complete understandings of
the NOFQD phenomena in the bosonic QAH systems.
We also make instructive mappings between the GL theory here and the Chern-Simon  GL theory
applied by one of the authors to study the exciton superfluid and quantum phase transitions (QPT)
in the bilayer quantum Hall systems at the total filling factor ν_T=1.

The scaling behaviours of various physical quantities such as the roton gaps, the specific heat, transverse magnetization and
the Lyapunov exponent/butterfly velocity

are derived.




In the conclusion section, in order to spell out the crucial connections and differences between the  OFQD and NOFQD. we classify the OFQD into two classes:
Type-I responses trivially to a deformation λ,  type-II does non-trivially and leads to the NOFQD.
Some interesting and insightful analogy with NAdS_2/NCFT_1 ( where N also means nearly )
correspondence in the context of Sachdev-Ye-Kitaev models are also hinted.

The ongoing experimental efforts to detect these novel NOFQD phenomena are critically examined.
Some future perspectives are outlined.

The rest of the paper is organized as into 4 parts.
 In part I consisting of Sec.2-5, we perform microscopic calculations at h=0 and strong h ≫ U at weak coupling limit.
 In part II consisting of Sec.6-9, to capture the NOFQD, we perform phenomenological effective action calculations at any h at weak coupling limit.
 Then in Part III consisting of Sec.10-11, we contrast the phenomenological effective action calculations in part II with the microscopic calculations in part I, therefore establish the intrinsic connections between the two complementary approaches.
 We discuss the experimental implications and  reach conclusions in Sec.12 and 13 respectively.






2.
The classical ground state at a high field h ≫ U and its possible breakdown in the small h ≲ U limit.

The experimentally realized quantum anomalous Hall model
of spinor bosons in a square lattice is described by
the Hamiltonian <cit.>

    ℋ=
    	   -t_0∑_⟨ ij⟩
    	a_i^†σ^z a_j
    	+it_so∑_⟨ ij⟩
    	a_i^†(𝐝_ij·σ⃗) a_j
       -h∑_i a_i^†σ^z a_i
    	+U/2∑_i n_i(n_i-1)
    	-μ∑_i n_i

where 
a_i=(a_i↑, a_i↓)^⊺
and ⟨ ij⟩ denotes a pair of nearest neighbor sites,
and 𝐝_ij=𝐫_i-𝐫_j.
U is the repulsive onsite interaction,
and μ is the chemical potential,





σ_x,y,z are three Pauli matrices,
and n_i=a_i^† a_i is the number of particles on site-i.




For the ^87Rb atoms used in the recent experiment <cit.>,
the two pseudo-spin components σ=↑,↓
denote the two hyperfine states |F=1,m_F=0⟩ and |F=1,m_F=-1⟩.





The Hamiltonian has the particle number conversation U_c(1) symmetry,
and a spin-orbital coupled [C_4× C_4]_D symmetry.
Under the mirror transformation
M =(-1)^i R_z(π)  T, h → -h, so h=0 has an enlarged Mirror symmetry.
So we only need to focus on h≥0 cases.
In this Letter,  we focus on the experimentally accessible regime 0≤ t_so≤ t_0 and U/t_0≪1.
The global phase diagram in most general parameter space
will be worked out in a much longer version <cit.>.

Any small h > 0 breaks the mirror symmetry and
split the two degenerate single-particle states
at momentum (0,0) and (π,π).
A straightforward mean field analysis predicts the condensation occurs at
either (0,0) or (π,π) depending on the sign of h.






Due to its spin alinement along the z axis, we name it Z-FM superfluid in Fig.1.

 When the spectrum minimum is located at (0,0) or (π,π), the non-interacting Hamiltonian takes a simple form

    H_0(k)=-[h+2t_0(cos k_x+cos k_y)]σ_z

   where the SOC t_so drops out.

   Let us consider h>0 case, thus the minimum is located at (0,0).
In the weak coupling limit U/t ≪ 1, by writing

    a_k↑→√(N_0)δ_k,0+ψ_k↑,  
    	a_k↓→ψ_k↓

  where N_0 is the number of condensate atoms,
  and N_s is the total number of lattice sites, thus n_0=N_0/N_s is the condensate fraction.




   We can perform the expansion ℋ=ℋ^(0)+ℋ^(1)+ℋ^(2)+⋯ where the superscript denotes the order in the quantum fluctuations.
    The zeroth order term ℋ^(0)=E_0= -1/2U n_0N_0  is the classical energy of the condensate.
  The vanishing of the linear term sets the value of the chemical potential μ=-h-4t_0+U n_0.
  Diagonizing  ℋ^(2) by a generalized 4 × 4 Bogliubov transformation leads to:

    H=E_0+ E^(2)_0 + ∑_k,s=±ω_s(k)( α_sk^†α_sk + 1/2)

where  E^(2)_0=-( h+4t_0+ 1/2 U n_0 )
and ω_±(k) are the two Bogliubov modes.

Since ω_+(𝐤)>ω_-(𝐤) always holds, so we focus on ω_-(𝐤)
which displays a gapless superfluid Goldstone mode near 𝐤=(0,0) and
a gapped roton mode at 𝐤=(π,π). In the long wavelength limit, we find
the superfluid Goldstone mode becomes linear
ω_-(k)=c|𝐤| where the SF velocity is

    c=√(2(4t_0^2-2t_so^2+th)Un_0/4t_0+h)

  which is shown in Fig.<ref>a. In the h → 0 limit, it reduces to c= √( n_0 U t (2- t^2_s/t^2) ).




In Fig.<ref>a, at a fixed t_so=1, as h decreases to 0 sitting on the boundary between  Z ↑ FM-SF at (0,0) and
  the Z ↓ FM-SF at (π,π), the roton mode gets lower and lower,
  then touches zero at Q⃗= (π,π). In the h → 0 limit, the roton mode can be
  extracted from ω_-(k) mode near k=(π,π)+q:

    ω_-(k)
    	=√(4h^2+v_4(q_x^2+q_y^2)^2)

where

    v_4=(2t_0^2-t_so^2)[t_0(4t_0+n_0U)-2t_so^2]/2t_0(4t_0+n_0U)

which depends on the the interaction n_0 U.

In the weak coupling limit n_0U≪ t_0, the v_4 can be simplified as
v_4=(2t_0^2-t_so^2)^2/(4t_0^2)=(v^2)^2 which remains finite even at the non-interacting limit n_0 U=0.
In the small h limit, the roton gap

    Δ^0_R(h)=2|h|

 which is independent of t_so and shown in the dashed line in Fig.4b.




As one decreases h from h ≫ U to h ≲ U, then to 0^+,
the roton mode gets lower and lower, then touches zero at h=0^+,
it reduces to the quadratic form ω_-(k)=√(v_4(q_x^2+q_y^2)^2)=v^2 q^2.
This behaviour may signify a possible first order transition at h_c=0^+.
This physical picture leads to the naive phase diagram shown in Fig.S1 in SM.
However, as to be shown below, this naive physical picture only holds when h is sufficiently large h > U,
but may break down <cit.> at a small h ≲ U.




3. The order from quantum disorder (OFQD) effects to determine the quantum ground state at h=0.

The single-particle energy spectrum U=0  of Eq.(<ref>) in the momentum space can be easily obtained:
ϵ_±,𝐤^h
=±√([h+2t_0(cos k_x+cos k_y)]^2
+4t_so^2(sin^2 k_x+sin^2 k_y)).
At h=0, its lower branch develops two degenerate minima at
momentum (0,0) and (π,π),
and two eigen-spinors are η_0=(1,0)^⊺
and η_π=(0,1)^⊺ respectively.
There is a classically degenerate family of states:

    Ψ_i=c_0η_0+(-1)^i_x+i_y c_πη_π

where the two complex number c_0 and c_π
satisfy normalization condition |c_0|^2+|c_π|^2=1.
This classically degenerate manifold is due to the spurious SU(2)_s symmetry at the mean field level.

By writing the spinor field as

Ψ= √(N_0)Ψ_0+ ψ, we can expand the Hamiltonian ℋ=ℋ^(0) + ℋ^(1)+ℋ^(2) +⋯.
The zeroth order term ℋ^(0) gives the classical ground state
energy E_0=-1/2Un_0N_0.
Setting the linear term ℋ^(1) vanish gives the value of the chemical potential μ=-4t+Un_0.
Diagonalizing the quadratic term  ℋ^(2) by a generalized 8 × 8  Bogliubov transformation <cit.> leads to:

    ℋ=E_0+E_ofd
    	+∑_l,𝐤∈ RBZω_l,𝐤α_l,𝐤^†α_l,𝐤

where E_ofd=-(4t_0+1/2n_0U)N_s
+1/2∑ω_l,𝐤
is quantum correction to the mean field ground-state energy,
and ω_l,𝐤 with l=1,..,4 represent 4 Bogoliubov spectrum.


It is convenient to parameterize c_0 and c_π as

    c_0=e^-iϕ/2cos(θ/2),  
        c_π=e^iϕ/2sin(θ/2)

then we can plot numerically  E_ofd(θ,ϕ)
as a function of θ and ϕ in Fig.2 where one can identify
the quantum ground states as θ=π/2 and ϕ_m=π/4+mπ/2 (m=0,1,2,3).
It  has a uniform density ⟨ n_i⟩=n_0 and
a XY-AFM ordered spin structure: ⟨ a_i^†σ⃗ a_i⟩
=(-1)^i_x+i_yn_0(cosϕ_m,sinϕ_m,0). It is 4 fold degenerate breaking the  joint C_4 symmetry.




After identifying the quantum ground-state as the N=2 XY-AFM SF state,
we can also evaluate all ω_1,2,3,4(k) in Eq.<ref>.
There are one linear ω_1 ∼ k SF Goldstone mode and
one ω_2 ∼ k^2 quadratic roton mode located at (0,0).
While ω_3(k) and ω_4(k) are two fully gapped higher energy modes.









4. The Roton mass gap generated from the OFQD mechanism in the N=2 XY-AFM state at h=0. 

Now we can have a better understanding of the spurious SU(2)_s symmetry and the origin of this spurious qudratic roton mode:
Indeed, it setting the two minima (0,0) and (π,π) into the kinetic energy, then the term involving t_so drops out,
so at the mean field level, the system has the SU(2)_s symmetry originating from the exact SU(2)_s symmetry at t_so=0.
While the N=2 XY-AFM state is a collinear state invariant under the rotation
e^i ασ⃗·n⃗, n⃗=(1,1,0), in fact, it is a XY-FM state in the rotated S̃Ũ(2) basis
𝐒̃_i=((-1)^i_x+i_yS_i^x,(-1)^i_x+i_yS_i^y,S_i^z) where the  SU(2)_s symmetry
is manifest ( See appendix B ). So its breaking to the FM state leads to the spurious quadratic FM roton mode which is remanent of the exact FM mode at t_so=0. However, the OFQD opens a gap to this spurious quadratic FM roton mode.


The OFQD mechanism picks up θ=π/2 and ϕ=π/4 as the quantum ground state.
Then we can expand the ground-state energy around the minimum as <cit.>:

    E_GS[θ,ϕ]=E_0+A/2δθ^2+B/2δϕ^2

where θ=π/2+δθ and ϕ=π/4+δϕ, the two coefficients A ∼ (n_0 U)^2/t and
B ∼ (n_0 U)^4/t can be extracted from Fig.M2.


  Using the commutation relations [1/2n_0δθ, δϕ ]= i ħ, one can see that
  the "quantum order from dis-order " mechanism generates a roton gap:

    Δ_R( h=0 )=2√(AB)/n_0

  which is the dashed line in Fig.M3b  ( also in the Fig.S1b in SM  ).
  The SF Goldstone mode remains unaffected due to its protection by the U(1) symmetry breaking.


5. The nearly order from quantum disorder (NOFQD) effects  at a small 0< h ≲ U:
 Mean field theory treatment in the microscopic calculation 



Any small h > 0 spoils the spurious SU(2)_s symmetry and also
splits the degenerate single-particle state
at momentum (0,0) and (π,π). So one can determine the ground state to be the classical Z-FM SF state
even at an infinitesimal small h.
However, when h ≲ U, the splitting is so small that any small interaction U
still mix the two states considerably. So even the spurious SU(2)_s symmetry disappears at any small h,
one must still take into account the OFQD effects at h=0
and study the delicate  competition between the OFQD generated effective potential in Fig.2 due to the
interaction U  with the Zeeman energy due to a small finite h.
We call this phenomena as NOFQD.





As shown in SM, the OFQD generated effective potential in Fig.2 can be cast into the form:

    E_ofd
    	=    E_ofd,0+N_s[1/4A(1+cos2θ)
       +1/64B(1-cos2θ)^2(1+cos4ϕ)+⋯].



  Appealingly, one can cast the energy density:
ℰ_ofd=E_ofd/N_s into a compact from

    ℰ_ofd
    	=
    	1/2A  (|c_0|^2-|c_π|^2)^2
    	+1/2B  [(c_0c_π^*)^2+(c_0^*c_π)^2]^2

where the coefficients A and B are positive numbers at least
in the weak coupling limit.







To get to finite temperatures, it is necessary to parameterize the quantum fluctuations in
the polar coordinate as:

    a_i(τ)
    	=√(n_i(τ))
    	e^iχ_i(τ)
    	[c_0,i(τ) η_0
    	+(-1)^ic_π,i(τ) η_π]

   and write the action as two parts
𝒮=∫ dτ(ℒ+δℒ_c).

    ℒ    =    ∑_𝐤a̅_𝐤
    	[∂_τ+ϵ_-,𝐤^h=0-μ]a_𝐤
    	+U/2∑_in_i(n_i-1)            
    δℒ_c
    	    =    ∑_i [ℰ_ofd(c_0,i,c_π,i)-nh(|c_0,i|^2-|c_π,i|^2)]

where  the second part contains the competition between the
order-from-disorder effects generated effective potential and the Zeeman effect


  With the parametrization Eq.(<ref>), we obtain

    ℒ_0
    	    =E_0
    	    +A/2cos^2θ
    	    +B/8sin^4θcos^22ϕ
    	    -n_0 hcosθ

where E_0 is independent of θ and ϕ.

The cos4ϕ term is a result of C_4 symmetry
and plays a similar role as the C_4 clock term.
Performing a minimization with respect to θ and ϕ leads to θ_0=arccos(h/h_c) and ϕ_0=π/4 with h_c=A/n_0.

In the h→ 0^+ limit, one recovers the N=2 XY-AFM SF state.

For 0<h<h_c, the spin-orbital  structure is
n_0((-1)^i_x+i_ysinθ_0cosϕ_0,
(-1)^i_x+i_ysinθ_0cosϕ_0,cosθ_0),
and named as a canted AFM (CAFM) SF state.
When h ≤ h_c, the spin structure becomes n_0(0,0,1) which is fully aligned to the z direction, so named the Z-FM state (Fig.3a).



   Unfortunately, in this microscopic NOFQD calculation, the BEC condensation density n_0 is taken as a constant
   at very beginning ( see above  Eq.<ref> ). The chemical potential μ=-h-4t + U n_0 also depends on n_0 and h.
   It is very difficult to capture the density fluctuations in a consistent way.
   To study the competition between the effective potential generated by the OFQD and the Zeeman field, then
   to capture the nature of the quantum phase transition near h ∼ h_c, one must take account the density fluctuations.
   To circumvent this difficult suffered in this this microscopic NOFQD calculation,
   one may turn to an independent phenomenological calculations based on the symmetries of the Hamiltonian.




6. Ginsburg-Landau type action from the symmetry principle to study the NOFQD phenomena: 

   To capture the density fluctuations in a consistent way, it is convenient to
   absorb the density/phase part into c_0 and c_π in Eq.<ref> and
   re-parameterize low energy quantum fluctuations at a small h ≪ U as

    a_iα=ψ_1η_0+(-1)^iψ_2η_π

   which includes both spin and charge fluctuations.

   Then one can find the density  n= |ψ_1|^2+|ψ_2|^2 and the spin:

    S^+    =(-1)^iψ_1^*ψ_2,   
    
        S^-   =(-1)^i ψ_2^*ψ_1,   
    
        S^z    =|ψ_1|^2-|ψ_2|^2


 Note that here we only expand the boson field in terms of the two minima of the kinetic term without assuming any symmetry breaking.
 In order to study the transition from the weak coupling SF to the strong
 coupling Mott phase an an integer filling, we construct the Ginsburg-Landau action consistent with all the symmetries
 of the Hamiltonain such as the translational symmetry, the U(1)_c symmetry, the [ C_4 × C_4 ]_D
 ( or [R_π/2× R_π/2]_D ) spin-orbital coupled symmetry
 and the  Mirror symmetry M at h=0:

    ℒ    = ψ_1^*∂_τψ_1+ψ_2^*∂_τψ_2
    	+v^2(|∇ψ_1|^2+|∇ψ_2|^2)    
       +r(|ψ_1|^2+|ψ_2|^2) +U/2(|ψ_1|^2+|ψ_2|^2)^2   
        +A/2(|ψ_1|^2-|ψ_2|^2)^2
    	+B/2[(ψ_1^*ψ_2)^2+(ψ_1ψ_2^*)^2]^2     
        -h(|ψ_1|^2-|ψ_2|^2) + ⋯

where the ⋯ means the higher order terms consistent with all the symmetries ( see SM Eq.S25 ),
r<0 ensures the system is always in a superfluid phase,
and v,U are all positive to make the system stable.
However, the sign of A and B can not be determined by the symmetry analysis, but can be
fixed by comparing with the microscopic calculations at the weak coupling limit.
The last h term is the Zeeman field coupled to S^z in Eq.<ref>.
It breaks the Mirror symmetry M.
The relations between these phenomenological parameters and the microscopic ones calculated in the previous sections
will be established in Sec.9.

Although this effective action was designed to study the transition from weak coupling to the strong coupling.
Here, we only use it to study the transition driven by the Zeeman field,
while leave the original goal to a future study <cit.>.



 It is convenient to write the two bosons in terms of its magnitude and phase
ψ_α=√(ρ_α)e^iθ_α,
then

    Ω    =r(ρ_1+ρ_2)+U/2(ρ_1+ρ_2)^2
    	+A/2(ρ_1-ρ_2)^2              
       + 2Bρ_1^2ρ_2^2cos^2[2(θ_1-θ_2)]
    	-h(ρ_1-ρ_2)

which leads to two sets of saddle point solutions.

For h< h_c=-rA/U=ρ_0 A:

    ρ̅_1  =(ρ_0+h/A)/2,  ρ̅_2=(ρ_0-h/A)/2,    
       θ_1-θ_2  =π/4+mπ/2,    m=0,1,2,3

 where ρ_0=ρ̅_1+ρ̅_2= -r/U is the total density of the spinor atoms <cit.>.
 Note that at these saddle points, the effective mass term r= - ρ_0 U is independent of h.
 From Eq.<ref>, S^z=ρ̅_1-ρ̅_2=h/A=ρ_0 cosθ_0, namely, cosθ_0=h/h_c.
 So the saddle point Eq.<ref> has the spin-orbital  structure
 ρ_0((-1)^i_x+i_ysinθ_0,(-1)^i_x+i_ysinθ_0,cosθ_0) named as a canted AFM (CAFM) SF state (Fig.3a).
 In the h→ 0^+ limit, it reduces to the the N=2 XY-AFM SF state with θ_0=π/2.
 When h → h^-_c, the spin structure becomes n_0(0,0,1) which is fully aligned to the z direction, so named the Z-FM state (Fig.3a).

 For h > h_c:

    ρ̅_1=ρ_0,  ρ_2=0

  which, according to Eq.<ref>, is the Z-FM state with θ_0=0 in Fig.3a.
  Note that at this saddle points, the effective mass term r=h-ρ_0 (A+U)=(h-h_c)-ρ_0 U depends on  h
  and lim_h→ h_c^+r=-ρ_0 U which approaches to the same value as that at h ≤ h_c.

  At the mean field level, we reproduced the two different ground states  h< h_c and h>h_c respectively achieved by the microscopic approach used in the previous sections. In the following, we will explore the quantum fluctuations, especially the Quantum phase transition between the two.



7. The canted anti-ferromagnetic (CAFM) SF phase at h<h_c. 

Below the critical point h<h_c, the bosonic system is in canted anti-ferromagnetic superfluid phase Eq.<ref>.
Without of loss generality, we choose θ̅_1=0 and θ̅_2=π/4
and write fluctuations as
ψ_α=√(ρ̅_α+δρ_α)e^i(θ̅_α+δθ_α),
for the notation simplicity, we drop the δ in δθ_α,
then

    ℒ_CA   =iδρ_1∂_τθ_1+iδρ_2∂_τθ_2
    	+v^2ρ̅_1(∇θ_1)^2 +v^2ρ̅_2(∇θ_2)^2     
       +U/2(δρ_1+δρ_2)^2
    	+A/2(δρ_1-δρ_2)^2
    	+8Bρ̅_1^2ρ̅_2^2(θ_1-θ_2)^2   
        +v^2/4ρ̅_1(∇δρ_1)^2+v^2/4ρ̅_2(∇δρ_2)^2

 where we still keep the gradient term on the density fluctuations
 in the 3rd line which are higher orders to the mass terms
 of the density fluctuations in the second line <cit.>.

 Following bilayer quantum Hall systems <cit.>,
 we introduce ρ̅_±=ρ̅_1±ρ̅_2,
 δρ_±= δρ_1±δρ_2 and  θ_±=θ_1±θ_2, then

    ℒ_CA   =i/2δρ_+∂_τθ_+
    	+i/2δρ_-∂_τθ_-
    	+v^2 ρ_0/4[(∇θ_+)^2+(∇θ_-)^2]  
       +U/2(δρ_+)^2
    	+A/2(δρ_-)^2
    	+B/2(ρ^2_0-ρ̅_-^2)^2( θ_-)^2     
        +v^2ρ_0/4(ρ^2_0-ρ̅_-^2)
    	[(∇δρ_+)^2+(∇δρ_-)^2]                     
        +v^2ρ̅_-/2∇θ_+·∇θ_-
        -v^2ρ̅_-/2(ρ^2_0-ρ̅_-^2)∇δρ_+·∇δρ_-

  where the B terms leads to the mass term of θ_-, the last line
  is the coupling between the + and - mode <cit.>.
  In the absence of the Zeeman field h=0, ρ̅_-=0, then the + and - modes decouple,
  this corresponds to the balanced BQHE.


  The general expressions for the two eigen-modes are complicated. So we only list its long-wavelength limit:

    ω^+_CA(k)    =√(2ρ_0Uv^2k^2)
    ω^-_CA(k)    =√(Δ^2_R + v^2_R k^2 )

where ρ_0=ρ̅_1+ρ̅_2=-r/U, θ_0=arccos(h/h_c)
and

    Δ_R      =      2√(B/A^3)(h^2_c-h^2),   
    
     v^2_R       =      2ρ_0(A+Bρ_0sin^2θ_0)v^2

The + mode is nothing but the superfluid Goldstone mode.
While the - mode is the roton mode whose energy takes a relativistic form with a roton gap:
 vanishes  as Δ_R ∼ h_c-h when h → h- h^-_c.
 Obviously, it is the B term which opens the gap to the roton mode.


8.  The Z-FM SF phase at  h ≥ h_c.

Above the critical point h > h_c,
it is more convenient to introduce a mixed coordinate system
with ψ_1=√(ρ̅_1+δρ_1)e^i θ_1   in the polar coordinate and ψ_2
in the Cardisian coordinate. So (δρ_1, θ_1 ) and ψ_2 are all small, the Lagrangian density becomes

    ℒ_FM   =iδρ_1∂_τθ_1
    	+v^2 ρ_0(∇θ_1)^2
    	+m^2_1/2(δρ_1)^2 +v^2/4ρ_0(∇δρ_1)^2 
       +ψ_2^*∂_τψ_2
    	+v^2|∇ψ_2|^2+r_2|ψ_2|^2
    	+1/2(U+A)|ψ_2|^4    
       +(U-A)δρ_1|ψ_2|^2 +B(ρ_0+δρ_1)^2|ψ_2|^4  
       +B/2(ρ_0+δρ_1)^2
    	[(e^-iθ_1ψ_2)^4+(e^i θ_1ψ_2^*)^4]

where m^2_1=U+A 4 and r_2=r+h+(U-A)ρ_0=2(h-Aρ_0)=2(h-h_c).


 To study the universality class of the transition, it is sufficient  only to keep
 the lowest order coupling terms:

    ℒ_FM   =ψ_2^*∂_τψ_2+v^2|∇ψ_2|^2
    	+r_2|ψ_2|^2+u_2 |ψ_2|^4    
       + i Γδρ_1∂_τθ_1
    	+v^2ρ̅_1(∇θ_1)^2
    	+1/2m^2_1(δρ_1)^2	  
       + γ_0 δρ_1|ψ_2|^2 +γ_1
    	[e^-i4 θ_1ψ^4_2+e^i4 θ_1ψ^*4_2 ]

  where we dropped (∇δρ_1)^2 term which is irrelevant anyway <cit.>,
  u_2=1/2(U+A+ 2 Bρ̅^2_1 ), γ_0=U-A, γ_1=ρ̅^2_1 B/2.
  Γ=1 at the bare level, but was introduced to keep track of its flow under RG.
  The first line is the quantum critical mode of ψ_2 in the universality class of zero density
  superfluid to Mott (SF-Mott) transition, therefore has the critical exponents z=2, ν=1/2, η=0.
  The second line stands for the gapless SF mode of (δρ_1, θ_1 ).
  The third line is the coupling between the two sectors with the coupling parameters γ_0,γ_1.

  From the quadratic parts of  Eq.<ref>, one can

   extract the two eigen-modes
  in the long wavelength limit:

    ω^+_FM   =√(2ρ_0(A+U)v^2k^2)
    ω^-_FM   =2(h-h_c)+v^2k^2

which has the dynamic exponent z=2.


   Contrasting Eq.<ref> with Eq.<ref>, one can see that
   due to the extra A term in the latter, the SF velocity increases  from h ≪ h_c to h ≫  h_c,
   while the difference ( v^+_FM )^2- (v^+_CA)^2= 2 ρ_0 A v^2  k^2 maybe a good measure of the A term.
   However, the SF density ρ_s = ρ_0 v^2 remains  the same  as shown in Fig.1b.
   The roton - mode's gap vanishes as Δ_R ∼ |h-h_c| on both sides,
   However, the roton gap at h > h_c seems independent of the B term,
   while that at h < h_c is a good measure of the B term.



9. The nature of the quantum critical point at h = h_c and new dangerously marginally irrelevant operators 


By keeping both SF  ( δρ _1, θ_1) mode
and the critical ψ_2 mode in Eq.<ref>, one can perform RG flow
on the two couplings γ_0, γ_1 between the two sectors.

     The RG analysis can be done by following the procedures developed in <cit.> which study
     the SF to metal transition in a high temperature superconductor.
     The first line in Eq.<ref> stands for the critical mode with the dynamic exponent z=2,
     one can scaling dimension of the critical field [ ψ_2( x⃗, τ) ]= 1,
     then [u_2]=0 is marginally irrelevant at the upper critical dimensiona d_u=2.
     The second line stands for the SF mode with its own dynamic exponent z_SF=1.
     This is contrast to the case in <cit.> where both the boson and the fermion have the same dynamic exponent z=2.
     However, to be consistent with the ψ_1 mode, one must still choose z=1.
     Then one can determine the scaling dimensions of the magnitude fluctuation [ δρ _1( x⃗, τ) ]= 2,
     and the phase mode [θ_1( k⃗,ω ) ]=-3.
     Then the Berry phase term [Γ]=-1, so it becomes irrelevant, then the magnitude fluctuation δρ _1
     and the phase mode are asymptotically decoupled.
     Similar phenomena happens in the two channel Kondo model <cit.> where the
     Berry term also becomes irrelevant and turns into a boundary condition for the Majorana fermions !
     The third line which stands for the coupling between the two sectors,
     so can determine [ γ_0]=0, so it is marginal.
     Similar, one can also determine [ γ_1]=0, so it is also marginal.


1. The effective action 

     To construct an effective action to study the nature of the QCP,
     one can simply integrate out  ( δρ _1, θ_1) in Eq.<ref>
     to study how the SF phase fluctuation affect the quantum critical behaviors of the ψ_2 mode.
     Because the δρ _1 is always massive with the mass m_1, so can be integrated out to generate
     a |ψ_2|^4 term, so can be absorbed into the coefficient u_2.
     One can also see ⟨ e^-i4 θ_1⟩= ⟨ e^ i4 θ_1⟩=e^-8 ⟨θ^2_1 ⟩
     is a finite constant at d=2+1 dimension which stands for the suppression of the BEC condensate at T=0.
     So integrating out ( δρ_1, θ_1 ) leads to the effective action describing the transition driven by  the roton dropping:

    ℒ_QC   =ψ_2^*∂_τψ_2+v^2|∇ψ_2|^2 +r_2|ψ_2|^2+u_2 |ψ_2|^4    
       +γ_1[ψ^4_2+ψ^*4_2 ] + ⋯

   where only the last term breaks the U(1) symmetry to the underlying [C_4 × C_4]_D symmetry.
   However, it can be shown this symmetry breaking γ_1 term is marginally irrelevant
   at the z=2, ν=1/2, η=0 universality class at d=2. But it is still dangerously irrelevant
   at h < h_c, leads to the CAFM state in Fig.3a.


    If the Berry phase term in Eq.<ref> had became a second derivative  |∂_τψ_2|^2,
    then the dynamic exponent would be z=1,  then it is in the 3d XY universality class.
    Then the γ_1  maybe shown to be irrelevant near this 3d XY
    QCP by ϵ= 3-d expansion <cit.>.
    But here it is a Gaussian fixed point with z=2, we expect its is marginally irrelevant.
    Because both u_2 and γ_1 are marginally irrelevant, doing RG calculations, one must treat
      u_2 and γ_1 on the same footing. Because the γ_1 term breaks the  U(1) symmetry,
      it may lead to logarithmic corrections <cit.> different from
      those just due to the u_2 term shown in Eq.<ref>.



In Eq.<ref>,at r_2 < 0,  ψ_2 has a condensation
ψ_2=√(ρ̅_2+δρ_2)e^i θ_2,
then to quadratic order, the Lagrangian becomes

    ℒ_R   =iδρ_2∂_τθ_2
    	+v^2[1/4ρ̅_2(∇δρ_2)^2+ρ̅_2(∇θ_2 )^2]
       +u_2(δρ_2)^2+16γ_1ρ̅_2^2( θ_2)^2

 where we have set the linear term  (r_2+2u_2ρ̅_2)δρ_2 vanishing  to find ρ̅_2=-r_2/(2u_2).
 Then one can extract the eigen mode
ω_R=√((v^2k^2+4ρ̅_2u_2)(v^2k^2+16γ_1ρ̅_2))≈√(64γ_1ρ̅_2^2u_2+ρ̅_2(4u_2+16γ_1)v^2k^2)
which recovers the ω_- mode in Eq.<ref>.
Note that it is the γ_1 term which leads to the mass of the phase θ_2.

2. Quantum critical scalings and quantum information scramblings in the QC regime near h_c. 





The order parameter is the magnetization in the XY plane ⟨ S^x ⟩.
It vanishes in the Z-FM SF phase, but non-zero in the CAFM SF phase.
Because the γ_1 term breaks the U(1) symmetry, but is marginally irrelevant at d=2.
There is an emergent U(1) symmetry at the QCP which indicates the CAFM SF to the Z-FM SF transition








is in the same universality class of the zero-density SF-mott transition with the
critical exponents z=2, ν=1/2, η=0  studied in Ref.<cit.>.
However, the two phases are very much different from the SF and Mott.
In the CAFM SF, there is cos 4 ϕ clock term in Eq.<ref> which is dangerously irrelevant near the QCP:
it is irrelevant near the QCP, but controls the quantum phase h < h_c. It leads to the CAFM SF phase,
breaks the presumably U(1) symmetry to the joint C_4 symmetry and also opens a gap Eq.<ref> to the roton mode.
There is a melting transition T_M above the CAFM SF phase ( Fig.1b ).
The universality class of melting process belongs to
2D q=4 state clock transition.


The standard scaling shows the roton gap should scale as Δ_R ∼ |h-h_c |^z ν∼ |h-h_c|
which is consistent with our specific calculations shown above.
The specific heat should scale as  C_R ∼ T^d/z∼ T.
Our specific calculations show that the specific heat due to the roton in the spin sector is
C_R=π^3/3v^2T.
Applying the scaling result for the U(1) conserved quantity in <cit.> to the order parameter
leads to :

    ⟨1/N_s∑_i[(S_i^x)^2+(S_i^y)^2]⟩
    	=2mT/4π1/{ln[Λ^2/(2mT)]}^4

where m=1/(2v^2)=t_0/(2t_0^2-t_so^2)
and Λ is a momentum upper cutoff <cit.>.


Obviously, due to the non-linearities in Eq.<ref> and <ref>, the system shows quantum chaos <cit.>
in the QC regime. The quantum information scrambling encoded in the spinor boson out of time ordered correlation function (OTOC):

    ⟨ψ^†_α(t, x⃗ ) ψ^†_α(0, 0) ψ_α(t, x⃗ )
     ψ_α(0, 0) ⟩∼ e^λ_L,α(t- x/v_B,α),

where α=1,2 stand for SF component and Roton component in Eq.<ref> respectively.
Near the lightcone x=v_B,α t ( Fig.<ref> ), the OTOC
is greatly enhanced in the QC regime in Fig.1.
The Lyapunov exponent  λ_L,α due to the critical roton mode reaches its maximal value λ_L,R∼ T,
while the butterfly velocity v_B,R∼ T^1-1/z. For z=2,  v_B,R∼√(T).
From the dimensional analysis <cit.>, we conclude  v_B,R∼  v √(T)<cit.>.




Note that the SF sector stays un-critical across the QCP.




We obtain the finite KT transition temperature ( Fig.1b ) T_KT=π/2ρ_s∼π n_0v^2.
The Lyapunov exponent due to the SF mode <cit.> is λ_L,sf∼ T^3/ρ^2_s and its
butterfly velocity v_B,sf∼ c= v √( 2 n_0 U) as can be seen from Eq.<ref>,<ref>.

From Eq.<ref>, using Eq.<ref> and shifting the orders of the operators in Eq.<ref>, 
one may also evaluate the transverse spin-spin correlation function,
⟨ S^+(t, x⃗ ) S^-(0, 0)  ⟩ = (-1)^i-j⟨ψ^†_1(t, x⃗ ) ψ_2(t, x⃗ ) ψ^†_2(0, 0) 
 ψ_1(0, 0) ⟩. Note that it can only be evaluated from the original action Eq.<ref>
 instead of the effective action Eq.<ref> where the SF model ψ_1 has been integrated out.





10. Contrast the GL approach with  with the microscopic calculations in studying the NOFQD phenomena 

The two approaches are completely independent and complementary to each other. Both have its own
advantages and limitations.
 It is very constructive to contrast phenomenological approach in Sec.5-8  with
 the microscopic calculations in Sec.2-4.
 We can see the following mappings between the parameters in the phenomenological GL action Eq.<ref> which holds at any h and those
 evaluated by the microscopic calculations at h=0 in Sec.5 and Sec.B in the SM, especially the effective potential generated by the OFQD:

    v^2   =(t_0^2-t_so^2/2)/t_0,  r=-(μ+4t)
          
    
     A   = A_m n^2,   B= B_m n^4,

  where the subscript m indicates the parameters evaluated by the microscopic calculations at h=0 in Sec.2-5.

  The interaction U and the Zeeman field h are the same in both approaches.
  Especially, it indicates that the A and B in the GL action are nothing but the effective potential generated by
  OFQD upto a density dependence. So both are positive at h=0 in the weak coupling limit.
  However, their signs may change as the interaction increases from the weak to strong couplings <cit.>.

  From Eq.<ref>, one can see the roton spectrum matches that achieved by microscopic calculation at its valid regime:
  high h limit. For example, the gap Δ_R= 2(h-h_c)  matches
  the microscopic calculations in the high h limit in Sec.2. This fact may also be related to the mapping shown
  in the first row of Eq.<ref> that the GL parameter r in Eq.<ref>
  is pinned at r=- ρ_0 U when h ≤ h_c, but increases linearly r=(h-h_c)- ρ_0 U matches
  the chemical potential -(μ+4t)=h-n_0 U in the microscopic calculations in Sec.2.


  The density dependence of A and B  in Eq.<ref> is important in the evaluations of
  excitation spectrum at h >h_c and h<h_c, especially the nature of the QCP. Had one ignored the quantum fluctuations
  in the density channel ( for example, by fixing it to be the BEC condensate n=n_0 ),
  one would have found the Zeeman field only couples to the spin sector, therefore leads to
  the un-physical physical picture: the SF sector is completely decoupled from the spin sector.
  Of course, A and B also depend on h, so the mapping Eq.<ref> only holds at h=0.


  Most importantly, for the very first time, we establish the connections between the GL action and the OFQD evaluated by
  microscopic calculations !
  This observation also endows the physical meaning of the OFQD from a complete new and deep  perspective.


11. Analogy and difference from the bilayer quantum Hall effects at the total filling factor ν_T=1. 

It is very instructive to compare the Bilayer quantum Hall effects (BLQH) <cit.>
at the total filling factor ν_T=1: Eq.<ref> is similar to the effective Chern-Simon Ginsburg-Landau (CSGL) composite boson
theory of the BLQH: in the absence of B term, the symmetry is U(1)_1 × U(1)_2 = U(1)_+ × U(1)_-,
while the B term breaks the U(1)_- symmetry and opens a gap to the putative Goldstone mode and transfers it into
a pseudo-Goldstone mode. The U term maps to the intra-layer Coulomb iteration,
the A term maps to the inter-layer capacitive term at a finite distance d.
However, there are also several  crucial differences:
(1) Here, both + and - sector are charge neutral. While in BLQH,  only - is the charge neutral,
there is a Chern-Simon term in the + charge sector
which opens a charge gap to the putative SF Goldstone mode, leads to the QH with the quantized Hall conductivity σ_H=1
in the charge sector and the topological charges of the merons.
(2) Here, as shown in Eq.<ref> at the microscopic level,  the A term and B term are the effective potentials  generated by OFQD
in Eq.<ref> which depend on the t/U. However, in BLQH,
the capacitive A term  is due to the difference between intra- and inter-layers Coulomb interactions tunable by the
distance d  of the two layers.
The B term is similar to the inter-layer tunneling term H_t=- t (ψ_1^*ψ_2 + h.c ) which can also
be manipulated by the in-plane magnetic field.
(3) Here the " tunneling term " between the two species ψ_1, ψ_2  is
 the 4th power dictated by the underlying [C_4 × C_4]_D  symmetry. However, it is just 1st power in BLQH.
 The power difference in the interlayer tunneling makes  dramatic differences in the two systems.
 Here, the 4th power was shown to be marginally and dangerously irrelevant near the QCP in Fig.1a, but opens the gap to the roton mode.
 However, the 1st power was known to be relevant in BLQH,, opens a gap to the interlayer SF Goldstone mode, essentially destroys the
 interlayer phase coherence.
(4) Here, there is a roton dropping transition tuned by the Zeeman field h when either ψ_1 or ψ_2
is depleted and the spin is polarized along the z direction in Fig.1a. In principle, in the absence of
the inter-layer tunneling term H_t, the ν_T=1 BLQH  is independent of the bilayer imbalance,
there is always a neutral gapless mode.
As the distance d increases, there maybe a roton dropping transition in the charge neutral sector
leading to a first order transition from the exciton superfluid (ESF) to pseudo-spin density wave (PSDW)<cit.>.
Both the roton mode at a finite wavevector |k⃗ |=k_0 and SF mode near the zero wavevector k⃗=0 happen in the - mode,
the origin of the roton mode is due to the inter-layer interaction, so very much different from the roton mode here tuned by the
Zeeman field in Fig.3. The main difference is that the BLQH is in a a continuous system, so  the roton is at a 1d roton ring
driving a 1st order transition, while here is on a square lattice, so the roton is at (π, π)  driving a 2nd order transition.
(5) Here, we did the microscopic calculations by the Bogoliubov calculations at any Zeeman field h
    which can be contrasted with the phenomenological GL action Eq.<ref>. While, in the BLQH, there is a also a microscopic
    calculation called Hartree-Fock (HF) projected within the Lowest Landau level (LLL) which was also compared to
    the phenomenological Chern-Simon GL effective action <cit.>.

















12. Detectable experimental implications in cold atoms.















In the ongoing experiment <cit.>,
the BEC has N∼ 3×10^5 atoms
and trapped within the diameter d = 80μm,
so one lattice site has about n = 10 atoms.
For typical experiment parameters,
optical lattice potential is 4.16E_r,
and Raman potential is 1.32E_r,
where E_r=375 nK denotes the recoil energy,
so the tight binding model parameter can be estimated
as t_0∼100 nK, t_so∼ 30 nK, so t_so/t_0 ∼ 1/3.
The short-range Hubbard like interaction
U=4πħ^2 a_s/m∫ d^2r|w(r)|^4,
where s-wave scattering length of the ^87Rb atoms a_s=103a_0
and a_0 is the Bohr radius and the mass of the bosons ^87Rb
lead to U∼ 10nK, so n_0 U ∼  100 nK.

Based on these experimental parameters,

one can estimate T_KT∼ 100 nK, the roton gap away from the critical point in Fig.3b Δ_R ∼ 1 nK,
so the meting transition T_M ∼ 1 nK and the critical Zeeman field h_c ∼ 1 mG.
So far, the experiments are operating at T ∼ 20nK.
As noted earlier,  there is also a specific heat contribution from the superfluid sector
C_sf=π/n_0Uv^2T^2 ∼ (T/n_0U) T/v^2 ≪  T/v^2 ∼ C_R
which is sub-leading to that from the critical roton mode in the spin sector.
Of course, the SF sector makes no contributions to the magnetization  Eq.<ref>.
In terms of quantum information scramblings, the Lyapunov exponent due to the SF mode λ_L,sf∼ T^3/ρ^2_s ∼  T (T/ρ_s)^2∼  T (T/T_KT)^2 ≪ T ∼λ_L,R, while its butterfly velocity v_B,sf∼ c= v √( 2 n_0 U)≫ v_B,R=v √(T).
This means the light cone of the SF mode x_sf= v_B,sf t is well within that of the roton mode x_R= v_B,R t ( Fig.<ref> ).
Due to the wide separation of the two light-cones, both light-cones could be detected in the  spinor boson ( see Eq.<ref> ), spin or density OTOC.

So we conclude that the current experimental operating temperatures  may not be able to get into the CAFM SF phase,
but are still low enough such that the dramatically quantum critical fluctuation effects in
the specific heat, Lyapunov exponent, butterfly velocity,
transverse magnetization, boson or spin correlation functions can be detected


by dynamic or elastic, energy or momentum resolved,
longitudinal or transverse Bragg spectroscopies
<cit.>,
specific heat measurements <cit.>
and in-situ measurements <cit.>.




13. Conclusions and discussions 



  In this work, we propose a new concept: NOFQD. We also perform both microscopic and
  effective GL action to implement this new concept in the specific QAH Hamiltonian Eq.1.
   We did microscopic calculations at both weak and
   strong couplings at both h=0 and h ≫  U. Then we we construct the effective GL action
   at any h  and contrast with the  microscopic calculations at both h=0 and very strong Zeeman field h ≫ U.
   However, only the GL action can capture the competition near h ∼ h_c.

In order to clarify the connections and relations between the conventional OFQD and the NOFQD studied here,
according to its response to a given deformation which breaks either an exact symmetry or a spurious symmetry or both,
we classify OFQD into two classes:

1. Type-I OFQD response trivially to a deformation: 

In type-I, at a high symmetry point, there is an exact symmetry and also a spurious continuous symmetry,
the OFQD breaks the spurious symmetry down to the exact symmetry and selects
the quantum ground state upto the exact symmetry.
Then a deformation λ away from the high symmetric point breaks the exact symmetry explicitly,
but may or may not break the spurious symmetry,
it just pick up a state from the exact degenerate manifold. If the spurious symmetry remains at λ≠ 0,
then the OFQD analysis at λ≠ 0 is very similar to that at λ=0 ( Fig.<ref>a ).



2.  Type-II OFQD response non-trivially to a deformation and leads to NOFQD 





In type-II, at a high symmetry point, there is an exact symmetry and also a spurious continuous symmetry.
The OFQD breaks the spurious symmetry down to the exact symmetry and selects
the quantum ground states upto the exact symmetry.
Then a deformation λ away from the high symmetric point keeps the exact symmetry, but breaks the spurious symmetry.
There is an competition between the effective potential generated by the OFQD at λ=0 and
the deformation λ as λ evolves,
which drives a QPT due to the spontaneously symmetry breaking of the exact symmetry.
The QPT happens when λ matches the scale of the effective potential.
This is what we called NOFQD: despite there is no spurious symmetry anymore when λ≠ 0,
but the effective potential generated at λ = 0 still has a dramatic effect when λ≠ 0  ( Fig.<ref>b ).



It maybe constructive to compare this classification with that of superconductors in the response to an orbital magnetic field:
in Type-I, there is  a direct transition from the normal state to a superconducting state at just one critical field h_c,
in Type-II, there is a mixed state intervening between the low h_c1 and high critical field h_c2.



3. Two examples of Type I OFQD under SOC 

In a recent work on complete/In-complete devil stair cases in the strong coupling limit of SOC <cit.>,
there is an exact [C_4 × C_4]_D symmetry and also a spurious U(1) symmetry,
along the diagonal line α=β near the Abelian point α=β in the SOC parameter space.
The OFQD analysis selects either Y-x or X-y state as the quantum ground state which rotates to each other under
the exact [C_4 × C_4]_D symmetry. However, slightly away from the diagonal line α≠β, then
it  breaks the exact [C_4 × C_4]_D symmetry explicitly,
the spurious U(1) symmetry also disappears. Then
a simple classical analysis picks up  Y-x or X-y state when α > β
or α < β respectively. Therefore, the OFQD at α =  β response trivially to the  deformation α≠β.


For another example, for the (α=π/2, β ) SOC case
which has an exact U(1)_soc symmetry at the spin anisotropic interaction λ=1.
There is also a spurious U(1) symmetry at λ=1.
At λ=1,  the OFQD at the weak coupling selects the PW-X phase upto
the exact U(1)_soc symmetry to be the quantum ground state.
For example, the Z-x phase and the PW-X transform to each other under the exact U(1)_soc symmetry.
When λ < 1, it explicitly breaks the exact U(1)_soc symmetry,
but the spurious U(1) symmetry  remains,
then a similar OFQD analysis selects the PW-X phase to be the quantum ground state.
When λ > 1, it also explicitly breaks the exact U(1)_soc symmetry, the spurious U(1) symmetry also disappears,
then a simple classical analysis alone picks up the Z-x state to be the ground state.
So the OFQD at λ=1 response trivially to the deformation λ≠ 1.

Both examples belong to Type-I, the deformation breaks the exact symmetry at the high symmetric point, but it may or may not break
the spurious symmetry. There is NO competition between the effective potential generated by the OFQD
at the high symmetric point and the deformation away from the high symmetric point, no associated QPT either,
therefore no NOFQD phenomena.

4. Two examples of Type II OFQD under SOC 

   For the QAH studied in this work,  at h=0, it has the M symmetry and the exact [C_4 × C_4]_D symmetry.
   At h=0 there is also a spurious SU(2)_s symmetry, an OFQD analysis selects the N=2 XY-AFM SF as the quantum ground state
   and also opens  a gap to the roton mode. This state has the d=4 degeneracy,  keeps the M symmetry, but breaks the
   the exact [C_4 × C_4]_D → 1 symmetry spontaneously.
   Any small h ≠ 0  breaks the M symmetry, but still keeps the exact [C_4 × C_4]_D symmetry.
   The spurious  SU(2)_s symmetry disappears, a simple classical analysis picks up the Z-FM SF which keeps the exact [C_4 × C_4]_D symmetry. Then there is an associated QPT from the N=2 XY-AFM SF state to the Z-FM as h increases ( Fig.1 ).
   The exact [C_4 × C_4]_D symmetry was spontaneously broken in the N=2 XY-AFM SF state, but restored in the Z-FM  state.
   Of course, one need to check if the 2nd QPT be pre-emptied by other competing phases.
   By both microscopic and effective GL action, we show that this does not happen .
   Despite there are many known examples of TYPE-I, we expect Type-II is also transformative to many other systems.
   For example, in a recent unpublished work <cit.>, we show that Type-II also appears in the interacting bosonic Haldane model in a honeycomb lattice subject
   to a sublattice potential.

4. Contrast to NAdS_2/NCFT_1 in SYK models 

It is very instructive to compare the NOFQD with the NAdS_2/NCFT_1 correspondence in the SYK model.
In the latter, one must move slightly away from the 1d CFT to take into accounts
the leading irrelevant operator -i ω to study the OTOC and  extract its quantum Lyapunov exponent.
Here, a small h is also slightly away from h=0 where the OFQD selects a quantum ground state such as the N=2 XY-AFM SF state.
In the RG sense, due to the small roton gap, a small h is  also irrelevant.
However, the main difference is that in the SYK models, the ground state is a gapless QSL  state which breaks the 1d CFT ( reparamatrization ) invariance spontaneously to SL(2,R) leading to a zero ( Goldstone ) mode, the irrelevant deformation
-i ω breaks the 1d CFT ( reparamatrization ) invariance explicitly, lifts the zero mode to a pseudo-Goldstone mode
leading to maximal chaos, but does not drive any QPT.
While here the N=2 XY-AFM SF state breaks the the exact [C_4 × C_4]_D symmetry spontaneously and has a roton gap,
the Zeeman field h keeps the exact [C_4 × C_4]_D symmetry, it reduces the roton gap and
just tilts the N=2 XY-AFM state into the CAFM state which has the same symmetry
as the N=2 XY-AFM SF state. However, a sufficiently large h does drive a transition to the Z-FM SF. Of course, NAdS_2/NCFT_1 correspondence is only for 1d, at higher dimension
such as NAdS_d+1/NCFT_d with d ≥ 2, there is no need
to move slightly away from CFT, so nearly is not necessary. Here, although we demonstrate it in only 2+1 dimension, the NOFQD happens in any dimension.

As stressed in Sec.5, the original goal of constructing the GL action is to study the transition from weak coupling to strong
coupling at an integer filling. Here, we use it to investigate the transition driven by the Zeeman field h
within the weak coupling regime, so the system
remains inside a SF phase. We leave our original goal to a future work <cit.>.
In a future work, directly inspired by the experimental realizations, we look at the deformation due to t_sx≠ t_sy which breaks
   the exact [C_4 × C_4]_D symmetry explicitly to [C_2 × C_2]_D.  We will do  microscopic calculations at both weak and
   strong couplings.





Acknowledgements 


We thank Bao-Zong Wang and Zhan Wu for helpful discussions on experimental detections.
We also thank B. Halperin for the critical discussions in <cit.> and some other helpful comments on the manuscript.
We acknowledge AFOSR FA9550-16-1-0412 for supports.


99

kane
M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators,  Rev. Mod. Phys. 82, 3045 (2010).

zhang
X. L. Qi and S. C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).

QAHthe Rui Yu, Wei Zhang, Hai-Jun Zhang, Shou-Cheng Zhang, Xi Dai1, Zhong Fang,
Quantized Anomalous Hall Effect in Magnetic Topological Insulators, Science 329, 61 (2010).

QAHexp Cui-Zu Chang, Jinsong Zhang, Xiao Feng, Jie Shen, Zuocheng Zhang, Minghua Guo1, Kang Li, Yunbo Ou, Pang Wei, Li-Li Wang, Zhong-Qing Ji, Yang Feng1, Shuaihua Ji, Xi Chen, Jinfeng Jia1, Xi Dai2, Zhong Fang2, Shou-Cheng Zhang3, Ke He2,?, Yayu Wang1,?, Li Lu2, Xu-Cun Ma2, Qi-Kun Xue1, Experimental Observation of the Quantum Anomalous Hall Effect
    in a Magnetic Topological Insulator, Science 340, 167 (2013).



expk40 Lianghui Huang, et.al,
Experimental realization of a two-dimensional synthetic spin-orbit coupling in ultracold Fermi gases, Nature Physics 12, 540-544 (2016).

expk40zeeman Zengming Meng, et.al,
Experimental observation of topological band gap opening in ultracold Fermi gases with two-dimensional spin-orbit coupling, arXiv:1511.08492.

2dsocbec
Zhan Wu, Long Zhang, Wei Sun, Xiao-Tian Xu, Bao-Zong Wang,
Si-Cong Ji, Youjin Deng, Shuai Chen, Xiong-Jun Liu, Jian-Wei Pan,
Realization of two-Dimensional spin-orbit coupling for Bose-Einstein condensates,
Science 354, 83 (2016).

origin  The end of the boundary sits at the origin t_so=0, h=0 which
	has the exact SU(2) symmetry, the other parts only have the emergent U(1) symmetry.
	The whole line has the dynamic exponent z=2. As demeonstrated in the text,
    it has the same universality class as the zero density SF-Mott transition subject to a new
    marginally and dengerously irrelevant anisotropic term.

our
Fadi Sun, Junsen Wang, Jinwu Ye and Youjin Deng,
Symmetry protected bosonic topological phase transitions: Quantum Anomalous Hall system of weakly interacting spinor bosons in a square lattice,
arXiv:1711.11580, substantially revised version as version 2.


RBZ This is because there are 2 valleies in Eq.<ref>, 2 spins and the particle-hole spaces leading to the 8 × 8
Bogoliubov transformation and the Reduced BZ in Eq.<ref>.

Stoof2001
D. van Oosten, P. van der Straten, and H. T. C. Stoof,
Quantum phases in an optical lattice,
Phys. Rev. A 63, 053601 (2001).

sm
See the supplementary materials for more details to establish the effective potential
Eq.<ref> and <ref> generated by the
order from quantum disorder.
The analytical approach is based on a perturbation theory in t_so/t_0,
and the numerical approach is a direct numerical evaluation which confirms the analytical calculations on
ω_l,𝐤 and E_ofd(θ,ϕ) defined in Eq.<ref>.











bert We thank  B. Halperin for alerting us this fact  which inspired us to discover
the nearly order from quantum disorder (NOFQD) phenomena and the associated quantum phase transition, then
develop both the microscopic and phenomelogical approaches to study the NOFQD phenomena.






Sachdev
Subir Sachdev,
Quantum Phase Transitions, 2nd ed.
(Cambridge University Press, Cambridge, 2011).


z2
Subir Sachdev, T. Senthil, and R. Shankar,
Finite-temperature properties of quantum antiferromagnets in a uniform magnetic field in one and two dimensions,
Phys. Rev. B 50, 258 (1994).



RGflow   For notational simplicity, dropping the subscript 2 and writing ψ in terms of its real and imaginary part
   ψ=ϕ_1+ i ϕ_2, then the interaction terms in Eq.<ref> can be written as
   (u-6γ) ( ϕ^2_1 + ϕ^2_2)^2 + 8 γ ( ϕ^4_1 + ϕ^4_2) where the second term
   is the cubic anisotropic term with N=2 real components. Its RG flow at ϵ=3-d > 0 was well known at any N component
   with 4 RG fixed points (1) Gaussian fixed point
   u^*=v^*=0 (2) Ising fixed point u^*=0, v^* > 0 (3) O(N) fixed point u^* > 0, v^* = 0
   and (4)  Cubic fixed point u^* > 0, v^*∼ (N-4) ϵ≠ 0.
   The fixed point structure changes from N > 4 where the cubic fixed point u^* > 0, v^* > 0 is the only stable fixed point
   to N < 4 where the O(N) fixed point is the only stable fixed point.
   So when N=2 < 4, it belongs to the latter case where the cubic anisotropy is irrelevant near the 3d XY fixed point.
   At the upper critical dimension ϵ =0, all the 4 fixed points  converge to the  Gaussian fixed point.
   It is interesting to repeat this calculation at d=2, z=2 where the upper critical dimension d=2 and ϵ=0.
   For example, we will  put u_2 and γ_1 in Eq.<ref> at the same footing to perform the RG analysis, also evaluate the new
   logarithmic corrections to Eq.<ref> due to the two marginally irrelevant couplings.
   In fact, just like the transverse magnitization in Eq.<ref>, the specific heat C_R, the Lyapunov
   exponent λ_L, the butterfly velocity v_B,R will all receive a logarithmic factor at a finite T.




devil Fadi Sun and Jinwu Ye, Fractals, complete and in-complete devil staircases
of strongly interacting spin-orbit coupled bosons in a lattice,  arXiv:1603.00451, substantially revised version.








SY S. Sachdev and J. Ye, Gapless spin liuid ground state in a random quantum Heisenberg magnet,"
Phys. Rev. Lett. 70, 3339 (1993).

Kit A. Y. Kitaev, Talks at KITP, University of California, Santa Barbara," Entanglement in Strongly-
Correlated Quantum Matter (2015).

sun   Fadi Sun and Jinwu Ye,  Periodic Table of SYK and supersymmetric SYK, Phys. Rev. Lett. 124, 244101 (2020).


on Debanjan Chowdhury, Brian Swingle, Onset of many-body chaos in the O(N) model, Phys. Rev. D 96, 065005 (2017).



becbragg  Si-Cong Ji, Long Zhang, Xiao-Tian Xu, Zhan Wu, Youjin Deng, Shuai Chen, and Jian-Wei Pan, Softening of Roton and Phonon Modes in Bose-Einstein Condensate with Spin-Orbit Coupling, Phys. Rev. Lett. 114, 105301
(2015).


bragg1  Jinwu Ye, J.M. Zhang, W.M. Liu, K.Y. Zhang, Yan Li, W.P. Zhang, Light scattering detection of various quantum phases of ultracold atoms in optical lattices, Phys. Rev. A 83, 051604 (R) (2011).


bragg2 Jinwu Ye, K.Y. Zhang, Yan Li, Yan Chen and W.P. Zhang, Optical Bragg, atom Bragg and  cavity QED detections of quantum phases and excitation spectra of ultracold atoms in bipartite and frustrated optical lattices,
   Ann. Phys. 328 (2013), 103-138.


heat1
Kinast, J. et al.
Heat Capacity of a Strongly Interacting Fermi Gas.
Science 307, 1296 (2005).

heat2
Ku, M. J. H. et al.
Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas.
Science 335, 563 (2012).

dosexp Gemelke, N., Zhang X., Huang C. L., and Chin, C.
In situ observation of incompressible Mott-insulating domains in ultracold atomic gases,
Nature (London) 460, 995 (2009).

blqh1 Jinwu Ye, Fractional charges and quantum phase transitions in imbalanced bilayer quantum Hall systems,
        Phys. Rev. Lett. 97, 236803 (2006).

blqh2  Jinwu Ye and Longhua Jiang, Excitonic  superfluid to pseudo-spin density wave transition in bilayer quantum Hall systems,
        Phys. Rev. Lett. 98, 236802 (2007).


blqh3  Jinwu Ye, Mutual Composite Fermion and Composite Boson approaches to balanced and im-balanced bilayer quantum Hall
       systems: an electronic analogy of Helium 4 system,  Annals of Physics, 323 (2008), 580-630.

strongcoupling Fadi Sun and Jinwu Ye,   Strong coupling limit of the bosonic quantum anomalous Hall, quantum compass model in a transverse field and weak coupling to strong coupling transition, preprint in preparation.

anotherwriting More straightforwardly and intuitively, just like in BLQH <cit.>, one can write the interaction term as
 U/2[ |ψ_1|^2+|ψ_2|^2 - ρ_0 ]^2 where ρ_0 is the average total atom density
 ( or electron density  in BLQH ). Expanding it out automatically leads to r=-ρ_0 U.

blqhdrop  In the BLQH <cit.>, the two terms in the gradient of  δρ_+ and δρ_- were dropped due to
their subleading to the corresponding mass terms.

pmmode In the absence of Zeeman field h=0, the ± modes in Eq.<ref> are decoupled,
so the ± mode stand for the gapless SF mode and gapped roton mode in ± channels respectively.
However, when h > 0, there is a mixing between the two channels whci lead to two new eigen-modes which are neither
+ or - mode.  But we still use the same symbols.

sfmetal J. Ye and S. Sachdev, Superconducting,metallic,and insulating phases in a model of CuO_2 layers,
  Phys.Rev.B 44, 10173 (1991).

kondoye123452
Jinwu Ye,
Abelian Bosonization approach to quantum impurity problems,
Phys. Rev. Lett. 79, 1385 (1997).

dimen In the whole manuscript, we take the lattice constant a=1, so ka and x/a are  dimensionless, then
the butterfly velocity v_B carry the dimension of energy.

haldane  Fadi Sun et.al , NOFQD in bosonic Haldane model under a sublattice potential and its experimental observation, unpublished.








Supplementary Materials for“Nearly order from quantum disorder phenomena:
its application and detection in the bosonic quantum anomalous Hall system”
















In the main text, we present an NOFQD analysis on
the weakly coupling bosonic Quantum Anomalous Hall model leading to the phase diagram in Fig.M1b.
In this Supplementary Materials (SM),
in the first section, we present a traditional approach which ignores the NOFQD phenomena  and leads to a first order transition at h=0.
The breakdown of this traditional approach at h ≲ U was used to motivate the importance of the NOFQD phenomena
and also inspire us to develop the
new symmetry based GL approach in the main text to study the  NOFQD phenomena in the most systematic and complete way.



In section two, we provide the explicit form of the
effective potential Eq.M15  generated by the OFQD analysis at h=0
 which was also predicted by the independent symmetry based GL approach.

,
 more details on deriving effective action in Sec. S2,
explanation of different coordinate systems in Sec. S3,
and discuss possible issues in Sec. S4.







§  A. THE CONVENTIONAL APPROACH WHICH IGNORES THE NOFQD PHENOMENA LEADS TO A FIRST ORDER TRANSITION AT H=0.  



In the weak coupling limit, the conventional approach is to treat the kinetic ( single particle ) term exactly, then
treat the weak interaction term by perturbation theory. In the opposite strong coupling limit, one may treat interaction exactly, then
treat the kinetic term by perturbation theory <cit.>. However, in the intermediate regime, usually no controlled analytical approach can be
found.

In this manuscript, we are focusing on the weak interaction U < t_0 regime, which
is also the experimentally relevant situation. So conventional approach can be used to study  the Z-FM in the h ≫ U limit.
However, when h becomes also weak, namely,
in the regime h ≲ U, it becomes quite tricky and difficult to study the delicate balance between the
weak interaction effects and the weak Zeeman effects.
Here, we show that if one still takes the conventional ( or traditional ) approach which treats the kinetic ( single particle ) term first,
then the interaction by perturbation, one reaches that there is a first order transition line
at h=0 ( Fig.S1a ). As shown in the main text, this conclusion is in-correct. This is because
the conventional approach fails to capture the NOFQD effect which
treats the kinetic term and interaction effects on the equal footing, so
completely breaks down in this regime. A new and systematic approach  is needed to capture the NOFQD effect.

Such a new approach is the GL action  developed in the main text.
As alerted in the first paragraph, in a general many body interacting system, it is very difficult to construct
a controlled analytic approach  when the two competing scales such as the kinetic energy and the interaction are comparable.
Here, we make two crucial observations (1) even in the presence of h, the quantum ground state can still be chosen from the manifold spanned by the family of states shown in Eq.M2 parameterized by Eq.M4.
(2) We managed to find an analytic formula for the effective potential generated by the order from disorder phenomena Eq.M5
which holds in the whole range ( 0< θ < π, 0 < ϕ < 2 π).
 ( For the explicit derivation of this formula, see the next section of the SM ).
 This complete formula is important, because the minimum (θ_0, ϕ_0 ) changes as h increases.
 In fact, this complete formula was also predicted  independently from the symmetry based GL action.

It is these two facts which enable us to democratically and systematically
capture the NOFQD effect which stands for the
delicate competition between the effective potential generated by the OFQD due to the weak interaction
and the Zeeman energy due to a small finite h, therefore map out the global phase diagram Fig.M1.
Especially, the first order transition at h=0 is transferred to the
two second order transitions at h= ± h_c,
then a "new" phase called canted AFM phase ( CAFM ) emerges in the narrow window of 0 < h < h_c
( it still has the same symmetry as the N=2 XY-AFM phase ).
Obviously, there is a qualitative change from first order transition reached by the conventional approach
whcih ignores the NOFQD effects and the
second order transition reached by the new approach which incorporates the NOFQD effects.
It is the second order transition associated with the  NOFQD  which leads
to all the sharp experimentally detectable  predictions presented in the main text.

Although this conventional approach breaks down and lead to in-correct physics at h ≲ U, especially  a first order transition at h=0,
it is still instructive to present this traditional approach in this SM, because
it was used to inspire the authors to discover the NOFQD phenomena and to develop the new and systematic GL action approach in the main text to explore the NOFQD phenomena <cit.>.
It may also help to  shed some lights on the physical implications of NOFQD phenomena on the phase diagram in Fig.M1.

(a) The conventional approach to find the Z-FM at any h >0. 

Any small h ≠ 0 breaks the mirror symmetry and
split the degenerate single-particle state
at momentum (0,0) and (π,π).
A straightforward mean field analysis predicts the condensation occurs at
either (0,0) or (π,π) depending on the sign of h.
the unique lowest energy single-particle state lead to

|Ψ_sf,↑⟩∼(∑_i a_i↑^†)^N|0⟩ if h>0
or |Ψ_sf,↓⟩∼[∑_i (-1)^i_x+i_ya_i↓^†]^N|0⟩ if h<0.
Due to their spin alinement along the z axis, we name it Z-FM superfluid.

Then after choosing  (0,0) as the unique BEC condensation point at h >0, we
apply  the Bogoliubov theory to study the effects of the interaction U at any h ≠ 0.
As shown in Sec.M2,
we obtain a linear gapless superfluid mode and a gapped roton mode with the gap at Q⃗= (π,π) listed in
Eq.<ref>

    Δ^0_R(h)=2|h|

 shown in Fig.S1b (  also shown in the dash line in Fig. M4b in the main text ).
 As one decreases h from h ≫ U to h ≲ U, then to 0^+,
the roton mode at Q⃗= (π,π) gets lower and lower, then touches zero at h=0^+.
This behaviour may signify a possible first order transition at h_c=0^+.
This physical picture leads to the naive phase diagram shown in Fig.S1a.

As shown in the main text, this conventional approach is valid at h > U, but breaks down when h ≲ U.
Then the next two steps

(b) The order from quantum disorder analysis to find the N=2 XY-AFM state at h=0  
and

(c) The Roton mass gap generated from the "order by disorder" mechanism in the N=2 XY-AFM state st h=0. 
were performed in Sec.3 and Sec.4 in the main text respectively.

(d) A naive ( also wrong ) Roton mass gap at h> 0 which ignores the NOFQD effects. 

  In the conventional approach taken in (a), we treat h exactly, then treat U by perturbation. So
  we may expect that the roton gap at any h ≠ 0 could be obtained simply by
  incorporating the roton gap Eq.<ref> at h=0  directly into the roton gap Eq.<ref> at h ≠ 0:

    Δ_R(h)=√(Δ^2_R( h=0 )+ 4h^2 )

  which is always above the dashed line and shown in the solid line in Fig.S1b.



  Eq.<ref> leads to the following scenario:
  The transition at h=0 is a first order phase transition in Fig.S1a driven by the roton dropping tuned by the Zeeman field h shown in the solid line in Fig.S1b.
  The roton still has a non-zero gap Δ_R( h=0 ) in Eq.<ref> before it sparks the first order transition.



 Now the contradiction comes: if it were indeed a first order transition at h=0, then the state would be a mixed state
  of the (0,0) spin up state and (π,π) spin down state
  with any ratio. However, by the OFQD phenomena at h=0, we found the state is a pure state: the N=2 XY-AFM state
  instead of a mixed state. This contradiction shows the conventional approach may break down when h → 0^+. 









Indeed, as shown in the main text, this naive scenario holds only when h is sufficiently large h > U,
but breaks down at a small h ≲ U. As explicitly demonstrated in the main text, to compute the roton spectrum at
h ≲ U, one must examine the NOFQD effects which
treat the effects of U and h on equally footing.
Namely, by treating the two nearly degenerate states at (0,0 ) and (π,π) on equal footing, one can
investigate the competitions between the effective potential  generated by the OFQD in Fig.M2 due to the
interaction U and the Zeeman energy due to a small finite h.
The systematic GL effective calculation developed in the main text precisely capture the NOFQD effects and
leads to the correct roton spectrum
in the second equation in Eq.M18 where one can extract the correct roton gap at any h< h_c:

    Δ^-_R=4ABsin^4θ_0 =2√(B/A^3)(A^2-n_0^2h^2)/n_0

 which is listed in Eq.<ref> and shown in the solid line in Fig.M4b.


 In the h→ h^-_c limit, the roton gap vanishes as Δ_R∼ h_c-h.
 Note that the new  and systematic GL action approach not only leads to the correct roton gap Eq.<ref>, but also
 the correct roton spectrum at any k in the long wavelength limit in
 the second equation in Eq.M18 ( the solid line in Fig.M3b) which is qualitatively different from
 the naively speculated roton spectrum Eq.<ref> ( the solid line in Fig.S1b ).


Therefore, the breakdown of the conventional approach at h ≲ U leads to a qualitative change of the physical picture:
The first order transition at h=0 splits into two second order transitions at
a finite h = ± h_c shown in Fig.M1a. The N=2 XY-AFM SF phase at h=0 turns into the canted AFM SF phase at h >0.
Of course, the two phases still have the same symmetry breaking patterns. It is the zero temperature 2nd order transition
associated with the NOFQD effects which leads to all the sharp predictions on
all the experimental measurable quantities analyzed in the main text.

In short, we demonstrate that the effective potential generated from the OFQD
at least has 3 important impacts in the bosonic QAH system:
(a) at h=0, it leads to a quantum ground state selection rule,
(b) it also opens a gap to the roton mode at h=0
(c) Most importantly, at small h ≲ U, it leads to the NOFQD effects.
Then  we develop a new  and systematic GL effective action scheme to capture the NOFQD effects which stands for the competition
between the effective potential generated by the "order from disorder" mechanism with the Zeeman energy due to a small finite h.
The competition leads to the corrected roton excitation spectra at any k in the long wavelength limit in the second equation in Eq.M18.
It is the roton gap closing at h=h_c which
drives a second order quantum phase transition (QPT) from the CAFM to the Z-FM in Fig.M1b.
The QPT at T=0 associated with the NOFQD effects leads to sharp and unique predictions  in several
physical quantities which can be easily detected in the current experiments.

















§  B. THE FUNCTIONAL FORM OF THE EFFECTIVE POTENTIAL GENERATED BY THE OFQD 


  Here, we will explicitly demonstrate the explicit functional form Eq.M5 which holds in the whole range 0< θ < π, 0 < ϕ < 2 π).
  This complete formula is important, because the minimum (θ_0, ϕ_0 ) changes as h increases.
It is this formula which enable us to democratically and systematically
capture the NOFQD effects, therefore map out the global phase diagram Fig.M1.
Most importantly,  this formula reached by the microscopic calculation at weak coupling precisely match
the most general form in the  phenomenological GL action
constructed just from the symmetry principle. This match establishes the deep connection between
the two independent, different, but complementary approaches.




 §.§ 1. Analytical approach: perturbation theory in t_so/t_0. 


In Fig.M1a, at small t_so/t_0 and h=0, we perform the perturbation calculation
in t_so/t_0 which is independent and complementary to the order from disorder analysis done in the main text.
This in-dependent approach can also be used as a check on the results achieved from the order from disorder analysis.



Before performing perturbation calculations,
it is convenient to apply transformation
a_i↓→ (-1)^ia_i↓
to the original Hamiltonian Eq.(1)
and obtain

    ℋ
    	=-t_0∑_⟨ ij⟩a_i^† a_j
    	+it_so∑_⟨ ij⟩(-1)^i
    	a_i^†(d_ij·σ⃗)a_j
    	-h∑_i a_i^†σ^z a_i
    	+U/2∑_i n_i(n_i-1)
    	-μ∑_in_i

Since the ground-state energy is a gauge invariant quantity,
we can use the Hamiltonian Eq.(S1) for the ground-state energy calculation.

In the absent of Zeeman field h=0,
we first separate Hamiltonian Eq.(S1) into two parts

    ℋ_0
    	=-t_0∑_⟨ ij⟩a_i^† a_j
    	+U/2∑_i n_i(n_i-1)
    	-μ∑_in_i,  ℋ_so
    	=it_so∑_⟨ ij⟩
    	(-1)^ia_i^†(d_ij·σ⃗)a_j



Becasue the Hamiltonian ℋ_0 has the spin SU(2)_s symmetry in the rotated S̃Ũ(2) basis
𝐒̃_i=((-1)^i_x+i_yS_i^x,(-1)^i_x+i_yS_i^y,S_i^z),  at the weak coupling U≪ t
one can condense the spinor bosons at k=0 with the spin symmetry breaking direction along (θ, ϕ).
So the ground-state is a superfluid with many-body wavefunction
|Ψ_sf⟩∼(∑Ψ_ia_i^†)^N|Vac⟩
where Ψ_i=c_0 η_0+c_πη_π
and |c_0|^2+|c_π|^2=1.
Notice we are useing the same notation η_0=([ 1; 0 ]) and
η_π=([ 0; 1 ]) as defined in main text.
The coefficients  parameterize the same  as Eq. M4  c_0=e^-iϕ/2cos(θ/2)
and c_π=e^+iϕ/2sin(θ/2).

Within Bogoliubov approximation,
the original boson operators is replaced by condensate part plus a fluctuations
a_i=√(N_0)Ψ_i+ψ_i,
thus one can expand ℋ_0
and keep terms up to second order in the fluctuations.
The quadratic theory of ℋ_0 after the  expansion is ℋ_0,Bog,
which can be diagonalized as

    ℋ_0,Bog
    	=const.+∑_𝐤(
    	    ω_1,𝐤α_1,𝐤^†α_1,𝐤
    	    +ω_2,𝐤α_2,𝐤^†α_2,𝐤)

where Bogoliubov transformation is given by

    ψ_𝐤↑
    	=-c̅_πα_1,𝐤
    	    +c_0(coshξ_𝐤α_2,𝐤
    		+sinhξ_𝐤α_2,-𝐤^†),  
        ψ_𝐤↓
    	=+c̅_0α_1,𝐤
    	    +c_π(coshξ_𝐤α_2,𝐤
    		+sinhξ_𝐤α_2,-𝐤^†)

and the excitation spectrum are ω_1,𝐤=4t_0-2t_0(cos k_x+cos k_y)
and ω_2,𝐤=√(ω_1,𝐤(ω_1,𝐤+2n_0U)).
The auxiliary angle is defined as ξ_𝐤 =-1/2arctanh[ω_1,𝐤/ω_1,𝐤+n_0U].
The ground-state |0⟩ of ℋ_0,Bog is  the vacuum of α_1,𝐤 and α_2,𝐤,
so that α_1,𝐤|0⟩=0 and α_2,𝐤|0⟩=0.



The SU(2)_s symmetry of ℋ_0 guaranteed that c_0 and c_π can take any values
as long as they satisfy normalization condition.  Namely, any spin symmetry breaking direction along (θ, ϕ) is equivalent.
However, ℋ_so explicitly breaks the SU(2)_s symmetry,
so the SU(2)_s symmetry becomes spurious and only exists at mean-field level.

Now the perturbation ℋ_so takes following form in the Bogoliubov approximation

    ℋ_so
    	=2t_so∑_𝐤
    	(γ_𝐤ψ_𝐤↑^†ψ_𝐤+𝐐↓+h.c.)

where γ_𝐤=sin k_x-isin k_y
and 𝐐=(π,π).

In order to apply the perturbation,
it is convenient to express ℋ_so
in terms of α_𝐤:

    ℋ_so   =t_so∑_𝐤
    	[f_1,𝐤α_1,𝐤^†α_1,𝐤+𝐐
    	+f_2,𝐤α_2,𝐤^†α_2,𝐤+𝐐
    	+f_3,𝐤α_1,𝐤^†α_2,𝐤+𝐐
    	+f_4,𝐤α_2,𝐤^†α_1,𝐤+𝐐
       +(f_5,𝐤α_1,𝐤^†α_2,-𝐤+𝐐^†
    	 +f_6,𝐤α_2,𝐤^†α_2,-𝐤+𝐐^†
    	 +h.c.)]

where the coefficients are

    f_1,𝐤   =-sinθ(e^iϕγ_𝐤-e^-iϕγ̅_𝐤)
    
        f_2,𝐤   =sinθ(coshξ_𝐤coshξ_𝐤+𝐐+sinhξ_𝐤sinhξ_𝐤+𝐐)
    	    (e^iϕγ_𝐤-e^-iϕγ̅_𝐤)
    
        f_3,𝐤   =-2coshξ_𝐤+𝐐
    	    [e^iϕsin^2(θ/2)γ_𝐤
    	    +e^-iϕcos^2(θ/2)γ̅_𝐤]
    
        f_4,𝐤   =2coshξ_𝐤
    	    [e^iϕcos^2(θ/2)γ_𝐤
    	    +e^-iϕsin^2(θ/2)γ̅_𝐤]
    
        f_5,𝐤   =-2sinhξ_𝐤+𝐐
    	    [e^iϕsin^2(θ/2)γ_𝐤
    	    +e^-iϕcos^2(θ/2)γ̅_𝐤]
    
        f_6,𝐤   =sinθcoshξ_𝐤sinhξ_𝐤+𝐐
    	(e^iϕγ_𝐤-e^-iϕγ̅_𝐤)




By treating t_so/t_0 as a small parameter,
one can apply perturbation theory to obtain

    E_ofd=δ E^(1)+δ E^(2)
    	+δ E^(3)+δ E^(4)+𝒪(t_so^5)

where we denote δ E^(n)∝ t_so^n.
Since ℋ_so has odd parity,
so all odd order terms vanish and

    δ E^(2)   =
    	⟨ 0|ℋ_so g ℋ_so|0⟩
    δ E^(4)   =
    	⟨ 0|ℋ_sog
    		  ℋ_sog
    		  ℋ_sog
    		  ℋ_so|0⟩
    	-⟨0|ℋ_sog
    		  ℋ_so|0⟩⟨0|ℋ_sog^2
    		  ℋ_so|0⟩

where g=∑_n≠0|n⟩⟨ n|/-E_n.

After a substitution of  ℋ_so in Eq.<ref>, we reach

    δ E^(2)   =-t_so^2∑_k[
    		|f_5,𝐤|^2/ω_1,𝐤+ω_2,𝐐-𝐤
    		+|f_6,𝐤|^2+f_6,𝐤f̅_6,𝐐-𝐤/ω_2,𝐤+ω_2,𝐐-𝐤]

and

    δ E^(4)
    	=-t_so^4∑_𝐤   [
    	|f_3,-𝐤f_5,𝐤|^2/ω_1,𝐤(ω_1,𝐤+ω_2,𝐐+𝐤)^2
    	+1/ω_2,𝐤|
    	f_4,𝐐+𝐤f_5,𝐐+𝐤/ω_1,𝐐+𝐤+ω_2,𝐤
    	+
    	f_2,-𝐤(f_6,𝐤+f_6,𝐐-𝐤)/ω_2,𝐤+ω_2,𝐐-𝐤|^2
    	. 
       .
    	+
    	1/ω_1,𝐤+ω_2,𝐤|
    	f_1,𝐤f_5,𝐐+𝐤/ω_1,𝐐+𝐤+ω_2,𝐤
    	+
    	f_2,-𝐤f_5,𝐐+𝐤/ω_1,𝐤+ω_2,𝐐-𝐤
    	+
    	f_3,𝐤f_6,-𝐤/ω_2,-𝐤+ω_2,𝐐+𝐤
    	+
    	f_3,𝐤f_6,𝐐+𝐤/ω_2,𝐐+𝐤+ω_2,𝐤|^2
    	]
    	+⋯

where ⋯ means ϕ-independent part.


After some manipulations, we found

    δ E^(2)   =t_so^2N_s[c_2+a(1+cos2θ)],
    δ E^(4)   =t_so^4N_s[c_4+bsin^4θ(1+cos4ϕ)+⋯],

where ⋯ also contains θ-only dependent terms like cos2θ and cos4θ which is subleading to δ E^(2) when t_so≪ t.
The coefficients are

    a   =1/N_s∑_𝐤|γ_𝐤|^2(
    	cosh^2ξ_𝐤sinh^2ξ_𝐐+𝐤
    	+coshξ_𝐤sinhξ_𝐤coshξ_𝐐+𝐤sinhξ_𝐐+𝐤/ω_2,𝐤+ω_2,𝐐+𝐤
    	-sinh^2ξ_𝐐+𝐤/ω_1,𝐤+ω_2,𝐐+𝐤)
    
    	b   =1/N_s∑_𝐤(γ_𝐤^4+γ̅_𝐤^4)[
    	-sinh^2(2ξ_𝐤+𝐐)/4ω_1,𝐤
    	-1/4ω_2,𝐤(sinh(2ξ_𝐤)/ω_1,𝐤+𝐐+ω_2,𝐤
    	-sinh(2ξ_𝐤+2ξ_𝐤+𝐐)/ω_2,𝐤+𝐐+ω_2,𝐤)^2 
             
    	+1/ω_1,𝐤+ω_2,𝐤(sinh(ξ_𝐤)/ω_1,𝐤+𝐐+ω_2,𝐤
    	+sinh(ξ_𝐤+𝐐)
    		cosh(ξ_𝐤+ξ_𝐤+𝐐)/ω_1,𝐤+ω_2,𝐤+𝐐
    	+cosh(ξ_𝐤+𝐐)
    		sinh(ξ_𝐤+ξ_𝐤+𝐐)/ω_2,𝐤+ω_2,𝐤+𝐐)^2]


For U>0 case, we  can prove that a,b>0.
So δ E^(2) and δ E^(4)
reach the minima at cos2θ=-1 and cos4ϕ=-1 respectively,
which tells the quantum ground-state sits at θ=π/2 and ϕ=π/4 which are nothing but the XY-FM state.
Transforming it back from the rotated S̃Ũ(2) basis to the original basis, it is nothing but the N=2 XY-AFM SF state.
So we reach the same ground state from the
two independent methods.

 When comparing Eq.<ref> with Eq.M5,

    E_ofd
    	=E_ofd,0
    	+N_s[A/4(1+cos2θ)
    	+B/16sin^4θ(1+cos4ϕ)] + ⋯

   where the ⋯ means the higher order terms respecting the same symmetry of the Hamiltonian.

Then we obtain the analytical expressions of A= a (2t_so)^2
and B= b (2t_so)^4 from the perturbation theory.
Then Eq.<ref> leads to Δ_R ∼ t^3_s at h=0 shown in Fig.S1b.
As shown in the main text, Eq.<ref> acn be written as the form in the GL action:

    A/2(|ψ_1|^2-|ψ_2|^2)^2+B/2[(ψ_1^*ψ_2)^2+(ψ_1ψ_2^*)^2]^2 +
    C(|ψ_1|^2-|ψ_2|^2)^4 + D [(ψ_1^*ψ_2)^2+(ψ_1ψ_2^*)^2]^4  + ⋯

  where ⋯ means the even higher order terms than C and D term
  respecting the [ C_4 × C_4 ]_D symmetry of the Hamiltonian. In fact, as shown below, the C and D terms can also be
  determined numerically by the microscopic calculations. Both found to be positive.



 §.§ 2. Numerical approach: compared with the analytical results


We will show that although the form Eq.(<ref>)
is obtained from the perturbation theory in a small t_so/t_0,
its form fits numerical data very well even upto t_so/t_0=1 .
The data is obtained from numerical evaluation
ω_l,𝐤 and
E_ofd(θ,ϕ) in Eq.M3.

    E_ofd(θ,ϕ)
    	=-(4t_0+1/2n_0U)N_s
    	+1/2∑_l,𝐤ω_l,𝐤

where l=1,..,4 and 𝐤 is defined in
reduced Brillouin zone RBZ={(k_x,k_y)||k_x|+|k_y|<π}


The comparison consist of two parts:
the global property and the local property.

For the global property,
we plot numerical results from Eq.M3
and the analytical result from Eq.(<ref>) in Fig.S2.
It is clear that the two results fit quite well in the whole range 0< θ < ϕ, 0 < ϕ< π/2.

On the other hand, if we only keep the form of the perturbation theory
and treat coefficients A and B as fitting parameters,
then the relative error can be controlled below 1%.
Figure S2 shows the form work very well
even for t_so is not small,
i.e. t_so/t=1.






    In fact, the higher order terms such as (|ψ_1|^2-|ψ_2|^2)^4 and [(ψ_1^*ψ_2)^2+(ψ_1ψ_2^*)^2]^4
    in Eq.(<ref>)  can also be determined from Fig.S2.
    The fact that the period of Δ E is only half of E can be cast into
	these high-order terms.
	For example, the 4th order perturbation also contain θ-only dependent terms
	like cos2θ and cos4θ, but they are proportional to t_so^4.
	The cos2θ terms do not change the form of Eq.(<ref>), they can be absorbed into the coefficient A anyway.
    But cos4θ does. So after the subtraction, only small sub-leading terms survive,
    then Δ E(θ,π/4)∼cos4θ can be cast into the form C (|ψ_1|^2-|ψ_2|^2)^4.
	Similar arguments hold also for ϕ and lead toΔ E(π/2,ϕ)∼cos8ϕ which can be cast into
    D [(ψ_1^*ψ_2)^2+(ψ_1ψ_2^*)^2]^4. Both C and D are found to be positive.





For the local property, we focus on the
derivative with respect to θ and ϕ:

    A'(θ_0)
    	=∂^2/∂θ^2
    	E_ofd(θ,ϕ)|_θ=θ_0,ϕ=π/4,
            
        B'(θ_0)
    	=∂^2/∂ϕ^2
    	E_ofd(θ,ϕ)|_θ=θ_0,ϕ=π/4

Analytical result Eq.(<ref>) predicts
A'_an=-Acos2θ_0 and B'_an=Bsin^4θ_0,
and numerical derivatives can be evaluated from the finite differences of Eq.M3.
We plot both results in Fig.S3(a) and (b),
and find the differences are quite small.
The insert of Fig.S3(b) also confirms  B'∼θ^4 behaviour.

As a final comparison,
we evaluate A' and B' as a function of t_so
and plot them in Fig.S3(c) and (d).
When extrapolating the data to the t_so→0 limit,
they reach agreement with the perturbation results with a high accuracy.
The insert of Fig.S3 (c) and (d) also shows that
the leading behaviours of A'∼ t_so^2 and B'∼ t_so^4,
in consistent with the perturbation theory listed below Eq.eq:pb.

In conclusion,
we show that the E_ofd obtained from the perturbation theory
fit very well the data from the direct numerical evaluations, even upto t_so/t=1.






§ C. SOME ADDITIONAL CALCULATIONS USING THE EFFECTIVE GL ACTION 


 In this appendix, we will use the effective action Eq.<ref> to perform alternative calculations to
 compute the excitation spectrum in both h>h_c and h < h_c. They reproduce the results Eq.<ref> and
 Eq.<ref> respectively achieved in the main text and also provide additional physical insights into the problem.




 §.§ 1. A rotated  polar coordinate calculation at h>h_c

When h>h_c, the system is in the Z-FM superfluid phase
with the saddle point solution Eq.<ref> ρ̅_1=ρ_0=-r/U, ρ̅_2=0, and no constraints on θ̅_1 and θ̅_2.
In order to avoid any singularities in using the polar coordinates, we introduce two rotated order parameters:

    ψ'_1=(ψ_1+ψ_2)/√(2),   ψ'_2=(ψ_1-ψ_2)/√(2)

then parameterize them as
ψ'_1=√(ρ'_1)e^iθ'_1 and ψ'_2=√(ρ'_2)e^iθ'_2.
It is easy to identify the saddle point solution
ρ̅_1=ρ_0, ρ̅_2=0, θ̅_1=θ̅_2=0 in old coordinates system
correspond to ρ̅'_1=ρ̅'_2=ρ_0/2, θ̅'_1=θ̅'_2=0 in new coordinate system.

For the notation simplicity, in the following, we drop the ' in ρ and θ.
After including fluctuations as
ψ_α=√(ρ̅_α+δρ_α)e^i(θ̅_α+δθ_α),
then

    ℒ_ZF
    	=   iδρ_1∂_τδθ_1+iδρ_2∂_τδθ_2
    	+v^2/2ρ_0[(∇δρ_1)^2+(∇δρ_2)^2]
    	+1/2v^2ρ_0[(∇δθ_1)^2+(∇δθ_2)^2]
       +U/2(δρ_1+δρ_2)^2
    	+A/2[4δρ_1δρ_2-ρ_0^2(δθ_1-δθ_2)^2]
    	
    	+h/2[1/ρ_0(δρ_1-δρ_2)^2+ρ_0(δθ_1-δθ_2)^2]

where the linear term drops out due to the saddle point condition
[r+(A+U)ρ_0-h](δρ_1+δρ_2)=0.
From Eq.(<ref>), one can easily extract two eigen-modes in the long wavelength limit:

    ω_-   =√(2ρ_0(A+U)v^2k^2)
    ω_+   =2(h-Aρ_0)+v^2k^2

  which are identical to Eq.<ref>.

 Similar to Eq.<ref>, one can also introduce ± in the new coordinates system:
  δρ_±=δρ_1±δρ_2
and δθ_±=δθ_1±δθ_2,
then

    ℒ_FM   =i/2δρ_+∂_τδθ_+
    	+v^2/4ρ_0[(∇δρ_+)^2+ρ_0^2(∇δθ_+)^2]
    	+A+U/2(δρ_+)^2
       +i/2δρ_-∂_τδθ_-
    	+v^2/4ρ_0[(∇δ n_-)^2+ρ_0[(∇δθ_-)^2]
    	+h-Aρ_0/2ρ_0[(δρ_-)^2+ρ_0^2(δθ_-)^2]

which lead to exactly the same two eigen-modes.




 §.§ 2. A mixed coordinate calculation at h<h_c

One may also use the mixed coordinate system representation to calculate the excitation spectrum
in the CAFM SF at h<h_c. Namely, using the saddle point in Eq.<ref> and substituting  ψ_1=√(ρ̅_1+δρ_1)e^i(θ̅_1+δθ_1)
and ψ_2=√(ρ̅_2)e^iθ̅_2+δψ_2 into Eq.<ref> leads to
the quadratic part of the Lagrangian:

    ℒ_CA
    	=   iδρ_1∂_τδθ_1
    	+δψ_2^*∂_τδψ_2
    	+v^2/4ρ̅_1(∇δρ_1)^2
    	+v^2ρ̅_1(∇δθ_1)^2
    	+v^2|∇δψ_2|^2+(r+h)|δψ_2|^2	
       +1/2(U+A)(δρ_1)^2+8Bρ_1^2ρ_2^2(δθ_1)^2
    	+[2(U+A)ρ_2+(U-A)ρ_1+4Bρ_1^2ρ_2]|δψ_2|^2	
       -i[1/2(U+A)-2Bρ_1^2]ρ_2(δψ_2)^2
    	+i[1/2(U+A)-2Bρ_1^2]ρ_2(δψ_2^*)^2		
       +(U-A)ρ_2^1/2δρ_1(e^-iπ/4δψ_2+e^+iπ/4δψ_2^*)
    	+8Bρ_1^2ρ_2^3/2δθ_1(e^+iπ/4δψ_2+e^-iπ/4δψ_2^*)

From Eq.(<ref>), one can extract the two eigen modes
which are exactly the same as Eq.<ref> obtained in the polar coordinates in the main text.

99





rh For the conventional strong-coupling approach in a SOC system, see Fadi Sun, Jinwu Ye, Wu-Ming Liu, Quantum magnetism of spinor bosons in optical lattices with synthetic non-Abelian gauge fields at zero and finite temperatures,
        Phys. Rev. A 92, 043609 (2015).

more To calculate the gap at h=0, this expansion of the effective potential near the minimum θ=π/2, ϕ=π/4 is enough.
However, to get the correct roton spectrum at any h >0  at any k in the long wavelength limit,
one need to find the complete form of the the effective potential  in all ranges of  0< θ < π, 0< ϕ < 2 π which is found to be Eq.M5 and re-written as Eq.M6 in the main text. Most importantly,
as shown in the main text, just from the symmetry point of view, the A term and B term
in the phenomelogical GL action Eq.<ref> precisely match Eq.M6 at the weak coupling limit. This endows the effective potential generated by OFQD
from a new and profound perspective.








