Luttinger hydrodynamics of confined one-dimensional Bose gases with dipolar interactions

Ultracold bosonic and fermionic quantum gases confined to quasi-one-dimensional (1D) geometry are promising candidates for probing fundamental concepts of Luttinger liquid (LL) physics. They can also be exploited for devising applications in quantum information processing and precision measurements. Here, we focus on 1D dipolar Bose gases, where evidence of super-strong coupling behavior has been demonstrated by analyzing the low-energy static and dynamical structures of the fluid at zero temperature by a combined reptation quantum Monte Carlo (RQMC) and bosonization approach. Fingerprints of LL behavior emerge in the whole crossover from the already strongly interacting Tonks–Girardeau at low density to a dipolar density wave regime at high density. We have also shown that a LL framework can be effectively set up and utilized to describe this strongly correlated crossover physics in the case of confined 1D geometries after using the results for the homogeneous system in LL hydrodynamic equations within a local density approximation. This leads to the prediction of observable quantities such as the frequencies of the collective modes of the trapped dipolar gas under the more realistic conditions that could be found in ongoing experiments. The present paper provides a description of the theoretical framework in which the above results have been worked out, making available all the detailed derivations of the hydrodynamic Luttinger equations for the inhomogeneous trapped gas and of the correlation functions for the homogeneous system.

2 be effectively set up and utilized to describe this strongly correlated crossover physics in the case of confined 1D geometries after using the results for the homogeneous system in LL hydrodynamic equations within a local density approximation. This leads to the prediction of observable quantities such as the frequencies of the collective modes of the trapped dipolar gas under the more realistic conditions that could be found in ongoing experiments. The present paper provides a description of the theoretical framework in which the above results have been worked out, making available all the detailed derivations of the hydrodynamic Luttinger equations for the inhomogeneous trapped gas and of the correlation functions for the homogeneous system.

Introduction
Ultracold bosonic and fermionic quantum gases are versatile and robust systems for probing fundamental condensed-matter, atomic and molecular physics problems and for devising useful applications for quantum atom optics, quantum information processing and precision measurements. Since the discovery of Bose-Einstein condensation (BEC) of Bose alkali atoms of 87 Rb and 23 Na in 1995 [1], followed by the realization of quantum degeneracy in Fermi gases [2] and the superfluid transition of 6 Li and 40 K in the presence of Fano-Feshbach resonances [3], the field of ultracold atoms has experienced tremendous developments, where experiment and theory have been stimulating each other in predicting and realizing novel quantum states. The reason why ultracold atoms are an ideal laboratory for fundamental and applied physics lies in the possibility of reaching extreme quantum regimes by tuning several 'knobs' in a clean and highly controllable way. Basic knobs are temperature, which can be reduced to a few nanokelvin, and the possibility of disentangling interaction and statistical effects since bosons, fermions and mixtures of them can be equally utilized. More interestingly, the short-range atomic interactions can be effectively tuned in strength and sign after variations of an external magnetic field by means of the Fano-Feshbach resonance mechanism [4].
More recently, it has been demonstrated that the range of the interactions can also be manipulated, and dipole interactions with a long-range-like anisotropic character have been observed in a BEC of 52 Cr atoms [5] after exploiting the large magnetic moments of this atomic species. A BEC containing up to 50 000 52 Cr atoms has then been realized below a transition temperature T c 700 nK [6,7] and its dynamical behavior is being investigated [8]. For dipolar gases subjected to anisotropic harmonic confinement, the stability diagram has been clarified to be governed by the trapping geometry [9,10] 7 , as corroborated also by path-integral quantum Monte Carlo (QMC) studies [11]. More recent diffusion Monte Carlo (DMC) simulations reach different conclusions as a result of a more accurate inclusion of the dipole-dependence of s-wave scattering length describing the short-range part of the interaction [12].
Tuning of the interactions can be combined with the possibility of storing ultracold atoms in one-, two-and three-dimensional (1D, 2D and 3D) optical lattices, thereby reducing their dimensionality and enhancing the role of quantum fluctuations [13]- [19] 8 . Under these conditions, ultracold atoms can in fact be utilized, e.g. for precision measurements [20] quantum computing [21] and the development of future atomtronic quantum devices. Dipolar atomic gases in reduced dimensionality have been investigated by means of simulational and theoretical methods. Theoretical studies at the mean-field level have predicted the rich quantum-phase diagram of dipolar bosons in optical lattices [22]. QMC studies of the ground state [23] together with perturbative calculations of the excitations [24] have been carried out in 2D and quasi-2D, respectively. In quasi-1D geometry, mean-field calculations have predicted the occurrence of a roton spectrum [25]. Recently, the physics of super-Tonks-Girardeau (TG) regime in trapped 1D dipolar gases has been discussed [26,27] and the influence of dipolar interactions on the confinement-induced scattering resonances has also been reported [26].
Special attention has been devoted to the strictly 1D case with repulsive interactions, in which the QMC determination of the equation of state [28] has suggested the existence of a crossover from a liquid-like, superfluid state with the characteristics of a TG gas [29] to 4 a quasi-ordered (particles localized at lattice sites), normal-fluid state with increasing linear density. Theoretical analysis of the simulational outcomes may become a useful tool for a theoretical comprehension of the underlying physics and become in turn a robust guide to the experiments. It is well known that the ground-state and the low-energy excitations of 1D liquids with sufficiently short-range interactions can be effectively described within the Luttinger-liquid (LL) Hamiltonian [30], yielding analytical expressions. For interacting fermionic particles in 1D, the LL Hamiltonian predicts the absence of the Fermi discontinuity at the Fermi surface and the separation between the charge (or density, for neutral fluids) and the spin degrees of freedom [30,31].
Along these lines, we have more recently provided the first evidence of two relevant aspects of the dipolar gas in 1D [32]. Firstly, from the methodological viewpoint, we have shown that the low-energy behavior of 1D dipolar atomic Bose gases can be effectively described within the LL framework, finding remarkably good agreement between the LL prediction for the static structure factor and the outcome of the reptation quantum Monte Carlo (RQMC) [33] simulations. Secondly, this combined LL and RQMC analysis has established that the 1D dipolar gas is always found in super-strongly coupled conditions (Luttinger exponent K < 1) at all densities in the crossover from the TG to an ordered dipolar-density-wave (DDW) regime. The dipolar tail of the interactions greatly enhances the correlation effects. As a result of such an analysis, validated analytical expressions have been provided for the density-dependence of all the relevant physical quantities of the homogeneous fluid, namely the ground-state energy per particle and the LL exponents. From the latter, further analytical expressions for the lowenergy observable properties of the fluid have been predicted. Among these quantities are the momentum distribution and the static structure factor, which can be measured in ongoing experiments.
As we have obtained from more recent simulational work on the dynamical density-density correlation function [34], the robustness of the LL behavior is also observed in the dynamical correlations. In accordance with two more refined LL predictions, no gap ever opens up at the first Brillouin zone and no long-range order can exist in the 1D gas of dipolar bosons.
All the above work is concerned with the homogeneous dipolar Bose gas. Under special conditions, where inhomogeneous effects play a minor role, it also provides quantitative predictions on observables which can be tested in experiments with e.g. polar molecules [32] where the whole crossover region could, in principle, be accessed. Atomic gases, however, live in confined geometries, usually given by harmonic traps. The question that arises immediately is then: Which are the signatures of such LL behavior in the crossover, and how can they be monitored within typical experimental set-ups? Moreover, the study of e.g. the collective modes under realistic trap conditions can be useful, in view of current proposals to implement a quantum memory, to an investigation of the decoherence mechanisms that may arise from the coupling of internal and external degrees of freedom [35]. To this aim, we have applied the local density approximation (LDA) to refined RQMC data for the equation of state of the homogeneous gas, to determine the frequencies of the collective trap modes in the crossover [36]. The use of two different methods, Luttinger hydrodynamics and a sum-rule (SR) approach, have led to the same results.
The present paper reviews the description of the theoretical framework in which the above results have been worked out, making available, for easiest use, the detailed derivations of the correlation functions for the inhomogeneous gas and of the hydrodynamic Luttinger equations.

5
The results of our work on the 1D dipolar Bose gases are discussed and placed in the realistic perspective of the experiments, which will be performed in traps.
The plan of the paper is as follows. In section 2.2, we give a brief overview of the LL theory for the homogeneous dipolar Bose gas from [32], providing additional details on the expressions relating the Luttinger exponent K and velocity of the excitations u to the compressibility, the static structure factor and the ground-state energy per particle. We complete our overview of the homogeneous gas by recalling how the Luttinger exponent K can be determined as a function of the density [32], by means of the RQMC simulation technique [37]. This is done in section 3, where a detailed comparison of the structure factor as obtained from RQMC and from the LL theory in the crossover is presented, together with a description of the ground-state properties. With the detailed knowledge of the homogeneous system at hand, we can turn our attention to the inhomogeneous LL. In section 4, the general hydrodynamic equations for the inhomogeneous LL are set up, providing the essential background on the use of the LDA as well as some useful analytic considerations. This theoretical framework is then applied to the special case of an external harmonic trapping in section 5, where the analytic expressions related to both the TG limit and the general interacting case are given. The use of the hydrodynamic equations allows the full determination of the collective modes in the trap. All the results of the numerical calculations which determine the behavior of these modes in the crossover are presented in section 6. In the first part, we focus on the breathing mode as derived from the SR approach and the LDA to the RQMC data for the ground-state energy. The excitation modes from the hydrodynamic equations of the LL are given in the second part, section 6.3. We finally draw conclusions in section 7.

1D models with 1/x β interactions
In this section and the forthcoming ones, we will be concerned with models of interacting spinless 1D particles with periodic boundary conditions (PBC). In the first quantization, the general Hamiltonians for such models read: Purely dipolar interactions correspond to V (x) = C/|x| 3 . Models of fermions with repulsive interactions falling off as 1/x β have been studied in the past by various authors [38]- [40]. The interactions were taken of the form V (x) = (a 2 + x 2 ) −β/2 to avoid strong short-range repulsion effects. For β < 1 it was shown that the long-range part of the interaction affected the dispersion relation of the plasmon modes, thus driving the system away from the LL fixed point. On the other hand, for β > 1 the plasmon dispersion remained of the form ω(q) = u|q| as in the case of short-range interactions. In the case of interest to us, namely the repulsive dipolar interaction with β = 3, these results show that the long-distance decay of the interaction potential is sufficiently fast to allow us to describe the fermionic system as a LL. We shall argue that this result remains true in the bosonic case, as we will see in section 2.2 that bosonization is also applicable to the description of the low-energy physics of interacting bosons in 1D.

A brief overview of the LL
The concept of LL [41,42] is central in the present understanding of strongly interacting low dimensional systems with both bosonic and fermionic degrees of freedom. For fermionic systems, LL behavior can be established by renormalization group methods, by Feynman diagram methods [43], by bosonized representation of Fermion operators [44] and by Ward identities [45]. Bosons having hard-core repulsions can be mapped on to fermions by the Jordan-Wigner transformation [46], and LL behavior can thus also be established for them. In the case of non-hard-core repulsion, LL behavior is more difficult to establish since, contrarily to the case of fermions, some non-perturbative methods are required. LL behavior has been established for 1D bosonic integrable systems [47]- [54] and also by numerical methods for various non-integrable systems [55]- [57].
In the LL theory, the system is described by a phase Hamiltonian [30,58] where the commutator [φ(x), (x )] = iδ(x − x ). The parameter K is called the Luttinger exponent. Its important property is that it decreases as the interactions become more repulsive.
In the case of a gas of hard-core bosons [59,60], such as the TG gas, the Luttinger exponent K = 1 as in the free-fermion case. With less strongly repulsive interactions, the exponent K > 1 like in the case, for instance, of the Lieb-Liniger gas with finite repulsion. The normal ordered density operator can be expressed as a function of the field φ as (see [30, p 72]):n where A m are non-universal dimensionless constants and n 0 is the average boson density. The fluctuations of the field φ are described by the Hamiltonian (2). The non-universal A m are known only in the case of a few integrable models [61] 9 . The expression (3) for the density is valid for both bosons and fermions. Using equation (3), we can show that if we perturb the original Hamiltonian with a two-body interaction v(x) such that dx|v(x)| < ∞, the new Hamiltonian describing the long-range physics will be of the form (2), with renormalized u and K . This explains why, in the case of fermionic systems, the LL behavior persists for sufficiently large β [39,40]. In the case of bosonic systems, we could imagine separating the interaction into finite-range part plus a long-range perturbation. If the finite-range perturbation is inducing a LL, then, provided that the long-range perturbation decays sufficiently rapidly, this will simply renormalize the velocity of the excitations and the Luttinger exponent. This line of reasoning supports our conjecture of section 2.1, namely that bosonic models with power-law repulsion decaying faster than 1/r at long distance possess a LL theory description at low energies. The boson-field operator ψ B can also be expressed in terms of the fields and φ that enter the Hamiltonian (2). In order to do that, one must introduce a field θ(x) defined by: 7 The expression of the field operator reads: where again the B m are non-universal constants. The physical meaning of equation (5) can be understood from the following heuristic argument [30]. Since the commutation relation of θ and φ is [ where (x) is Heaviside's unit step function, for a state with given φ(x), e iθ (x ) acts by transforming φ(x) → φ(x) − π (x − x ). According to equation (3), the expectation value of the density operator in this state is then shifted by δ(x − x ), i.e. e iθ (x ) creates a particle at position x . Since the operators e iθ(x) and e iθ (x ) commute for x = x , it follows that the operators e iθ(x) behave as boson-field operators.
An interesting property of the Hamiltonian (2) that can be derived from the commutation relation for θ (x) and φ(x ) is that of duality. Indeed, differentiating their commutation relation with respect to x shows that θ (x) and −∂ x φ(x)/π are fields canonically conjugate to each other. Therefore, by exchanging the fields θ and φ, one can map a LL of exponent K on to a LL of exponent 1/K . By straightforward calculations it can be shown that the correlation functions of e inθ and e imφ are algebraically decaying with exponents that depend on the Luttinger parameter. Then duality allows one to show that the product of the critical exponents of e inθ and of e imφ has to be a constant. As a consequence of that, in the K → 0 limit the density-density correlations are decaying slowly, indicating a tendency to crystallization, whereas the single-particle density matrix ψ B (x)ψ † B (x ) is decaying rapidly, indicating that there is no tendency to Bose condensation. When K → ∞ the situation is reversed; there is a strong tendency to Bose condensation, with rapid decay of density-density correlations, indicating no tendency to form a crystal.
The Hamiltonian (2) can also be described from the point of view of a phonon picture [38]. If we start from the Hamiltonian (1) with V some repulsive potential, in the classical limit M = ∞, we can minimize the ground-state energy by choosing a state such that x n = na, with a lattice spacing, a. For M < ∞, after expansion to second order around the classical minimum, with x n = na + u n and with [u n , p m ] = iδ n,m , we find a phonon Hamiltonian of the form: where V is the second derivative of V . The continuum limit of this Hamiltonian, with (x = na) = p n /a and φ(x = na) = u n , leads to the LL Hamiltonian (2) provided that the potential is sufficiently short ranged. The equation for the density (3) is also recovered in this limit. Therefore, a LL can be viewed as a 1D solid melt by zero point phonon fluctuations. The stability of the solid requires a sufficiently long-range interaction to induce a vanishing phonon density of state for q → 0. The case of the Wigner crystal with 1/r forces is marginal [38], with only logarithmic correction to the phonon dispersion. We note that the present strong repulsion limit of the LL is actually independent of the statistics of the particles in (1), thus justifying that in the strong repulsion case, bosons as well as fermions can form a LL state as we have already discussed in section 2.1. Here, we focus on the case of long-range dipolar interactions where the tendency to form a crystal can be well established. In the following sections, we present calculations of the Luttinger exponent, of the velocity of the excitations, and of the structure of the liquid.

The Luttinger exponent K and the velocity of excitations u
At the end of section 2.1 and in section 2.2, we have argued that the low-energy, long-wavelength behavior model of 1D dipolar bosons can be described within the framework of LL theory. Here, we detail the relations connecting the Luttinger exponent and the velocity of the excitations with the density and the compressibility of the microscopic model [42]. These relations are used in section 3 to provide a non-perturbative definition of the Luttinger exponent in terms of the ground-state energy density computed via RQMC. The knowledge of the Luttinger exponent gives us access to the structure factor, the computation of which is recalled in section 2.4. This LL structure factor will then be compared to the one obtained from RQMC in section 3.2. (1) is characterized by Galilean invariance under the boosts p → p + Mv and x → x + vt. This property allows us to express u K as a function of the density [42]. Indeed, we first notice that when the center of mass of the system is moving with a velocity v, the total momentum of the system is P = N mv and the total energy is E = E GS + 1 2 N Mv 2 , where E GS is the ground-state energy measured in the center-of-mass reference frame. Secondly, under a Galilean boost the field operator transforms as [62]

Consequence of Galilean invariance. Hamiltonian
. This implies that the field (x) transforms as π (x) → π (x) + Mv/h. Reporting this latter expression in the bosonized Hamiltonian (2), we find the following expression of the energy: Equating the excess energy (7) to the total kinetic energy N Mv 2 /2 with the total number of particles, N , we obtain This result can be stated as follows: in a Galilean-invariant system, the product u K depends only on the density. Let us note that the present derivation holds both for fermions and bosons and does not require the Luttinger theorem to be applicable. Another useful consequence of Galilean invariance is that the total momentum is given by: 2.3.2. Compressibility. The compressibility 1/χ of a 1D interacting system can be obtained from the energy per unit length e(n 0 ) as The expression of χ as a function of u and K can be obtained by a bosonization procedure [30,42]. To this aim, we first notice that the density of particles at position x is related to φ via ∂ x φ = −πn. If we imagine a uniform compression, then ∂ x φ = −πδn. The energy cost of this compression is straightforwardly obtained from the effective Luttinger exponent as Lh 2π As a consequence, the second derivative of the energy per unit length is simply πu/K . The compressibility is thus related to n 0 u and K by Since the value of K decreases as long as the interactions become more repulsive, we see that repulsive interactions reduce the compressibility of the system, as expected.
Combining the relations (8) and (12), we finally obtain The first expression yields the sound velocity [63]. Since in a non-interacting Fermi gas, (d 2 e(n 0 )/dn 2 0 ) = π 2h2 n 0 /M, the second expression leads to the conclusion that K = 1 in a non-interacting Fermi gas. Therefore, a calculation of the second derivative of the energy per unit length in a Galilean-invariant 1D fluid is the needed ingredient to calculate the Luttinger exponent K .
An approximated expression for K can be derived also by means of a classical argument. Let us consider a system of particles interacting through the potential V (x) = C/|x| β . A classical treatment of the interaction would suggest that these particles will form a 1D crystal, with interparticle spacing a. The particle linear density would be n 0 = N /L = a −1 , whereas the energy per unit length would be: (ζ being the Riemann zeta function). If the particles are confined in a region of order a, the kinetic energy associated with their zero-point motion is overcome by the potential energy (15). The classical potential energy thus yields a good first-order approximation to e(n 0 ). The approximated expression for K using (14) is We observe that at high density K → 0, i.e. the fluid's behavior progressively resembles that of a classical system [27]. This same result would be recovered by considering the phonons of the 1D crystal.

Static structure factor and momentum distribution function
The static structure factor S(q) of the fluid is defined as: where L is the system size, and we are assuming PBC so that correlation functions are translationally invariant.
Using equation (3), we obtain the density-density correlation function as a sum of correlation functions e 2imφ(x) e −2imφ(0) . The following expression for the correlation function is obtained by standard calculations [30] and yields: where a 0 is a short distance cutoff parameter. To obtain the expression of the structure factor S(q), we have to find the Fourier transform of this expression. We note that since our system is periodic, one only needs to consider q = 2π n L with n an integer. Indeed, one first needs to use the generalized binomial expansion (equation (3.6.8)) [64] to obtain a series expansion for the Fourier transform of (18). This is resummed using the definition of the Gauss hypergeometric function 2 F 1 in [64] to find the Fourier transform of (18) in closed form: This implies that: We have to consider two cases, referring to whether 1 − 2m 2 K > 0 or 1 − 2m 2 K < 0. Firstly, if 1 − 2m 2 K > 0, one can safely take the limit a 0 /L → 0 in the hypergeometric function [64] to obtain: where B is Euler's Beta function. As L → ∞, m diverges as m (q = 0) ∼ L 1−2m 2 K . This implies the presence of quasi-Bragg peaks in S(Q) for |Q| = 2πmn 0 , whose intensity is diverging with system size. As K gets smaller, additional quasi-Bragg peaks appear. Secondly, in the case 1 − 2m 2 K < 0, we apply equation (15.3.3) of [64] followed by equation (15.1.20), thereby yielding: Therefore, no divergent Bragg peaks occur in this case. The expression of S 0 (Q) can be easily obtained and results in: where K is the Luttinger exponent.
Within the LL theory we can also calculate the momentum distribution function n(k), which is accessible in experiments via e.g. the analysis of time-of-flight images. The bosoncreation operator in the continuum may be represented [30] as in equation (5). Using the expression for the correlator of the field θ(x), we finally obtain that n(k) ∝ k 1/2K −1 in the L → ∞ limit. The divergence of the momentum distribution at k = 0 is expected to be progressively reduced while K decreases, until it disappears for K < 1/2.
Whether the description of the 1D dipolar Bose gas is consistent with that of a LL, can be checked by determining the Luttinger exponent K from the ground-state energy per particle e(n) and the structure factor S(q) computed by means of QMC techniques which, for a bosonic system, yield exact numerical results for these properties.
In particular, we resort to the RQMC [37]. While we refer to the lecture notes of Baroni and Moroni for a detailed description of the method [65], here we sketch just a few basic concepts. In RQMC, the probability distribution trial [66] method samples the mixed probability distribution trial , being the ground-state wavefunction. In RQMC, for a sufficiently large number of time slices M and sufficiently small time steps in the middle of the path, the ground-state probability distribution is sampled instead. This allows a direct and simple way to obtain pure estimations of ground-state space diagonal operators, such as the static structure factor S(q). In the following, we focus on the purely dipolar bosonic system, namely the α = 3 case of the more general Hamiltonian (1).
The 1D dipolar Bose gas can be described as a system of N atoms, with linear density n, mass M, and with permanent dipole moments arranged along and orthogonal to a line. We consider here the case in which all the dipoles are parallel to each other, resulting in purely repulsive interactions. In particular, we consider the situation of the Stuttgart experiment [5], where an external polarizing magnetic field is present and only one spin component is effectively entering the dynamics 10 .
The system is characterized by the strength of the interactions C dd , resulting from either magnetic C dd = µ 0 µ 2 d or electric C dd = d 2 / 0 dipoles 11 . An effective Bohr radius can be defined from C dd as r 0 ≡ MC dd /(2πh 2 ) . The Hamiltonian in effective Rydberg units This formulation clearly shows that the model is entirely specified by the dimensionless coupling parameter nr 0 , so that in the high-density limit the system becomes strongly correlated and a quasi-ordered state occurs, where the potential energy dominates. The liquid state is instead reached in the low-density limit when the system eventually can be mapped on to a fluid of non-interacting fermions, the TG gas. We study the system Hamiltonian as a function of density in order to follow the evolution from the quasi-ordered state to the TG limit. Our trial wavefunction is made up of products of two-body Jastrow factors trial = i< j e u(|z i −z j |) . Here, we choose to focus on the long-range behavior of the system and use 10 Interesting analogies with the Einstein-de Haas effect have been pointed out [92,93], these effects are intimately related to the fact that the dipolar interaction does not allow the conservation of the spin and the orbital angular momentum independently, as for the usual contact interaction, but the total one L + S. 11 µ d and d are the magnetic and electric dipole moments and µ 0 and 0 are the vacuum permittivities.  (26) fitting the RQMC data, once the time step and the size effects have been removed. Long-dashed and dotted lines: high-and low-density energy limits, respectively. the trial wavefunction: which for K = 1 (that is, in the very-low-density limit) recovers the spinless free fermionic wavefunction. Different choices of wavefunctions, as for example the product of Gaussians centered on the 'crystal sites', relevant at high density, have led to different time step extrapolations for ground-state energy and yet did not lead to significant differences in the computation of the static structure factors. We performed simulations for N = 40, 60, 80, 100 up to, in selected cases, N = 200 bosons placed in a square box with PBC. To reduce the inevitable size-effects, the interaction potential V PBC = L V (x + L) has been resummed over ten simulation boxes of length L, instead of resorting to the more common truncation potential method [67]. Finally, the energies have been extrapolated to their thermodynamic limit after using the relation E N = E ∞ + b/N a with a = 2 or 3, once their dependence on the time step has been removed.

Ground-state energies
The energy per particle e p obtained from RQMC simulations is shown in figure 1. It reproduces the energy per particle of a free spinless Fermi gas ζ (3)(nr 0 ) 3 with ζ (3) = 1.202 05 at low densities, and the energy density of a classical lattice of dipoles in the high-density limit [28,68].  (14) after differentiating twice expression equation (26) fitted from the RQMC data. Squared symbols: K (n) obtained from the slope of the structure factor (q → 0). The upper (dotted) curve represents the classical limit where the kinetic energy for zero-point motion of the particles is neglected in the calculation of the total energy per unit length. The thermodynamic energy per particle in Rydberg units can be represented as an analytical function of nr 0 : Expression (26) can be useful for further theoretical investigations. As a first application, differentiation of the energy per unit length e(nr 0 ) = nr 0 e p (nr 0 ) leads to the determination of K from equation (14), as a function of nr 0 . The resulting K (n) is displayed in figure 2. K (n) fulfills the TG nr 0 1 limit K (nr 0 1) → 1. For comparison, the classical prediction is also shown and represented by the upper curve. The velocity of the excitations as obtained from equations (14) and (8) is instead displayed in figure 3, again together with the classical and the TG limit.

Static structure factor from RQMC
As already briefly stated above, we resort to RQMC because it allows simple access to pure estimator of the static structure factor S(q). The low-energy behavior provides an independent way to extract the Luttinger exponent by means of equation (23). Figure 4 shows S(q → 0) for different sizes corresponding to N = 40, 60, 80 and 100 particles, and in the case nr 0 = 100. Note that from now on wavevectors q will be intended in units of n. The RQMC data support the LL prediction of a linear dispersion for S(q → 0) = K q/(2π) (dashed line) with a Luttinger exponent K = 0.11. The K -values fitted at any given density from S(q → 0) are displayed as squares in figure 2, falling in excellent agreement with the results obtained from the second derivative of the energy. This is of course a consequence of the consistency of both the RQMC simulation and the LL theory, since the compressibility obtained from long-wavelength behavior of the fluid structure and from the ground-state energy give the same result.  If LL theory were indeed to describe the low-energy behavior of the dipolar Bose gas in the whole crossover, we also expect to find in the RQMC data for S(q) the same quasi-Bragg peaks at q = 2mπ which are predicted by equations (20) and (21). The heights of the peaks should also scale as AN 1−2m 2 K , with K being the exponent already fitted at the given density from either S(q → 0) or d 2 e(n)/dn 2 (notice that the constant A cannot be determined within the LL theory). Finally, we expect that, according to equation (21), all the curves representing S(q)/N 1−2m 2 K versus N (q − 2mπ )/40 for different N values should collapse in a single curve. The RQMC data are found to meet all these expectations, as they are displayed for two selected densities nr 0 = 1000 in figure 5 and nr 0 = 100 in figure 6. In the case of nr 0 = 1000, the results are detailed for the three peaks at q = 2π m with m = 1, 2 and 3. In the case of nr 0 = 100, the results are detailed for the two peaks at q = 2πm with m = 1 and 2. S(q)/N 1−2m 2 K is displayed versus N (q − 2mπ )/40.  Figure 6. The same as in figure 5, but for nr 0 = 100. Here also, S(q)/N 1−2m 2 K versus N (q − 2mπ )/40 to evidence the scaling behavior. From left to right, the cases with m = 1 and 2 are displayed, corresponding to the peaks at q = 2π and 4π . Note that q is in units of n.
As can be seen, deviations of the LL prediction from the RQMC data are observed since, at variance with the LL curves, the RQMC peak is asymmetric. A possible explanation for this discrepancy is that the peak asymmetry is a consequence of the parabolic dispersion of the interacting bosons. In bosonization indeed, the spectrum is linearized. Band-curvature terms become nonlinear in or φ and have, in principle, to be added to the Hamiltonian equation (2). These terms are irrelevant in the RG sense though, and are therefore neglected. The calculation of the real-space correlation functions with the fixed point Hamiltonian yields the expression (18), which is even in x and real. This produces an approximate structure factor that is symmetric around q = 2mπ . However, we know that in the case of the free fermions, where the structure factor can be computed exactly, there is no such symmetry. Asymmetric structure factors are also obtained within exact calculations on fermionic Calogero-Sutherland models [69] which indeed possess a nonlinear spectrum. This suggests that the neglect of the nonlinearity of the dispersion in bosonization results is at the origin of the symmetry of the structure peaks. The terms which are nonlinear in φ or added to the Hamiltonian (2) can break the symmetry φ → −φ and render the correlator in equation (18) complex. This in turn, leads to a structure factor which is no longer symmetric around 2mπ. Using the phonon picture of the LL discussed in section 2.2, this peak asymmetry could be alternatively viewed as the occurrence of anharmonic phonon effects.

The inhomogeneous LL: LDA and hydrodynamic approach
In the previous sections, we have established the LL behavior of a 1D gas of dipolar bosons with PBC. In this section, we investigate the signature of the LL effects of the 1D gas of interacting particles in the trapped conditions which are closer to the experimental realizations. In section 5 we will consider in particular the most common case of a confinement due to a harmonic potential V (x) = Mω 2 0 x 2 /2 in the axial direction. If the trapping in the transverse plane is sufficiently tight with respect to that in the longitudinal direction, the transverse excitations are frozen. In such a regime, the system is quasi-1D. If the longitudinal potential is sufficiently shallow, the dynamical behavior of the confined gas can be determined by resorting to an LDA to the equation of state, in which the energy of the inhomogeneous system is a functional of the local density n(x), i.e. it is expressed as the integral over x of the energy density of the homogeneous gas with density n(x). By construction, the validity of such an approach is limited to dynamical effects in which the typical range of variation for n(x) is much larger than the average interparticle distance in the longitudinal direction. In section 4, we show how LDA can be exploited to yield an effective inhomogeneous LL theory. In section 4.2, we will provide the equations that are necessary to compute the spectrum and correlation functions of the inhomogeneous LL. These results will be applied to the dipolar gas in section 6.

LDA
Formally, the LDA can be derived from a density functional E[n(x)] which gives the groundstate energy for a given average density n(x). The functional that has to be minimized is then: where V (x) is the trapping potential, and λ is a Lagrange multiplier that fixes the particle number. In LDA, the functional E[n] is assumed to be a local functional E[n] = dxe(n(x)).
Varying the functional with respect to the particle density n(x), one obtains the equation: where the chemical potential appears as the functional derivative: If an analytic expression of µ as a function of n is given, equation (29) would allow to find n(x) by inverting the relation equation (28). A simple physical interpretation of equation (28) is as follows. If we assume that the density n(x) is a slowly varying function of the position, we can imagine the system as being cut into segments of size x in which the density is almost constant. In each of these subsystems, quasi-equilibrium conditions set in, with average density n(x). This implies that the density-(and thus space-) dependent chemical potential is adjusted so that no current is flowing from one of the subsystems to its neighbors.
Turning to the LL exponents, if we assume that quasi-equilibrium conditions hold in each subsystem, then K and u can be determined as functions of the local density by the same relation valid at equilibrium [30,42]. Within LDA, we thus obtain: while the system is described by the following Hamiltonian: Let us make a few remarks. Firstly, if we have a smooth variation of the potential δV (x), after completing the square in equation (32) we can obtain the variation of the density in linear response as: After comparing with equation (28), equation (31) is recovered. Till now, we have not discussed the boundary conditions imposed on the fields φ and at the edges of the system. To derive these conditions, let us note that sincen = −∂ x φ/π , the continuity equation ∂ tn + ∂ x j = 0 implies j = ∂ t φ/π [30]. Imposing that no current is going in or out at the edges of the system, therefore j (±R, t) = 0, ∀t amounts to requiring [70]: Let us note that from the equations of motion we have j (x) = u(x)K (x)π = π 2 n(x)/m , and that at the edge lim x→R n(x) = 0. Thus, these boundary conditions impose no constraint on , contrarily to the case of the homogeneous gas.

Equation of motions and correlation functions 4.2.1. Equations of motion and eigenmodes.
Using the equations of motion derived from the Hamiltonian (32), we find that: Combining these two equations and using the boundary condition (34), we have: We note that if we use the relation ∂ x φ = K (x)/u(x) in equation (38), we obtain a stationary solution. The Fourier decomposition of (38) yields: with boundary conditions ϕ n (±R) = 0, and the orthogonality relation: The above equations are related to those derived in [71] using the hydrodynamic approach. The solution of the eigenvalue equation (39) gives access to the eigenmodes of the trapped gas. We also note that the eigenvalues and eigenfunctions of equation (39) satisfy a variational principle [72]. Namely, if f (x) is a function satisfying the boundary conditions f (±R) = 0 and orthogonal to the first k eigenfunctions of (39) with respect to the scalar product (40), then the inequality holds:

Expansion of the fields.
In terms of the eigenfunctions {ϕ n (x)}, which are also eigenfunctions of the Sturm-Liouville operator, we have the following expansion for the fields: where [a n , a † m ] = δ n,m and the second term in (43) comes from the addition of N particles to the system. Indeed, with this form of φ and , we have the commutator: where the last line is a consequence of the completeness of the basis [72,73] formed by the eigenfunctions of the Sturm-Liouville operator {ϕ n (x)}. The Hamiltonian, expressed in terms of the boson creation and annihilation operator reads: A decomposition similar to (43) can also be introduced for the field θ(x), π (x) = ∂ x θ : where θ 0 is such that e iθ 0 N e −iθ 0 = N + 1, and: Indeed, the above equation implies that: with: i.e. the functions {ϑ n (x)} form an orthogonal family of functions. A variational principle can also be derived from equation (51) namely: for a function g orthogonal to the first k functions ϑ n with respect to the inner product of equation (52). It is interesting to note that the functions ϑ n satisfy the same differential equation and the same orthogonality condition as φ n after the transformation K (x) → K −1 (x). Thus, the duality θ(x) → φ(x) K → K −1 that exists in the homogeneous case (see section 2.2) is changed in the inhomogeneous case into θ (x) → φ(x) and K (x) → K −1 (x). However, as in the case of the homogeneous LL with open boundary conditions, the boundary conditions satisfied by ϑ n are not the same as those satisfied by ϕ n . As a result, the duality transformation cannot be used to obtain the correlation functions of the field θ from those of the field φ or vice versa. However, the variational principle (53) is immediately deduced from (42), by means of the duality transformation.

Correlation functions.
The decomposition of the fields equations (43)-(48) permits us to recover easily the linear response result equation (33) and the relations (30) and (31). To this aim, we calculate the zero-temperature Matsubara correlation function: which in terms of the functions {ϕ n (x)} is: Knowing thatn(x) = − 1 π ∂ x φ(x) + · · ·, the density-density response function is readily obtained as: where we can neglect the higher-order oscillating terms, since they are expected to oscillate rapidly over the typical length of variation of the applied potential. Using the relation (49) and the closure relation (52), we find: so that (33) is recovered.

The harmonic trapping potential: non-interacting versus interacting case
In this section, we apply the results of the preceding section 4 to the case of harmonic trapping. We will consider a solvable model which is convenient to obtain the eigenmodes of the trapped gas in the two extreme limits of low and high density. But before proceeding further, let us show that within the LDA, we can recover the exact frequency of the center of mass oscillations of a system with Galilean invariant interaction trapped in a harmonic potential. We start from the differential equation (38) which, using (31), can be rewritten as: If we make φ(x) = λn(x), we find that ∂ x φ∂ n µ = −mω 2 0 λx and thus the differential equation is satisfied with φ = λn(x) and ω = ω 0 . Since n(±R) = 0, φ automatically satisfies the boundary conditions (34). The interpretation of this equation is quite straightforward. If we write φ(x, t) = φ(x − Acos(ω 0 t), 0) which describes an oscillation of the center of mass of frequency ω 0 , and take the derivative with respect to the position, we recover precisely the solution above. This establishes that the LDA correctly reproduces the center of mass oscillations.
As a particular case let us take [74]: This form of the Luttinger parameter is obtained for instance in the case of µ(n) = gn γ . Indeed, using equations (30), (31) and (28), we find u(x) and K (x) in the form of equation (59), with: As the chemical potential has to be an increasing function of particle density, we must have α > −1/2. For γ = 2, we recover the trapped free fermions or the trapped TG gas case. For γ = 3, we obtain the dipole gas in the high-density limit. We will first discuss the free fermion or TG gas case, owing to its simplicity, and then we will turn to the more general case.

Non-interacting fermions or TG bosons case
In the case of non-interacting fermions or TG bosons in a harmonic trap, we have α = 0 and K = 1. The eigenvalue equations (39) and (51) then become identical, and can be solved by a simple change of variables [75], s(x) = x 0 (dx /u(x )). The most general expression of the eigenfunction can be obtained [75] in terms of the Chebyshev polynomials of the first kind T n and second kind U n as: with ω n = nω 0 , and n 1. This leads to the following expressions for the fields: Using these expressions, it is a straightforward task to compute the correlation functions [75]. We note that computing the current, j (x) = u(x) (x), shows that the current is proportional to (1 − x 2 /R 2 ) 1/2 and thus vanishing on the edges as expected.

General case
Here, we study the general case where α = 0. In this case, an approach based on Sturm-Liouville equations is needed. Upon substitution of K (x) = K 0 (1 − x 2 /R 2 ) α , the Sturm-Liouville equation (39) becomes: whose solutions ϕ n (x) can be written in terms of Gegenbauer (also called ultraspherical) polynomials [64,76]. More precisely, we have from equation (22.6.6) in [64]: Note that in the free fermion case α = 0, using equation (22.5.34) of [64], the Gegenbauer polynomials reduce to the Chebyshev polynomials U n (x/R) previously obtained in equation (64). The expansion (68) was obtained previously in [74]. The factor A n comes from the normalization integral of Gegenbauer polynomials, equation (22.2.3) in [64] and is given by: Using the completeness relation of the Gegenbauer polynomials, it is straightforward to check that the solutions (68) form a complete system on [ − R, R]. Therefore, there are no other eigenfunctions than the ones of equation (68). In terms of the parameter γ , the eigenvalues are: In section 6 we will compare those eigenvalues with those derived from the variational approximation. The first eigenvalue n = 0 is ω 0 , which is exactly the trapping frequency. This result is expected since it corresponds simply to harmonic oscillation of the center of mass of the trapped cloud. The second eigenvalue is ω 0 (2 + γ ) 1/2 which will be again in agreement with the variational method. In particular, in the case of γ = 3, we obtain the asymptotic limit of the dipolar gas. The resulting expression of φ(x) is simply obtained by substituting (68) in (43). Interesting physical information comes from the expression of cos 2φ(x) . This can be derived as a series of squares of Gegenbauer polynomials. In particular, we find: n n!(n + α + 1) (n + 2α + 2) 1 √ (n + 1)(n + 2α + 1) Using a formula for the product of Gegenbauer polynomials [77], the sum (72) can be expressed as an integral. The details can be found in appendix B where the expression for the general correlation function is given. The main result of the calculation is that φ 2 (x) is a function of x/R which is minimal near the edges. As a result, Friedel-type oscillations [78]- [80] will be observed in a 1D gas with open boundary conditions. These oscillations will decrease when going to the center of the trapped gas. In the case of α = 0, φ 2 (x) will display logarithmic behavior. Instead, we now turn to the family of functions {ϑ n }. They satisfy the Sturm-Liouville equation (51): whose solution is [64]: Finally, the normalization of B n is found from (22.2.3) in [64] in the form: From there, the expression of θ (x) is found as: In the limit α → 0, using equations (22.5.4) and (22.5.28) from [64], we recover the expression of the θ field in terms of Chebyshev polynomials of the first kind [75].

Collective excitations of the dipolar 1D gas
In this section, we discuss the collective excitations of the dipolar 1D gas by combining the LDA and SR approaches, as corroborated by RQMC. Using the RQMC data, we first discuss the breathing mode as obtained from an SR approach with the use of the LDA density profiles.

Equation of state
An analytic expression of the density profile in a given potential can be obtained by inverting the equation of state (29). To this aim, we first resort to the RQMC data, and obtain the dependence of the energy per particle on the density using the fitting (26) in section 3.1. The chemical potential can be obtained from the fitted expression as µ RQMC (n) = (1 + n(∂/∂ n ))ε(n). Using this expression and the relations (14), we obtain easily the eigenvalue equations (39) that correspond to a given density.
Under these conditions, the ground-state density profile n(x) of the trapped dipolar gas is obtained by plugging µ RQMC (n) into the equation (28)

Breathing mode from an SR approach
We first determine the effect of the interactions on the frequency of the lowest compressional (breathing) mode of the trapped gas by an SR approach. This mode is coupled to the ground state by the operatorX 2 = N i=1 x 2 i . An upper bound for the frequency ω B can be derived by the SR and satisfies the inequality [81]- [83]: where m 1 is the first moment of theX 2 m 1 ( In order to determine m −1 (X 2 ), we follow the useful procedure of [83] and notice that the perturbation −σX 2 , which couples to the density, has the overall effect of shifting the strength of the external trapping potential by the perturbation strength σ , namely Mω 2 0 /2 → Mω 2 0 /2 − σ . Thus, the static polarizability χX 2 (0) reads χX = N −1 x 2 n(x)dx. Let us note that the same inequality (78) is valid also within the LL hydrodynamic approach. Using the operatorX 2 = dx g(x)∂ x φ(x), one obtains:  (79), while the symbols represent the results of LL hydrodynamics equation (39).
and the variational principle of equation (53) is recovered. For µ(n) = λn γ , a closed form expression of B in equation (79) is obtained by the following scaling argument. Firstly, using equation (28) to compute the particle density, and noting that the particle number N is fixed, we have that R ∼ ω −2/(γ +2) 0 . This gives us . Using the logarithmic derivative to express (79), we find B = ω 0 (2 + γ ) 1/2 . In particular, in the low-density limit, we have γ = 2 and thus B = 2ω 0 , while in the large-density limit γ = 3 and B √ 5ω 0 . For intermediate densities, we have to resort to a numerical estimation of equation (79) using the LDA density profile derived from equations (28) and (15). The result is represented by the dashed line in figure 8, showing the smooth evolution of the breathing mode frequency ω B from the low-density to the high-density regime. Figure 8 enables e.g. the identification of the interaction regime of the trapped dipolar gas, by means of one of the best handled experimental probes available with cold atomic (molecular) quantum gases [84,85]. We remark that, in principle, equation (79) only yields an upper bound for the frequency of the lowest compressional mode. Moreover, computing the frequencies of the higher modes of the trapped gas by the SR approach is less straightforward than in the case of the compressional mode. In section 6.3, we switch to the hydrodynamic LL theory which simplifies greatly the computation of the frequencies of the higher modes. (ω n /ω 0 ) 2 versus N (r 0 /a ho ) 2 from equation (39). From bottom to top: the modes with n = 3, 4 and 5 as in the legend. The lines are a guide to the eye.

Collective modes from the hydrodynamic Luttinger equations
In this section, we apply the LDA treatment developed in sections 4.1 and 4.2 to the dipolar Bose gas. Using the results of section 4.2.1, the eigenfrequencies are exactly known in the two asymptotic limits of low and high density. Indeed, insertion of µ(n) ∝ n γ in equations (28) and (39) yields solutions of the form ϕ n (x) = A n (1 − x 2 /R 2 ) 1/γ C (1/γ +1/2) n (x/R), the associated eigenvalues being ω 2 n = ω 2 0 n[2 + (n − 1)γ ]/2 [36,83], where C (α) n are Gegenbauer polynomials and A n normalization factors [74]. Thus, in the two opposite low-density (γ = 2) and highdensity (γ = 3) limits, one obtains, respectively, ω 2 n = n 2 ω 2 0 and ω 2 n = n(3n − 1)ω 2 0 /2. For intermediate densities, equation (39) has to be solved numerically inserting as ingredients the LDA density profiles and the analytical expression for ∂ n µ obtained from equation (28) and the RQMC fit equation (15). This computation was performed with the help of the Sledge program available online [86]. In the numerical solution, we have taken special care of the finite mesh-size effects for best accuracy. The values of the breathing frequency ω B /ω 0 obtained from equation (39) are represented in figure 8 by the filled symbols, and agree up to the second digit with the SR result. The eigenfrequencies of the higher modes with n = 3, 4 and 5 are displayed in figure 9 as functions of N (r 0 /a ho ) 2 , showing the same smooth crossover behavior between the two opposite TG and DDW regimes. The exact frequencies in these asymptotic regimes are recovered by the numerical calculations.
Thus, a comparison of the measured excitation frequencies with those predicted from (39) with our RQMC data provides a test of LL behavior in trapped 1D dipolar Bose gases within the validity of LDA.

Conclusions
Bosons arranged in a 1D geometry in the presence of repulsive long-range interactions realize an unusually super-strongly correlated quantum system, which is useful for devising novel quantum phases as well as technological applications. We have devoted the present paper to provide a detailed overview of a number of our more recent findings, which put 1D dipolar Bose gases under a new perspective, from both physical and methodological points of view.
From the viewpoint of the underlying physics, we can conclude that the ground-state and the low-energy excitation behavior of 1D dipolar Bose gases evolves from a TG gas at low densities to a quasi-ordered state at high densities, which we have called DDW. The unifying theoretical framework describing this crossover is the LL theory, which also provides analytical expressions for the correlation functions once the Luttinger parameters K and u are determined. The Luttinger parameters have been extracted from RQMC calculations on a system with PBC. The correlation functions obtained from LL theory have then been systematically compared with RQMC data, quantitatively demonstrating that the system is always in a very strongly correlated LL regime with K < 1 at all densities.
Turning to the comparison with experiments, we have to keep in mind that actual experimental systems will be experiencing a confining potential. One can imagine experimental situations where the effect of the confinement can be neglected and thus our predictions for the structure factor or the momentum distribution of the homogeneous LL can be tested by e.g. Bragg spectroscopy [87]. However, in this paper, we have focused our attention on the simplest experimental probe that can be imagined, namely a measurement of the eigenfrequencies of the collective trap modes across the whole crossover. We have shown that for a system trapped in a shallow potential, an effective LL hydrodynamic theory can be set up. We have discussed in detail the general equations and their range of validity, pushing the derivation of analytical solutions as far as possible. As often occurs with cold atomic gases, the use of the LDA can be considered a first step towards a quantitative prediction. Within this approximation, we have determined the lowest collective modes of the harmonically trapped dipolar gas in the crossover. For the breathing mode, we have found remarkable agreement with the results of an SR calculation.
In [32], we have provided estimates for the system parameters, suggesting that the 1D limit considered in the present article is within the reach of current experimental techniques. In particular, while in the case of Cr atoms possible ranges of nr 0 would remain limited to the TG regime, SrO molecular gases [88] are among promising candidates possessing sufficiently large electric moments. Under these conditions, values of nr 0 spanning from the TG to the DDW regime seem to be realistic and compatible with the realization of an effective quasi-1D condition. Interesting applications have already been proposed, such as the design of quantum memories [35], which indicate the possibility of realizing these systems.
Concerning the possible theoretical extensions of the present work, a first direction is the study of arrays of coupled 1D dipolar gases [89], which is relevant for dipolar bosons trapped in 2D optical lattices. Solid-like phases, as well as more exotic sliding density waves can be obtained in that case as a result of interchain dipolar interaction. Another related direction of work is the study of a quasi-1D situation in which the transverse trapping is not tight enough to allow only one transverse level to be accessible. This case can be viewed as the realization of a multicomponent dipolar gas, and can be treated using a formalism similar to that of [89]. This study would be interesting to understand the limit where many transverse modes are occupied, as it occurs e.g. in the case of dipolar bosons in a cigar-shaped trap. A further possible case study is concerned with non-polarized dipolar gases, in which some kind of spin-density separation can be expected [31]. Finally, the out-of-equilibrium dynamics of dipolar gases is an open issue, which could be induced by suddenly changing the orientation of the field polarizing the dipoles.