A two-dimensional algebraic quantum liquid produced by an atomic simulator of the quantum Lifshitz model

Bosons have a natural instinct to condense at zero temperature. It is a long-standing challenge to create a high-dimensional quantum liquid that does not exhibit long-range order at the ground state, as either extreme experimental parameters or sophisticated designs of microscopic Hamiltonians are required for suppressing the condensation. Here we show that synthetic gauge fields for ultracold atoms, using either the Raman scheme or shaken lattices, provide physicists a simple and practical scheme to produce a two-dimensional algebraic quantum liquid at the ground state. This quantum liquid arises at a critical Lifshitz point, where a two-dimensional quartic dispersion emerges in the momentum space, and many fundamental properties of two-dimensional bosons are changed in its proximity. Such an ideal simulator of the quantum Lifshitz model allows experimentalists to directly visualize and explore the deconfinement transition of topological excitations, an intriguing phenomenon that is difficult to access in other systems.


Supplementary Note 1
Quartic dispersions in shaken lattices By shaking a lattice, an effective spin-orbit coupling can be produced, where band indices play the role of spin [1,2]. If one shakes a one-dimensional lattice, it has been shown that a quartic dispersion could be produced [1][2][3]. This method can be straightforwardly generalized to two dimensions.
For a shaken square lattice with the lattice potential V (x, y) = V 0 cos 2π a (x + f cos ωt) + cos 2π a (y + f cos ωt) , the Hamiltonian is given in equation (4) of the main text. The interband coupling C can be written as (1) where J 1 (πf /a) is the Bessel function, and W s (x, y), W px (x, y) and W py (x, y) are the Wannier wave functions of the three bands. In the small coupling limit |C| δ L , A slight difference with the Raman scheme is that the lattice potential reduces the symmetry to a four-fold one, instead of a full rotation symmetry in the momentum space. The quadratic term vanishes at the critical values δ L± c , At δ L+ c , Away from the critical point, however, the dispersion becomes quadratic, as show in Supplementary Fig. 1.

Ω ≥ Ω c
Expanding the Lagrangian about the mean-field solution, up to quadratic order of the fluctuations we have where is defined in such a way that it is dimensionless and equals to 1 at Ω = Ω c .
The fields δρ, δφ and δχ are all gapped and can be integrated out, giving the effective Lagrangian for θ In particular, for the isotropic case with η = 1 we obtain Alternatively, at the critical field Ω = Ω c , we have in which the ∼ ∂ 2 x term is suppressed, leading to the vanishing of T BKT .
At the quantum critical point Ω = Ω c and η = 1, we have Z 0 = 1 and the usual dominant spatial derivative ∼ ∇ 2 is suppressed, giving the Quantum Lifshitz model where the term ∼ ∂ 2 τ ∇ 2 can be dropped as it now enters as a higher order correction in the absence of the usual ∼ ∇ 2 term.

Ω < Ω c
In the same manner the Lagrangian is expanded as whereg 0 = g 0 + g s cos 2 φ 0 /2. Note the additional dependence on φ 0 as it is now a function of Ω.
In integrating out each of the massive fields δρ, δφ and δχ, the coefficients in the effective Lagrangian receive corrections and the expressions are substantially more complicated than those for Ω ≥ Ω c . The corrections are expressed in terms of symbols Z s and Y s, which are all defined in such a way that they are dimensionless and equal to 1 at Ω = Ω c (η). The effective Lagrangian for θ is and the lengthy expressions for Z and Y in terms of the system parameters are listed in the last section of this supplementary note. Note also that the coefficients of the effective Lagrangian are continuous at Ω = Ω c .

Condensate fraction and correlation functions
The effective Lagrangian was computed via derivative expansion and terms up to 4 th order in derivatives are kept. The most general effective Lagrangian considered can be parameterize as where terms like ∼ (∂ τ θ)(∂ y θ), which can contribute to the observables to the same order of our calculations, are not generated in integrating out the heavy modes and are therefore absent.
For the isotropic case with η = 1 and Ω > Ω c , we haveα x =α y =α andβ x =β y =β xy /2 =β. The characteristic momentum scale is given by q * = α/β and the integral is given by whereq andr are both dimensionless (measured in units of q * and (q * ) −1 ). For |r| (q * ) −1 , the integral is dominated by the small momentum contribution and so giving the standard long-range correlation in the condensate.
Alternatively, for intermediate values of |r| such that ξ |r| (q * ) −1 , the large momentum contribution of the integral dominates and it gives a power-law correlation function where K = 4π α τβ is the effective Luttinger liquid parameter. As such the correlation function has a crossover behavior set by the characteristic length scale (q * ) −1 . In the limit Ω → Ω c , (q * ) −1 diverges and therefore the correlation function is dominated by the power-law behavior.
Expressions for Z's and Y 's in the effective Lagrangian for Ω < Ω c

Supplementary Note 3 Effective theory for shaken lattices
We analyze the long-wavelength limit of the four-band model described in the main text and Supplementary Note 1, for which the free single-particle Hamiltonian is given by H 0 (k) − µ, with µ denoting the chemical potential. In particular, we focus on the deep lattice limit, for which we can keep only the on-site interaction termŝ U = x,σ=s,px,py,dxy Note that we have dropped terms like , which vanish under time-averaging. Deep in the superfluid phase, ψ σ (x) = b xσ is smooth on lattice scale and we consider the energy functional of the four-component Schrödinger field Ψ L that is given by equation (21) of the main text. To simplify expressions, we define¯ In addition, we also absorb the phase of C into the definitions of φ and ζ and henceforth we use C to denote its magnitude. In these variables, we have E =E KE − µρ +¯ ρ cos 2 χ + δ ρ cos 2 χ cos(2ν) + p ρ sin 2 χ + √ 2Cρ sin(2χ) [cos ϕ cos ξ (cos ν cos φ + sin ν cos ζ) + sin ϕ sin ξ (cos ν sin φ − sin ν sin ζ)] where E KE is the kinetic energy which can be straightforwardly evaluated from equation (4) and (21) of the main text.
The mean-field solution is analyzed by taking the mean-field values of the massive fields ρ 0 , χ 0 , ν 0 , ϕ 0 , φ 0 , ξ 0 and ζ 0 to be spatially homogeneous, and θ = κ(x + y). The mean-field equations are then obtained by requiring ∂ F0 E = 0 for F being any of the massive fields or κ. In particular, we have φ 0 = ξ 0 = ζ 0 = 0 for weak interactions. Note that when U pp ρ 0 a 2 dominates over C, there are additional non-trivial solutions with ξ 0 → ±π/4. These solutions correspond to spin-canting resulting from strong repulsive inter-band interaction between bosons in the p x and p y bands, but they are unattainable when 2 √ 2U pp ρ 0 a 2 < C cot χ 0 csc 2 χ 0 , which is satisfied in the weakly interacting regime we are focusing on.
The phase-boundary between the zero-and finite-momentum phases is determined by the mean field equation In the zero-momentum phase, κ = 0 so Eq.(26) is satisfied trivially and due to the symmetry betwene p x and p y bands we have ϕ 0 = 0. In the finite momentum phase, however, κ = 0 and (ϕ 0 , ν 0 , χ 0 ) satisfy Eq.(26) non-trivially. The phase boundary is then determined by requiring (ϕ c 0 = 0, ν c 0 , χ c 0 ) to satisfy the two sets of mean-field equations on both sides of the transition.
The coefficients for the quadratic and quartic terms are found in terms of the microscopic parameters as where we have introduced the rescaled parameters K 0 = K 0 cos 2 χ 0 cos 2 ν 0 ;K 1,2 = K 1,2 sin 2 χ0 2 ;K 3 = K 3 cos 2 χ 0 sin 2 ν 0 ;C = C sin χ0 cos χ0 At the critical point, α L c = 0 and the long-wavelength behavior of the system is controlled by the quartic terms. In particular, the system is stable at the critical point since We have also verified that the excitation spectrum derived from L L does recover the single-particle energy at the critical point in the non-interacting limit.
Numerical study on the shaken lattice model To be more concrete, we study the shaken lattice model numerically. The on-site interaction parameters U are given by where W σ (x) denotes the 2D Wannier functions (WFs) for the σ-band, and g is the contact interaction strength of the bosons. Since we are considering separable static lattice potential V s (x, y) = V 0 (cos(k r x) + cos(k r y)), the 2D WFs can be factorized into products of 1D ones. In particular, as we are only interested in the deep-lattice limit it is justified to approximate the 1D WFs by normalized Gaussians (32) where ξ 0,1 are dimensionless measure of the localization of the WFs, and the 2D WFs are given by W s (x, y) = W 0 (x)W 0 (y); W px (x, y) = W 1 (x)W 0 (y); W py (x, y) = W 0 (x)W 1 (y); W dxy (x, y) = W 1 (x)W 1 (y) We consider a lattice with V 0 = 16E r , for which all microscopic parameters were found numerically. We take ρ 0 = 3a −2 for the numerical calculations and the mean-field phase-diagrams are shown in Supplementary Fig. 2.
The correlation function, plotted in Supplementary Fig. 3, was then numerically found for µ g = gρ 0 = 0.02E r and πf /a = 0.015. At the critical point with δ L c = 0.351E r , the correlation function decays algebraically, similar to the case in the Raman scheme.
A recent work [4] has also studied the phase diagram of a shaken square lattice. Without a microscopic calculation, the effective Lagrangian provided in Ref. [4] does not include the cross term (∂ x ∂ y θ) 2 . Such a term in general exists in the presence of interaction as it renders the system non-separable. Our theory including only onsite interactions does not lead to the term |Φ∂Φ| 2 in Ref. [4], which is irrelevant to our discussions here. Nevertheless, their conclusion that the phase boundary f c increases with increasing interaction agrees with ours.