Interacting Bose gases in quasi-one dimensional optical lattices

We study a two-dimensional array of coupled one-dimensional (1D) tubes of interacting bosons. Such systems can be produced by loading ultra-cold atoms in anisotropic optical lattices. We investigate the effects of coupling the tubes via hopping of the bosons (i.e. Josephson coupling). In the absence of a periodic potential along the tubes, or when such potential is incommensurate with the boson density, the system undergoes a transition from an array of incoherent Tomonaga-Luttinger liquids at high temperature to an anisotropic Bose-Einstein condensate (BEC), at low temperature. We determine the transition temperature and long wave-length excitations of the BEC. In addition to the usual gapless (Goldstone) mode found in standard superfluids, we also find a gapped mode associated with fluctuations of the amplitude of the order parameter. When a commensurate periodic potential is applied along the tubes, they can become 1D Mott insulators. Intertube hopping leads to a deconfinement quantum phase transition between the 1D Mott insulators and the anisotropic BEC. We also take into account the finite size of the gas tubes as realized in actual experiments. We map out the phase diagram of the quasi-1D lattice and compare our results with the existing experiments on such systems.


Introduction
One dimension (1D) provides fertile grounds for studying the physics of strongly correlated quantum many-body systems. It is a well established theoretical result that, because of the enhanced role of quantum fluctuations in low dimensionality, there is no broken continuous symmetry and therefore long range order, even at zero temperature, in a 1D quantum many body system. Instead, in a large class of systems possessing a gapless spectrum, power laws characterize the various correlation functions. This class of systems is known as a Tomonaga-Luttinger liquid (TLL).
One of the simplest examples of a TLL is the Lieb-Liniger model [1], which describes a system of non-relativistic bosons interacting via a contact (Dirac delta-function) potential. This model has been solved exaclty using the Bethe ansatz [1], which has provided valuable insights into its thermodynamic properties and excitation spectrum. However, calculation of the asymptotics of correlation functions is not an easy task using this method, and in this regard it is nicely complemented by the harmonic-fluid approach introduced by Haldane [2]. The combination of these two methods has allowed us to understand that, as the the interaction is increased, the system exhibits no quantum phase transition. Instead, it smoothly crosses over from a weakly interacting (often called quasi-condensate [3,4,5]) regime, where phase fluctuations decay very slowly, to a strongly interacting regime, the Tonks regime [6]. In the latter, repulsion between the bosons is so strong that it leads to an effective exclusion principle in position space, and makes the system resemble a gas of non-interacting spinless fermions in many properties. Thus it is sometimes stated that the bosons in the Tonks regime "fermionize". This very simple model solved by Lieb and Liniger in 1963 has recently received a great deal of attention thanks to spectacular advances in confining quantum degenerate gases of alkali atoms in low dimensional traps [7]. While at the beginning only weakly interacting 1D Bose gases were available [8,9], the application of a deep optical potential parallel to the 1D axis to arrays of 1D Bose gases has recently allowed the group at ETH (Zürich) [10,11] to demonstrate a quantum phase transition from a superfluid TLL to a 1D Mott insulator, which is the low dimensional analog of the superfluid to Mott insulator transition first demonstrated in the 3D optical lattice by Greiner and coworkers [12]. Furthermore, Kinoshita et al. [13] have succeeded in reaching the Tonks regime by creating an array of tightly confined 1D Bose gases. A different path [14,15] to realize a strongly interacting 1D Bose gases has been followed by Paredes et al. [16], who applied a deep periodic potential to an array of (otherwise weakly interacting) 1D Bose gases. Thus they were able to dramatically increase the ratio of the interaction to the kinetic energy of the atoms. In both cases [13,16], the observation of fermionization signatures was reported.
In the past, much experimental effort focusing on low dimensional systems has been spent on studying electronic systems such as spin chains and quasi-1D metals (for a review see [17]). Theoretically at least, electrons in 1D become a TLL (provided there are no gap-opening perturbations) as soon as interactions are taken into account in a non-perturbative way. However, for the electronic systems comparison between theory and experiment has proved problematic partly because of material-specific issues such as the existence of many competing phenomena at similar energy scales (e.g. in Bechgaard salts [18,19]) or the always-present disorder in solid state systems [17]. Furthermore, one major complication stems from the fact that these materials are indeed three dimensional (3D), as there is always a small amplitude for electrons to hop from chain to chain, whereas Coulomb interactions (often poorly screened because of bad metallicity) also enhance 3D ordering tendencies. These inter-chain couplings cannot be neglected at low temperatures, which leads to the interesting physics of dimensional crossover : at temperatures high compared to the energy scale determined by these inter-chain couplings, the system effectively behaves as a collection of isolated 1D chains. However, when the temperature is lowered below this scale, 3D coherence develops between the chains and, at zero temperature, a complete 3D description of the system is required. Much exciting physics is expected to occur in this crossover from the exotic 1D TLL state to a conventional 3D (but strongly anisotropic) Fermi liquid, or an ordered 3D phase [17]. Indeed, these phenomena may well have relevance to the pseudogap (and other anomalous) behavior of the quasi-2D high temperature superconductor materials. For instance, one of the many interesting questions that arise in this context is the following: how [19] do the exotic low-energy excitations (bosonic collective modes) of the 1D system transform into the Landau quasi-particles characteristic of Fermi liquids and which resemble more the constituent electrons?
The dimensional crossover phenomena described above has also a counterpart in Bose gases loaded in quasi-1D optical lattices. From a theoretical point of view, the fact that bosons can undergo Bose-Einstein condensation and, therefore, can be collectively described using a macroscopic quantum field makes more amenable these systems to a theoretical analysis relying on a mean field approximation. Such a mean-field analysis is undertaken in this work: here we treat the intertube hopping via a mean field approximation, while quantum fluctuations within a tube are treated accurately at low energies using the harmonic-fluid approach.
From the experimental point of view, the recent availability of gases of ultracold bosonic atoms in optical lattices has revived the interest in understanding the properties of coupled systems of bosons in 1D. Furthermore, the already demonstrated high degree of tunability of these systems can provide a clean setup to study the physics of dimensional crossover. Thus, as we show in section 3 of this work, the phenomenology of Bose-Einstein condensation in the most anisotropic version of these lattices is very different from the one exhibited by Bose gases condensing in the absence of a lattice or even in shallow isotropic optical lattices.
In addition to studying the properties of Bose condensed systems in lattices, some experimental groups have applied a tunable optical potential along the longitudinal direction of an array of effectively uncoupled 1D Bose gases. It has been thus possible [10], as mentioned above, to drive a quantum phase transition from the 1D superfluid (TLL) state to a 1D Mott insulator. Under these circumstances, the effect of a small hopping amplitude (i.e. Josephson coupling) between neighboring tubes in the array becomes most interesting: as we show in section 4, for a finite value of the hopping amplitude, the system undergoes a new type of quantum phase transition, known as "deconfinement" transition, from the 1D Mott insulator to an anisotropic Bose-Einstein condensate. Thus the tendency of bosons to develop phase coherence overcomes the localization effect of the optical potential. As explained in section 6, we believe that some features of this deconfinement transition may have been observed in the experiments of Ref. [10].
This paper is organized as follows: section 2 defines our model of a 2D array of 1D tubes of bosons coupled by weak intertube tunneling, and describes the harmonic-fluid approach that allows us to derive a low-energy effective Hamiltonian. Then we discuss separately the phenomenon of Bose-Einstein condensation and the properties of the condensate in an anisotropic lattice (section 3), and the deconfinement quantum phase transition (section 4). In the former case, we obtain the zero-temperature condensate fraction, the condensation temperature, and the excitations of the condensate using the mean field approximation and including gaussian fluctuations about the mean field state. We also show that some of the results obtained by these methods are recovered from the self-consistent harmonic approximation (section 3.2). These calculations are described in pedagogical detail in Appendix A, Appendix B, and Appendix C. Moreover, in section 3.1.3 we discuss the excitation spectrum of the Bose-condensed phase, and in particular the existence of a gapped mode related to oscillations of the amplitude of the condensate fraction, in the light of a phenomenological Ginzburg-Landau theory. Differences between this theory and the more conventional Gross-Pitaevskii theory, where the gapped mode is absent, are exhibited. The transition from the (anisotropic) regime described by the Ginzburg-Landau theory to the one where Gross-Pitaeskii theory holds is also outlined there.
The deconfinement transition is described in section 4. There we show how the renormalization group can be used to understand the competition between the localization in the optical potential and a small tunneling amplitude between neighboring 1D Bose systems. To assess the nature of the deconfinement transition, we map (see section 4.2) the mean-field Hamiltonian, at a specific value of the system's parameters, to a spin-chain model. Finally, section 5 discusses the effect of the finite size of the 1D system that are realized in actual experiments, and how it modifies the phase diagram obtained from the renormalization-group analysis. We argue that, for an array of finite 1D Bose gases, there is a extra energy scale below which the tunneling of atoms between tubes is effectively blocked. Finally, in section 6 we present a discussion of the experimental observations and how they compare with the predictions of our theory. A brief account of some of our results has appeared elsewhere [20].

Models and Bosonization
In this section we introduce the models for an interacting Bose gas in a quasi-1D optical lattice that we study in the rest of the paper. We also briefly describe the harmonicfluid approach (henceforth referred to as "bosonization"), which allows us to obtain a convenient effective low-temperature description of these models.

Models of a Bose gas in a quasi-1D optical lattice
In this work we study the low-temperature equilibrium properties of systems recently realized by several experimental groups [7,13,10,16,21,22,23] which have succeeded in loading ultracold atomic gases in highly anisotropic (quasi-1D) optical lattices. In the experiments carried out by these groups, a Bose-Einstein condensate of ultracold alkali atoms (in particular, 87 Rb) is placed in a region of space where two pairs of non-coherent and mutually orthogonal counter-propagating lasers intersect. The intensity of the lasers is increased adiabatically, and because of the AC stark effect, the atoms experience a periodic potential [24] V ⊥ (y, z) = V 0y sin 2 (2πy/λ) + V 0z sin 2 (2πz/λ), where usually V 0y = V 0z = V 0⊥ , with V 0⊥ ∝ E 2 0 , E 0 being the maximum intensity of the electric field of a laser whose wavelength is λ. The strength of the optical potential, V 0⊥ , is frequently given in units of the recoil energy, E R =h 2 π 2 /(2Ma 2 ), where a = λ/2 is the period of V ⊥ (y, z) and M the atom mass ‡. For V 0⊥ ≫ µ 3D , where µ 3D ≈ 4πh 2 a s ρ 3D 0 /M [25] is the chemical potential of the BEC (a s being the s-wave scattering length of the atomatom interaction potential, and ρ 3D 0 the mean atom density), the gas atoms are mainly confined to the minima of V ⊥ (x, y), thus forming long, effectively 1D, gas tubes as shown in figure 1 (there are additional harmonic potentials along the x, y, z-axes that trap the atoms). Under these conditions, we can project the boson field operator Ψ(x, y, z) onto the lowest (Bloch) band of the periodic potential, where R = (m y , m z ) a, with integers m i = −M i /2, . . . , +M i /2 (i = y, z), is the location of the different potential minima; W R (x, y) is the Wannier function of the lowest Bloch band centered around R. The system just described is called a 2D (square) optical lattice. In typical experiments M y ∼ M z ∼ 10 2 , and thus the number of gas tubes is of the order of a few thousand separated half a laser of wave-length, λ/2 ≃ 0.4 µm.
In addition, in several of the experiments [10,22,16] an extra pair of lasers was applied to define a longitudinal periodic potential V || (x) = V 0|| sin 2 (2πx/λ) = 1 2 V 0|| + u(x), and u(x) = u 0 cos(Gx), where u 0 = −V 0|| /2 and G = 4π/λ = 2π/a. Since, for a commensurate potential, i.e. when the period a matches the mean interparticle distance, the variation of the strength of u(x) is responsible for the transition [2,26,17] to a Mott insulating state where the bosons are localized, we shall henceforth refer to ‡ For 87 Rb, E R /h = α/λ 2 , where α ≃ 14.424 µm 2 kHz and λ is the laser wavelength in µm. Thus, for λ = 0.826 µm used in the experiments of Ref. [10], E R /h = 21.14 kHz L R' R Figure 1. Sketch of a two-dimensional array of coupled tubes (also called quasi-1D optical lattice in this work). Bosons can hop from tube to tube with amplitude J, which can be controlled by the height of an optical potential transverse to the tube axis. Along the tube, a periodic optical potential (for clarity, shown here for only one of the tubes) can be also applied. This potential can drive a quantum phase transition where each tube transforms into a 1D Mott insulator.
u(x) as the "Mott potential". In this work we are interested in the extreme anisotropic case where the strength of this potential |u 0 | ∼ V 0|| ≪ V 0⊥ . In experiments, this was achieved by increasing V 0⊥ so that the motion of the atoms effectively becomes confined to a one-dimensional tube, where they undergo only zero-point motion in the transverse directions. This picture is equivalent to the projection onto the lowest Bloch band described above. Just like in three dimensional space, the atomic interactions in the tube are still well described by a Dirac-delta pseudo-potential, v a (x) = g 1D δ(x), where the effective coupling g 1D has been obtained in Ref. [27]: The constant C ≃ 1.4603, and the oscillator length in the transverse direction, ℓ ⊥ = (h/(Mω ⊥ )) 1/2 (wherehω ⊥ /E R ≃ 2(V 0⊥ /E R ) 1/2 for deep lattices [12], is the transverse oscillation frequency).
Taking into account the above considerations, we are led to consider the following Hamiltonian: where , and commutes otherwise as corresponds to bosons. The last term in equation (3) describes the hopping of bosons between two neighboring tubes. This term is also known as Josephson coupling. The hopping amplitude [28]: for V 0⊥ ≫ E R . We shall assume a very deep transverse potential so that J ≪ µ 1D , where µ 1D = µ 1D (ρ 0 ) is the chemical potential in each tube, as can be obtained from the exact solution of Lieb and Liniger [1] (or measured experimentally). Since experimentally µ 1D < 0.1hω ⊥ [7], assuming J ≪ µ 1D requires at least that V 0⊥ > 10 E R . Such limit allows us to regard the system as an array of weakly coupled 1D interacting Bose gases.
In the experimental systems, the Mott potential u(x), besides the periodic part described above, also contains a slowly varying harmonic term that confines the atoms longitudinally. Additionally, there are also harmonic traps in the transverse directions. This makes the array, not only finite, but also inhomogeneous. Whereas we shall discuss finite-size effects further below (see section 5), it has proven so far intractable in 1D (except for a few limiting cases [29]) to deal with the harmonic confinement analytically and only numerical results are available [30,31,32,33,34,35]. Thus, with the exception of section 7 where we discuss the experimental consequences of our work, we shall neglect the effect of the harmonic potential, and treat the system as an array of gas tubes of length L and hard-wall box (open) boundary conditions [36,37,38], each one containing N 0 = ρ 0 L bosons (in typical experiments, N 0 ranges from a few tens to a few thousand particles). This makes sense, however, as much of the interest of this work is focused on computing the phase diagram as well as some of the thermodynamic properties of the different phases, for which one is required to take the thermodynamic limit L → ∞ at constant density ρ 0 = N 0 /L.
Finally, when the axial potential u(x) is made so deep that |u 0 | ≫ µ 1D (but still |u 0 | ≪ V 0⊥ ), it becomes convenient to perform a further projection of Ψ R (x) onto the lowest Bloch band of the Mott potential u(x). Thus, where [b mR , b m ′ R ] = δ R,R ′ δ m,m ′ , commuting otherwise, and w m (x) is the Wannier orbital centered around site m. This projection leads to the so-called Bose-Hubbard model [39]: where n mR = b † mR b mR is the boson number operator at site (m, R). Note that |u 0 | ≪ V 0⊥ implies that J x ≫ J, and therefore we are dealing with a very anisotropic (quasi-1D) version of the Bose-Hubbard model.

Bosonization and low-energy effective theory
In order to treat the interactions in a non-perturbative way, it is convenient to rewrite (3) in terms of collective variables that describe phase and density fluctuations. The method of bosonization [2,17,40,37] allows us to obtain such a description. The description is accurate as long as one is interested in the low-temperature and low-frequency properties of the system. In this subsection we shall briefly outline this method, but we refer the reader to the vast literature on the subject [2,17,40,37] for more technical details.
To study the model of equation (3), we use the density-phase representation of the field operator. To leading order, where is A B is a model dependent (i.e. non-universal) prefactor (see below). The density operator, The fields θ R (x) and φ R (x) describe collective fluctuations (phonons) of the phase and the density, respectively §. Since the density and the phase are canonically conjugate variables, it is possible to show [37,17] that θ R (x) and φ R (x) obey the commutation relation: is small compared to the equilibrium density ρ 0 . Note that this approach does not assume the existence of a condensate and thus properly takes into account the quantum and thermal fluctuations, which are dominant in 1D. Indeed, it is the particle density and not a "condensate" density that appears in the density-phase representation of equations (9,10) (unlike usual mean field approaches). Furthermore, the terms of equation (10) with m = 0 account for the discrete particle nature of the bosons, which is crucial for obtaining a correct description for the transition from the 1D superfluid phase to the 1D Mott insulator [2,26,17,41,42].
Using equations (9) and (10), the Hamiltonian in (3) becomes where δ = (G − 2πρ 0 )/π = 2(a −1 − ρ 0 ) is the mismatch between the density ρ 0 and the periodicity of the Mott potential. The two last terms are related to the Mott potential and the Josephson coupling, respectively. We use the following dimensionless couplings to characterize their strength: where a 0 ≈hv s /µ 1D is a short-distance cut-off. The dimensionless parameter K and the sound velocity v s depend both on the density, ρ 0 , and on the strength of the boson-boson interaction. Indeed, for the Lieb-Liniger model, they depend on a single dimensionless parameter, the gas parameter [1] γ = Mg 1D /h 2 ρ 0 . Analytical results for these parameters are only available in the small and large γ limit (see ref. [37] for results in the intermediate regime) : whereas Galilean invariance fixes the product [2,37]: Therefore, 1 ≤ K < +∞ for the Lieb-Liniger model, K = +∞ corresponding to free bosons and K = 1 to the Tonks limit. For the Lieb-Liniger model, the formula [37]: provides a reasonably good interpolation (accurate to within 10% of the exact results) for the non-universal prefactor of the field operator. It is worth stressing that in (11) we have retained only terms which can become dominant at low temperatures. In the language of the rernormalization group (RG) these correspond to marginal or relevant operators. Whenever, for reasons to be given below, any of the terms retained in equation (11) becomes irrelevant in the RG sense, we shall drop it and consider the remaining terms only. We emphasize that the validity of (11) is restricted to energy scales smaller than min{T, µ 1D }, and therefore its analysis will allow us to determine the nature of the ground state and low-energy excitations of the system. Interestingly, the effective Hamiltonian of equation (11) also describes the low energy properties of the anisotropic Bose-Hubbard model of (6). This result can be derived by a lattice version of the bosonization procedure [17], followed by a passage to the continuum limit. However, the relationship of the effective parameters K, v s , g u , and g J to microscopic parameters J x , J, and U is no longer easy to obtain analytically (except for U small [17] and certain limits in 1D [15]). Thus, in this case K, v s , g u , and g J must be regarded as phenomenological parameters, which must be extracted either from the experiment or from a numerical calculation.

Bose-Einstein Condensation
Let us first discuss the phenomenon of Bose-Einstein condensation in anisotropic (quasi-1D) optical lattices. This is the phase transition from a normal Bose gas to a Bose-Einstein condensate (BEC), which exhibits some degree of spatial coherence throughout the lattice. Coherence arises thanks to the Josephson term, which couples the 1D interacting Bose gases together. In the absence of the Josephson coupling and the Mott potential (i.e. for J = 0 and u(x) = 0), the system behaves as an array of independent TLL's [2,17,40,37]: the excitations are 1D sound waves (phonons) and, at zero temperature and for L → ∞, all correlations decay asymptotically as powerlaws (e.g. phase correlations , which implies the absence of long-range order. However, as soon as the Josephson coupling is turned on, bosons can gain kinetic energy in the transverse directions, and therefore become delocalized in more than one tube. This process helps to build phase coherence throughout the lattice and leads to the formation of a BEC. Nevertheless, the building of coherence just described can be suppressed by two kinds of phenomena: one is quantum and thermal fluctuations, which prevent the bosons that hop between different gas tubes from remaining phase coherent. The other is the existence of a Mott potential that is commensurate with the boson density. For sufficient repulsion between the bosons, this leads to the stabilization of a 1D Mott insulator, where bosons become localized. Whereas we shall deal with the latter phenomenon in section 4, we deal with the fluctuations in this section. In this section we have in mind two different experimental situations: one is a 2D optical lattice where the Mott potential is (experimentally) absent. The other is when the Mott potential is irrelevant in the RG sense. The precise conditions for this to hold will be given in section 4 but, by inspection of (11), we can anticipate that one particular such case is when the Mott potential is incommensurate with the boson density (i.e. δ = 0). Under these circumstances, the cosine term proportional to g u oscillates rapidly in space and therefore its effect on the ground state and the long wave-length excitations is averaged out.

Mean field theory
To reduce the Hamiltonian (3) to a tractable problem, we treat the Josephson coupling in a mean-field (MF) approximation where the boson field operator Ψ R (x) is replaced by ψ c + δΨ R (x), where ψ c = Ψ R (x) and δΨ R (x) = Ψ R (x) − ψ c . We then neglect the term ∝ δΨ R (x)δΨ R ′ (x), which describes the effect of fluctuations about the mean field state characterized by the order parameter ψ c . This approximation decouples the 1D systems at lattice sites R = R ′ so that the problem reduces to a single "effective" 1D system in the presence of a (Weiss) field proportional to |ψ c |. The bosonized form of the Hamiltonian for this system is (we drop the lattice index R): where z C is the coordination number of the lattice (z C = 4 for the square lattices realized thus far in the experiments). The phase of the order parameter θ c = arg ψ c has been eliminated by performing a (global) gauge transformation θ(x) → θ(x) + θ c , and henceforth we shall take ψ c = |ψ c | to be real. We note that, up to the last constant term, H MF eff defines an exactly solvable model, the sine-Gordon model, of which a great deal is known from the theory of integrable systems as well as from various approximation schemes [40,17,43]. The mean-field Hamiltonian must be supplemented by a self-consistency condition, which relates the effective 1D system to the original lattice problem: Note that the average over the cosine must be computed using the eigenstates of H MF eff .
3.1.1. Condensate fraction at T = 0: To obtain the condensate fraction at zero temperature, we first compute the ground-state energy of H MF eff using the following result for the ground state energy density of the sine-Gordon (sG) model [43]: where p = β 2 /(8π − β 2 ) = 1/(8K − 1), with β 2 = π/K. The excitations of the sG model are gapped solitons and anti-solitons (and bound states of them for β 2 < 4π) [40,17]. The soliton energy gap (also called soliton "mass") is [43]: and where Γ(x) is the gamma function. Hence, the ground state energy of H MF eff is given by: The self-consistent condition (18) is equivalent to minimizing the mean-field free energy with respect to the order parameter. At T = 0, the entropy is zero and the condensate fraction is thus obtained by minimizing E MF 0 (|ψ c |) with respect to ψ c = |ψ c |, which yields: In the above expression, we have introduced the dimensionless parameter η = A B (K)(ρ 0 a 0 ) 1/2K . We note that the condensate fraction at T = 0 has a power-law dependence on the amplitude of the Josephson coupling, i.e.
The dependence of the exponent on the parameter K, which varies continuously with the strength of the interactions between the bosons in the tubes, is a consequence of the fact that the appearance of a BEC at T = 0, as soon as J is made different from zero, is affected by the characteristic 1D quantum fluctuations of the TLL's. This behavior is markedly different from what can be obtained by a Gross-Pitaevskii meanfield treatment of the system that overlooks the importance of 1D quantum fluctuations. Further below we shall also encounter other instances where the present results strongly deviate from Gross-Pitaevskii theory.
To show the effectiveness of the fluctuations to cause depletion of the BEC, let us estimate the fraction of bosons in the condensate at T = 0 using equation (24). For K = 1 (i.e. the Tonks limit, where phase fluctuations are the strongest), we take a 0 ρ 0 = 1 and assume that µ 1D = 0.1hω ⊥ [7]. Therefore, for V 0⊥ = 20 E R , the hopping amplitude J ≃ 3 × 10 −3 µ 1D , and hence, using equation (24), ρ −1 0 ψ 2 c ≈ 20 %. We also note that for weakly interacting bosons, i.e. as K → +∞, ψ 2 c → ρ 0 , as expected. In the experiments of Ref. [10], very small coherence fractions were observed. The experimental coherence fraction is defined as the fraction of the total number of atoms that contributes to the peak around k = 0 in the momentum distribution, as measured by time of flight. The fraction was quantified by fitting a gaussian distribution to the peak [10]. We note that such a procedure yields a non-vanishing coherence fraction even for the case of uncoupled 1D tubes, where the momentum distribution exhibits a peak with a power-law tail ∼ |k| 1 2K −1 and a width ∼ max{L,hv s /T } [37,36]. Thus, identifying the condensate fraction with the experimental coherence fraction may not be appropriate as there is no condensate in such a 1D limit. Furthermore, as we show below, in the anisotropic BEC phase discussed above, this experimental procedure would also pick up a non-condensate contribution to the coherence fraction. Nevertheless, provided heating effects and thermal depletion of the condensate are not too strong, one can regard the experimental coherence fraction as an upper bound to the condensate fraction, and since a strong decrease of the coherence fraction was reported in the experiments of Ref. [10], our results showing a strong depletion of the condensate are consistent with the experimental findings.
3.1.2. Condensation temperature: At finite temperatures, thermal (besides quantum) fluctuations cause depletion of the condensate. Ultimately, at a critical temperature T = T c , the condensate is destroyed by the fluctuations. To compute the critical temperature we follow Efetov and Larkin [44] and, taking into account that near T c the order parameter ψ c (T ) is small, consider a perturbative expansion of the self-consistency condition (18) in powers of ψ c . To lowest order, where Ψ(x, τ ) = (ρ 0 A B ) 1/2 e iθ(x,τ ) , τ being the imaginary time, T the ordering symbol in (imaginary) time, and the Weiss field, The correlation function: where the average . . . 0 is performed using the quadratic part of H MF eff . Hence, analytically continuing to real time, equation (26) reduces to: which determines T c ; g R 1 (q, ω) is Fourier transform of the retarded version of (28). This function is computed in the Appendix B (cf. equation (B.15)), and leads to the result: where with B(x, y) = Γ(x)Γ(y)/Γ(x, y) the beta function. Thus, the dependence of T c on J is a power-law whose exponent is determined by K. It is interesting to analyze the behavior of T c in the limit of weak interactions, i.e. for K ≫ 1. Using that F (K ≫ 1) ∼ K, along with Galilean invariance, equation (15), v s = v F /K, we find that where ǫ F =hπρ 0 v F /2 = (hπρ 0 ) 2 /2M. Thus, we see that as K → +∞, T c ∼ J 1/2 . This result comes with two caveats. The first one is that, as we turn off the interactions (i.e. for g 1D → 0) so that K can diverge, the chemical potential µ 1D ≃ g 1D ρ 0 → 0. Therefore since the effective Hamiltonian, equation (11), is valid only for energy scales below µ 1D , the range of validity of the result is going to lower and lower J in this limit. Indeed one cannot simply let K → ∞ and keep J fixed. One would get On the other hand, we know that as g 1D → 0 T c must remain finite because a non-interacting Bose gas in a 2D optical lattice indeed undergoes Bose-Einstein condensation (see below). Therefore, if T c /µ 1D diverges as µ 1D → 0 for K ≫ 1, it is because of this fact. The second caveat is that the behavior T c ∼ J 1/2 implied by (32) is not the correct one in the non-interacting limit. To see this, consider the condition for the Bose-Einstein condensation of a non-interacting Bose gas in a 2D optical lattice (we take L → ∞ but keep the number of lattice sites where ǫ(q) and ǫ ⊥ (Q) = −2J (cos Q y a + cos Q z a − 2) are the longitudinal and transverse dispersion of the bosons, respectively. In the latter case, we have chosen the zero-point energy (i.e. the bottom of the band) to be zero. In the J → 0 limit the 1D systems become decoupled and T c vanishes (i.e. there is no BEC in 1D). Therefore, we expect T c to be small for small J. It is thus justified to expand about Q = (0, 0) the transverse dispersion ǫ ⊥ (Q) ≃ Ja 2 Q 2 . Furthermore, ǫ(q) =h 2 q 2 /2M, for free bosons, and thus (33) leads to the law T c ∼ J 2/3 , which is different from (32). The difference stems from the assumption of a linear longitudinal dispersion of the excitations implied by the quadratic part of (11) and (17). Had we insisted in assuming a linear longitudinal dispersion of the bosons ǫ(q) =hv 0 |q|, where v 0 has dimensions of velocity, we would have found that T c ∼ J 1/2 as implied by (32). This argument implicitly assumes that in the weakly interacting limit, the difference between particles and excitations vanishes. This seems quite reasonable, as the spectrum of excitations of the non-interacting Bose gas coincides with the single-particle spectrum. Nevertheless, we can expect equation (30) to provide a good estimate of the condensation temperature as long as J/µ 1D ≪ 1. This is of course particularly true in the strongly interacting limit where K ∼ 1. The interpolation formula for the prefactor A B , equation (16), allows to be more quantitative than Efetov and Larkin and to obtain an estimate of T c for K = 2, for instance. For the Lieb-Liniger model this value of K corresponds to γ ≃ 3.5 [37] and to (30) yields T c /µ 1D ≃ 0.6 and therefore T c /h ≃ 10 kHz or T c ≃ 70 nK. We note this temperature seems to be well above the experimental estimate for the temperature of the cloud, T < 10 −3h ω ⊥ ≈ 10 −2 µ 1D [7]. Thus, we expect the experimentally realized systems to be in the anisotropic BEC phase described in this section, and the optical lattice to exhibit 3D coherence.

Excitations of the BEC phase:
Let us next compute the excitation spectrum of the BEC. Due to the existence of coherence between the 1D systems of the array, this turns out to be very different from the "normal" Bose gas (TLL) phase, where the gas tubes are effectively decoupled (i.e. ψ c = 0) because fluctuations destroy the long range order. The most direct way of finding the excitation spectrum is to consider configurations where the order parameter varies slowly with x and R and time. Thus the order parameter becomes a quantum field, Ψ c (x.R, t). One then asks what is the effective Lagragian, or Ginzburg-Landau (GL), functional, L GL (x, R, t), that describes the dynamics of Ψ c (x, R, t). There are at least two ways of finding this functional: one is to perform a Hubbard-Stratonovich transformation and to integrate out the original bosons to find an effective action for the fluctuations of the order parameter about a uniform configuration. To gaussian order (random phase approximation, RPA) this is carried out in Appendix A. This method has the advantage of allowing us to make contact with the microscopic theory (the Hamiltonian (11) in our case) and thus to relate the parameters that determine the dispersion of the excitation modes to those of the microscopic theory. Another method, which we shall follow in this section, is to guess L GL using symmetry principles only. To this end, we first write Using global gauge invariance, to the lowest non-trivial order in Ψ c , we can write Since we are interested in the ordered phase where Ψ c = ψ c , we must have a < 0. Thus, we can write V as follows: where λ is a phenomenological parameter. On the other hand, in order to guess the form of the kinetic term T global gauge invariance alone is not sufficient. Further insight can be obtained by close inspection of the effective low-energy theory described by (11). If we write it down in Lagragian form: we see that it is invariant under Lorentz transformations of the coordinates (x, t) where: being β a real number. Thus, T can only contain Lorentz invariant combinations which are the lowest order derivative terms allowed by global gauge invariance. This leads to where v || = v s and Z 2 and v ⊥ are again phenomenological parameters. Higher order terms and derivatives are also allowed, but the above terms are the ones that dominate at long wave-lengths and low frequencies. Therefore, the GL functional reads: We find the excitation spectrum by setting Ψ c (x, R, t) = [ψ c + η c (x, R, t)] e iθc(x,R,t) and keeping only terms up to quadratic order, which yields: Hence 2 gives the dispersion of the Goldstone mode θ c , and ω 2 Hence we conclude that the BEC phase has two excitation branches: a gapless mode, which is expected from Goldstone's theorem, and a gapped mode, the longitudinal mode, which describes oscillations of the magnitude of the order parameter about its ground state value ψ c . It is worth comparing these phenomenological results with those of the analysis of gaussian fluctuations (RPA + SMA) about the mean-field state described in Appendix A. There we find, in the small (q, Q) limit: where v is the soliton gap. One noticeable difference of these results from those of the GL theory is that the values of v ⊥ predicted by the RPA+SMA for the longitudinal and Goldstone modes differ, i.e. v (+) ⊥ . This is due to a poor treatment of the Landau-Ginburg theory that neglected terms of L GL of order higher than quadratic in the fields η c and θ c . These terms represent interactions between the Goldstone and longitudinal field, as well as self-interactions of the longitudinal field. On the other hand, the RPA treatment of Appendix A takes into account some of these interactions and thus yields a correction to the dispersion of the Goldstone and longitudinal modes, which leads to different values of v ⊥ .
It is interesting to consider the behavior of the system in the limit of strong Josephson coupling where J > µ 1D . Under these conditions the motion of the atoms is not mainly confined to 1D and thus they can more easily avoid collisions. Therefore, the effect of fluctuations is expected to be less dramatic and assuming the existence of a 3D condensate is a good starting point. This leads to Gross-Pitaevskii (GP) theory, which, assuming a slowly varying configuration of the order parameter, has the following Lagrangian: is the Wannier function at site R). The condensate fraction ψ c is the solution to the (time-independent) GP equation: where µ is the chemical potential. Formally, and ignoring the different values of the couplings, L GP looks almost identical to L GL except for the fact that However, this change has a dramatic effect on the excitation spectrum as we show in what follows.
and carrying out the same steps as above, we arrive at: From the above expression, we see that the amplitude η c has no dynamics, i.e. there is not term involving ∂ t η c that is not a total derivative . The term −2ψ c η c ∂ t θ c , however, couples the dynamics of the phase θ c and the amplitude η c . This is very different from the situation with the GL theory encountered above. We can integrate out η c , e.g. by using the equations of motion for η c , which yield η c (x, Thus, the lagrangian becomes: being v 2 || =h|ψ c | 2 λ/(2MZ 1 ) and v 2 ⊥ = y ⊥ λ|ψ c | 2 /(2Z 1 ). This means that a BEC described by the GP theory has a single low-energy mode, whose dispersion, ω(q, Q) = v 2 || q 2 + v 2 ⊥ Q 2 , vanishes in the long wave-length limit. This is precisely the Goldstone mode, which corresponds to the one that we have already encountered above when analyzing the excitation spectrum of the GL theory that describes the extremely anisotropic limit where J ≪ µ 1D . The presence of the Goldstone mode as the lowest energy excitation in both theories is a direct consequence of the spontaneous symmetry breaking of global gauge invariance. Thus both kinds of superfluids are adiabatically connected as J is increased. On the other hand, the latter theory also exhibits a mode, the longitudinal mode, that is gapped for (q, Q) = (0, 0). However, the fate of the longitudinal mode as J increases is not clear from the above discussion. To gain some insight into this issue, let us reconsider the derivation of the GL Lagrangian in the light of the previous discussion of GP theory. Apparently, the Lorentz invariance discussed above forbids the existence of a term ihZ 1 Ψ * c ∂ t Ψ c (x, R, t). Nevertheless, we have to remember the emergence of such an invariance in equation (35) is a consequence of the linear dispersion of the TLL modes of the uncoupled tubes. The dispersion is linear for energieshω ≪ µ 1D , but for J ∼ µ 1D non-linear corrections must be taken into account. These break the effective Lorentz invariance of (35) and thus a term of the total derivatives have been dropped from both L GP and L GL .
form ihZ 1 Ψ * c ∂ t Ψ c is allowed. Phenomenologically, we can consider a modified GL theory of the form: where Z 1 is very small for J/µ 1D ≪ 1 but grows with J and should grow large for J/µ 1D > 1. We can perform the same analysis as above to find the excitation spectrum, which leads to the following quartic equation for the dispersion of the modes: where ω 2 ± (q, Q) are the dispersion of the Goldstone and longitudinal modes for Z 1 = 0, which have been found above. We see that for small q and Q there is a solution such that ω 2 = Aq 2 + BQ 2 . Furthermore, there is also a gapped mode, for which ω 2 → ∆ 2 + + (Z 1 /Z 2 ) 2 as q and Q) tend towards zero. Thus we see that the gap of the longitudinal mode depends on the ratio Z 1 /Z 2 . For Z 1 ≪ Z 2 , the dispersion of the longitudinal mode is essentially the same as that obtained from the GL theory. However, as Z 1 /Z 2 grows with J, the gap increases and therefore the exciting the longitudinal mode becomes more costly energetically. At this point is useful to recall that the RPA analysis in Appendix A indicates the existence of a continuum of soliton-antisoliton excitations above a thresholdhω = 2∆ s . The stability of the longitudinal mode is ensured if ∆ + < 2∆ s , which seems to be the case provided one trusts the single-mode approximation employed to estimate ∆ + . However, as the gap is pushed up towards higher energies by an increasing Z 1 /Z 2 , the longitudinal mode is likely to acquire a large finite lifetime by coupling to this continuum. This will eventually lead to its disappearance from the spectrum.

Self-consistent harmonic approximation (SCHA)
An alternative approach to obtain the properties of the BEC phase is the selfconsistent Harmonic approximation (SCHA) (see e.g. Ref. [17]). The idea behind this approximation is to find an optimal (in the sense to be specified below) quadratic Hamiltonian that approximates the actual Hamiltonian of the system. The method is more suitably formulated in the path-integral language, where the role of the Hamiltonian is played by the euclidean action. The latter is obtained by performing, in the Lagrangian of equation (35), a Wick rotation to imaginary time, t → −iτ , and defining the action as The partition function reads: In the SCHA, S[θ] is approximated by a gaussian (i.e. quadratic) action: where θ Q (q, ω n ) is the Fourier transform of θ R (x, τ ). This gaussian form is motivated by considering the action, equation (50), in the limit of large J. We then expect that the replacement cos ( This yields a quadratic action and a phonon propagator the form of (58). To find the optimal G v (q, Q, ω n ), we use Feynman's variational principle [45]: where . . . v stands for the trace over configurations of θ R (x, τ ) using e −Sv [θ] and Therefore, by extremizing we find the optimal G v : where F (Q) = a −2 2 t 1 − e iQ·t , with t the vectors that connect a lattice point to its nearest neighbors, and G −1 0 (q, ω n ) = K πv || ω 2 n + v 2 || q 2 the free phonon propagator (v || = v s ). This self-consistency equation (56) can be solved by rewriting it into two equations as follows: defining we arrive at: The set of equations (57) and (58) are solved in Appendix C), for the variational parameter v 2 ⊥ which is the transverse phonon velocity of the 3D superfluid. At T = 0, the equation for v 0 ⊥ = v ⊥ (T = 0) can be solved analytically. We merely state here the result (the details of the calculation can be found in Appendix C): where B 0 ≃ 1.2971 for a 2D square latttice. Note that the variational propagator G v (ω n , q, Q) has a pole, which for small q and Q, is located at for a square 2D transverse lattice. This is the Goldstone mode with v ⊥ ∼ (J/µ 1D ) 2K/(4K−1) , a power law with an exponent that agrees with the RPA+SMA result of Appendix A, where it is shown that v . However, the SCHA fails to capture the existence of the longitudinal mode. This is because, as described above, the SCHA is tantamount to expanding the cosine in the Josephson coupling term. Such an approximation neglects the fact that the phase θ R (x, τ ) is 2πperiodic, and that as a consequence the system supports vortex excitations. On the other hand, this is captured by RPA because it does not disregard the periodic aspect of the phase since the mean-field state is described by the sine-Gordon model. The latter has soliton, anti-soliton, and breather excitations, which are a direct consequence of the periodicity of θ R .
In the weakly interacting limit (i.e. for K ≫ 1), equation (59) becomes: where we have used that At T > 0, the variational equations must solved numerically (see Appendix C). The absence of a solution to equation (57) such that for T > T c v ⊥ = 0 allows us to also determine the condensation temperature. Thus for weakly interacting Bose gas (K ≫ 1 or γ small), we find: This is the same form as that obtained for large K from the RPA analysis of section 3.1.2, equation (32). Thus, the same caveats as those discussed there apply to it.
This is indeed an artifact of the SCHA, which incorrectly captures the second-order character of the transition from the BEC to the "normal" (TLL) phase.

Momentum distribution function
Besides the excitation spectrum, the SCHA also allows us to obtain the momentum distribution of the system in the BEC phase. The latter can be measured in a time-offlight experiment [25], and it is the Fourier transform of the one-body density matrix: , where r ⊥ = (y, z), n(k, K) can be written as: where W (K) is the Fourier transform of the Wannier orbital W 0 (r ⊥ ). Using the bosonization formula (9) and within the SCHA where G v (q, Q, ω n ) is given by equation (58). The result on the last line is the asymptotic behavior for large x and R. Thus we see that for large distances G v (x, R, t) decays very rapidly and therefore it is justified to expand e Gv(x,R,t=0 Keeping only the lowest order term, we obtain the following expression for (k, K) near zero momentum: where Z 0 = |ψ c | 2 π/2K, and the (modulus squared of the) condensate fraction, where C 0 ≃ 0.2365 for a 2D square lattice. Thus the SCHA approximation also yields a condensate fraction that behaves as a power law of J/µ 1D . The exponent is the same as the one obtained in section 3.1.1 from mean-field theory, equation (25). The momentum distribution can be also obtained at finite temperatures from: Thus for k and K near zero, and T > 0 we find with Z 0 (T ) defined as above with φ c replaced by ψ c (T ).
There are several points that are worth discussion about the results of equations (68) and (72). The first is that these form are anisotropic generalizations of the result that can be obtained from Bogoliubov theory for a BEC [25]. However, for sufficiently large q we expect a crossover to the same power-law behavior of a single 1D Bose-gas ∼ |k| (for T = 0 and L → ∞). This is nothing but a reflection of the dimensional crossover physics: bosons with kinetic energy ǫ(k) =h 2 k 2 /2M ≫ J (but such that ǫ(k) ≪ µ 1D for the bosonization to remain valid) cannot "feel" the BEC and the intertube coherence. However, they still "feel" the characteristic quantum fluctuations of the TLL state. The second observation is related to our discussion of section 3.1.1 on the experimental "coherence fraction": it is clear by fitting a gaussian to n(k, K) about zero momentum, one does not only pick a contribution from |ψ c (T )| 2 but also from the tail of the noncondensate fraction. The presence of a trap (and the inhomogeneity of the system caused it) should not dramatically modify this conclusion. However, as both contributions are positive, we expect that the experimentally measured coherence fraction is indeed an upper bound for |ψ c (T )| 2 .

Deconfinement
In this section we shall obtain the phase diagram of the system for the case where the Mott potential is relevant. One necessary condition for this is that δ = 0 and u 0 = 0 so that g u = 0 in equation (11). As it has been mentioned above, under these conditions, and provided repulsion between the bosons is strong enough so that the parameter K < 2, the system undergoes a transition to a Mott insulating state for J = 0 and arbitrarily small u 0 [2,26,17,41]. In the Mott insulating state bosons are localized about the minima of the potential. The question thus is what happens as soon as a small Josephson coupling between neighboring tubes is present. The latter favors delocalization of the bosons leading to the appearance of a BEC, and therefore competes with the localization tendencies favored by the Mott potential. The competition between these two opposite tendencies leads to a new quantum phase transition known as deconfinement. For small values of |u 0 | and J, it can be studied by means of the renormalization group (RG). This is a much better approach than comparing ground state energies of mean-field theories for each phase, as it takes into account the fluctuations, which are entirely disregarded within any mean-field approximation. However, although the RG approach allows us to estimate the position of the boundaries between the Mott insulating and BEC phases, in the present case it cannot provide any insight into the nature of the deconfinement transition. For this, we shall rely on the mean-field approximation and map the mean field hamiltonian at the special value of K = 1/2 to a spin-chain model.

Renormalization Group calculation
Let us return to the low-energy effective field theory defined by the Hamiltonian in equation (11). Since it is a description valid for low-energies and temperatures, the couplings K, g u , g J , etc. all depend on the cut-off energy scale, which is typically set by the temperature for Λ = T < µ 1D . Thus, as the temperature is decreased, the couplings get renormalized due to virtual transitions to states with energy > Λ = T , which now are excluded from the description. The nature of the ground state is thus determined by the dominant coupling as T = Λ → 0.
The change ('flow') for the couplings can be described by a set of differential equations. For small values of g J and g u the RG flow equations can be obtained perturbatively to second order in the couplings. The details of such a calculation are given in Appendix D. The result reads: where ℓ ≈ ln µ/T . The coupling g F is generated by the RG and describes an interaction between bosons in neighboring tubes. However, since initially there is not such an interaction (i.e. g F (0) = 0) and as dg F /dℓ = O(g 2 J ) its effect on the flow of the other couplings is not very important.
The phase boundary between the BEC and 1D Mott insulator phases can be obtained by studying the asymptotic behavior of the solutions of equations (73) to (76). As the temperature is lowered, i.e. as ℓ → +∞, in the regime 1 4 ≤ K < 2 ¶ the effective couplings g u (ℓ) and g J (ℓ) grow as both the Mott potential and the Josephson coupling are relevant perturbations in the RG sense. If g u (ℓ) grows larger, the localization tendencies of the Mott insulator phase set in, whereas if g J (ℓ) grows larger, delocalization tendencies of the BEC phase described in section 3 dominate in the ground state. Both phases correspond to two different strong coupling fixed points of the RG. When both perturbations are relevant, one has a crossing of the two strong coupling fixed points, and since we obtained the above equations assuming that both g u and g J are small, the RG is not perturbatively controlled when g u and g J grow large. Note that, given this limitation, we cannot entirely exclude that intermediate phases might exist. We will come back to that particular point in the next section.
The relative magnitude of the initial (i.e. bare) values of g J and g u determines which coupling grows faster to a value of order one, where the above RG equations cease to be valid. The faster growth of one coupling inhibits the other's growth via the renormalization of the parameter K: If g u (ℓ) grows faster, K flows towards 0, leading to the 1D Mott insulator. This reflects the fact that, for large enough g u , the Mott potential drives the system towards an incompressible phase, where boson number fluctuations are strongly suppressed by the localization of the particles. On the other hand, if g J (ℓ) grows faster, K also grows and leads to the 3D BEC phase. A rather crude estimate of the boundary can be obtained by ignoring the renormalization of K. In such a case, dividing (74) by (75), keeping only the leading order terms, and integrating this equation up to the scale ℓ * , where both g J (ℓ * ) = g u (ℓ * ) = s, where s ∼ 1 we get that the critical ¶ The regime where K < 1 describes a 1D system of bosons interacting via long-range forces [17]. Thus it is not physical for a system of cold atoms interacting via short-range forces. However, the discussion in this section will be done for a general value of K. value of g c J = g J (0) for a given value of g u (0) scales as: In the next section we will see that this power-law is also recovered in the mean-field treatment of the following section. However, for a more precise estimate of the boundary one needs to integrate numerically the RG-flow equations and take into account the renormalization of K, which can be important especially near K = 2. This is shown in figure 2, where we compare the estimate of (77) with the phase boundary obtained from a numerical integration of the RG equations. The RG equations also allow to extract the gap of the Mott insulator as well as the soliton energy gap that characterizes the decay of the one-body correlation functions (i.e. the healing length) in the BEC phase. This is achieved by integrating the equations until the fastest grown coupling is of order one. One can then use that, at this scale ℓ * , the gap ∆(ℓ) is also of order one, and since gaps renormalize like energies, i.e. ∆(ℓ) = ∆(ℓ = 0)e −ℓ . For example, in the BEC phase, using equation (74), keeping only the leading term, and integrating this equation up to a scale ℓ * ≈ log(Λ/∆ s ) where g J (ℓ * ) ∼ 1 gives ∆ s ∼ (J/µ 1D ) 1/(2−1/2K) , which has the same power-law behavior as the one obtained from the mean-field theory, equation (42). For 1 4 < K < 2, a better estimate (numerical) can be obtained by accounting for the renormalization of K using the full set of RG equations. However, for K > 2, the renormalization of K is negligible and the previous estimate of ∆ s should be accurate. Another instance where it can be necessary to take into account the renormalization of K is the evaluation of the condensation temperature T c in the presence of the Mott potential. In particular, for K < 2, although the initial conditions are such that the coupling g J is the dominant one, K may undergo an important renormalization if g u is non-zero. The renormalized value of K could be then used in combination with the RPA formula, equation (30), for T c .

Mean field theory of the deconfinement transition
The RG approach explained above does not provide any insight into the nature of the deconfinement transition, e.g. it cannot answer questions such as whether the transition is continuous (i.e. second order) or not, or what universality class it belongs to. In order to address, at least partially, some of these issues we will treat again the interchain coupling in a mean field approximation. We begin by rewriting the mean-field Hamiltonian of equation (17) in the presence of the Mott potential: The new mean-field Hamiltonian is now a double sine-Gordon theory, which is not exactly solvable for general values of the parameter K [46]. However, for K = 1 2 the solution takes a particularly simple form as one can directly relate it to the results obtained in section 3.1 for the sine-Gordon model. As it has been pointed out above, the regime K < 1 is not physical for models like the Lieb-Liniger model or the 1D Bose-Hubbard model. However, we expect that the solution at this point captures some of the essential features of the deconfinement transition for K ≥ 1. Indeed in particular in the RG flow nothing special occurs at K = 1 so we can expect the two limits to be smoothly connected. The key observation behind the solution at K = 1 2 is that, as shown in Appendix E, at this point (78) is the effective low-energy description of a Heisenberg anti-ferromagnetic spin chain in the presence of a staggered magnetic field: , where a 0 is the short-distance cut-off. We can perform a rotation about the OY axis, so that the spin variables transform as: where tan ϕ = h z /h x , and the mean-field Hamiltonian becomes: Applying the the spinsS m , the bosonization rules the spins given in Appendix E, the above Hamiltonian becomes: This is another sG model, similar to the one that we encountered in section 3.1. One can again make use of equations (19,20) and (21) to obtain the ground-state energy density, The behavior of the order parameter ψ c follows from the extremum condition: The algebra is essentially very similar to that in section 3.1.1, and we just quote the results here. There exists a critical tunneling J c , when the condensate fraction ψ c → 0, i.e. h x → 0. We find at K = 1/2, Thus the condensate fraction grows continuously from zero at the deconfinement transition from the Mott Insulator at J < J c to the 3D anisotropic superfluid at J > J c according to: where κ(1/3) has been defined in (22) and at K = 1/2 the dimensionless ratio η = ρ 0 a 0 A B (1/2). Note that the scaling J c ∝ u 2/3 0 in (86) agrees with that deduced from the RG approach of equation (77) upon setting K = 1/2.
It is instructive to compare this deconfinement scenario in 3D with the case where only two 1D systems are coupled together (via the same Josephson coupling term), i.e. a bosonic two-leg ladder, which has been studied in [47]. In both the quasi-1D optical lattice and the ladder, there is the competition between the Josephson coupling that favors superfluidity, and the localizing tendency of the commensurate periodic potential on top of which the interacting bosons hop. This is manifested in a very similar structure of the bosonized Hamiltonians. For the ladder case, it is convenient use a symmetric and anti-symmetric combination of the fields of the two chains, which leads to slightly different forms for the RG equations. Nevertheless, the mentioned competition occurs for the same range of K: 1/4 < K < 2, where, in the ladder case, K characterizes to the correlations of antisymmetric fields. Qualitatively, the critical intertube coupling for the deconfinement transition has the same power-law dependence on g u as in equation (77). However, as the ladder is still effectively 1D, the nature of the deconfinement transition is expected to be very different from the one studied here, where a fully 3D (anisotropic) BEC results. In Ref. [47], the deconfinement transition in the ladder was found to be in the Berezinskii-Kosterlitz-Thouless universality class, at least when (for fixed interaction strength U) the intertube hopping is larger than the intratube hopping (note that this limit has not been studied in this work as it would correspond to coupled 2D systems, with very different physics compared to the coupled 1D systems.) Note that in the entire phase diagram (in the thermodynamic limit, see figure 2), we have found evidence for the BEC and the 1D Mott insulating phases only. We have not found any evidence (or hints of evidence) for more exotic states such like the sliding Luttinger liquid phase [48,49,50] or some kind of supersolid phase. The sliding Luttinger liquid [48,49,50] is a perfect insulator in the transverse direction (like our 2D Mott insulator for finite size systems); along the tube direction, the system is still a Tomonaga-Luttinger liquid, but the parameters K and v s are now functions that depend on the transverse momentum. This state may be stabilized by a strong interaction of the long-wave length part of density fluctuations between tubes [49]. We note that, although the RG procedure described above generates such a density-density interaction between the tubes, in the bare Hamiltonian there are no such interactions to start with (atoms in different tubes have negligible interactions) and therefore, the couplingg F is likely to remain small as the RG flow proceeds to lower energies, at least within weak coupling perturbative RG we use here. Thus it seems unlikely that the system can enter the regime of a sliding Luttinger liquid [49,50] phase. An even more exotic possibility would be a phase where the tubes are 1D Mott insulators (in the presence of a commensurate periodic potential and strong enough interaction in-tube), while at the same time, there would be coherence between the tubes in the transverse direction. In this scenario, presumably both the couplings g u (Mott potential) and J (Josephson coupling) would be relevant. We cannot rule out such a possibility in the strong coupling limit. However, we have no evidence for it in our weak coupling analysis. Nevertheless, in a quasi-1D system and on physical grounds, it seems hard to imagine a phase where the boson density is kept commensurate within each tube, as required by the existence of the Mott insulator, while at the same time bosons are able to freely hop to establish phase coherence between tubes (but not within tubes!).
As far as the competition between the Josephson coupling and the Mott potential is concerned, mean field theory is a useful tool (in higher dimensions as well as when the tubes are coupled) for providing insights into which possible phases and what their properties are. However, it is usually inappropriate for predicting the correct critical behavior. For cold atomic systems, because the systems are finite and of the existence of the harmonic trap that makes the samples inhomogeneous, it is not very realistic to study critical behavior as the latter is usually strongly perturbed [34]. Mean field theory seems thus perfectly adapted to the present study. Let us however more generally comment on the universality at the transition. Because the field cos(2φ R ) is a vortex (instanton) creation operator [17] for the field θ R , the Hamiltonian for each tube R can be faithfully represented by a classical 2D XY Hamiltonian involving It thus becomes apparent that the coupling between the tubes leads to an anisotropic version of the 3 + 1 dimensional XY model. The transition will thus be in the universality class of the classical D = d + z = 4 XY model in agreement with general considerations on superfluid to insulator transitions [39]. In the present case, the transition is tuned by varying the anisotropy of the model δ ∼ J − J c . For such a transition, the superfluid density of the classical model is expected to behave as δ ν(d+z−2) , where for the commensurate Mott insulator to superfluid transition z = 1 and, at d = 3, which is the upper critical dimension, the critical exponent ν = 1/2 is mean-field-like. One would thus find that the superfluid density behaves as δ, in agreement with our result for K = 1/2 of equation (87), upon identifying the superfluid density of the classical XY model with the square BEC fraction of the quantum model, ψ 2 c . Note that in the absence of the self-consistency condition (18), the mean field Hamiltonian is equivalent to an XY model in presence of a symmetry breaking field [51,52]. It has been recently argued that such a model can have a sequence of two phase transitions [53,54]. Since in the present case the coefficient of the Mott potential (∼ cos(2φ R )) term is small, the system would be in the regime where it undergoes a single phase transition even for the model with the symmetry breaking field. But, more importantly, the model that we have studied has to be supplemented by the selfconsistency condition, equation (18), which stems from our mean-field treatment of the original Josephson coupling cos(θ R − θ R ′ ). This condition will modify the properties of the phase where this coupling is relevant, and for the reasons explained above, should bring back the model in the universality class of the 3 + 1-dimensional XY model.

Finite size effects
Actual experimental systems consists of finite 2D arrays of finite 1D systems (tubes). The finite size of the sample has some well known consequences in condensed matter physics such as inducing the quantization of the excitation modes and turning phase transitions into more or less sharp crossovers. The latter can be understood using the renormalization-group analysis of section 4.1. The finite size of the tubes introduces a length scale into the problem, namely the size of the tube L, which effectively cutoffs the RG flow at a length scale of the order of L (or, for that matter, the smallest length scale characterizing the size of the sample). Thus, flows cannot always proceed to strong coupling where the system acquires the properties of one of the thermodynamic phases discussed above. Furthermore, as we show below, the finite size L leads to the appearance of a new phase.

Finite size tubes and 2D Mott transition
Typically in the experiments [9,7] a number N 0 of bosons ranging from a few a tens to a few hundreds are confined in a tube of ∼ 10µm long. This leads to energy quantization of the sound modes in longitudinal direction, but also to a finite energy cost for adding (or removing) one atom to the tube. This energy scale is E C =hπv s /KL [2,37]. The tube can be thus regarded as a small atomic "quantum dot" with a characteristic "charging energy" E C . If E C is sufficiently large, it can suppress hopping in an analogous way as the phenomenon of Coulomb blockade suppressing tunneling through electron quantum dots and other mesoscopic systems. In order to know when this will happen E C must be balanced against the hopping energy E J . The latter is not just JN 0 as one would naïvely expect for non-interacting bosons. Interactions and the fluctuations that they induce in the longitudinal direction alter the dependence of E J on N 0 as we show in Appendix F.
Let us consider a system of finite tubes in the BEC phase (i.e. the Mott potential is irrelevant or absent). In the limit of small J (J/µ 1D ≪ 1) of interest here, the longitudinal fluctuations with wavenumber q = 0 can be integrated out as discussed in Appendix F. The resulting effective Hamiltonian is a quantum-phase model identical to the one used to describe a 2D array of Josephson junctions: where N R is the particle-number operator of tube at site R, and θ 0R the canonically conjugate phase operator: The model (88) has been extensively studied in the literature (see e.g. Ref. [55] and references therein). It exhibits a quantum phase transition transition between a two-dimensional superfluid (2D SF) phase (a BEC at T = 0) and a two-dimensional Mott insulator (2D MI) at a critical value of the ratio E J /E C . For commensurate filling (integer N 0 ) using Monte Carlo (MC) van Otterlo et al. [55] found (E J /E C ) c ≃ 0.15. Using the forms for E C and E J derived above, and assuming that the tubes are described by the Lieb-Lininger model, the MC result reduces to J c /µ 1D ≃ 0.3N −3/2 0 in the Tonks limit (K = 1) and to J c /µ 1D ≃ 0.15N −2 0 for weakly interacting bosons (K ≫ 1) [20]. A similar result has been more recently obtained in Ref. [56] using the random phase approximation (RPA). Our result, however, includes critical fluctuations beyond the RPA as it relies on the numerically exact results from MC simulations. In figure 3 we show the complete zero-temperature phase diagram including the possibility of a phase where the tubes are effectively decoupled in a 2D MI phase (labeled Decoupled tubes in the diagram). Note that in the 2D MI, only the phase coherence between different tubes of the 2D lattice is lost. However, apart from the finite size-gap ∼hπv s /L, there is not gap for excitations in the longitudinal direction. This makes this phase different from the 1D MI described above, where there is a gap (depending very weakly on the size L) for longitudinal excitations. If the Mott potential is applied to the system in this phase, it will cross over to a 1D Mott insulator for a finite, but small value, of the dimensionless strength g u . This crossover should take place for K < 2 [2,26,17,41] and is thus represented by a vertical dashed line in figure 3.
In the superfluid phase, the lattice will also exhibit a crossover as a function of temperature: If the condensation temperature T c found in section 3.1.2 is larger than the finite-size gaphv s /L, the gas should behave as a 3D BEC and cross over to a 2D SF for T ∼hv s /L. However, if T c ≪hv s /L the system should behave as a 2D quasi-condensate with no long range phase coherence (and therefore ψ c = 0) at finite temperature. Sincehv s /L ∼hω 0 , ω 0 being the longitudinal trapping frequency, and experimentally [7] ω 0 ∼ 0.1 kHz whereas T c /h ∼ 10 kHz, we expect to be in a regime where the BEC is always observed. The presence of harmonic confinement means that at least some of the tubes will have a filling deviating from commensurate filling. The critical value of J/µ 1D is reduced relative to the half-filled case, as shown in figure 4.
Finally, let us estimate the typical value of the critical ratio J c /µ 1D . Using the above results and assuming that µ 1D ≈ 0.1hω ⊥ [7] and N 0 = 100, we conclude that J c /hω ⊥ ≃ 10 −4 for K = 1 whereas J c /hω c ≃ 10 −5 for K ≫ 1. The deepest lattices used in the experiments of Refs. [7,10] had V ⊥ = 30 E R , which corresponds to J/hω ⊥ = 5 × 10 −5 . This is of the same order of magnitude as J c /hω ⊥ for K ≫ 1, although a larger number N 0 will further reduce this estimate, and therefore for this lattice the system may be well in the BEC phase (see figure 3). However, in the experiments [7,10] no coherence was apparently observed between the tubes. The explanation for this contradiction is that the value of tunneling timeh/J must be compared with the duration of the experiment t exp and it turns out that for a lattice depth of V 0⊥ = 30 E R this time is indeed longer than t exp . Therefore, over the duration of the experiment almost no tunneling process occurs and the build-up of intertube coherence can thus hardly take place.

Comparison to experiments
In this section we discuss and summarize our results in connection with the recent experiments on anisotropic optical lattices [21,7,10,22,16,13]. The complete phase diagram, taking into account the finiteness of the 1D gas tubes, is shown in figure 3. At zero temperature, we predict three phases: an anisotropic BEC, a 1D Mott insulator, and a phase where the gas tubes are decoupled (i.e. phase incoherent), namely a 2D Mott insulator.
However, as mentioned in section 2, the experimental systems are indeed harmonically confined. Unfortunately, so far it has proven very difficult to deal analytically and simultaneously with strong correlations and the harmonic trap. On the other hand, numerical methods can tackle this situation [30,31,32,33,34,35]. In the harmonic trap, the cloud becomes inhomogeneous, and thus as demonstrated in a number of numerical simulations [33,34,35,57,58], for a wide range of system's parameters, several phases can coexist in the same trap. These results have confirmed the experimental observations [12] as well as previous theoretical expectations based on mean-field calculations [24]. The numerical simulations have been mostly restricted to 1D systems [57,33] and isotropic lattices [34], we expect most of their conclusions to also apply to the anisotropic optical lattices that interest us here. Thus, we also expect coexistence of the phases described above when the system is confined in a harmonic trap. Qualitatively, this indeed can be understood using the local density approximation, assuming that the trap potential gives rise to a local chemical potential which varies from point to point. Thus, from strong interactions we expect the formation of domains of the 1D Mott insulating phase near the center of the trap, whereas the superfluid phase should be located in the outer part of the cloud. If the transverse lattice is very deep (i.e. if V 0⊥ is large) there also can be regions where the number of boson per tube is small and the the phase coherence is lost because the charging energy E C becomes locally larger than E J . Nevertheless, confirmation of this qualitative picture will require further experimental and numerical investigation, which we hope this work helps to motivate. The coexistence of different phases may be revealed by analyzing the visibility of the interference pattern observed in an expansion experiment [59,58]. Let us finally mention that the existence of the trap can also lead to new and poorly understood phenomena such as suppression of quantum criticality [34].
Taking into account the above, we only aim at a qualitative description of some experimental observations. Thus, we have shown in 3.1.1 that the strong quantum fluctuations characteristic of 1D systems strongly deplete the condensate fraction in the BEC phase that forms for arbitrarily small intertube Josephson coupling. The condensate fraction can be as low as ≈ 10 % for typical experimental parameters (V 0⊥ = 20 E R ). We have also argued that the experimentally measured coherence fractions may also include some of the non-condensate fraction, and therefore, provided heating effects of the sample in the measurement can be disregarded, the coherence fraction should be an upper bound to the BEC fraction. Quantum and thermal fluctuations are also responsible for the power-law behavior of the condensation temperature with J/µ 1D . Estimates of the latter, when compared to the experimental estimates of the temperature of the lattice [7] and the finite-size excitation gap means that for V 0⊥ = 20E R the system should be in the BEC phase described in section 3 provided the Mott potential is sufficiently weak. In this phase, the system exhibits two kinds of low energy modes in the long wave-length limit: a gapless Goldstone mode, and a gapped longitudinal mode. Whereas the former corresponds to fluctuations of the phase of the order parameter, the latter is associated with fluctuations of its amplitude, namely the BEC density.
If the strength of the Mott potential is increased, the system (or at least part of it near the center of the trap) should undergo a deconfinement quantum phase transition (or crossover, given the finiteness of the sample) to a 1D Mott insulating phase. This phase is characterized by the existence of a gap in the excitation spectrum which is largely independent of the size of the tube. The hopping between the tubes is suppressed for temperatures (or energies) below this gap, and it is incoherent above it. For a deeper optical lattices, we predict that hopping between tubes can also be suppressed by the "charging energy" of the (finite-sized) tubes. We also notice that in some experiments the duration time of the experiment is also a limiting factor for the establishment of phase coherence even if hopping is not effectively suppressed by the charging energy.
Stöferle et al. [10] have studied an optical lattice system for several values of the depth of the transverse optical lattice potential (hence J). By analyzing the width of the momentum peak near k = 0 in the interference pattern after free expansion, they were able to deduce the energy absorption in response to a time-dependent modulation of the Mott potential. From this analysis, they could tell when the system behaves as a superfluid or as a Mott insulator. These experiments were analyzed recently [42,60,61,62] and we will not repeat the conclusions of the analysis here. In the same series of experiments, they also studied the dependence of the width of the momentum peak near k = 0 in the ground state on the strength of the Mott potential, for several values of V 0⊥ . For a BEC or a 1D superfluid, the inverse of the width is set by a geometrical factor and it is proportional to the size of the size of the cloud (or, in the 1D case, the thermal lengthhv s /T , whichever the shortest). For a Mott insulator, however, the inverse of the width is set by the correlation length ξ c ≃hv s /∆ c , where ∆ c is the Mott gap. Thus, an increase of the width as a function of the Mott potential reveal the opening of a gap in the spectrum of the system [33]. The experimental curves of the width as a function of the Mott potential present upturns around the points where the system should enter the Mott insulating phase in the thermodynamic limit. Our expectation, based on the phase diagram of figure 3, is that the deconfinement phase transition should take place for values of the Mott potential slightly larger than in the pure 1D case (represented by the vertical line at K = 2 in figure 3). It is necessary to recall that, in the Bose-Hubbard model (cf. equation (6)) that describes the systems of Ref. [10], the Mott potential controls not only the dimensionless parameter g u that enters in the RG equations of section 4.1, but also the ratio of the interaction U to the kinetic energy J x , and therefore the Luttinger parameter K as well. The latter decreases towards the Tonks limit (K = 1) with increasing U/J x . The theoretical expectation is in good qualitative agreement with the experimental observation.
The finite extent of trapped cloud (either longitudinally or transversally) limits the minimum momenta of the modes. Thus the energy of the lowest modes in the 3D SF phase can be directly estimated from our results by using the minimum available momentum in (40,41). For instance, for a finite 2D lattice containing M y × M z tubes (i.e. an atom cloud of size L×M y a×M z a) the lowest available momentum is ∼ π/(M i a). Putting this value into (40,41) shows that the frequency of the lowest transverse modes decreases with decreasing J. This analysis neglects the possibility of phase coexistence described above, which calls for further investigation of these issues.
We have been not able to calculate the finite temperature excitation spectrum across the dimensional crossover, but it is interesting to speculate what happens to the Tomonaga-Luttinger liquid (TLL) spectrum, especially when the longitudinal momentum q is near 2πρ 0 . On general grounds, at temperatures large compared to the dimensional crossover scale ∼ T c , the system should exhibit essentially the 1D Longitudinal mode ω q Goldstone mode "particle−hole" continuum Figure 5. Schematic excitation spectrum of the BEC phase. Shaded lines and curves represent different types of excitations. The frequency of the excitation is ω and q is momentum in the longitudinal direction (i.e. parallel to the tube axis). The continuum of particle and hole excitations stems from the spectrum a 1D system, as found by Lieb in Ref. [63]. At lower energies, however, the system develops the phase coherence characteristic of a BEC and therefore must obey Landau's criterion of superfluidity.
In particular, this means that the system must exhibit a gap at q = 2πρ 0 , where ρ 0 is the linear density of the 1D system. In the long wave-length limit, we find the system exhibits two kinds of excitations: a gapless (Goldstone) mode, and a gapped longitudinal mode. There is also a continuum of excitations above the longitudinal mode (see discussion at the end of section 3.1.3).
spectrum of a TLL, with its characteristic continuum of (particle-hole) excitations extending roughly between q ≈ 0 and q ≈ 2πρ 0 [63]. As the temperature or the excitation frequency get below the crossover scale, the low-energy excitations around 2πρ 0 must disappear: being the BEC a 3D superfluid, by Landau's criterion, there can be no low energy excitation that contributes to dissipation of superfluid motion. The gap about q = 2πρ 0 in the spectrum of the anisotropic BEC phase would be thus reminiscent of the roton gap in liquid Helium, which may be regarded as the remnant of the crystallization tendency at the momentum of where the roton minimum occurs [64]. In our case, because of the underlying 1D physics, the roton gap should have a continuum of excitations above it, which is the remnant of the 1D continuum of excitation in the TLL (see figure 5). This may account for the "relatively broad" energy absorption spectrum observed in the "1D to 3D crossover" regime in the experiment of Stöferle et al. [10,42]

Conclusions and discussion
In summary, through a variety of methods, we have arrived at the following picture For the case where the Mott potential is either absent, incommensurate, or generally, irrelevant in the RG , the existence of an arbitrarily weak intertube tunneling immediately turns the 1D arrays into a 3D BEC at zero temperature, albeit with anisotropic properties. At finite temperature, thermal fluctuations cause a phase transition at T = T c where the 3D coherence (and the BEC) is lost. Since T c ≪ µ 1D in the quasi-1D optical lattices studied here, for T c < T ≪ µ 1D the system behaves as a 2D array of phase-incoherent Tomonaga-Luttinger liquids. We have obtained explicit expressions for T c , the condensate fraction at T = 0, the excitation spectrum, and momentum distribution function, which could be compared to measurements in the experimental systems. In particular, the power-law behaviors of J/µ found for T c and condensate fraction are consequences of the strong fluctuations that dominate the properties of the uncoupled 1D systems. In the presence of a relevant (in the RG sense) Mott potential, a finite intertube tunneling (J c ) is required to overcome the localization tendencies of the interacting bosons. This process takes place in the form of a deconfinement quantum phase transition. Our main result in this case is the phase diagram shown in figure 3 as a function of J/µ 1D and the parameter K characterizing the interactions in the tube. For the deconfinement transition, we also find a power-law dependence of J c on the strength of the Mott potential u 0 , and a power-law dependence of the condensate fraction near the transition (indeed, these two results are strictly valid for K = 1/2, which is outside the range of values for bosons with short-range interactions, but we speculate that these results are adiabatically connected to the physical (for bosons) regime of K > 1 as we have found no evidence of any quantum phase transitions for 1/2 ≥ K ≥ 1 in our RG calculation.) Finally, for comparisons to experiments, we have looked into trap effects by studying tubes of finite length and described the transition where the system goes from a system of uncoupled 1D Luttinger liquids to the 3D superfluid as (renormalized) J overcomes the "charging energy" of the finite length tube. As discussed in details in the previous section we can compare our predictions with the experiments of Stöferle et al. [10]. Some of our predictions, in particular the large reduction of the coherent fraction is in qualitative agreement with the experiment. Clearly, further experimental and theoretical studies are needed to provide a more thorough and quantitative picture of the dimensional crossover and the deconfinement transition. From the theoretical point of view, taking into account the harmonic trap in a more quantitative way would be a step towards a more quantitative description of the experiments. From the experimental point of view, investigation of larger systems with more atoms per tube and more tubes, and in longer experiments that allow to probe smaller values of the Josephson coupling J, would permit a more complete exploration of the phase diagram of the quasi-1D lattices. Also, most of the current experiments in these lattices have been carried out using rather deep Mott potentials in order to reduce the ratio of kinetic to interaction energy (this is the regime where the Bose-Hubbard model of equation (6) is applicable). However, It would be very interesting to conduct experiments where the interaction energy between the atoms is tuned independently of the longitudinally applied (Mott) potential. Increasing the interaction strength in the absence of Mott potential was recently achieved by Kinoshita et al. [13] by confining the atoms in tighter 1D traps. Under these conditions, and for interaction strength γ ≈ 3.5 [2,41] in 1D, the application of the Mott potential of strength much smaller than the chemical potential µ 1D should suffice to cause localization of the bosons in the 1D Mott insulating phase. To the best of our knowledge, this "weak lattice" Mott insultator has not yet really been realized. In this phase quantum effects are stronger, and it would be then very interesting to drive the system across the deconfinement transition by changing J.

Acknowledgments
We thank T. Esslinger, A. Nersesyan, and M. Gunn for useful conversations and/or correspondence. MAC and AFH thank the hospitality of the University of Geneva where this work was started. AFH is funded by EPSRC (UK), MAC by Gipuzkoako Foru Aldundia and MEC (Spain) under Grant FIS2004-06490-C03-01. This work was supported in part by the Swiss National Science Foundation through MaNEP and division II.

Appendix A. Mean-field theory (MFT) and gaussian fluctuations (RPA)
In this section we shall employ the functional integral to derive the mean-field Hamiltonian theory of section 3.1, and to go further by taking into account gaussian fluctuations about the mean-field state. The method is known as the random-phase approximation (RPA). Our treatment follows Ref. [65].
Let us consider the model defined by equation (11). Using the coherent-state pathintegral, the partition function can be written as follows: where the (euclidean) action, H(ψ * R , ψ R ) being the Hamiltonian density corresponding to equation (11). Since the Mott potential plays no role in the BEC phase (i.e. it is irrelevant in the RG sense), we drop it from the Hamiltonian H(ψ * R , ψ R ). However, the action S [ψ * , ψ] still contains the non-trivial Josephson coupling between neighboring 1D systems.
In the above expression we have introduced the following notation: where t runs over set of vectors that connect a lattice point to its nearest neighbors. We perform a Hubbard-Stratonovich (HS) decoupling of S J [ψ * , ψ]: To make sense of J −1 R−R ′ the Fourier transform of J R,R ′ is needed, (A.7) The last result applies to a square lattice (in the OYZ plane) with lattice parameter a. Hence, Thus, after performing the HS decoupling, the action becomes: To perform the mean field approximation we formally integrate out the boson degrees of freedom described by ψ R (x, τ ) and ψ * R (x, τ ), and obtain the effective action for the HS fields: We have introduced the notation S GL for the effective action to indicate that it is indeed the GL functional.
To obtain the mean-field equations we perform a saddle-point expansion. The saddle-point is determined by the following equations: We look for uniform solutions to the above equations, i.e. ξ R (x, τ ) = ξ c and ξ * R (x, τ ) = ξ * c . Furthermore, because of global gauge invariance, the particular choice for the phase of ξ 0 should not matter, that is, if (ξ * c , ξ c ) is a saddle point then so is (e −iθc ξ * c , e iθc ξ c ). Thus S ψ [ξ * c , ξ c ] = S ψ (|ξ c | 2 ). This implies that the above equations reduce to: where the prime stands for the derivative with respect to ρ ξ = |ξ| 2 . Furthermore, using equation (A.13), equation (A.14) can also be written as: where it is understood that the expectation value of ψ R (x, τ ) is taken by tracing out the boson degrees of freedom in the presence of the constant HS (ξ * c , ξ c ), i.e. . using S[ξ * c , ξ c , ψ * , ψ] (cf. equation (A.9)). Therefore, (A.15) must be solved self-consistently for (ξ * c , ξ c ): This is the self-consistency condition mentioned in section 3.1. Upon multiplying equation (A.17) by J R ′′ −R from the left and summing over R, we arrive at: where z C is the number of nearest neighbors (z C = 4 for a 2D square lattice). The above system of equations defines the mean field theory, which we have solved, in Hamiltonian form, in section 3.1. We next take into account fluctuations by considering the gaussian corrections to the saddle point (ξ * c , ξ c ). Let To proceed further we need to integrate out the boson degrees of freedom and obtain S ψ [ξ * , ξ] explicitly. Since we are interested in the properties at long wave-lengths and low frequencies and temperatures, we first apply bosonization and obtain: In the above expression, In this (bosonized) form, integrating out the boson degrees of freedom amounts to tracing out the phase field θ R (x, τ ). This will carried out in what follows using the cumulant expansion: where . . . sG denotes the trace over (smooth) configurations of θ R (x, τ ) using the weight e −S[ξ * c ,ξc,θ] (cf. equation ( A.21), note that this action defines a collection of independent sine-Gordon models, which explains the notation). The RPA corresponds to keeping the terms of the above expansion up to quadratic order in δξ and δξ * . Thus, using that we obtain: We have introduced the following notation: Since in this action the different lattice sites are uncoupled, the correlation functions are diagonal in R and R ′ . Thus we can write the quadratic action that describes the (gaussian) fluctuations around the saddle point as follows: Note that the terms in S RPA GL linear in Ξ and Ξ † vanish because we expand around the saddle point. We have also dropped the constant term that is proportional to the mean field free energy. In the above expression the following spinors have been introduced: where we have introduced the sources q * R (x, τ ) and q R , x, τ ), which in the original action are coupled to ψ R (x, τ ) − ψ c and ψ * R (x, τ ) − ψ * c , respectively. The matrices: , (A.34) with 1 = (R, x, τ ) and 2 = (R ′ , x ′ , τ ′ ), and (p, q = ±): To perform the averages using S sG a gauge transformation is applied θ R (x, τ ) − θ c → θ R (x, τ ) so that ξ c = −z C Jψ c can be treated as a real number.
Let us now pause to consider the consequences of global gauge invariance on the correlation functions of (A.35). As we have found above, the GL functional is a function of |ξ c | 2 only. Thus, if we expand S ψ (|ξ| 2 ) around the saddle point solution: Thus, by setting δξ * , δξ = const. and taking β → ∞ in equation (A. 25), and comparing with equation (A.36), the following results can be obtained: where we have used that lim β→+∞ S ψ (|ξ c | 2 )/(M 0 Lhβ) = E sG (|ξ c | 2 ) is the energy density of the sine-Gordon model, equation (19). We shall use the above identities to calculate of the energy dispersion of the modes in the BEC phase. This phase exists below the critical temperature (T c to be determined below), and is characterized by having ψ c ∝ ξ c = 0 so that global gauge symmetry is spontaneously broken. To compute the excitation spectrum in this phase, we diagonalize the matrix J −1 + G (henceforth we drop the index R − R ′ ). J −1 is already diagonal, whereas the matrix G is rendered diagonal by: (A.41) Thus, whereg 1 (q, ω n ) is the Fourier transform anomalous correlator g ++ 1 (x, τ ) = g −− 1 (x, τ ). In what follows we denote the eigenvalues of G(q, ω n ) as c ± (q, ω n ) = g 1 (q, ω n ) ±g 1 (q, ω n ). Next we express S RPA GL in terms of the eigenvectors of G(q, ω) and to subsequently integrate out the transformed HS fields. In the following we use a matrix notation where the summations over R, R ′ as well as over q, ω n and multiplication byh −1 are implicitly understood. Hence, Let us make the replacement G = UCU † , so that We shall obtain the generating functional in terms of these new sources. We can now integrate over the auxiliary fields X † , X. Since they are related to Ξ † , Ξ by a global unitary transformation the measure of the functional integral is not affected. Thus we are free to shift where N ′ is an unimportant constant. Although we have been careful to treat the matrices in the last expression as non-commuting, all of them are diagonal and the last expression can be simplified to: Thus, the system has two modes, whose propagation is described by ∆ = C(1 + J C) −1 . When expressed in Fourier components this propagator reads: The excitations are poles of the propagator: to know more about the correlation functions c ± (k, ω) for general k, ω. Let us for a bit return to real space; we first notice that To make contact with the notation of Ref. [66], we introduce Φ = θ/β, β 2 = π/K, x = (v s τ, x), and µ = √ A B ρ 0 ξ c /hv s = √ A B ρ 0 Jz C ψ c /hv s . The excitations of such a sine-Gordon model are the solitons and anti-solitons. Furthermore, β 2 < 4π (i.e. K > 1 4 ) the model also has breather excitations [40,17], which are soliton-anti-soliton bound states, and whose zero-momentum energies ("mass" gaps) are [40,17]:  (20)). The operators sin βΦ(x) and cos βΦ(x) have zero conformal spin and therefore can only create breathers (if they exist) or pairs of solitons and anti-solitons. Furthermore, since under charge conjugation [40], (A. 53) it follows that Thus the cosine operator can only have poles corresponding to breather states that are even under conjugation, whereas the sine can only have poles for breathers that are odd under conjugation. Taking all these facts into account we can write the spectral function of the correlators of these two operators [40]: where ρ sin,cos ss (q, ω) describes the continuum of soliton and anti-soliton excitations. We have also used that the solitons with even quantum number n are even and those with odd n are odd under C [40]. To make further progress, we assume that the behavior of the above correlation functions is dominated by the lowest energy pole (single-mode approximation, SMA), i.e. the lowest breather ∆ 1 = 2∆ s sin(pπ/2). Hence, the retarded correlation function: j=y,z (1 − cos Q j a) /h 2 ≃ (v s q) 2 + (∆ 1 a/2h) 2 Q 2 , the last expression being valid for |Q| ≪ a −1 . However, we emphasize that note that as K → +∞ (i.e. the bosons become weakly interacting) the number of breather poles proliferate and the SMA may break down. Going beyond the SMA requires some knowledge of the spectral weights w n . Instead, we show the same form of the dispersion for the lowest-energy Goldstone mode can be also obtained from the self-consistent Harmonic approximation described in section 3.2.
Using similar methods we can obtain the energy dispersion of the longitudinal mode, which is the solution of the equation: We again resort to the SMA, but taking into account that this time the pole is given by the lowest energy breather that is even under C, i.e. ∆ 2 = 2∆ s sin(πp). We first fix the pole reside z + by studying the limiting behavior of c(q, ω) as q, ω → 0. This can be obtained from (A.36,A.38), along with (A.39,A.40), which imply: where we have used gauge invariance to set θ c = 0 so that |ξ c | = ξ c , together with which can be directly obtained from the expressions for E sG and ψ c in section 3.1. Hence, assuming that (SMA): ]. If, within the SMA, we solve equation (A.60), we arrive at: where ∆ 2 + = ∆ 2 2 (8K − 2)/(8K − 1) is the energy gap of the longitudinal mode at (q, Q) = (0, 0), and F (Q) = a −2 t 1 − e iQ·t /2. At the transition temperature, T = T c , the order parameter ψ c = ξ c = 0 and the matrix J −1 R−R ′ + G R−R ′ becomes singular. For ξ c = 0, the anomalous correlators g ++ 1 = g −− 1 = 0, and g +− 1 = g −+ 1 = g 1 , which leads to: Upon multiplying by J R ′′ −R , summing over R, and integrating over x and τ , we obtain: This is precisely the same condition for the BEC temperature derived in section 3.1.2. The expression for g R 1 (q, ω), which is needed to obtain T c explicitly is evaluated in the appendix that follows this one.
Finally, to make connection with the discussion of section 3.1.3, we explain that the identification of the mode labeled as '−' with the Goldstone mode and the one labeled as '+' as the longitudinal mode is indeed quite natural. If we perform a gauge transformation such that θ c = 0, the matrix U of equation (A.41) produces the transformation: Let us consider (small) fluctuations of the order parameter of the form considered in section 3.1.3: . Thus we see that δξ

Appendix B. Finite-temperature phase susceptibility
In this appendix we compute the Fourier transform of the retarded phase-susceptibility of a 1D system of interacting bosons. This correlation function is to leading order the one-body boson Green's function: In the above expression, we have replaced the boson field by its bosonized form to leading order: where A B = A B (K) is a non-universal prefactor. To compute this correlation function, the (gaussian) action (A.22) must be used. Either by direct computation using functional integrals [17], or by means of a conformal transformation of the correlation function at T = 0 [37], one obtains the following result: In order to obtain the retarded version of this correlation function, we shall use the relationship (see e.g. [17]): where g F 1 (x, t, T ) is the time-ordered (à la Feynman) propagator. The latter can be obtained from (B.1) by analytical continuation: τ = it + ǫ(t), where ǫ(t) = sgn(t)ǫ (ǫ → 0 + ). Hence, where t ǫ = t − iǫ(t). We next exponentiate the above result and use that which follows upon expanding in powers of ǫ, as well as the identity ln(a + iǫ) = ln |a| + iπθ(−a) (the branch cut of the logarithm is put on the negative real axis). Thus, the imaginary part of g F 1 (x, t, T ) does not vanish provided that Re R(x, t) < 0, that is, if where D(K) = ρ 0 A B (K) sin π 4K /h. Since we are interested in the Fourier transform of g R 1 (x, t, T ), it is convenient to express it in terms of the imaginary part of g F 1 (x, t, T ): where we have used that dt dx = dξ + dξ − /2v s and have introduced k ± = ω/v s ± k. In the third step we employed that Im g F 1 (x, t, T ) is only non-zero for −v s t < x < v s t. Finally, we have expressed the integral in terms of "light-cone" coordinates ξ ± to be able to benefit from the separability of Im g F 1 in terms of these coordinates. Hence, where (see e.g. [17]): , (B.14) B(x, y) being the beta function, B(x, y) = Γ(x)Γ(y)/Γ(x + y). In the calculation of the critical temperature we need the value of g R 1 (k, ω, T ) for q = 0 and ω = 0: Appendix C. Self-consistent Harmonic Approximation calculation In this appendix, we fill in some of the details of the SCHA calculation. First, note that in equation (53), we can rewrite Then, the variational approximation to the free energy becomes and In evaluating S int v , we can take advantage of point group symmetry of the lattice in the perpendicular directions and replace G v (0, R−R ′ , 0)−G v (0, 0, 0) inside the exponential by the more symmetric expression 1 , where t are the unit lattice vectors in the perpendicular directions. Then, the outer RR ′ just gives a factor of z c . Taking the variational derivative of F , we get: with the free phonon propagator: This single self-consistent equation for G v (q, Q, ω n ) can be solved by rewriting it into two equations: Substituting equation (C.7) into (C.6), we get the equation for the variational parameter, the perpendicular phonon velocity v ⊥ : . (C.8) The Matsubara sum can be done to give: We first study the variational equation (C.9) at zero temperature. At T = 0, coth(x/2T ) → sign(x). To perform the summations over momentum, we take the continuum approximation for both the parallel (q) and perpendicular (Q) momenta, with a cut-off for the parallel momentum to be Λ = µ 1D /v || (µ 1D is the 1D chemical potential of an isolated tube). Thus, with 1 M 0 Q −→ a 2 π/a −π/a dQydQz (2π) 2 , and 1 L q −→ Λ −Λ dq 2π , we get: where v 0 ⊥ = v ⊥ (T = 0) and B = a 2 2 π −π dQ x dQ y (2π) 2 F (Q) ln a 2 F (Q) ≈ 0.836477. (C.12) Putting equations (C.11) and (C.9) together then gives equation (59) of Section 3.2.
We have assumed here v 0 ⊥ /a ≪ µ. Otherwise, our starting point using the Haldane harmonic fluid construction is not valid in principle.
When K → ∞ (towards the non-interacting Bose gas limit), A B (K) ≈ 1, and using the Lieb and Liniger solution [1] we can show that µ 1D ≈ v || ρ 0 π/K. Thus, asymptotically for large K, we get: Thus, v 0 ⊥ /a < µ 1D if and only if z c J < µ 1D . But since for the nearly free Bose gas, µ ∝ 1/K 2 is small, this demands a very hopping amplitude J.
In equation (59), there is a singularity as K → 1/4, which is a signature that for K < 1/4 the Josephson coupling becomes irrelevant in the RG sense. As this regime is beyond the reach for bosons with a Dirac-delta interaction, we will not consider it any further.
We now look at the variational equation (C.9) at finite temperature. First, defining v ⊥ (T c ) = v c ⊥ , we rewrite (C.9) as: (C.14) Following Donohue [65], we approximate crudely coth(x) ≈ 1/x for |x| ≤ 1, and coth(x) ≈ 1 otherwise. Next, we develop a series expansion in v c ⊥ /aT c < 1 for the integrand in I(T c , v c ⊥ ) − I(0, v c ⊥ ). It will turn out that at the transition We have added here a term proportional tog F , which is not present initially, but which is generated at second order ing J after integrating out short-distance degrees of freedom (see below). It represents a density interaction between nearest neighbor 1D systems. We next introduce the chiral vertex operators V + β (R,z) = : e iβφ R+ (z) : and V − β (R, z) = : e iβφ R− (z) :, and define z = v s τ + ix andz = v s τ − ix, and d 2 z = v s dτ dx = dzdz/2i. To explicitly display the scaling dimensions of the various operators in H eff , we introduce the dimensionless couplings: Thus, the interactions are described by The RG is performed in real space on the perturbative expansion of the partition function Z(a 0 ) = Z 0 T exp [−S int ] . We first consider an infinitesimal change of the short distance cut-off a 0 → a ′ 0 = (1 + δℓ)a 0 , where 0 < δℓ ≪ 1. Subsequently, we integrate out short-distance degrees of freedom in the range a 0 < |z| < a 0 (1 + δℓ)a 0 . Following Ref. [67], we use operator-product expansions (OPE) to evaluate expectation values of products of operators at two nearby points in space and time. For z → w (respectively,z →w): : + · · · (D.6) ∂φ R+ (z)∂φ R+ (w) = − 1 (z −w) 2 + : ∂ φ R+ z +w 2 2 : + · · · (D.7) V + β (R,z)V + −β (R,w) = 1 (z −w) β 2 1 + iβ(z −z)∂φ R+ z +w 2 + · · · , (D.8) z + w 2 ) + · · · (D.10) Therefore, we begin with by considering the perturbative expansion of the partition function at the new scale a ′ 0 = (1 + δℓ)a 0 , Z((1 + δℓ)a 0 , {g i (ℓ + δℓ)}) = 1 + Z (1) ((1 + δℓ)a 0 , {g i (ℓ + δℓ)}) + Z (2) ((1 + δℓ) At second order, there are ten terms in the expansion of Z (2) [(1 + δℓ)a 0 , {g i (ℓ + δℓ)}]. As the manipulations are standard but rather long, we will just look at two of them in detail to illustrate the procedure. To see what kind of terms are generated from the product of two operators, we shall use "center-of-mass" coordinates: u = z − w (respectively, u =z −w) and v = (z + w)/2 (respectively,v = (z +w)/2), and we split the integrals of the operator products appearing at second order as follows: The terms generated at second order by the RG transformation stem form the second integral, where the relative coordinate is restricted to the infinitesimal range a 0 < |u| < a 0 (1 + δℓ). To O(g F g J ), there are two possible operator pairings because the hopping term couples operators from neighboring chains: Focusing on the terms where R = T and R ′ = T ′ , or R = T ′ and R ′ = T and using the OPE's (D.10,D.11) with β = ±1/2 √ K, we arrive at the following expression: where the first factor of two stems from the two terms of O(g F g J ) and the second from the two possible pairings of R, R ′ and T, T ′ . Thus we see that to second order in the couplings the coupling g J renormalizes as follows: )g J (ℓ) + g F (ℓ)g J (ℓ) 4K(ℓ) δℓ.
(D. 21) In the limit δℓ → 0, this leads to the following differential equation: (D. 22) Other terms are dealt with in a similar fashion. It is also worth noticing that terms of O(g 2 K , g F g K , g K g J , g K g J ) need not be taken into account because at each RG step we set g K = 0 by properly renormalizing K. Thus g K (ℓ + δℓ) = 0, but new terms of the form of the operator proportional to g K are generated at O(g 2 u , g 2 J ). In order to get rid of them, we perform an infinitesimal canonical transformation. To understand this, let us consider a simplified version of the above model, where the chain index R is dropped, and the Hamiltonian reads: We assume the dimensionless coupling δg ≪ 1. The precise functional forms of H ± are of no importance in what follows. Next consider the following infinitesimal transformation: To O(δg) one can show that this is a canonical transformation which, to O(δg 2 ), brings the first term of the Hamiltonian to the diagonal form: 26) conditions [36,37]. In such a finite systems, the phase field operator consists two terms: θ R (x, τ ) = θ 0R (τ ) + Θ R (x, τ ), (F.1) where Θ R (x, τ ) describes phase fluctuations (phonons) of wave number q = 0. The operator θ 0R is canonically conjugate to the number operator N R , i.e. : [θ 0R , N R ′ ] = iδ R,R ′ . Introducing (F.1) into (49), the quadratic part of S reduces to The hopping term (50), however, couples θ 0R (τ ) and Θ R (x, τ ) in a non-linear fashion: In order to obtain an effective action in terms of θ 0R (τ ) only we need to integrate out Θ R (x, τ ). We shall do this perturbatively assuming that J/µ 1D is small. Thus, let us define: where we have used the cumulant expansion and employed the following notation: . The lowest order contribution to S eff J [θ 0 ] comes from the term S J Θ , which yields: Using that Θ 2 R (0) Θ = 1 2 K −1 ln(L/a 0 ), where a 0 is the short-distance cut-off, the above expression leads to: Assuming that each 1D tube is described by the Lieb-Liniger model, A B = (a 0 ρ 0 ) 1/2K ≃ (K/π) 1/2K [37], where a 0 ≃hv s /µ (i.e. approximately the healing length, ξ). Hence, the effective hopping energy of the quantum-phase model is E J = E J (N 0 ) ≃ J(N 0 ) 1−1/2K , where we have used that ρ 0 L = N 0 . The complete effective action reads: where the charging energy E C =hπv s /KL. This euclidean action corresponds to the following Hamiltonian (we restore the chemical potential term): This is the quantum-phase model used in section 5.1.