Critical behavior of a chiral superfluid in a bipartite square lattice

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 [1,2,3,4] or spinor Bose gases [5], 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 [6,7,8,9,10,11,12,13,14,15,16,17,18,19] and from experimentalists [20,21,22,23,24]. According to Feynman's "no-node" theorem [25], such states cannot be a ground state of a conventional bosonic Hamiltonian with shortrange 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 [12]. 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 [11,21]. 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 ± ip y order) that breaks TR symmetry has been realized by carefully tuning the lattice parameters [21,23,24]. Since such a state hosts spatially staggered orbital currents, it is dubbed a chiral superfluid. Similar p x ± ip y paring has been proposed for the A phase of superfluid 3 H [26,27] and for Sr 2 RuO 4 [28].
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 section 2 we develop the field theoretical description of the mixed orbital model in a bipartite optical square lattice in the low energy limit. Section 3 is devoted to a mean-field analysis of the effective model. In section 4, we study the critical behavior of the model in the framework of a one-loop RG calculation. In section 5, 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 figure 1(a) and described by the Hamiltonian, with the tight-binding model of a bipartite optical square lattice Here r = √ 2a(n x , n y ) with n x,y ∈ Z, 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 where k = (k x , k y ) and we set a = 1 in the following calculations. We illustrate the lowest band in momentum space in figure 1(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 Z 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 where 4π 2 e 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 equation (1) by an effective Hamiltonian with two components, where R = (r, z), µ j being the chemical potential of the j-th component and the effective mass being m * = (m 2 xy m 0 ) 1/3 . As an example, we describe the experimental parameters of [21]; for 87 Rb atoms and for typical laser intensities that were used, we have m 0 0.2m xy and m xy = 2 √ 2 2 / (λ 2 L J) with the laser wavelength λ L = 1064nm 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 Λ 2 q /2m * . We further consider the on-site interaction; see [23], which gives the following terms, where being an angular momentum operator. U s is the on-site interaction among s-orbitals. U p and U p are intra-and interorbital 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 [19]; 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 (4), we represent the effective interaction as whereg j (j = 1, 2) is the intra-component interaction andg 12 is the inter-component one. In this expression, there is an additional term with the interaction strengthg u , which is an Umklapp term. This additional scattering process is not present in the previously studied triangular lattice system [29,30,31], 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 [31]. 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: The full symmetry of the effective action is U (1) × Z 2 × Θ, where the U (1) symmetry corresponds to the invariance of the model under the global phase shift for both components, Z 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, 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 (7) as sin 4 θ +g 12 cos 2 θ sin 2 θ +g u cos 2 θ sin 2 θ cos 2φ .
First we note that for the frustrated triangular optical lattice, the Umklapp interaction does not occur (g u = 0), and we haveg 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 equation (10) occurs at θ = 0 or π/2. From equation (9), this means that bosons will condense in one of the energetic minima to break the Z 2 symmetry [31]. However, in the square bipartite lattice that we study in this paper, the Umklapp interactiong 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 equation (10); this corresponds to a chiral superfluid state given by a complex coherent superposition of two single particle states. This state breaks the time-reversal symmetry Θ, i.e., the chiral Z 2 symmetry, in addition to the U (1) continuous symmetry of the phase (The situation is similar to the fully frustrated XY models [32]). Comparing this to the single condensate at θ = 0, we find that the chiral superfluid state occurs when where we setg 1 =g 2 =g 0 . The above condition is also discussed in [19]. In terms of the interaction parameters in equation (6), we findg 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 [31] and ignoring quantum fluctuations, we calculate the RG equations at one-loop order as [33], 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, , and F (µ Λ ) = 1/(1 − µ Λ ). If g u = 0, the flow equations become the ones of the twocomponent φ 4 -theory [31]. In the critical regime, µ Λ 1, where we can approximate F (µ Λ ) ≈ 1, ‡ the RG equations exhibit four fixed points, except the trivial one, (µ * Λ , g * 0 , g * 12 , g * u ) = (0, 0, 0, 0), ‡ This corresponds to the lowest order of the -expansion [31,34,35,36]. The effect of higher order contributions to the fixed points is discussed in appendix.
All these fixed points are unstable, indicating that the system undergoes a first-order transition. In figure 2(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 equation (12d) is which explicitly shows that g u does not change signs under the RG. In figure 2(b), we show the flow diagram on the g u = 0 surface, which is also introduced in [31]. 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 We emphasize that G 1 > 0 coincides with the mean-field stable condition for a chiral superfluid, (11), 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.
Finally the third separatrix surface is G 2 ≡ g 0 − g 12 − g u = 0, and its RG equation The flow equation is similar to equation (15), and it guarantees that no flow passes through the G 2 = 0 plane. As illustrated in figure 2(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 figure 3. For small positive initial interactions, a typical flow has monotonically increasing g u , which can also be seen from equation (14). For the evolution of g 0 and g 12 , there are three regimes that we can characterize: (i) 0 < l < l 1 (the initial regime): the linear terms in RG equations are dominant, and g 0 and g 12 gradually increase.
(ii) 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. (iii) 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ḡ i /(1−kl), where k is the inverse length scale at which these coupling constants diverge [37]. 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 This implies that G 1 steadily increases and thus stabilizes the chiral superfluid order, (11). 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 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 [31]. At oneloop order, one contribution to the RG equation for g 6 is dg 6 /dl −24T Λ 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. q /(2π 2 Λ T Λ ). The phases are determined by the sign of the chemical potential, when the RG flow is terminated at g u = 2. Colors represent the final value of µ. The blue region is the condensed phase separated from the thermal phase with a first-order (discontinuous) phase transition. The pink area represents the thermal phase. As µ Λ increases or T Λ decreases, the condensed region becomes larger.

Phase diagram
The RG analysis can be also used to study the critical behaviors of the phase transition between the two phases. In figure 4, 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 ± ip 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 (17), the mean-field free energy can be written by two order parameters P ± ≡ ψ 1 ± iψ 2 as whereḡ > 0 is the single coefficient characterizing the asymptotic flow (17), 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 [38,39,40]. 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 [30,40].
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 [41,42]. 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 [43,44,45,46]. 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.
Appendix A.1. Derivation of Hubbard parameters We start from a following potential This has two local minima in a unit cell (see figure A1(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 ( figure 1(b)). The energy difference between local minima at A sites and B sites is . 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 we obtain a band dispersion as figure A1(b). We fit the obtained dispersion by the following tight-binding model where e 1 = ( √ 2a, 0) and e 2 = (0, √ 2a). The first term is already given in equation (2). The degeneracy at the Γ point is achieved by setting B = A + 4J . The fitted band dispersion is plotted in figure A1(b) as dots. The fitted parameters as functions of the potential depth V 0 is given in figure A2. 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 [47]. These Wannier functions give the bare interactions of a Bose-Hubbard model as 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 figure A3. We find that U s ∼ U p ∼ U p /3 for moderately strong potential depth.

Appendix A.2. Effective interactions in field-theory approximations
In this subsection, we show that the condition for the chiral superfluidity, (11), is satisfied in general based on the Hubbard parameters determined above. With a general With equation (8), we can show that thẽ 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.