Extended fermionization of 1D bosons in optical lattices

We present a model that generalizes the Bose–Fermi mapping for strongly correlated one-dimensional (1D) bosons trapped in combined optical lattice plus parabolic potentials, to cases in which the average number of atoms per site is larger than one. Using a decomposition in arrays of disjoint strongly interacting gases, this model gives an accurate account of equilibrium properties of such systems, in parameter regimes relevant to current experiments. The application of this model to non-equilibrium phenomena is explored by a study of the dynamics of an atom cloud subject to a sudden displacement of the confining potential. Excellent agreement is found with results of recent experiments, without the use of any adjustable parameters. The simplicity and intuitive appeal of this model make it attractive as a general tool for understanding bosonic systems in the strongly correlated regime.


Introduction
Cold bosonic atoms in optical lattices have recently been used to create quasi-one dimensional (1D) systems [1]- [6]. In such experiments, arrays of 1D tubes are realized by first magnetically trapping a Bose-Einstein condensate (BEC) in a parabolic potential, and then imposing upon it a deep 2D optical lattice, which restricts atomic motions to 1D. These defect-free highly controllable atomic systems offer an excellent opportunity to directly study strongly correlated regimes.
For low densities or large interaction strengths, a 1D gas of ultracold bosons behaves as a gas of impenetrable particles, known as a Tonks-Girardeau (TG) gas. Reference [7] shows that there is a one to one mapping between the eigenenergies and eigenfunctions of TG bosons and those of non-interacting fermions. This correspondence is known as fermionization. Two recent experiments [4,5] successfully reached this parameter regime. In [5] the TG regime was achieved by adding an optical lattice in the tubes' direction which increases the effective mass and therefore the ratio between interaction and kinetic energy. When the lattice is present and for large enough interactions a commensurate homogenous system not only fermionizes but also undergoes the superfluid to Mott insulator (MI) transition [8]. For an inhomogeneous system, the commensurability criterion changes; in particular, for the experimentally accessible case of confinement by a parabolic potential, the MI can be realized for any number of particles [9,10].
In the TG regime with at most one particle per lattice site, single-particle solutions of the periodic plus parabolic potential have been successfully used to describe equilibrium and non-equilibrium properties of the system, both in the presence and in the absence of the MI [9,10]. However, 1D experiments have been realized in a parameter regime where on-site particle densities are larger than one and standard fermionization treatments are inapplicable [3,6]. Here, we present a novel theoretical model for treating this strongly correlated regime, previously studied only by means of numerical simulations [11]- [14]. We show that the gas develops a layered domain structure, which can be approximated by a set of independent TG gases. Quantum phase transitions are found to occur in a layer-by-layer basis. This extended 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT fermionization (EF) model has broad implications for the structure of 1D quantum gases. Comparison with exact Monte Carlo simulations shows that, in the appropriate parameter regimes, the EF approach accurately describes equilibrium properties such as the density profile, the momentum distribution and the ground-state energy. Moreover, we show that a time-dependent version of the EF model reproduces the dipole oscillations of atoms subject to a suddenly displaced potential, in excellent quantitative agreement with recent experiments, without the use of adjustable parameters [3,6]. Some striking features of the experimental results are therefore explained for the first time. The insight and simplicity of the EF method make it a useful new tool for understanding strongly correlated bosonic systems.

The Bose-Hubbard (BH) Hamiltonian
The BH Hamiltonian describes the system's dynamics when the lattice is loaded so that only the lowest vibrational level of each lattice site is occupied and tunnelling occurs only between nearest-neighbours [15] Hereâ j is the bosonic annihilation operator of a particle at site j,n j =â † jâ j , and the sum i, j is over nearest neighbours. In equation (1) the hopping parameter decreases for sinusoidal lattices with the axial depth V o as [9] is the lattice recoil energy, m the mass of the atoms, λ/2 the lattice spacing, and β a dimensionless constant proportional to the atomic s-wave scatting length a s . In current experiments, where the 1D geometry is obtained by tightly confining the atoms in two dimensions, β equals β = 4/( √ 2π)(a s /λ)(V ⊥ /E R ) 1/2 . Here V ⊥ is the lattice depth in the transverse directions. The energy is proportional to the curvature of the parabolic potential.

Homogeneous lattice
For a homogeneous system ( = 0) with N atoms and M sites, the spectrum is fully characterized by the ratio, γ = U/J , 4 and the filling factor, N/M = (n − 1) + N/M. Here n is the smallest integer larger than N/M and N < M. If the lattice is commensurately filled ( N = 0), the ground state is an MI for γ > γ c (n − 1), where γ c ≈ 4 [17]. For the incommensurate case, N > 0, if γ/n γ c , the population of states with more than n atoms per site is suppressed by a factor of ≈J/U, so the BH system may be treated as one with two states per lattice site (containing n − 1 and n atoms respectively). In this case, it is well known [18] that the BH model maps to the XX quantum spin model. Thus, the standard Bose-Fermi mapping techniques derived for systems with n = 1 can be applied, if J is replaced by nJ. Here we present a generalization of standard fermionization ideas, to treat the systems that are of greatest interest in current experiments: those for which n is greater than 1 and varies across the system.

Low density regime: standard fermionization
When the parabolic trap is present ( > 0), the density profile of the atomic cloud is determined by an interplay of U, J, and N. The system is fermionized (with at most one atom per site) if [9] γ γ c , U > ( The first inequality is the same as for a homogeneous lattice, but the second is specific to trapped systems: it suppresses multiple occupancy of single sites. In this fermionized regime, the density at the trap centre increases as J decreases and, when the condition 2J ((N − 1)/2) 2 is satisfied, central trap sites begin to have unit filling [9,10]. If the inequality is also satisfied, the ground state is a unit-filled Mott state in all N sites. In this case, number fluctuations occur mainly at the edge of the density distribution, due to tunnelling of atoms to empty sites. Such fluctuations are proportional to N, which is the trap gradient at the site (N − 1)/2. Number fluctuations at the trap centre are due to a small admixture of particle-hole excitations in the ground state for finite J [19]. By using perturbation theory, it can be shown that this admixture is proportional to 2 √ 2(J/U).

High density regime: the EF model
For ((N − 1)/2) 2 U it is energetically favourable for atoms to pile up at the trap centre. In fact, superfluid and MI phases with different filling factors can coexist, due to the interplay between on-site interactions and the external potential [11,15,20]. In the trivial case J = 0, Fock states with a definite number of atoms in each site are eigenstates of the Hamiltonian. The density profile becomes a 'cake' structure with maximal occupation n max at the trap centre if where [20]. In fact, this profile is most conveniently viewed as a 'layer cake' structure of n max vertically stacked horizontal layers (see figure 1(a)). The number of atoms in the nth horizontal layer, N n , can be shown to satisfy subject to n max n=1 N n = N. In the limiting case of J = 0, atoms are frozen within each layer. Each layer can then be thought of as an independent Mott state with unit filling. For J > 0, atoms are no longer frozen, and it is not obvious that the layers should remain distinct. However, if the number fluctuations in adjacent horizontal layers do not overlap in space, all layers can still be thought of as independent. This situation pertains to many cases of experimental interest in the strong coupling regime, in particular when   In this case, all layers except the top can be viewed as Mott states with unit filling factor. In analogy to equation (2), the first inequality insures that within those layers, the average kinetic energy required for one atom to hop onto an already occupied site is insufficient to overcome the potential energy cost; thus particle-hole excitations, that is couplings between the layers, are suppressed by a factor of ≈nJ/U. The second inequality guarantees that in those layers number fluctuations are localized at the outer edges of each layer, leading to the presence of alternating regions of superfluid and insulating phases in the total density profile. Equation (6) guarantees the validity of using fermionization techniques layer by layer; in fact it insures that atoms in a given layer n sit on a well-defined plateau of n − 1 atoms frozen in a Mott state (which in turn is conveniently decomposed in terms of n − 1 vertically stacked unit-filled Mott states). Then, in all layers except the top, excitations take place only at the sites at the outer edges of the layer and involve states with n and n − 1 atoms per site only.
then, to a very good approximation, atoms in all layers can be treated as TG bosons with an effective hopping energy nJ. The top layer is governed by Jn max , N n max , and ; these determine whether it is a Mott state or not. Atoms in this layer do not need to be localized in a Mott state, since they do not have to form any plateaux for atoms in an upper (n max + 1) layer. Equation (7) insures the suppression of the admixture of states with n max + 1 atoms. Thus, excitations in the top layer involve states with n max and n max − 1 atoms per site only. Under these conditions, singleparticle solutions can be successfully used to obtain expressions for all many-body observables. We refer to this generalization of the Bose-Fermi mapping to spatially varying density distributions as EF.

Computation of physical observables
Before proceeding to apply the EF model to the system of [6], we shortly review how to compute physical observables by means of standard fermionization techniques, and their extension to the EF model. Fermionization of TG bosons in an optical lattice is realized by mapping equation 1 into a non-interacting fermionic Hamiltonian by using the Jordan-Wigner transformation (JWT). The latter expresses the bosonic operators {â j } in terms of fermionic operators {ĉ j } a j =ĉ j exp iπ Once the non-interacting fermionic ground-state, | F , is known the different many-body observables of the bosonic system are calculated by applying the JWT and Wick's theorem. Local correlations such as density n j and site number fluctuations n 2 j are quite simple: Here, j is the lattice site index. In the EF model, fermionization is applied separately to individual layers, and results for the various layers are then combined. For the particular case when there are only two layers in the system, as in [6] (see below), the total density is given by Here | (n=1,2) F are the many-body ground states of a system of N n=1,2 non-interacting fermions with hopping energies nJ. Number fluctuations for each layer can be also calculated using | (n=1,2) F in equation (9).

7
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Application of the EF model: explanation of the NIST 'damping' experiment
In the following we show the success and the limitations of the EF approach, by applying it to a system with the same parameters as the ones used in an experiment recently performed at NIST [6].
In the NIST experiment [6] approximately N T = 1.4 × 10 5 87 Rb atoms were trapped in an array of 1D tubes with N 80 atoms in the central tube. An additional periodic potential was added along the direction of the tubes and its depth V o was varied for different experiments. The frequency of the parabolic potential in the tubes' direction was such that = 7.4 × 10 −4 E R , with E R the photon recoil energy. Here we focus on the central tube only and consider V o > 2E R , where the tight-binding Hamiltonian equation (1) is expected to be valid. For these parameters, equation (4) yields n max = 2.

Static properties
In figure 1 we show comparisons between the density and number fluctuations for the central tube calculated by using the EF approximation and exact quantum Monte-Carlo numerical simulations based on the Worm algorithm [21], for lattice depths V o /E R = 11 (a), 9 (b), 7 (c), 5 (d), 3 (e), 2 (f). In the numerical simulations the temperature is 0.01J. In panel (a), the atoms belonging to the second layer have been shaded to explicitly show what we mean by horizontal layers. In all the plots the dotted (black) and dash-dotted (black) lines correspond to the density and number fluctuations as numerically calculated by using the Monte-Carlo code, respectively. The dashed (red), the solid-black (blue) and solid-grey (green) lines correspond to the density n j , and to the number fluctuations for the atoms in the first n (1) j and second layer n (2) j , as calculated with the EF model, respectively. The conditions for EF to be applicable, equations (6) and (7), are strictly valid for V o 5E R . Consistently, figures 1(a)-(c) show that for V o /E R = 11, 9, 7 the density profile and number fluctuations are very well reproduced by the model, except for the fluctuations of order 2 √ 2J/U in the flat region of the density profile. These small fluctuations are due to the particle-hole excitations which are neglected in the model, because they couple the two layers (of course, they can be added by hand, since they are well understood [19]). The model predicts that some sites at the trap centre have exactly a filling factor of two when 2(2J) ((N 2 − 1)/2) 2 . This condition is fulfilled for V o = 11E R , and in fact a flat density distribution with two atoms per site is observed in figure 1(a) at the trap centre, both in the analytical and numerical results. This confirms the validity of the idea of thinking of the atoms in the second layer as TG-bosons with effective hopping energy 2J.
For V o = 5E R , γ is barely 2γ c , and N 1 < J so that fluctuations at the edge of the first layer extend far enough to overlap with the ones of the second layer. We only expect EF to give qualitative predictions in this regime. For V o = 3, 2E R the conditions equations (6) and (7) are not satisfied, and the model fails to reproduce the exact results. Nevertheless, we notice that the EF model predicts the formation of a Mott state in the lowest layer for V o 3E R because γ ≈ γ c and (N 1 − 1) 2 /4 > 2J at 3E R . Accordingly, the numerically obtained fluctuations show the appearance of a flat region at the cloud's edge for V o = 3E R . Such flat region signals the formation of a Mott state as it evolves for deeper lattices into the observed dip in the fluctuations and disappears for shallower lattices ( figure 1(e)).

Centre of mass oscillations
In the NIST experiment [6], centre of mass dipole oscillations were induced by a sudden displacement of the harmonic potential by δ = 8 lattice sites. In a pure harmonic trap (without the lattice potential), the dipole mode is always stable and the cloud of atoms oscillates back and forth without altering its shape regardless of the magnitude of the initial displacement, the quantum statistics of the atoms and the strength of the interactions. On the other hand, theoretical studies done at mean field level using the so-called Gross-Pitaevskii equation [22]- [24] demonstrated that in the presence of an optical lattice, interacting bosons exhibit dynamical instabilities that may result in a very large damping of the dipole oscillations if the quasimomentum of the cloud of atoms is larger or of the order of πh/λ.
The initial displacement used in the NIST experiment could not accelerate the condensate to a quasi-momentum larger than about 0.1πh/λ and nevertheless a surprisingly large damping of the dipole oscillations was observed even for very shallow optical lattices. The damping increased dramatically with increasing lattice depth and at V o ≈ 3E R a sudden increase of about two orders of magnitude in the viscosity was observed (see figure 2). For V o 3E R , the gas became nearly immobile for times an order of magnitude longer than the single-particle tunneling time. We refer to this regime as the overdamped regime. The viscosity b was obtained in the experiment by fitting the centre of mass position to the formula m * ẍ = −bẋ − mω 2 T x, where m and m * are the atomic and effective masses and ω T is the magnetic trapping frequency.
Since these results were first presented a number of theoretical treatments have been put forward that, directly or indirectly, address various relevant aspects of the underlying dynamics, from different perspectives. It has been shown, for instance [25], that the momentum cutoff for the dynamical instability may be substantially lowered for commensurate lattices, and, probably more relevant for the experimental situation, that the boundary between regular and irregular motion becomes 'smeared out' due to quantum fluctuations. Numerical calculations based on a truncated Wigner representation [26] for shallow lattice depths have also shown that 9 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT quantum fluctuations enhanced by the lattice and the magnetic confinement generated a small, but non-negligible, atom population in the unstable velocity region of the corresponding classical dynamics causing damping of the centre of mass motion of the whole system. In [27] several of us have showed that time-dependent Hartree-Fock-Bogoliubov (HFB) calculations can describe the main features of the damping mechanism for V o < 2E R . There, a formula for the damping of the fluctuation-dissipation type was derived. This formula is based on a picture where the coherent motion of the condensed atoms is disrupted as these atoms try to flow through a random local potential created by the irregular motion of the noncondensed atoms.
In particular, the results of [26,27] provided results for the damping of the oscillations which quantitative agree with the experimental results in the shallow lattice regime. However the techniques used in these papers are not applicable to deeper lattices. Attempts to describe the overdamped dynamics in the strong correlated regime were presented in [9,28] by using standard fermionization techniques, valid for deep lattice, nevertheless these do not work to describe the NIST experiment as there the central lattice sites were populated with more than one atom. In this section we show that a time-dependent version of the EF model is capable to reproduce the dynamics in the strongly correlated regime without the use of adjustable parameters and in excellent quantitative agreement with the experimental results.
In the EF model, the centre of mass position of the atoms in the central tube is given by where Here, n = 1, 2; N is the number of atoms in the central tube, d is the lattice spacing and | (n) F (0) is the displaced ground-state of a system of N n non-interacting fermions with hopping energies nJ, respectively. In the experiment the measured viscosity was an average over all the tubes. The average can be calculated assuming that all tubes evolve independently and thus is the average centre of mass position. Here, P(N ) is the probability of having a tube with N atoms. Assuming an initial Thomas-Fermi distribution of the 3D system, P(N ) When there is at most one atom per site, atoms in the Mott state have been shown to be responsible for the overdamped dynamics [9,28]. Here, the two-layer model predicts that a Mott state is created in the lowest layer at V o ≈ 3E R , figure 1, and the formation of this Mott state explains the striking increase of the viscosity around V o 3E R (see also [3]). Atoms in the second layer are not required to be localized and in general their dynamics increases the mobility of the cloud. Nevertheless because there are always more atoms in the lowest layer than in the upper one the overall dynamics is overdamped above 3E R in accordance with the experimental observations (again, no adjustable parameters have been used in our theory). Notice that at 3E R the overall density profile is still an inverted-parabola (figure 1(e)) and does not exhibit any signature of a Mott phase. This is why the overdamped dynamics is unexpected in this regime. The EF model, based on the idea that phase transitions take place on a layer-by-layer basis, provides an explanation of the many-body dynamics.
In figure 2 the predictions of the EF model (red circles and blue triangles) for the viscosity coefficient are compared with the experimental data (black squares). Strictly speaking, the fluctuations in the two layers do not overlap during the dynamics only for V o > 7E R , thus we expect the model to be valid for the deepest lattice depths only. Nevertheless, because the dynamics is dominated by the overdamped behaviour of atoms in the lowest layer, figure 2 shows that the theoretical calculations are within the experimental error bars already for V o 4. In this case, dynamical properties are not strongly affected by the partial overlap between the layers. Figure 2 also shows that the viscosity for the central tube and the average are very similar and this is because tubes with about N = 80 atoms have the largest weight. The average shows a larger damping because it takes into account contributions from tubes that do not have extra atoms in the second layer. In the inset we also compare the experimental centre of mass position of the atomic cloud after 90 ms with the model solutions. The agreement between experiment and theory is consistent with the one found for the viscosity.

Summary
We have developed a simple model that generalizes the Bose-Fermi mapping to deal with spatially varying atomic densities in parameter regimes of direct relevance to current experiments. The model is relevant for 1D trapped gases where the co-existence of Mott-insulating regions with different occupation numbers is permitted. We presented the necessary conditions for the model to be valid and showed the usefulness and limitations of the method by comparing its predictions for some physical observables with numerical Monte Carlo simulations. Very good agreement between the model and the numerical solutions was found in the parameter regime where the model is valid. Finally we used the EF model to study the overdamped dynamics of the centre of mass after a sudden displacement of the trapping potential, and found good agreement with recent experiments. In particular, the overdamped motion was linked to the presence of a Mott state in the lowest atomic layer. In conclusion, the EF model provides intuitive insight and quantitative utility for the understanding of strongly correlated inhomogeneous 1D systems.