Critical behavior of a chiral superfluid in a bipartite square lattice
Junichi Okamoto Wen-Min Huang Robert Höppner Ludwig Mathey


Critical behavior of a chiral superfluid in a bipartite square lattice]Critical behavior of a chiral superfluid in a bipartite square lattice



^1 Zentrum für Optische Quantentechnologien and Institut für Laserphysik, Universität Hamburg, 22761 Hamburg, Germany
^2 The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg, Germany
^3 Department of Physics, National Chung-Hsing University, Taichung 40227, Taiwan
^*wenmin@phys.nchu.edu.tw


We study the critical behavior of Bose-Einstein condensation in the second band of a bipartite optical square lattice in a renormalization group framework at one-loop order. Within our field theoretical representation of the system, we approximate the system as a two-component Bose gas in three dimensions. We demonstrate that the system is in a different universality class than the previously studied condensation in a frustrated triangular lattice due to an additional Umklapp scattering term, which stabilizes the chiral superfluid order at low temperatures. We derive the renormalization group flow of the system and show that this order persists in the low energy limit. Furthermore, the renormalization flow suggests that the phase transition from the thermal phase to the chiral superfluid state is first order. 




















§ INTRODUCTION

Unconventional Bose-Einstein condensates (BECs) whose order parameter space is not simply the usual U(1) symmetry have been extensively studied. Examples from the field of ultracold atoms include Floquet engineered Bose gases <cit.> or spinor Bose gases <cit.>, where the order parameter space has an additional Ising component or even more complex symmetry groups. The experimental realization of such systems in ultracold atomic systems are ideal to investigate phase transitions of those complex orders due to well-defined and tunable nature of these systems. 

Recently, BECs that break time-reversal (TR) symmetry have attracted increased attention from theorists <cit.> and from experimentalists <cit.>. According to Feynman's “no-node" theorem <cit.>, such states cannot be a ground state of a conventional bosonic Hamiltonian with short-range interactions, since breaking TR symmetry inevitably leads to a wave function with a node in real space. Therefore, to create BECs without TR symmetry, the assumptions of the no-node theorem have to be circumvented. One approach uses a long-lived metastable state of ultracold bosons in bands of higher orbitals <cit.>. Experimentally, a BEC in a p-band has been realized by first populating particles in a staggered pattern in a checkerboard lattice and then suddenly changing the potential shape. Due to the large anharmonicity in the energy spectrum, the life time of the metastable BEC is longer than 100ms <cit.>. Since this BEC is not a ground state, breaking TR symmetry is not in contradiction with the no-node theorem. Indeed, a complex coherent superposition of the two p-band condensates (i.e., p_x ± i p_y order) that breaks TR symmetry has been realized by carefully tuning the lattice parameters <cit.>. Since such a state hosts spatially staggered orbital currents, it is dubbed a chiral superfluid. Similar p_x ± i p_y paring has been proposed for the A phase of superfluid ^3H <cit.> and for Sr_2RuO_4 <cit.>.

In this paper, we investigate the stability of the chiral condensate and its critical behaviors by a renormalization group (RG) analysis. The analysis addresses the competition of chiral and non-chiral condensation, and the critical behavior. While the chiral superfluid has been confirmed experimentally, it is still important to know how stable and general the state is. In particular, since the energies of the competing non-chiral BEC and of the chiral BEC are close at the mean-field level, it is not trivial which of the two BECs becomes dominant at low temperatures. We find that the stable condition of the chiral BEC is always preserved at low energy scales, and the transition is expected to be first-order. 

The paper is organized as follows. In EFT we develop the field theoretical description of the mixed orbital model in a bipartite optical square lattice in the low energy limit. MF is devoted to a mean-field analysis of the effective model. In rg, we study the critical behavior of the model in the framework of a one-loop RG calculation. In conc, we conclude. The details of calculations not covered in the main texts are summarized in appendix.





§ EFFECTIVE FIELD THEORY

  


The system that we consider here is sketched in lattice(a) and described by the Hamiltonian, 

    H_0= ∫ dz
    {∑_r,i b^†_i(r,z)(-ħ^2∂^2_z/2m_0-μ_ 3D) b_i(r,z)+H^ xy_0},

with the tight-binding model of a bipartite optical square lattice

    H^ xy_0   =   J∑_r[b_1^†(r,z)b_2(r+d_1,z)+b_1^†(r,z)b_3(r+d_2,z)
          
    -b_1^†(r,z)b_2(r-d_1,z)-b_1^†(r,z)b_3(r-d_2,z)+ h.c.].

Here r = √(2) a (n_x, n_y) with n_x,y∈ℤ, d_1/2=(a/√(2),± a/√(2)). b_i(r,z) with i=1,2,3 represent the annihilation operators of bosons at the s, p_x and p_y orbitals respectively. We assume that bosons move freely along the z-direction. The hopping amplitudes between neighboring p-orbitals, J_∥ and J_⊥, are set to be zero for simplicity (see appendix for a more general discussion). Converting the orbital representation into a band representation, we obtain three bands. The metastable BEC in experiments is loaded in the lowest band, whose dispersion is given by

    ϵ(k,k_z)=-2J√(1-cos(√(2)k_x)cos(√(2)k_y))+ħ^2 k_z^2/2m_0,

where k=(k_x,k_y) and we set a=1 in the following calculations. We illustrate the lowest band in momentum space in lattice(c). We note that there are two energetic minima at k_1=(π/√(2),0) and k_2= (0,π/√(2)). These energetic minima are degenerate, thus giving rise to the ℤ_2 symmetry of the noninteracting Hamiltonian. 

At low temperatures, bosons predominantly occupy momentum states near the two minima, and then condense below a critical temperature. To describe the critical behavior, we expand the bosonic operators near the two minima as

    b_α(r,z)   =   1/√(N)∑_ke^ik·ru_α(k)ϕ(k,z)
       ≃   2a^2/√(N)∑_j=1,2e^ik_j·ru_α j∫_|q_j|<Λ_qd^2q_j/4π^2 e^iq_j·rϕ_j(q_j,z)
       ≡   √(2)a∑_j=1,2ψ_j(r,z)u_α j e^ik_j·r,

where ψ_j(r,z)≡√(2)a/√(N)∫_|q_j|<Λ_qd^2q_j/4π^2e^iq_j·rϕ_j(q_j,z) with ϕ(k_j+q_j,z)≡ϕ_j(q_j,z). Λ_q is the momentum cut-off, and N is the number of the unit cells. The kernel u_α(k_j)=u_α j represents the projection of the wave function of orbital α on the wave function of the lowest band in the vicinity of the minimum k_j. These are given by (u_11,u_12)=(-i/√(2),i/√(2)), (u_21,u_22)=(1/2, -1/2) and (u_31,u_32)=(1/2,1/2). We use the field decomposition to approximate the full Hamiltonian of (<ref>) by an effective Hamiltonian with two components, 

    H_0^ eff=∑_j=1,2∫ d^3R  ψ^†_j(R)[-ħ^2/2m^*∇_ R^2-μ_j]ψ_j(R),

where R=(r,z), μ_j being the chemical potential of the j-th component and the effective mass being m^*=(m_ xy^2m_0)^1/3. As an example, we describe the experimental parameters of  <cit.>; for ^87Rb atoms and for typical laser intensities that were used, we have m_0≃0.2m_ xy and m_ xy=2√(2)ħ^2/(λ_L^2J) with the laser wavelength λ_L=1064 nm and J/E_ rec=0.13. We note that to simplify the RG analysis we use an isotropic effective model, and the momentum cut-off in the field decomposition sets the energy cut-off of the effective Hamiltonian as ϵ_Λ=ħ^2Λ_q^2/2m^*.

We further consider the on-site interaction; see <cit.>, which gives the following terms,

    H_I   =   ∫ dz∑_r{U_s/2n_s(R)[n_s(R)-1]
          +U_p/2n_p(R)[n_p(R)-1 ]-U'_p/2[ L_z^2(R) - n_p(R) ] },

where n_s =b_1^† b_1, n_p=b_2^†b_2+b_3^†b_3, and L_z=i(b^†_2b_3-b^†_3b_2) being an angular momentum operator. U_s is the on-site interaction among s-orbitals. U_p and U'_p are intra- and inter-orbital interactions among p-orbitals. In the tight-binding approximation, the strength of the on-site interactions can be calculated from the contact interaction and the Wannier functions <cit.>; for the details, see appendix. In the standard harmonic approximation, the on-site interactions follow U'_p= U_p/3. The precise ratio between U_s and U_p depends on the depth of the optical potential since the harmonic frequency of the s-orbital sites are different from the one of the p-orbital sites. For a moderately deep potential, we find U_s ∼ U_p. 

Within the field-theory approximation field, we represent the effective interaction as 

    H^ eff_I=∫ d^3R{g̃_1/2ψ_1^†(R)ψ_1^†(R)ψ_1(R)ψ_1(R)+g̃_2/2ψ_2^†(R)ψ_2^†(R)ψ_2(R)ψ_2(R)
           +g̃_12ψ_1^†(R)ψ_2^†(R)ψ_2(R)ψ_1(R)+g̃_u/2[ψ_1^†(R)ψ_1^†(R)ψ_2(R)ψ_2(R)+ H.c.]},

where g̃_j (j=1,2) is the intra-component interaction and  g̃_12 is the inter-component one. In this expression, there is an additional term with the interaction strength g̃_u, which is an Umklapp term. This additional scattering process is not present in the previously studied triangular lattice system <cit.>, which demonstrates that these two systems are in different universality classes. Our model is more general in the sense that three coupling constants flow independently under the renormalization equations. The Umklapp interaction allows interchange of bosons between the two components by lattice assisted collisions. In other words, the effective interaction only enforces conservation of the total boson number, in stead of the boson number of each component as, for instance, in the frustrated triangular optical lattice <cit.>. Equivalently, we only have one global U(1) symmetry, instead of two U(1) symmetries for each component. The bare values of the coupling constants in terms of U_s, U_p and U'_p are:

    g̃_1/2 =  2 U_s+U_p + 3 U_p'/8,  g̃_12= 2U_s+U_p-U_p'/4,  g̃_u = 1/2g̃_12.


The full symmetry of the effective action is U(1)×ℤ_2 ×Θ, where the U(1) symmetry corresponds to the invariance of the model under the global phase shift for both components, ℤ_2 is the exchange of the two components, and Θ is the time-reversal symmetry.












§ MEAN-FIELD THEORY

Before proceeding to the RG analysis, we study the ground state for the bare interactions within a zero-temperature approach. We assume that the bosons perfectly condense at the two energetic minima so that a many-body trial wave function is represented as,

    |Ψ⟩_θ,ϕ=1/√(N!)[cosθψ_1^†+e^iϕsinθψ_2^†]^N|0,0⟩,

where N is the total number of bosons, and |m,n⟩ stands for bosons' occupation numbers m(n) at momentum k_1(k_2)  respectively. We will use the angles θ and ϕ as variational parameters. θ determines the relative population of the two minima and ϕ denotes the relative phase of the two-component condensates. Using the trial wave function, we compute the energy of the interacting effective Hamiltonian Hieff as

    ⟨ H^ eff_I⟩_θ,ϕ=N(N-1){g̃_1/2cos^4θ+g̃_2/2sin^4θ+g̃_12cos^2θsin^2θ
          +g̃_ucos^2θsin^2θcos2ϕ}.

First we note that for the frustrated triangular optical lattice, the Umklapp interaction does not occur (g̃_u = 0), and we have g̃_1=g̃_2 ≈ 2g̃_12, which follows from the common origin of these terms, i.e., the contact interaction between the atoms. In this case, the minimum of (<ref>) occurs at θ=0 or π/2. From (<ref>), this means that bosons will condense in one of the energetic minima to break the ℤ_2 symmetry <cit.>. However, in the square bipartite lattice that we study in this paper, the Umklapp interaction g̃_u>0 exists due to the bare on-site repulsive interactions U_s ∼ U_p ≫ U'_p >0. In this case, another energetic minimum may appear at (ϕ, θ)=(±π/2, π/4) in (<ref>); this corresponds to a chiral superfluid state |Ψ⟩= [ψ_1^†± i ψ_2^†]^N|0,0⟩ /√(2N!) given by a complex coherent superposition of two single particle states. This state breaks the time-reversal symmetry Θ, i.e., the chiral ℤ_2 symmetry, in addition to the U(1) continuous symmetry of the phase (The situation is similar to the fully frustrated XY models <cit.>). Comparing this to the single condensate at θ = 0, we find that the chiral superfluid state occurs when

    G_1 ≡g̃_0-g̃_12+g̃_u>0,

where we set g̃_1=g̃_2=g̃_0. The above condition is also discussed in <cit.>. In terms of the interaction parameters in (<ref>), we find g̃_0-g̃_12+g̃_u = U'_p/2 >0, and thus a chiral superfluid order occurs for any repulsive interaction. Even if we include non-zero values of J_⊥ and J_∥, we find that the condition is still satisfied as long as the energetic minima are located at k_1=(π/√(2),0) and k_2= (0,π/√(2)) (see appendix). However, the energy difference between the normal and chiral condensates is of the order of ∼ U_p' at the mean-field level, and at low temperatures the coupling constants get renormalized under the RG flow. Then it is nontrivial which superfluid order emerges at low temperatures. To study this problem more systematically, and to study the critical behavior in the low-energy limit, we apply the renormalization group method in the next section.   




§ ONE-LOOP RENORMALIZATION GROUP METHOD

Following the RG method employed in <cit.> and ignoring quantum fluctuations, we calculate the RG equations at one-loop order as <cit.>, 

    dμ_Λ/dl   =   2μ_Λ-T_Λ F(μ_Λ)(2 g_0 + g_12),
    dg_0/dl   =   ϵ g_0-T_ΛF(μ_Λ)^2 ( 5g_0^2 + g_12^2 + g_u^2 ), 
    dg_12/dl   =   ϵ g_12-T_ΛF(μ_Λ)^2 ( 4g_0g_12+2g_12^2+4g_u^2), 
    dg_u/dl   =   ϵ g_u-T_ΛF(μ_Λ)^2 (2g_0g_u+4g_12g_u),

where l= ln(Λ_q/Λ_b) is the logarithm of the ratio between the bare momentum cutoff Λ_q and the running cutoff Λ_b. Here ϵ = 4-d with d being the spatial dimension of the system, and ϵ = 1 for our three-dimensional model. We have also defined dimensionless parameters,  μ_Λ=μ_1/ϵ_Λ=μ_2/ϵ_Λ, T_Λ=k_BT/ϵ_Λ, g_i=g̃_i Λ_q^3 /(2π^2ϵ_Λ), and F(μ_Λ) = 1/(1-μ_Λ). If g_u=0, the flow equations become the ones of the two-component ϕ^4-theory <cit.>. In the critical regime, μ_Λ≪ 1, where we can approximate F(μ_Λ) ≈ 1, [This corresponds to the lowest order of the ϵ-expansion <cit.>. The effect of higher order contributions to the fixed points is discussed in appendix.] the RG equations exhibit four fixed points, except the trivial one, (μ^*_Λ,g^*_0,g^*_12,g^*_u)=(0,0,0,0),

    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (1/5,1/5T_Λ,0,0),
    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (1/4,1/6T_Λ,1/6T_Λ,0),
    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (1/5,1/10T_Λ,1/5T_Λ,1/10T_Λ),
    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (1/5,1/10T_Λ,1/5T_Λ,-1/10T_Λ).


All these fixed points are unstable, indicating that the system undergoes a first-order transition. In flowD(a), we show the fixed points as red points for g_u≥0. 



 §.§ Basic structures of RG equations

While the full RG equations are complicated and can only be solved numerically,  the basic structure of the equations give useful insight in the RG flow. In particular, we find three separatrix surfaces of the RG equations. RG flow cannot pass through these surfaces, and therefore the asymptotic behavior of the RG flow is severely constrained by the initial condition. 

The first separatrix is g_u=0; the formal solution of (<ref>) is

    g_u(l) = g_u (0)  exp{∫_0^l dl' [ 1-T_ΛF(μ_Λ)^2 (2 g_0 + 4  g_12) ] },

which explicitly shows that g_u does not change signs under the RG. In flowD(b), we show the flow diagram on the g_u=0 surface, which is also introduced in <cit.>. The two fixed points (g^*_0,g^*_12,g^*_u)=(0,0,0) and (1/(5T_Λ),0,0) are unstable, and the one at (1/(6T_Λ),1/(6T_Λ),0) is marginally unstable.  



The second separatrix surface is G_1≡ g_0-g_12+g_u = 0. The RG equation for  G_1 is 

    dG_1/dl=G_1-T_ΛF(μ_Λ)^2G_1(5g_0+g_12-3g_u).

We emphasize that G_1 > 0 coincides with the mean-field stable condition for a chiral superfluid, chiralC, and therefore, for general repulsive interactions that give G_1 >0 as an initial condition, a chiral superfluid order always persists in the low-energy limit. In flowD(c), we illustrate the flow diagrams on the G_1=0 plane for g_u>0. There are two unstable fixed points, fp2 and fp3, on the plane, in addition to the trivial one (g^*_0,g^*_12,g^*_u)=(0,0,0). In particular the ray from (g^*_0,g^*_12,g^*_u)=(0,0,0) to (1/10T_Λ,1/5T_Λ,1/10T_Λ) separates the RG flows into two parts. For initial conditions g_0>g_u on the G_1=0 plane, we find that system eventually flows into the fixed point (g^*_0,g^*_12,g^*_u)=(1/6T_Λ,1/6T_Λ,0), where the system is reduced to the standard two-component ϕ^4 theory. However, if the initial conditions deviate from the G_1=0 plane by an arbitrarily small amount, g_u will grow to a large positive value, and thus the fixed point fp2 is actually unstable. For initial conditions g_0 < g_u on the G_1=0 plane, the flow first approaches (1/10T_Λ,1/5T_Λ,1/10T_Λ) and then runs away to larger positive values of g_u.







Finally the third separatrix surface is G_2 ≡ g_0-g_12-g_u =0, and its RG equation is

    dG_2/dl=G_2-T_ΛF(μ_Λ)^2G_2(5g_0+g_12+3g_u).

The flow equation is similar to (<ref>), and it guarantees that no flow passes through the G_2=0 plane. As illustrated in flowD(d), we find that the Umklapp interaction g_u always grows up on the G_2=0 plane. The increase of g_u is also found when the initial couplings deviate from the G_2 plane as we discuss below. 





 §.§ RG flows for the effective model

Now let us analyze the RG equations for the relevant parameter regime of our model. For our effective model, we have U_s ∼ U_p ∼ U'_p/3 >0, and the RG flow is constrained to the space of g_u>0, G_1 >0 and G_2 <0. Due to this constraint, as we will show below, the possible phases of our effective model are either the thermal gas phase (μ_Λ→ - ∞) or the chiral superfluid phase (μ_Λ→ + 1). The phases are determined by the non-universal nature of the RG flows and initial conditions. 

First, when μ_Λ remains positive, a typical flow of coupling constants behaves as in RG_flow. For small positive initial interactions, a typical flow has monotonically increasing g_u, which can also be seen from (<ref>). For the evolution of g_0 and g_12, there are three regimes that we can characterize:


  * 0 < l < l_1 (the initial regime): the linear terms in RG equations are dominant, and g_0 and g_12 gradually increase.

  * l_1 < l < l_2 (the intermediate regime): the quadratic terms become more important with increasing g_u, which eventually make g_0 and g_12 negative. 

  * l_2 < l< l_3 (the asymptotic regime): the quadratic terms give asymptotically diverging behaviors, while the coupling constants are still smaller than the unity. 

Of course, we should stop the RG flow before any of the coupling constant becomes order of unity near l_3, above which the cubic terms become dominant. In the asymptotic regime, μ_Λ approaches one, and the quadratic terms in the RG equations become dominant. Ignoring the linear terms, the flow can be analyzed by an ansatz g_i (l) = k g̅_i/(1- kl), where k is the inverse length scale at which these coupling constants diverge <cit.>. We find that the only asymptotic flow constrained in the space of g_u>0, G_1 >0 and G_2 <0 is given by 

    (g̅_0, g̅_12, g̅_u) = 1/T_ΛF(μ_Λ)^2( -1/10, -1/5, 1/10).

This implies that G_1 steadily increases and thus stabilizes the chiral superfluid order, chiralC. At the same time, the quartic interactions g_0 and g_12 are renormalized to negative values. This indicates the breakdown of the quartic effective field theory, and we need to include higher order interactions such as

    g_6 ∫ d^3R{ |ψ_1 (R)|^6 + |ψ_2 (R)|^6 },

which is generated from three quartic vertices after one-loop renormalization. We note that at tree-level the g_6 contribution is marginal, and stabilizes the system <cit.>. At one-loop order, one contribution to the RG equation for g_6 is d g_6/ dl ≃  - 24 T_Λ F(μ_Λ)^2 g_0 g_6. Since g_0 flows to negative values, the above contribution further stabilizes the system. 

Second, when the chemical potential becomes negative, F(μ_Λ)=1/(1-μ_Λ) gets suppressed, making the linear terms in the RG equations more dominant. Therefore, typical RG flows give a simple scaling behavior, (μ_Λ,g_i ) ∼(- e^2l,  e^l ). This corresponds to the thermal gas phase without condensation.









 §.§ Phase diagram


The RG analysis can be also used to study the critical behaviors of the phase transition between the two phases. In phaseD, we illustrate the phase diagram as a function of U_s and U_p. The pink region represents the thermal gas phase, in which μ_Λ flow to negative values under RG transformations. Th blue region is the chiral p_x± i p_y superfluid order with time-reversal symmetry breaking, where μ_Λ approaches one. As μ_Λ(0) increases or T_Λ decreases, the phase boundary is shifted to enlarge the superfluid region. We note that along the asymptotic flow asymptotic, the mean-field free energy can be written by two order parameters P_±≡ψ_1 ± i ψ_2 as

    F^Λ_ MF∼μ_Λ (|P_-|^2 + |P_+|^2) - g̅ (|P_-|^4 + |P_+|^4) + 𝒪 (P_-^6, P_+^6),

where g̅ >0 is the single coefficient characterizing the asymptotic flow asymptotic, and the last term is the higher order correction stabilizing the system. This form suggests that the transition is first order. Another indirect support for this scenario is obtained by considering the strong coupling limit for g_0 and g_u, while g_12 is set to 0. In this case, the model has effectively two XY spins on each space point that are orthogonal to each other due to the strong g_u interaction. Such a model is known as Stiefel's V_2,2 model, and Monte Carlo studies show that it undergoes a first order transition <cit.>. For moderate interaction strength, we speculate that the first order nature becomes weaker. Clarifying if the transition remains first order for weak interactions seems to require an extensive numerical simulations, which are beyond the scope of the paper. [The difficulty of determining the order of transitions is a common problem in the frustrated spin systems, whose effective Ginzburg-Landau models are similar to ours <cit.>.]

To detect this first order transition in experiments, we suggest two possibilities. The first one consists of measuring the condensate fraction as a function of temperature, preferably in a box potential <cit.>. The measured temperature dependence will approach a non-zero jump at the condensation temperature for large systems. As a second approach we suggest to measure the spatial evolution of a phonon pulse in a condensate in a smoothly varying trap <cit.>. Here, the phonon velocity will vary as the phonon pulse approaches the interface of the condensate and the thermal gas. For a second order transition the pulse velocity will smoothly approach zero, whereas for a first order transition the velocity will approach a non-zero value before the pulse is reflected, which gives a clear indication of a first order transition.







§ CONCLUSIONS


In this paper we have investigated the critical behavior of unconventional Bose-Einstein condensates in the second band of an optical lattice. We have demonstrated that an Umklapp process between the two minima of the dispersion stabilizes a chiral superfluid state that breaks time reversal symmetry, first at the mean-field level and then within a renormalization group calculation. The latter shows that this stability is always persistent at low energy scales after integrating out thermal fluctuations. We obtain this result by identifying three separatrix planes in the RG flow, which constrain the low energy behavior to a stable regime. Furthermore, the RG flow suggests that the phase transition of the chiral superfluid state to the thermal state is of first order, in contrast to the usual second order transition of a conventional condensate. 




We thank A. Hemmerich for helpful discussions on the experimental aspects of the system. J.O., R.H., and L.M. acknowledge the support from the Deutsche Forschungsgemeinschaft (through SFB 925 and EXC 1074) and from the Landesexzellenzinitiative Hamburg, which is supported by the Joachim Herz Stiftung. W.M.H. especially acknowledges the support from Ministry of Science and Technology, Taiwan through Grant No. MOST 104-2112-M-005-006-MY3. 





§ EFFECTIVE INTERACTIONS FOR GENERAL HOPPING AMPLITUDES AND INTERACTIONS

In this appendix, we start from a single particle picture in a bipartite optical lattice to derive the hopping amplitudes and interaction parameters of a Bose-Hubbard model. We then show that the condition for the chiral superfluid state G_1 = g̃_0 - g̃_12 + g̃_u >0 in (<ref>) is preserved as long as the band has two minima at k_1=(π/√(2),0) and k_2= (0,π/√(2)).



 §.§ Derivation of Hubbard parameters

We start from a following potential

    V(𝐫) = - V_0| cos[  k_0 (x+y)/√(2)]+ e^iβcos[  k_0 (x-y)/√(2)]  |^2.

This has two local minima in a unit cell (see app1(a)): one is a shallow local minimum hosting a s-orbital like Wannier state (denoted as A sites), and the other is a deep minimum hosting two p-orbital like Wannier sates (denoted as B sites). k_0 = π/a with a being the distance between neighboring A and B sites (lattice(b)). The energy difference between local minima at A sites and B sites is Δ V = E_A^0 - E_B^0 = - 4 V_0cos(β). We tune β so that the doubly degenerate first excited states in site B is close to the ground state in site A; in particular, we choose β so that the three bands are exactly degenerate at the Γ point in the following.

 

After solving the single particle Schrödinger equation 

    [ ħ^2 ∇^2/2 m  + V(𝐫) ] ψ^n_𝐤 (𝐫) = E^n_𝐤ψ^n_𝐤 (𝐫),

we obtain a band dispersion as app1(b). We fit the obtained dispersion by the following tight-binding model

    H^ xy_0    = J ∑_𝐫∈ A[b_1^†(r,z)b_2(r+d_1,z)+b_1^†(r,z)b_3(r+d_2,z)
       -b_1^†(r,z)b_2(r-d_1,z)-b_1^†(r,z)b_3(r-d_2,z)+ h.c.] 
       - J_⊥∑_𝐫∈ B, ν = x,y[ b_2^†(r,z)b_3(r+e_ν,z)+b_3^†(r,z)b_2(r+e_ν,z)+ h.c.] 
       - J_∥∑_𝐫∈ B, ν = x,y∑_i=2,3[ b_i^†(r,z)b_i(r+e_ν,z)+ h.c.] 
       + ϵ_A ∑_𝐫∈ A b_1^†(r,z) b_1(r,z) +  ϵ_B ∑_𝐫∈ B∑_i=2,3 b_i^†(r,z) b_i(r,z),

where e_1 = (√(2) a,0) and e_2 = (0, √(2)a). The first term is already given in (<ref>). The degeneracy at the Γ point is achieved by setting ϵ_B = ϵ_A + 4 J_∥. The fitted band dispersion is plotted in app1(b) as dots. The fitted parameters as functions of the potential depth V_0 is given in app2. We note that when the potential is relatively deep, the hopping between p-orbitals (J_⊥ and J_∥) is much smaller than the one between s- and p-orbitals, J. Thus, in practice, we can ignore J_⊥ and J_∥.

 

Now we use the Bloch wave functions for the obtained band dispersions to construct localized Wannier functions. Here we employ a simple projection approach <cit.>. These Wannier functions give the bare interactions of a Bose-Hubbard model as

    U_s   =   g ∫ dz |w_z(z)|^4∫ d^2 r|w_1( r)|^4,
    
    U_p   =   g∫ dz |w_z(z)|^4∫ d^2 r|w_2/3( r)|^4, 
    
    U'_p   =   g∫ dz |w_z(z)|^4∫ d^2 r|w_2( r)|^2|w_3( r)|^2,

where w_z(z) is the Wannier function of a harmonic trap along the z-axis, and w_i( r), i=1,2,3 are the Wannier functions of the s-, p_x and p_y-orbitals on the xy-plane respectively. g is the contact interaction strength. The obtained values are plotted in app3. We find that U_s ∼ U_p ∼ U'_p/3 for moderately strong potential depth. 



 §.§ Effective interactions in field-theory approximations

In this subsection, we show that the condition for the chiral superfluidity, chiralC, is satisfied in general based on the Hubbard parameters determined above. With a general dispersion in (<ref>) and ϵ_B = ϵ_A + 4 J_∥, the projection of the Wannier orbitals to the two minima in the lowest band becomes

    (u_11, u_12)    =  1/√(2+|λ|^2)(λ, λ^* ),
    
    (u_21, u_22)    =  1/√(2+|λ|^2)(1, -1) ,
    
    (u_31, u_32)    =  1/√(2+|λ|^2)(1,  1).

with 

    λ = i Δ_J - √(2 J^2 + Δ_J^2)/J and Δ_J = J_⊥- J_∥.

With (<ref>), we can show that the 

    g̃_0-g̃_12+g̃_u  = 4  U'_p J^4  [2 J^2+Δ_J^2] [ J^2- Δ_J (√(2 J^2+Δ_J^2)-Δ_J )]/3 (Δ_J (√(2 J^2+Δ_J^2)-Δ_J )-2 J^2)^4.

For J ≫ J_∥, J_⊥, the above quantity is always positive. Therefore, even when the system deviates from the simple limit of J_⊥ = J_∥ = 0, the condition of the chiral superfluidity is still satisfied.


 



§ THE Ε-EXPANSION ANALYSIS OF THE FIXED POINTS

In rg, we have shown the fixed points that correspond to the lowest order of the ϵ-expansion <cit.>. To estimate the high-order effect in the ϵ-expansion, we investigate the fixed points of the RG equations without expanding F (μ_Λ). Finding the zeros of the right-hand sides of muFlow-guFlow leads to four fixed points, except the trivial one, (μ^*_Λ,g^*_0,g^*_12,g^*_u)=(0,0,0,0),

    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (ϵ/5+ϵ,5 ϵ/(5+ϵ)^2 T_Λ,0,0),
    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (ϵ/4+ϵ,8 ϵ/3(4+ϵ)^2 T_Λ,8 ϵ/3(4+ϵ)^2 T_Λ,0),
    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (ϵ/5+ϵ,5 ϵ/2(5+ϵ)^2 T_Λ,5 ϵ/(5+ϵ)^2 T_Λ,5 ϵ/2(5+ϵ)^2 T_Λ),
    (μ^*_Λ,g^*_0,g^*_12,g^*_u)   =   (ϵ/5+ϵ,5 ϵ/2(5+ϵ)^2 T_Λ,5 ϵ/(5+ϵ)^2 T_Λ,-5 ϵ/2(5+ϵ)^2 T_Λ),

Taking the lowest order in ϵ and setting ϵ = 1 recovers the fixed points in fp1-fp4. We emphasize that the above fixed points still lie on the three separatrix planes g_u=0, G_1=0 and G_2=0 in the g_0-g_12-g_u space. This fact guarantees the persistence of the chiral superfluid in the low energy limit even when we include the higher order terms in ϵ in our one-loop RG study. 





10
url<#>1#1urlprefixURL 



Struck2011
Struck J, Ölschläger C, Targat R L, Soltan-Panahi P, Eckardt A,
  Lewenstein M, Windpassinger P and Sengstock K 2011 Science 333
  996

Struck2013
Struck J, Weinberg M, Ölschläger C, Windpassinger P, Simonet J,
  Sengstock K, Höppner R, Hauke P, Eckardt A, Lewenstein M and Mathey L
  2013 Nat. Phys. 9 738

Parker2013
Parker C, Ha L and Chin C 2013 Nat. Phys. 9 769

Clark2016
Clark L W, Feng L and Chin C 2016 Science 354 606

Kawaguchi2012
Kawaguchi Y and Uedaa M 2012 Phys. Rep. 520 253

Isacsson2005
Isacsson A and Girvin S M 2005 Phys. Rev. A 72 053604

Liu2006
Liu W and Wu C 2006 Phys. Rev. A 74 013607

Kuklov2006
Kuklov A B 2006 Phys. Rev. Lett. 97 110405

Wu2006
Wu C, Liu W, Moore J and Sarma S 2006 Phys. Rev. Lett. 97 190406

Lim2008
Lim L K, Smith C M and Hemmerich A 2008 Phys. Rev. Lett. 100
  130402

Stojanovic2008
Stojanovic V M, Wu C, Liu W V and Das Sarma S 2008 Phys. Rev. Lett.
  101 125301

Wu2009
Wu C 2009 Mod. Phys. Lett. B 23 1

Lewenstein2011
Lewenstein M and Liu W V 2011 Nat. Phys. 7 101

Cai2011
Cai Z and Wu C 2011 Phys. Rev. A 84 033635

Li2012
Li X, Zhang Z and Liu W V 2012 Phys. Rev. Lett. 108 175302

Martikainen2012
Martikainen J P and Larson J 2012 Phys. Rev. A 86 023611

Cai2012
Cai Z, Duan L M and Wu C 2012 Phys. Rev. A 86 051601

Liu2013
Liu B, Yu X L and Liu W M 2013 Phys. Rev. A 88 063605

Li2016
Li X and Liu W V 2016 Reports Prog. Phys. 79 116401

Olschlager2011
Ölschläger M, Wirth G and Hemmerich A 2011 Phys. Rev. Lett.
  106 015302

Wirth2011
Wirth G, Ölschläger M and Hemmerich A 2011 Nat. Phys. 7 147

Soltan-Panahi2012
Soltan-Panahi P, Lühmann D S, Struck J, Windpassinger P and Sengstock K
  2012 Nat. Phys. 8 71

Olschlager2013
Ölschläger M, Kock T, Wirth G, Ewerbeck A, Morais Smith C and
  Hemmerich A 2013 New J. Phys. 15 083041

Kock2015
Kock T, Ölschläger M, Ewerbeck A, Huang W M, Mathey L and Hemmerich
  A 2015 Phys. Rev. Lett. 114 115301

feynman1998statistical
Feynman R P 1998 Statistical Mechanics: A Set of Lectures 
  (New York: Avalon Publishing)

Leggett1975
Leggett A J 1975 Rev. Mod. Phys. 47 331

dobbs2000helium
Dobbs R 2000 Helium Three (Oxford: Oxford University Press)

Mackenzie2003
Mackenzie A P and Maeno Y 2003 Rev. Mod. Phys. 75 657

Antonenko1994
Antonenko S A and Sokolov A I 1994 Phys. Rev. B 49 15901

Kawamura1998
Kawamura H 1998 J. Phys. Condens. Matter 10 4707

Janzen2016
Janzen P, Huang W M and Mathey L 2016 Phys. Rev. A 94 063614

Villain1977
Villain J 1977 J. Phys. C Solid State Phys. 10 1717

stoof2008ultracold
Stoof H T C, Dickerscheid D B M and Gubbels K 2009 Ultracold Quantum
  Fields (Dordrecht: Springer Netherlands)

cardy1996scaling
Cardy J 1996 Scaling and Renormalization in Statistical Physics (Cambridge: Cambridge University Press)

Wilson1972
Wilson K G and Fisher M E 1972 Phys. Rev. Lett. 28 240

Domany1977
Domany E, Mukamel D and Fisher M E 1977 Phys. Rev. B 15 5432


Balents1996
Balents L and Fisher M P A 1996 Phys. Rev. B 53 12133

Kunz1993
Kunz H and Zumbach G 1993 J. Phys. A Math. Gen 26 3121

Loison1998
Loison D and Schotte K D 1998 Eur. Phys. J. B 743 735

Itakura2003
Itakura M 2003 J. Phys. Soc. Japan 72 74

Gaunt2013
Gaunt A L, Schmidutz T F, Gotlibovych I, Smith R P and Hadzibabic Z 2013 Phys. Rev. Lett. 110 200406

Corman2014
Corman L, Chomaz L, Bienaimé T, Desbuquois R, Weitenberg C,
  Nascimbène S, Dalibard J and Beugnon J 2014 Phys. Rev. Lett.
  113 135302

Tey2013
Tey M K, Sidorenkov L A, Guajardo E R S, Grimm R, Ku M J H, Zwierlein M W, Hou
  Y H, Pitaevskii L and Stringari S 2013 Phys. Rev. Lett. 110
  055303

Sidorenkov2013
Sidorenkov L A, Tey M K, Grimm R, Hou Y H, Pitaevskii L and Stringari S 2013
  Nature 498 78

Weimer2015
Weimer W, Morgener K, Singh V P, Siegl J, Hueck K, Luick N, Mathey L and Moritz
  H 2015 Phys. Rev. Lett. 114 095301

Singh2016
Singh V P, Weimer W, Morgener K, Siegl J, Hueck K, Luick N, Moritz H and Mathey
  L 2016 Phys. Rev. A 93 023634

Marzari2012
Marzari N, Mostofi A A, Yates J R, Souza I and Vanderbilt D 2012 Rev. Mod.
  Phys. 84 1419




