Strongly interacting spin-orbit coupled Bose-Einstein condensates in one dimension

We theoretically study dilute superfluidity of spin-1 bosons with antiferromagnetic interactions and synthetic spin-orbit coupling (SOC) in a one-dimensional lattice. Employing a combination of density matrix renormalization group and quantum field theoretical techniques we demonstrate the appearance of a robust superfluid spin-liquid phase in which the spin-sector of this spinor Bose-Einstein condensate remains quantum disordered even after introducing quadratic Zeeman and helical magnetic fields. Despite remaining disordered, the presence of these symmetry breaking fields lifts the perfect spin-charge separation and thus the nematic correlators obey power-law behavior. We demonstrate that, at strong coupling, the SOC induces a charge density wave state that is not accessible in the presence of linear and quadratic Zeeman fields alone. In addition, the SOC induces oscillations in the spin and nematic expectation values as well as the bosonic Green’s function. These non-trivial effects of a SOC are suppressed under the application of a large quadratic Zeeman field. We discuss how our results could be observed in experiments on ultracold gases of Na in an optical lattice.


I. INTRODUCTION
Ultracold gases of spin-1 bosons offer an exciting platform to understand the interplay of superfluidity and magnetism [1,2]. Depending on the species of atom, the spin dependent interactions can either be ferromagnetic (as in the case of 87 Rb) or antiferromagnetic (as in the cases of 23 Na), which induces ferromagnetic and polar superfluidity, respectively [3,4]. As a result, in addition to the condensate breaking the U(1) charge symmetry of the system, the SU(2) spin symmetry can also be broken due to the spinful hyperfine interactions. With the development of artificial gauge fields in ultracold atoms, it is now possible to couple the internal hyperfine spin states to their momentum through an engineered spinorbit coupling (SOC) [5,6]. This has now been realized in gases of fermions [7] or bosons [8][9][10][11][12] with a SOC in one and two dimensions [13][14][15][16]. In bosonic gases this induces superfluid order at a non-zero momentum that is dictated by the SOC [17][18][19][20][21][22][23][24][25][26], while symmetry protected topological phases are possible [27,28] and have been observed in Fermi gases [16]. This opens an interesting avenue to explore intertwined order between superfluidity, magnetism, and topology in ultracold gases.
By introducing an optical lattice, both dimensionality and the strength of correlations can be controlled with great accuracy. This allows experimentally realizing phenomena in one-dimension (1D) where strong correlations are significantly enhanced. For example, Luttinger liquid physics has been observed in strongly interacting fermionic gases [29,30]. In addition, 1D SOC is much easier to realize experimentally as compared to its 2D * Electronic address: junhlee@umd.edu † Electronic address: jed.pixley@physics.rutgers.edu analog. This experimental prospect therefore requires a detailed theoretical understanding of the strongly correlated problem. Fortunately in 1D, the existence of powerful analytical and numerical techniques make this understanding possible.
Recently, significant theoretical progress has been made in understanding ultracold gases with a SOC using mean field theories [18][19][20][21][22][23]25] and variational wavefunctions [17,24,31]. One can also view the SOC induced "hopping" between different hyperfine states as a "synthetic dimension" that carry topological edge currents in the superfluid regime [10,25,[32][33][34][35][36], which has also been explored in 1D ladder models [37][38][39][40][41]. In the strongly correlated regime [42,43], the SOC spin-1 Bose-Hubbard model at the odd integer filled Mott lobes can be mapped to an insulating quantum spin-1 magnet in a helical magnetic field which tunes a quantum phase transition [44,45]. However, on the other hand, the strongly correlated superfluid regime of dilute bosons in 1D in the presence of a SOC has not yet received much attention, despite the potentially rich magnetic phenomena due to the spin-1 nature of the problem, which extends beyond the spin-1/2 case [32,46].
Here, we present a comprehensive study of strongly interacting SOC-ed polar superfluidity in 1D. Using a field theoretic framework, we develop the theory for strongly interacting superfluidity in the spin-1 Bose Hubbard model as the complexity of the problem is increased to include a quadratic Zeeman field, a transverse magnetic field, and then finally a SOC. In each case, we verify our theoretical predictions using precise density matrix renormalization group (DMRG) calculations. As a result, we are able to isolate and determine the effect of each of these perturbations on the polar superfluid behavior of spin-1 bosons in a 1D optical lattice. As we show below, in the absence of any perturbing fields the spin-1 Bose-Hubbard model at a fixed dilute filling dis-plays a transition in the excitation spectrum; at sufficiently large interactions the system remains gapless but forms a molecular superfluid phase as the single particle excitations gap out and the two-particle excitations become gapless. We choose to avoid this extreme interaction limit and focus on the experimentally relevant regime with gapless single-particle excitations. First, introducing a quadratic Zeeman field and a transverse magnetic field (which can be considered as a SOC with zero wave vector) in this regime we interestingly find that the spin sector remains quantum disordered in a remarkably robust spin-liquid phase. We determine the Luttinger parameter in the charge sector as well as nematic correlations, which inherit the density-density response due to the gapped spin-liquid sector. In the presence of a full SOC, the ground state displays superfluidity at zero and non-zero momenta concomitant with the existence of a strong coupling charge density wave oscillating at the bosonic particle density. This imprints strong density oscillations in the nematic correlation function and the von Neumann entanglement entropy. Lastly, we discuss how these phases can be observed in experiments on ultracold gases of 23 Na in an optical lattice.
The remainder of the paper is organized as follows: In Sec. II we discuss the model and the DMRG approach we have used. In Sec. III we present the results of our field theoretical analysis and in Sec. IV we verify the physical predictions of the field theory using DMRG. We discuss the implications of our results and their experimental realization in Sec. V. The detailed derivation of the field theory is exposed in Appendix A (effective field theory Hamiltonian), Appendix B (contribution of phase slips), Appendix C (reduction of Luttinger parameter).

II. MODEL AND METHODS
We focus on the spin-1 Bose Hubbard model in the presence of a quadratic Zeeman field and a SOC. This is given by the following lattice Hamiltonian in 1D: where H loc = j g 0 2 :n 2 j : + Here we introduced three component bosonic onsite creation and annihilation operators b † j = (b † x,j , b † y,j , b † z,j ), b j = (b x,j , b y,j , b z,j ) T , the number operatorn j = b † j b j , and the spin operatorŜ where S = (S x , S y , S z ) and (S a ) bc = −i abc are the FIG. 1: Single particle dispersion of spin-1 bosons in the corotating frame for Θ = π/4 and q = 0.2t. A Zeeman field of h = 0.2t opens gaps at the level crossings between different mz states.
spin-1 matrices. This basis is related to the hyperfine eigenbasis by b The parameters of the Hamiltonian contain the hopping strength t that we take as the unit of energy, on site repulsion g 0 > 0, antiferromagnetic spin exchange interaction g 2 > 0, chemical potential µ, quadratic Zeeman field q, and a helical Zeeman field that gives rise to a SOC: Here h is the strength of the helical field and Θ is the pitch of the spiral, i.e., the SOC wave vector. Unless it is explicitly restored, we set the lattice constant a to one. We consider both Θ = 0 where h j induces SOC, and Θ = 0 where h j is merely a transverse field. Upon transforming into a co-rotating frame all terms are invariant except for the helical Zeeman field h j → (h, 0, 0) and the kinetic term The spin-orbit coupling and Θ being its wave vector is manifest in this co-rotating frame. The corresponding single particle dispersion is shifted depending on the hyperfine eigenstate, see Fig. 1. We also define the nematic operatorN ab,j = b j † N ab b j with N ab = δ ab 1 − {S a , S b }/2 to probe nematic order, which shows non-trivial behavior in the polar superfluid phase [47,48].
We solve the Hamiltonian in Eq. (1) with two complementary methods, namely analytical field theoretical calculations and DMRG simulations. The results of the field theoretic analysis are presented below and the detailed derivations are provided in the Appendices. The DMRG directly simulated Eq. (1) on a 1D lattice with open boundary conditions. We focus on a dilute filling of ρ 0 = 1/5 and work in the strong interaction regime so that we can truncate the local bosonic Hilbert space. For all of the results presented here we consider a truncated bosonic Hilbert space to at most two bosons per site. We have checked that in this strongly correlated dilute regime the particle number fluctuations are always small making this approximation very accurate. Note that we also verify that the truncation of the local Hilbert space to two bosons per site in the regime of strong coupling is valid analytically in the derivation of the low energy field theory given below. In the numerical calculations, we used a field strength of h = 0.1t and a SOC wave vector Θ = π/10 on a system of size L = 200 unless stated otherwise. We monitor the convergence of the DMRG by specifying a truncation error of 10 −10 , and a maximum number of 1000 states were kept to obtain the ground state within the truncation error. Lastly, the DMRG calculations are performed using the ITensor library [49].

III. EFFECTIVE FIELD THEORY
In this paper we focus on the strong coupling regime g 0 , g 2 t, where the interaction energies parametrically exceed the bandwidth. For clarity and completeness, we also discuss analytical results in the opposite limit [50] to provide a complete understanding of the problem. The analytical strong coupling calculations [51] are derived in the dilute limit of small superfluid density corresponding to 0 < µ + 2t t and perturbatively in δH [Eq. (1d)].
Both in the strong and weak coupling limits, Eq. (1) in the lab frame maps to the following Hamiltonian density for the bosonic three spinor fields ψ(x) in the continuum The parameters of this theory depend non-trivially on the microscopic parameters of Eq. (1). The weak and strong coupling asymptotes of this functional dependence are compared in Table I and a derivation of Eq. (5) for the strong coupling limit is given in Appendix A. We highlight that onsite eigenstates of Eq. (1c) with up to only two bosons are involved in the derivation; three or more boson eigenstates enter the derivation only for higher order terms. This analytically demonstrates that in the dilute limit the constraint we have imposed on the local Hilbert space in our DMRG studies is a very accurate approximation for Eq. (1). quantity weak coupling strong coupling (1), in the weak coupling limit g0,2 t and strong coupling limit 0 < µ + 2t t g0,2. Within our perturbative calculations h(x) and q are unchanged. Here, we have restored the lattice constant (denoted as a) in order to make both units of energy and length manifest.
A. Low-energy field theory Without symmetry breaking terms (i.e. h = q = 0) the theory displays a [U (1) . This symmetry is spontaneously broken to O(2) on the mean field level, where the field takes the form ψ MF = √ ρ 0 e iϑn (n ∈ S 2 and ρ 0 =μ/g 0 ). The low energy field theory of Goldstone modes is obtained from Eq. (5) in the rotating frame by Gaussian integration of longitudinal, massive fluctuations about the mean field solution. There are two kinds of longitudinal fluctuations: the total density δρ(x, τ ) = ρ(x, τ ) − ρ 0 with gap Λ c = ρ 0g0 and massive spin fluctuations with gap Λ s = ρ 0g2 (see Ref. [50] for details).
Integrating out the massive spin excitations leads to an effective Lagrangian density in the lab frame (LF) It is customary [52] to relabel field integration variables δρ → −φ /π to make the Luttinger liquid nature of the first three terms in Eq. (7a) apparent. The last two terms in Eq. (7a) correspond to a non-linear sigma model (NLσM) in the spin sector. We define the dimensionless stiffness and velocity in charge and spin sector as K c,s = π ρ 0 /[mg 0,2 ] and v c,s = ρ 0g0,2 /m, respectively. In the weak coupling regime, the stiffnesses K c,s are both large, while in the strong coupling regime they can be small (see also Fig. 7 below). For example, as g 0 → ∞ the Luttinger parameter K c → π (ρ 0 a)/2 + (ρ 0 a) 2 − 4(ρ 0 a) 3 . Finally, the second line, Eq. (7b) contains the leading perturbative corrections due to Eq. (1d). The result in the rotating frame (RF), is obtained via Eq. (6) that amounts to the replacements h(x) → h(1, 0, 0), (8) The integral over the superfluid phase ϑ = ϑ smooth + ϑ vortex incorporates both smooth fluctuations and phase slips (i.e. space-time vortices). While the smooth part enters in the form of Eq. (7), the summation over vortex configuration leads to an additional term of the form (see Appendix B) Here, y is the fugacity (Boltzmann weight) of the vortex which is typically ln(y) ∼ −K c . This term is highly oscillatory and produces a contribution to the action that averages to zero unless ρ 0 is an integer. For non-integer ρ 0 , while this term is not relevant in the renormalization group sense, it is responsible for imprinting density oscillations in various observables that are significantly enhanced by a SOC as we demonstrate below.
B. Unperturbed theory: h, q = 0 Before analyzing the implications of the symmetry breaking terms due to the quadratic Zeeman field and the SOC, we briefly discuss the "unperturbed theory," i.e., the spin-1 Bose-Hubbard model defined in Eqs. (1b),(1c) which lead to the Lagrangian Eq. (7a), (11) [51,53]. In the dilute limit 0 < ρ 0 1 of major interest in this paper, the charge sector is a 1D superfluid, i.e., K c > 1. In this case the cosine in Eq. (11) wildly oscillates in real space and is ineffective. Contrary, at integer filling, e.g., ρ 0 ∈ Z, no such oscillations occur and the system undergoes a superfluid to Mott insulating transition as K c drops below 2.
Unlike the various scenarios in the charge sector, the spin sector is always quantum disordered in the absence of symmetry breaking terms (Mermin-Wagner theorem) [53]. Since the primary order parameter ψ does not display off-diagonal long range order but ψ T ψ does, this is an example of quantum vestigial order [54]. The spin-liquid gap ∆ SL in the sigma model part of Eq. (7a) is of order Λ s e −2Ks in the weak coupling regime and of the order Λ s /K s at strong coupling.
Employing the strong coupling parameters of Tab. I, and fixed superfluid density ρ 0 =μ/g 0 , the chemical potential µ(ρ 0 ) drops below the lower band edge −2t when Note that this condition is ρ 0 independent. When g 2 exceeds this line, on-site "molecules" of two bosons with S = 0 form. This can be viewed as the strong coupling . 2: (a-c) One-and two-particle gaps extracted from the ground state energy calculations of the finite size DMRG. The energy gap for the thermodynamic limit is obtained by the extrapolation to 1/L → 0. (a) One-particle gap versus 1/L for g2 = 0.5t, g0 = 5.0t, (b) One-particle gap versus 1/L for g2 = 2.0t, g0 = 1.0t, and (c) Two-particle gap for g2 = 2.0t, g0 = 1.0t. (d) Phase diagram of the unperturbed model in the plane of density-density interaction constant g0/t and spinspin interaction parameter g2/t. The color code represents the numerically obtained single particle gap µc. Below the orange line defined by Eq. (12) µc vanishes: providedn establishes long range-order the system is a nematic superfluid. Above the line, a molecular phase occurs, here the spin-1 bosons form bound states on each site. The green dots represent the position of actual data points and the plot is constructed by interpolation of the data. The naturally occuring ratio of coupling constants for 23 Na (g2/g0 ≈ 1/32) is well in the phase of vanishing single particle gap.
limit of the aforementioned vestigial order. At large g 0 /t, Eq. (12) reproduces the simple condition 2g 2 = g 0 + 2t at which the local 2-boson S = 0 configuration becomes energetically advantageous to the single boson state, cf. Eq. (1c) with µ = −2t (for a table of the eigenstates and energies, see Tab. II in App. A.).
We now verify these expectations using DMRG. In our DMRG calculations we always have a finite-size gap that we can use to determine the nature of one-and twoparticle excitations. To compute the charge gap in the thermodynamic limit and observe the "molecular phase" transition, we calculate the length dependence of the finite-size charge excitation gap in the n-particle sector following Ref. [55], focusing on n = 1 and 2. The finitesize chemical potential to add or remove n particles via the difference in ground state energy is: where E(N, L) is the energy of the system of size L with N particles. The nature of the gaps in the single (n = 1) and double (n = 2) particle sectors follow from the dependence of µ n± (N, L) as a function of L while fixing the density ρ ≡ N/L, as shown in Fig. 2  to µ n± ∼ a n± + b n± /L with a n+ = a n− and have µ n+ − µ n− → 0 as L → ∞. Whereas in the presence of a finite charge gap, we fit the chemical potential [55] to µ n± ∼ a n± + b n± /L 2 with a n+ > a n− and find that µ n+ − µ n− > 0 in the thermodynamic limit. We find that for g 2 g 0 [ Fig. 2 (a)] the system is in a robust superfluid phase with a gapless single particle sector. The single particle excitations become gapped in the opposite limit of g 2 g 0 [ Fig. 2 (b), (c)], however the two-particle (i.e. molecular) excitations remain gapless which is in excellent agreement with the field theoretic analysis. The numerically calculated single particle gap together with the analytical phase boundary are summarized in the g 2 − g 0 phase diagram in Fig. 2 (d).

C. Review: Weak coupling theory
As we explained in the previous section, the weak coupling regime implies an exponentially small gap (∆ SL ∼ Λ s e −2Ks ) in the quantum disordered spin sector. Therefore, the system is very susceptible to symmetry breaking perturbations and a moderate spiral Zeeman field ∆ SL is sufficient to drive the spin sector into an easy plane (spin flop) phase [50]. This is nicely illustrated within weak coupling renormalization group (perturbative in 1/K s and ∆ h ) with flow equations: Here, l is the running logarithmic scale. Depending on the relative magnitude of ∆ h and Λ s e −2Ks the system either flows into a spin liquid phase (for small ∆ h ) or into an easy plane wheren ⊥ (1, 0, 0) [see Fig. 3]. The easy plane model at sufficiently large h contains two different nematic phases and a spin disordered phase in the parameter space spanned by K s and . When the renormalized K s | l=ln(Λs/∆ h ) is larger than 2, there is a direct transition between the two nematic states at = 0. The critical theory is characterized by a spin-charge separated line of pairs of Luttinger liquid fixed points. On the contrary, the transition is indirect as a function of at smaller K s with an intermediate spin liquid state. In the regime where the transition is split, the critical state between nematically ordered and disordered states is a rather exotic c = 3/2 conformal field theory [56][57][58]. It consists of a spin-charge locked pair of a Luttinger liquids in charge space and a Majorana (Ising) critical state in spin space where v c = v s at the critical fixed point. Such a field theory attracted substantial attention recently since it represents a rather simple example of supersymmetric field theories [59] and is related to topological superconductivity [60,61].

D. Strong coupling limit
In contrast to the weak coupling case, K s ∼ 1 at strong coupling. According to the RG estimate from Eq. (14) ∆ h Λ s /K s would be needed to drive the system into the easy plane. However, in this regime the low-energy many body theory [Eq. (7), (11)] is not applicable. Physically, when h is that large, the single particle spin polarizing term in Eq. (1d) is larger than the many-body interaction terms, Eq. (1c). The system is then close to the conventional BEC ground state of fully polarized (i.e., essentially spinless) bosons, instead of being in the vicinity of the spin-nematic BEC.
Since we are interested in the non-trivial regime when many-body effects dominate over h, we assume ∆ h < Λ s /K s and the spin sector is always spin disordered in the remainder of the paper. At time scales beyond 1/∆ SL , it is justified to integrate outn from the low-energy many body theory to obtain an effective Luttinger liquid action of the charge excitations, which is valid at largest length/time scales: Somewhat counterintuitively, spin-orbit coupling (the quadratic Zeeman splitting) has an indirect impact on the charge sector, as it reduces (enhances) the Luttinger parameter K c . In App. C we derive the correction to v c /K c due to the fifth (sixth) term in Eq. (7), proportional to ρ|n | 2 (qρnS 2 zn ). Using a discretization of the field theory on the scale of the coherence length ξ s ∼ v s /∆ SL we integrate gapped fluctuations in the spin sector and obtain where and α, β, γ, δ are non-universal numerical coefficients. Note that the Θ induced suppression of K c can be suppressed when q/∆ SL g 2 /(g 0 K 2 c K s ) (which is still much smaller than unity). We test this prediction in Sec. IV B 3 numerically and find that the SOC wave vector drives a charge density wave by making K c (Θ) < 1.

IV. OBSERVABLES
In this section we determine the consequences of our field theoretic results on physical observables such as the nematicity, entanglement, and correlation functions. We verify this by the finite size DMRG calculations on the SOC S = 1 Bose-Hubbard model in the lab frame [Eq. (1)], which shows excellent agreement with the field theory results. We reiterate the parameter regime of the numerical calculation which is in the strong coupling limit [t g 0,2 ], dilute filling [ρ = 1/5], and we use h = 0.1t, Θ = π/10 on a L = 200 lattice.
A. Effects of homogeneous fields q, h = 0 and Θ = 0 To understand the effect of the symmetry breaking field and the SOC separately, we begin by analyzing the situation without the SOC, i.e. Θ = 0, but with nonzero fields [q, h = 0]. Note that especially h = 0 but Θ = 0 leads to a homogeneous transverse magnetic field [Eq. (2)], and this allows us to build up our intuition for this case before moving onto the effect of a full SOC. The main result for this is that the model remains "stuck" in the spin liquid phase despite tuning the degeneracy lifting quadratic Zeeman and transverse fields, if we stay in the non-trivial regime at which many-body effects are dominant. To demonstrate this we first analyze the nematic order parameter N zz − N yy , which should vanish linearly as q → 0 in the spin liquid phase [50]. In addition, we use the entanglement entropy to determine the number of gapless modes and show that it is independent of the fields. This also substantiates the evidence for the gapped spin liquid phase since the only gapless excitations result from the charge sector of the theory.
The expectation value of the nematicity tensor Nzz − Nyy as function of q at g0 = 5.0t and Θ = 0. The inset is the same plot in a log-log scale to show the power law behavior, which holds for a wide range of g2. The gray dashed line is a guide to the eye which has a slope of 1, indicating Nzz − Nyy ∼ q as q → 0.
We furthermore study various correlators and the Luttinger parameter (of the charge sector) which provides a comprehensive understanding of the model.

Nematicity tensor
While all spin and nematic correlators are short ranged, the presence of a (quadratic) Zeeman field induces a finite expectation value of the nematicity tensor. Even in the spin disordered phase, the linear field h impliesn ⊥ (1, 0, 0) locally. Therefore, the only non-trivial expectation value of the nematicity tensor is N zz −N yy , with N zz + N yy = ρ 0 being fixed by our normalization convention.
In the quantum disordered spin liquid phase, the field theory expectation [50] is that N zz − N yy ∼ q, since in any (quantum or thermally) disordered phase the expectation value of the order parameter vanishes linearly as a function of its conjugate variable. In Fig. 4 we numerically demonstrate this behavior for a number of parameters quite clearly. This serves as a strong numerical evidence for the system robustly remaining a spin-liquid in the presence of the fields.

Entanglement entropy
Another evidence for the gapped spin-liquid would be added if we can observe the nonzero spin gap. An indirect method to detect the gap is by counting the gapless modes, or calculating the central charge, of the system. For the current case of Θ = 0 we have an algebraically ordered superfluid in the charge sector, i.e., a Luttinger liquid with K c > 1, which is known to contribute central charge c = 1 to the system. Considering The entanglement entropy as a function of the bipartite position x, for g0 = 5.0t, g2 = 1.0t, and Θ = 0. The inset is the same data with the horizontal axis as 1 6 ln L π sin πx L , the slope determines the central charge per Eq. (19). The gray dashed lines are guide to the eye which corresponds to c = 1 result while explicit linear fits gave c = 1.002, 1.003, and 1.004 for q = 0.0t, 0.1t, and 1.0t respectively. the spin sector, the system is clearly in a gapped spinliquid phase in the unperturbed regime (h, q = 0) with the gap ∆ SL ∼ Λ s /K s . Due to this spin-liquid gap, the spins do not contribute to the central charge and thus the total central charge will be c = 1. And if the system remains in this gapped spin-liquid after turning on the fields (h, q = 0), the central charge will as well remain c = 1.
To extract this numerically, we analyze the dependence of von Neumann entanglement entropy S(x) as a function of the position x of the bipartition. We fit S(x) to the well known result from conformal field theory with open boundary conditions [62][63][64] where c is the central charge and d is a nonuniversal constant. We calculate this for a number of parameters in Fig. 5 together with the fit to the form of Eq. (19).
We consistently obtain c ≈ 1 up to q in the order of t which adds strong evidence that the ground state of the model across this parameter regime has a spin-sector that remains in a gapped spin-liquid phase.

Bosonic correlators and Luttinger parameter
We now investigate the correlation functions at Θ = 0 and study the behavior of the Luttinger parameter. We begin by studying bosonic correlators G 0 (x) = b z,j+x b † z,j of the spin-0 projection of the bosonic field (note that S z (0, 0, b z,j ) T = 0) as well as the total Green's function G tot (x) = tr[ b j+x b † j ]. Since our numerics have open boundary conditions, we calculate the correlation functions at the center of the chain and set j = L/2 to avoid boundary effects as much as possible. Since the spin sector is gapped, the correlators are dictated by the Luttinger liquid charge sector and behave as: for both α = 0, tot. This power-law behavior is demonstrated in Fig. 6 for various parameters. From the power-law fit of the correlations we extract the Luttinger parameter K c for various values of g 0 , g 2 , and q. We find K c has a very weak dependence on q, whereas the g 0 and g 2 dependence is prominent. This is in agreement with the analytical results, according to which the q dependence enters only via weak fluctuation corrections, Eq. (16). To understand this we make a comparison between the numerics and the analytical expectation K c = π ρ 0 /[mg 0 ] for q = 0, based on Tab. I. The result is presented in Fig. 7, which demonstrates good qualitative agreement between the two. Moreover, this shows K c > 1 in the wide parameter regime of strong coupling.

Nematic correlators
We now turn to the behavior of the nematic correlation function, and compute the connected correlation function that is defined as C zz (x) = : N zz (x) :: Within the field theory description, we can understand the individual charge and spin contributions to C zz (x) by introducing source fields q → q + δq(x) in Eq. (7), (11) and appropriately differentiating with respect to δq(x), before taking the limit δq(x) → 0 at the end. This generates a vertex ρ(x)n(x)S 2 zn (x) whose relation to C zz (x) is given by: Thus, the nematic correlation function receives a contribution from the spinn(x) and the charge ρ(x) sectors of the field theory.
Since the spin sector is gapped, integrating outn we obtain: where δρ(x) = ρ(x) − ρ 0 and C is a non-universal constant C ∝ y 2 [see Eq. (11) and below for definition of y]. This result has important implications from the value of the Luttinger parameter. When the Luttinger parameter obeys K c > 1 as in Θ = 0 [ Fig. 7], the asymptotic power law regime for x 1 is dominated by the 1/x 2 contribution, while the second term stemming from Eq. (11) is subdominant and only generates weak oscillations in the amplitude. On the other hand, K c < 1 implies the oscillatory second term dominates C zz and thus the ground state will be in a charge density wave state with a wave The Luttinger parameter Kc in the charge sector for g2 = 0.5t and g2 = 1.0t obtained from numerical calculations, while g0 = 5.0t is fixed. The dashed curves denote the plot for q = 0 and Kc in the strong coupling limit from analytical calculations. Note that Kc > 1 for a wide range of parameters in this limit.
vector Q CDW = 2πρ 0 . For Θ = 0, we calculate C zz in Fig. 8 which shows a 1/x 2 power-law decay, consistent with K c > 1 from the previous section. A q-independent weak oscillation with a wave vector 2πρ 0 is also apparent from the data.
B. Effect of the spin-orbit coupling: h, q, Θ = 0 We move on from the spatial uniform transverse magnetic field and now consider the effect of a SOC on the strong coupling superfluidity of polar spin-1 bosons, by considering Θ = 0. As we consider a regime with a robust gapped spin-liquid phase, the physics with the SOC is very rich and a correlated charge density wave state also appears. In the strong coupling regime under investigation, the behavior is entirely due to the powerlaw decay of density-density correlations, see Eq. (22). In the present case of ρ0 = 1/5, Kc is larger than 1 and the long range asymptotics is given by Czz x −2 , which is in good agreement with the numerical calculation.

Spin and nematic texture
As in the transverse field case, we start our analysis with the spin and nematic expectation values. In the lab frame, the average spin component will try to locally anti-align with the magnetic field along the chain. As a result of the finite SOC (with wave vector Θ = π/10 in the DMRG), the spin expectation values S x (x) and S y (x) follow the pattern of the helical magnetic field. Explicit forms are given by: A x and A y are the amplitudes for each spin expectations. This functional form can be understood analyti- Spin expectation values (a) Sx(x) and (b) Sy(x) for g0 = 5.0t, g2 = 0.5t, Θ = π/10, and two values of q = 0.0t, 1.5t. Note the oscillation wavelength 2π/Θ is indicated as an arrow between two maxima, and the nonzero q suppresses the oscillation. cally by considering the transformation of theŜ x andŜ y operators from the lab frame to the rotating frame using Eq. (6). On the other hand, the spin component perpendicular to the field is suppressed due to g 2 > 0 and we find S z (x) ≈ 0. In Fig. 9 we show plots of S x (x) and S y (x) for two different values of q, showing oscillations at the wavelength of the SOC. The oscillation is suppressed by the quadratic Zeeman field as expected. Upon the unitary transformation from the lab to rotating frame of the bosonic operators, the nematicity tensors of components N xx (x) and N yy (x) pick up a contribution from the spatially dependent phase factor that is not present for Θ = 0. The functional forms are obtained as with amplitudes A, B, and phase φ. On the other hand, N zz (x) remains invariant under the unitary transformation to the rotating frame and does not acquire any oscillatory behavior due to the SOC. Rather, the oscillations occur from the charge density modulation with a wave vector of 2πρ 0 : This can be understood by considering the N zz (x) being generated through the vertex ρ(x)nS 2 zn [see Eq. (21)]. In Fig. 10  N yy (x) and N zz (x) . We observe oscillations, which are suppressed with the quadratic Zeeman field, with different wavelengths originating from the SOC and charge density, respectively.
To determine whether the model remains in the spinliquid phase we use these functional forms to extract an estimate of the difference of the nematic expectation values N zz − N yy . However, since they both oscillate at different periods we first fit the data to the functional forms given in Eq. (24), and then determine N zz − N yy via the following procedure: we evaluate N zz by averaging N zz (x) over the lattice (we exclude some sites at the boundary during averaging to avoid boundary effects), we extract A yy from the fit of N yy (x) to the functional form above and use N zz − A yy as a proxy for N zz − N yy . We expect that N zz − N yy vanish linearly in the spin liquid regime like N zz − N yy ∼ [50], where = q + Θ 2 /2m. As shown in Fig. 11 we find good agreement with this vanishing linearly with , however due to the oscillation periods being distinct this leads to a non-perfect estimate of N zz − N yy and shifts the zero away from = 0.

Entanglement entropy
We again look at the entanglement entropy and calculate the central charge for additional evidence of the spin gap. As shown in Fig. 12, we find that entanglement entropy is very weakly affected by a quadratic Zeeman field and obtain a central charge c ≈ 1 from the linear fit of S(x) versus log L π sin πx L [see Eq . 19], which is in excellent agreement with the expectation that the spin sector remains gapped and the only gapless modes are due to the superfluidity in the charge sector. This also is in agreement with our results for Θ = 0 [Sec. IV A 2], thus we conclude the model remains in the spin-liquid FIG. 12: The entanglement entropy as a function of the bipartite position x for g0 = 5.0t, g2 = 0.5t, and Θ = π/10. The inset is the same data with the scaled horizontal axis to obtain the central charge from the slope (Eq. (19)). The gray dashed lines are guide to the eye which correspond to c = 1, while explicit linear fits gave c = 0.995, 1.002, and 1.004 for q = 0.0t, 0.1t, and 1.5t respectively.
phase even in the presence of a full SOC.
However, comparing with the case of Θ = 0, we find that the oscillations in the entanglement entropy are much larger in the case of nonzero SOC. These oscillations occur with a period given by 1/ρ 0 and are thus due to the oscillation in the charge density. As we demonstrate below, the SOC induces a charge density wave of period 1/ρ 0 due to the Luttinger liquid in the charge sector having K c < 1. [See also Eq. (22) and the discussion below]

Bosonic correlators and Luttinger parameter
We now turn to the bosonic Green function of each spin state, see Fig. 13. For the spin-0 component these are given by G 0 (x) = b † 0 (x)b 0 (0) and the spin-(±1) component of the bosonic correlator is G ±1 (x) = b † ±1 (x)b ±1 (0) . Applying the transformation from the rotating frame to the lab frame allows us to deduce the functional form of G α (x). Since the spin-0 component is unaffected by this transformation, the form remains as in Eq. (20): The DMRG results for G 0 (x) are presented in Fig. 13(a) and we extract the Luttinger parameter K c from a fit to the power-law form. Interestingly, distinct from the case with Θ = 0, we now find that K c strongly depends on the quadratic Zeeman field. As shown in Fig. 14, our data fits remarkably well to a simplified variant of the field theoretical result Eq. (16) with two fitting parameters A sin 2 (Θξ s ) and B.
In contrast to the spin-0 Green function, the spin-(±1) components do alter as we transform to the lab frame This suggests that the power-law form is identical to the spin-0 case but it acquires an oscillatory component due to the SOC, consistent with the data shown in Fig. 13 (b). To demonstrate this, we first extract K c from G 0 (x) using Eq. (26) and then plot x 1/(2Kc) G +1 (x) in the inset, which does not decay and oscillates with a period 2π/Θ thus confirming the functional form in Eq (28). Lastly, the positive quadratic Zeeman field strongly suppresses G ±1 (x) as expected.
The extracted Luttinger parameter in the charge sector K c for a finite Θ as a function of the quadratic Zeeman field is given in Fig. 14 for various values of g 2 . This demonstrates that the finite SOC leads to K c < 1 in small q, which induces a charge density wave state due to the functional form of the charge correlation function [see Eq. (22) and the discussion below]. By applying a large quadratic Zeeman field, the effect of SOC and thus the charge density wave is suppressed inducing a crossover from K c < 1 to K c > 1. The proximate charge density wave regime is the reason that the Luttinger parameter is so sensitive to tuning q in contrast to the limit of Θ = 0. The charge density wave can be clearly seen in the nematic correlation function, which we now turn to.

Nematic correlators
As a result of the SOC driving K c < 1, we expect that the nematic correlation function in Eq. (22) is dominated by the oscillating term with a power law given by 2K c . We demonstrate this by plotting the nematic correlator C zz (x) for a number of different values of q in the presence of SOC in Fig. 15(a). For q = 1.0, we can check from Fig. 14 that K c > 1 and C zz (x) show similar behavior as in Θ = 0 case. However, as we decrease q to the regime where K c < 1 in Fig. 14, we find that oscillations enhance as well as the power of the decay changes. If we use the K c value extracted from G 0 (x) [Fig. 14] to Eq. (22), we find excellent agreement between the numerics and the functional form, which is shown in Fig. 15(b). This also confirms the emergence of a charge density wave from SOC with the wave vector Q CDW = 2πρ 0 . Thus, we reach one of our main conclusions: In the presence of large interactions a SOC induces a strong coupling charge density wave phase in dilute polar superfluids.

V. DISCUSSION
In summary, we have presented a combined numerical and analytical study of polar spin-1 lattice bosons at non-integer filling in one dimension under the influence of spin-orbit coupling and quadratic Zeeman field. Complementary to the previous study at weak coupling [50], we here concentrated on the limit when interaction effects are stronger than the kinetic energy. Our main finding, which is supported by the excellent agreement between analytics and numerics, is that in this regime the spinliquid gap is substantial and therefore the perturbative inclusion of symmetry breaking terms is insufficient to restore the algebraic nematic order. At the same time, the robustness of the spin-liquid phase does not render the spin sector entirely innocuous: we have demonstrated that spin-orbit coupling is capable of tuning the charge sector into a charge density wave by reducing the Luttinger parameter K c below unity. A qualitative explanation of this reduction of K c may be understood in the limit of large helical background magnetization h and negligible quadratic Zeeman field q. We first discuss this limit in the case of vanishing spinorbit wave vector Θ = 0. Then, only the b x boson is of importance and our model displays conventional BEC of spinless bosons. We repeat that the superfluid stiffness is K c ∼ √ tρ 0 a/ √μ . The first factor accounts for the intuitive increase in stiffness with increasing hopping strength while the remaining factors stem from the on-site mean field solution and are independent of the kinetics. Now we return to Θ = 0, in the presence of the such a SOC the BEC has a rotating on-site polarization. Therefore, the overlap of neighboring single-particle wavefunctions of adjacent sites is substantially weakened due to the spin dependent hopping and the numerator in K c is reduced.
Lastly, we conclude with the experimental realization of our strong coupling theory using ultracold gases of the polar spin-1 boson 23 Na. A natural generalization of the experimental setup in Ref. [65] by including a onedimensional optical lattice, should be able to straightforwardly realize the spin liquid phase we have uncovered in the limit of no spin orbit coupling in Sec. IV A. The ability to tune the quadratic Zeeman field across the nematic transition in the weak coupling limit implies such a transition can also be studied here. A clear cut signature of the spin liquid regime would be given by the difference in nematic expectation values [66] vanishing linearly with decreasing quadratic Zeeman field (as in Fig. 4). The realization of our newly discovered strong coupling charge density wave phase that is induced by spin orbit coupling is in principle also possible within existing experimental setups. However, it requires long coherence times for 23 Na atoms in miscible F = 1 hyperfine states -a requirement which so far has been challenging due to strong magnetic noise. We are hopeful that the most recent experimental breakthrough in shielding techniques has overcome this bottleneck [67]. Thus, we expect that a spin-orbit coupling can be induced in polar spin-1 bosons in the near future and the non-trivial predictions of our theory can be tested. In particular, the strong coupling charge density wave can be observed either directly, through measuring the charge response via single-site imaging techniques [68,69] and Bragg scattering [70,71], or indirectly, using nematic tensor components [66]. a. Solution of local problem and molecular phase The weak coupling limit of Eq. (5) follows trivially from the continuum limit of Eq. (1). Therefore, this section focuses on the strong coupling limit, where we perturb about local eigenstates, Tab. II. To determine the latter, note that :Ŝ 2 :=Ŝ 2 − 2n and :n 2 :=n(n − 1). Eigenvalues follow fromn → n andŜ 2 → S(S + 1) for conserved quantum numbers n and S. The structure of eigenstates follows from b. Derivation of effective continuum field theory in the strong coupling limit As a first step, we decouple the hopping term This also demonstrates that the matrix t is positive definite in the infrared limit of interest [72].
The overall strategy is to derive an effective action for ψ = Ψ/(E 1 √ a). To this end, we express the non-local The matrix elements of δH = b † δhb, of Eq. (1c) in the single particle sector are obviously given by the matrix form of δh.

c. Effective Action
We begin with the derivation of the effective action S[ψ] by focusing only on quadratic terms of the kind In the limit E 1 /T → +∞, this leads to In addition to the conventional time derivative termψψ there is a term with two derivatives. However, in the interesting regime of time scales τ Λ c 1 it is suppressed and henceforth omitted.
We now determine all other static terms in Eq. (5). To this end, it is sufficient to consider time independent field configurations. The bare partition function is Z = 1 + 3e −βE1 + e −βE2,0 + 5e −βE2,2 . (A9) We will consider sufficiently large E 1 /T, E 2,0 /T, E 2,2 /T → ∞, and will only keep the contribution of occupied states if the contribution of empty states vanishes. We will further use the following identities: We obtain the following perturbative correction to the ground state energy Restoring slow time dependence of fields and dτ E = δS leads to the remaining terms in Eq. (5). Collecting all terms and rescaling Ψ → ψ leads to the identification of parameters of the field theorỹ which leads to the expressions in Tab. I of the main text (the leading order in µ/t+2 1 is presented there).

Appendix B: Phase slips
In this appendix we derive the effective action of phase phase slips, Eq. (11). It is sufficient to consider the first three terms of Eq. (7a) for the sake of this derivation. As mentioned in the main text, we introduce the field φ by means of δρ = −φ /π. Amongst all boundary conditions of the fields, the important one is ϑ(x, β) = ϑ(x, 0) + 2πf (x) where f (x) ∈ Z∀x is a piecewise constant function. In order to introduce vortices in ϑ we split ∂ µ ϑ = ∂ µ ϑ reg + A µ , where the gauge potential accounts for vortices µν ∂ µ A ν = 2π i n i δ( x − x i ). It is convenient to choose a "Landau" gauge in which A µ = (0, 2π i n i δ(τ − τ i )θ(x − x i )). Note that, contrary to usual Berezinskii-Kosterlitz-Thouless physics, also non-neutral configurations i n i = 0 are consistent with the periodic boundary conditions and kept. To avoid double counting, we keep only vortices of n i = ±1 but allow them to sit on top of each other (i.e. effectively creating double vortices). Furthermore, avoiding double counting also implies that we do not count permutations of equivalent sets of vortex positions { x i } twice.
We then obtain the Lagrangian as L reg + δL where We have dropped the reg subscript in ϑ reg . The total amplitude in a sector of a total of N = n +n vortices, where n (n) is the number of vortices with positive (negative) winding is The combinatorial factor is the number of possibilities (n+n)! n!n! to arrange n vortices with positive winding if there are n +n vortices in total divided by the number of configurations with equivalent spatial ordering (n +n)!. The Boltzmann weight of a vortex is denoted y. Summation over N leads to a Dirac function δ(θ), so that the overall theory is given by the effective Lagrangian, Eq. (11).
For the calculation of density correlators perturbatively in y we may integrate ϑ and obtain In this part of the appendix, the index c in K c , v c is suppressed.
x → x /∆x and This implies for the density correlator [52] ρ(x)ρ(0) = ρ 2 0 + K 2π 2 x 2 + const. In this appendix we determine the SOC induced corrections to K c /v c in the strong coupling limit, Eq. (16) and (17). Since the NLσM sector is gapped one may integrate out the spin sector, and the term ρ|n | 2 /2m leads to additional terms (φ ) 2 i.e. to a renormalization of K c . Here we estimate these terms by evaluation of |n (x, τ )| 2 |n (x , τ )| 2 .
In view of the short range correlations in spin space,n decays on the scale ξ s and we discretize the field theory in segments of length ξ s . The spin sector of the Goldstone theory, Eqs. (7a), (7b), is then In Hamiltonian formulation [72], the time derivative term becomes H top = 1 2I i L 2 i , i.e. it is a sum over tops with moment of inertia I = ξ s /g 2 . The energy levels l(l + 1)/2I have eigenstates given by spherical harmonics Ω|l, m = Y m l (Ω), where Ω is the solid angle parametrizing the target manifold of the sigma model.