Density distribution of a trapped two-dimensional strongly interacting Fermi gas

We calculate and measure the density distribution and cloud size of a trapped two-dimensional $^{6}$Li Fermi gas near a Feshbach resonance at low temperatures. Density distributions and cloud sizes are calculated for a wide range of interaction parameters using a local density approximation (LDA) and a zero-temperature equation of state obtained from quantum Monte Carlo simulations reported by G. Bertaina and S. Giorgini, Phys. Rev. Lett. \textbf{106}, 110403 (2011). We find that LDA predictions agree well with experimental measurements across a Feshbach resonance. Theoretical results for Tan's contact parameter in a trapped gas are reported along with predictions for static structure factor at large momentum which could be measured in future Bragg spectroscopy experiments on two-dimensional Fermi gases.


I. INTRODUCTION
The striking ability to manipulate and control ultracold atomic 6 Li and 40 K Fermi gases has allowed the experimental investigation of strongly interacting twocomponent Fermi gases [1,2]. Using highly anisotropic pancake-shaped potentials to confine atoms in the lowest axial mode [3][4][5], it has become possible to realize experimentally two-dimensional (2D) Fermi systems. At low temperatures these can exhibit exotic properties such as Berezinskii-Kosterlitz-Thouless (BKT) [6][7][8] or inhomogeneous Fulde-Ferrell-Larkin-Ovchinikov (FFLO) superfluidity [9][10][11][12]. It has also been proposed that a 2D interacting Fermi gas may provide useful insights into hightemperature superconductivity [13,14] and itinerant ferromagnetism [15,16]. To date, a weakly interacting 2D Fermi gas has been imaged in situ [3] and used to characterize the crossover from two to three-dimensions [5]. The observation of 2D confinement induced resonances and measurement of the molecular binding energy using rf-spectroscopy of a strongly interacting Fermi gas was recently reported [4]. However, the thermodynamic properties of strongly interacting Fermi gases and fermionic mixtures [17] in reduced dimensions are yet to be fully explored.
On the theoretical side, numerous studies of 2D Fermi gases have been presented, addressing superfluid transitions [7,18], the effects of harmonic trapping [19], and population and mass imbalance [8,20,21]. Of particular importance, the equation of state of a 2D uniform Fermi gas was recently obtained through quantum Monte-Carlo (QMC) simulations [22] through the crossover from a Bardeen-Cooper-Schrieffer (BCS) superfluid of Cooper pairs to a Bose-Einstein condensate (BEC) of tightlybounded molecules [2]. Tan's universal many-body contact parameter was also found in 2D using the adiabatic relation [22,23]. While the most theoretical studies have relied on a perturbative or mean-field approach, the ab-initio QMC results at zero temperature [22] should provide a quantitative description of the many-body ground state of a strongly interacting 2D Fermi gas.
In this work, we measure the density distribution and cloud size of a two-dimensional trapped 6 Li Fermi gas at low temperatures in the BEC-BCS crossover and compare the data with theoretical predictions. The theoretical density distribution is calculated using a local density approximation (LDA) [24], based on the zero temperature QMC equation of state [22]. We find good qualitative agreement between experiment and theory in the strongly interacting regime near Feshbach resonance. When the atom number becomes large, we observe also substantial deviation from the 2D equation of state when new transverse vibrational modes are populated.
We also give theoretical predictions for Tan's contact parameter in a trapped gas and the static structure factor at large momentum, which could be measured in future Bragg spectroscopy experiments on 2D Fermi gases. Nontrivial pair correlations are reflected in the manybody part of the contact parameter. The many-body contribution to the contact exhibits a maximum near the 3D Feshbach resonance. In the deep BEC limit, however, the two-body contribution to the contact arising from the molecular state dominates. This paper is structured as follows: In the next section, we introduce the LDA and the QMC results for the zero-temperature equation of state for a uniform strongly interacting Fermi gas. In Sec. III, we discuss the production of a strongly interacting trapped 2D Fermi gas of 6 Li atoms and how to calculate the density distribution and cloud size from the QMC equation of state within LDA. The experimental procedure for the density measurements is briefly summarized. In Sec. IV, we present the theoretical density distributions and sizes and compare these with the experimental measurements. In Sec. V. we find Tan's contact and the large-momentum static structure factor of a 2D trapped Fermi gas. Conclusions and future perspectives are given in Sec. VI.

II. LDA AND 2D UNIFORM FERMI GAS IN THE BEC-BCS CROSSOVER
The equation of state of a strongly interacting Fermi gas in homogeneous space provides a convenient way to calculate the density distribution in a harmonic trap using the local density approximation (LDA) [24]. The basic idea is that for a sufficiently large number of particles in a slowly varying trapping potential V ext (r) we may treat the trapped Fermi system as a collection of many independent units that behave locally as a uniform Fermi gas. The correlation between different units, for example, the surface energy of each unit, is assumed to be negligibly small. Therefore, the local chemical potential of a unit at position r may be written as, where n (r) is the local density and µ is the chemical potential at the trap center. At zero temperature the local chemical potential µ [n (r)] of the locally uniform unit depends on the local density n (r) only. Hence, given the local equation of state µ[n(r)] at position r, we could solve inversely the density n (r) = n(µ − V ext (r)). The chemical potential at the trap center µ is set by the normalization condition´drn (r) = N , where N is the total number of particles. The LDA has been shown to work well in a wide range of situations [24]. It is valid for either non-interacting or strongly interacting Fermi gases in different geometries from 3D to 1D. The essential ingredient of the LDA is the local uniform equation of state µ(n). For a non-interacting twocomponent (spin-1/2) 2D Fermi gas at zero temperature, the chemical potential is simply the Fermi energy µ = E F = 2 k 2 F /(2m), where m is the mass of fermions and the 2D Fermi wave-vector is given by k F = (2πn) 1/2 . Therefore, the non-interacting chemical potential is proportional to the density, µ = π 2 n/m. The mean energy per particle is E/N = 2 k 2 F /(4m) = E F /2. For an interacting 2D Fermi gas in the BEC-BCS crossover, the energy per particle E/N has been calculated by Bertaina and Giorgini as a function of the interaction strength, by using the fixed-node diffusion Monte Carlo method [22]. Here, we extract the chemical potential from the QMC data of the energy per particle, since µ = ∂E/∂N . A suitable parameterization of the QMC data for E/N is therefore needed, as we discuss in detail below.
In the 2D BEC-BCS crossover, a peculiar feature of the contact interactions (between two fermions with unlike spins) is that any attraction, however small, will support a two-particle bound state with energy ǫ B = −4 2 /[exp (2γ) ma 2 2D ], where γ ≃ 0.577216 is the Euler's constant and a 2D is the 2D s-wave scattering length [14,25]. This is in sharp contrast with the 3D BEC-BCS crossover, where a two-body bound state appears on only one side of the Feshbach resonance where the 3D scattering length is positive [2]. The scattering length in 2D a 2D is always positive due to the existence of the bound state. The unitarity limit with an infinitely large scattering length (a 2D → +∞) is in fact trivial: it corresponds simply to the non-interacting (BCS) limit. In the opposite (BEC) limit with infinitely small scattering length (a 2D → 0 + ), where the energy of the bound state is infinitely large, two fermions are tightly bound to form a composite molecule. There will be a repulsive interaction between two composite molecules, characterized by an effective scattering length a d > 0.
The interaction strength in 2D may be expressed as a dimensionless interaction parameter η = ln(k F a 2D ). The weakly interacting BEC and BCS limits correspond to η → −∞ and η → +∞, respectively. The strongly interacting crossover regime occurs at about η = 0, where a 2D ∼ k −1 F . We interpolate the QMC data for the 2D equation of state (E/N − ǫ B /2)/E F G in the BCS-BEC crossover with a smooth, continuous analytical function f (η) that consists of three parts. On the BEC side in the range η < −1/2 we use the equation of state for a gas of composite molecules with a molecular scattering length a d [22]: where m d = 2m and N d = N/2 are respectively the mass and number of molecules, and n d = n/2 is their density. By assuming that a d /a 2D = α m ≈ 0.6 [22], Eq. (2) turns into (3) On the BCS side, for η > 2.72 we consider a Padé-type approximate function where the Padé coefficients a 1 = 0.164106, a 2 = 0.702385 and b 1 = 1.16411, b 2 = 2.40527 are obtained by minimizing the standard deviation between the QMC data and the values of the function (4) under the constraint of f BCS (η → +∞) → 1 − 1/η, which is required in a weakly interacting normal Fermi liquid in 2D. In the crossover regime between the BEC and BCS limits we use a sixth order polynomial function The coefficients are selected to provide the fit best to the the QMC data and to ensure continuity of the equation of state function and of its first and second derivatives at the two connection points η = −0.5 and η = 2.72. We find that, c . The blue line on the BCS side is the Padé-type approximation function of Eq. (4) that minimizes the standard deviation with respect to the QMC data [22]. The dot-dashed line in the crossover regime is a polynomial fit ensuring the continuity of the function and its first and second derivatives at the two connection points, see Eq. (5). The black circles show the QMC data of the equation of state. The chemical potential derived from the interpolating energy function is plotted by the thin black line.
As shown in Fig. 1, the interpolating function f (η) for the equation of state provides and excellent fit to the QMC data. By using µ = ∂E/∂N , we find that, The dimensionless chemical potential f µ (η) is shown in Fig. 1 by the thin black line.

III. A 2D STRONGLY INTERACTING FERMI GAS IN HARMONIC TRAPS
A. Achieving the 2D regime For an atomic Fermi gas in a harmonic trap, the 2D regime is achieved if the chemical potential µ and temperature T are sufficiently small compared to the excitation energy in one dimension (z). We consider a spin-1/2 Fermi gas of N 6 Li atoms with equal spin-populations in a highly oblate harmonic trap, where ω x ≈ ω y = ω ⊥ and ω z are the trapping frequencies in the radial (x, y) and axial (z) directions, respectively. The trap aspect ratio λ = ω z /ω ⊥ ≫ 1. The basic requirement for achieving a 2D trapped Fermi gas may be estimated by considering first the zero-temperature and non-interacting limits, in which the properties in a 2D harmonic trap, can be conveniently understood using LDA. As µ = π 2 n/m in an ideal 2D uniform Fermi gas, we obtain that . In other words, we expect a Thomas-Fermi (TF) distribution, where n T F is the TF peak density and ρ T F is the TF radius. Once the distance ρ is larger than the TF radius, the density is necessarily zero. The TF peak density and radius are related to the chemical potential by µ = π 2 n T F /m and µ = mω 2 ⊥ ρ 2 T F /2, respectively. Using the normalization condition´drn (r) = N , it is straightforward to obtain that, n T F = N 1/2 /(πa 2 ⊥ ) and ⊥ . As the lowest excitation energy in the z direction is ω z , one finds that the 2D regime can be reached if µ, E F < ω z and T < ω z /k B . The former condition requires that the total number of atoms N must be less than a 2D critical number, N 2D , equal to the number of single particle states with energy less than the lowest lying state with one transverse excitation. It is straightforward to show that N 2D = λ 2 . In our experiment with 6 Li atoms, the trapping frequencies are ω z ≃ 2π × 2800 Hz and ω ⊥ ≃ 2π × 47 Hz, leading to λ ≈ 60, a ⊥ ≈ 6.0 µm and N 2D ≈ 3600.

B. Achieving the strongly interacting regime
Experimentally, the strongly interacting regime is reached by tuning an external magnetic field B near a Feshbach resonance (B 0 = 834 G) for 6 Li atoms, for which the s-wave scattering length can be changed precisely to arbitrary values [26]. Here, a bg = −1405a 0 with a 0 ≈ 0.529 × 10 −10 m is the background scattering length, ∆B = 300 G is the width of Feshbach resonance and α = 0.0004 G −1 . In our highly oblate geometry, the tight-confinement in the z-direction induces a bound state in the 2D x − y plane. Therefore, one can express the 2D scattering length in terms of the 3D scattering length [1,28], where a z ≡ [ /(mω z )] 1/2 and b ≈ 0.915. In Fig. 2, we plot the dimensionless interaction parameter η = ln(k F a 2D ) as a function of the magnetic field at the 2D critical number of atoms, N = N 2D = 3600. As indicated by the dotted line, at the location of the 3D Feshbach resonance (B = B 0 ), the interaction parameter η ≈ 1. ⊥ . The 2D scattering length is calculated using Eq. (11) with az ≈ 770 nm for ωz = 2π × 2800 Hz.

C. Theoretical density distributions
Let us now consider the density distribution of an interacting 2D Fermi gas in the BEC-BCS crossover within LDA. The simple relation n(r) = mµ(r)/(π 2 ), useful for the ideal gas, is no longer applicable. We have to obtain numerically the local density from the local chemical potential by using the f µ -function, defined in Eq. (6). That is, we need to solve the following equation to find the density n(r) from, In the inversion procedure, our analytic interpolating f µfunction appears to be very convenient. The density distribution n(r) is calculated for an initially chosen value of chemical potential µ. We then adjust µ to satisfy the number equation´drn (r) = N . At the final stage, we quantify the cloud size using root mean square (rms) radius, In Fig. 3, we plot the theoretical density distributions at three values of the interaction parameter η, which correspond to the weakly interacting BEC and BCS sides, and the strongly interacting crossover regime. The density and distance from the center of the cloud are plotted in units of the TF density n T F and the TF radius ρ T F , respectively. One finds that with decreasing the interaction parameter η from the BCS side to the BEC side, the 2D cloud becomes denser and narrower in size, as anticipated. On the BEC side, our analytic interpolating f µ -function leads to the following asymptotic behavior, where α m = a d /a 2D ≈ 0.6 is the ratio between 2D molecular and atomic scattering length. Hence, with decreasing η (→ −∞) the peak density increases as [−4η] 1/2 and the radius of the cloud decreases as [−4η] 1/4 . On the other hand, on the BCS side the density distribution converges to the ideal Fermi gas result of Eq. (9), n (ρ) In accord with these asymptotic density distributions, the cloud sizes are given by, and in the BEC and BCS limits, respectively. Here, ρ 2 IG = ρ T F / √ 3 is the rms cloud size of an ideal 2D Fermi gas. We have checked that these analytic results agree well with our numerical calculations in the appropriate limits. We have also checked the sensitivity of the calculated density profiles to the form of the equation of state. The widths of the theoretical density profiles do not vary by more than 5% which is small on the scale of the width changes as the interaction strength is varied over the range considered here.

D. Experimental measurements
To measure the density distribution of a strongly interacting 2D Fermi gas, we use an experimental setup similar to the one used in our previous work [5]. In brief, we start with a cloud of approximately 1 × 10 5 6 Li atoms in the two hyperfine states |F = 1/2, m F = ±1/2 in a far detuned 3D optical dipole trap. The cloud is evaporatively cooled to the lowest possible temperature near the Feshbach resonance. At this stage, the number of atoms is controlled by further lowering the trap depth to spill atoms out of the dipole trap. We then ramp on a 2D optical trap in 200 ms to create a highly oblate trap with trapping frequencies ω z ≃ 2π×2800 Hz and ω ⊥ ≃ 2π×47 Hz, which gives an aspect ratio of approximately 60. Finally, we tune the interaction strength by adiabatically ramping the external magnetic field to 810 G, 834 G and 992 G, where the cloud is held and imaged. The dimensionless interaction parameters ln(k F a 2D ) at these fields are about −0.5, +0.6, and +5, respectively.
The critical number of atoms for reaching 2D regime is N 2D ≈ 3600. Depending on the depth of the dipole trap, the final number of atoms in the cloud can be varied in the range of 500 to 5000 atoms. The final temperature of these small 2D and quasi-2D clouds is difficult to determine when interactions are present. However, we anticipate it to be approximately 0.1 T F based on applying the same preparation procedure to clouds with a larger atom number.
Before the imaging, we allow a short time of flight (500 µs). This allows us to resolve the density distribution in the tightly-confined z-direction, since this time scale is long compared to 1/ω z ≈ 57 µs. It is however much short compared to 1/ω ⊥ ≈ 3.4 ms, and therefore, the cloud distribution in the radial direction is essentially equivalent to the in situ profile. The imaging beam propagates roughly along the radial x-direction, which means that there is an automatic integration over the x-direction for the total density distribution. We then integrate these distributions over the z-direction to generate a double-integrated column densityñ(y). The rms cloud size is calculated from the first moment of the onedimensional profile ρ 2 =´dyñ(y)y 2 /´dyñ(y). Theoretically, we perform the same integration procedure in the x-direction for the 2D density distribution. This should lead to the same distribution as the experimentally measured double-integrated column densityñ(y).

IV. COMPARISON BETWEEN THEORY AND EXPERIMENT
In Fig. 4, we compare the LDA column densityñ(y) (lines) with the experimental measurements (solid circles) at three magnetic fields. To reduce the experimental noise, we average the experimental density distributions over many images with a range of atom numbers that are all well below N 2D . Accordingly, the theoretical lines correspond to the average number of atoms, while the standard deviation in the atom number is illustrated by shaded region. We observe good qualitative agreement between theory and experiment with no adjustable parameters. The distributions become wider with increasing ln(k F a 2D ) however, there are notable discrepancies between theory and experiment, particularly in the wings of the clouds. This is due to a number of effects including the finite imaging resolution and recoil induced blurring during the imaging pulse. The combination of these two artifacts is to lower the effective resolution of the imaging system to approximately 6 µm. This alone however, is not enough to fully account for the observed discrepancies. The remaining differences are most likely due to the finite temperature of the clouds which will show up most in the wings of the distribution.
In Fig. 5, we plot the rms cloud size of the 2D Fermi gas in units of the TF radius as a function of the interaction parameter η = ln k F a 2D . The lines and symbols show, the LDA predictions and experimental data, respectively. In the strongly interacting regime at the magnetic fields B = 810 G and 834 G, the LDA predictions agree quantitatively well with the experimental data (solid red circles and green squares). At the high field B = 992 G, the gas is more weakly interacting and the experimental data (blue diamonds) lie slightly above the theory curve. In the inset of Fig. 5, we present the actual cloud size as a function of the number of atoms. The agreement between theory and experiment near the Feshbach resonance (810 G and 834 G) becomes more apparent.
In all three sets of measurements, the radial cloud size drops below the true 2D prediction for the largest atom numbers. This happens when the first transverse excited  state becomes energetically accessible and leads to a drop in the growth rate of the radial cloud size. In this quasi-2D regime, shell structure associated with resolving the discrete transverse states can dramatically affect the density profiles [5]. The atom numbers for which the radial width departs from the 2D theory are slightly below what we would expect for an ideal Fermi gas. This could be due to interactions or scattering resonances [27][28][29][30][31] but also may arise through thermal excitations. For a strongly interacting Fermi gas with contact interactions, the asymptotic behavior of various physical quantities in the limit of short-distance or largemomentum is governed by a single parameter, called the contact, which measures the density of fermionic pairs within a short distance. This was first discussed by Tan in 2008, when he derived a set of exact universal relations for strongly interacting Fermi gases [23]. Being an important many-body parameter, Tan's contact is also related to the thermodynamics via the adiabatic relation [23]. The contact in a homogeneous 2D Fermi gas was calculated by Bertaina and Giorgini [22]. Here we use these results to find the contact and static structure factor in a trapped system. In 2D, the adiabatic relation takes the form [32], where the derivative is taken at constant entropy. The calculation of Tan's contact parameter for a 3D strongly interacting Fermi gas was performed recently [33], by using the similar adiabatic relation.
The zero-temperature contact of a homogeneous 2D Fermi gas can be calculated by substituting our interpolated energy per particle E/N = ǫ B /2 + E F G f (ln[k F a 2D ]) into Eq. (19). We find that I = I 2b + I mb , where is the contribution from the two-body bound state and is the contribution from the many-body correlations, respectively. The two-body contact I 2b increases monotonically from the BCS to the BEC limit. In contrast, the many-body contact I mb should exhibit a maximum at the crossover regime, according to the behavior of the energy f -function (see Fig. 1). As shown in Fig. 6 by a thin line, the maximum of the many-body contact occurs at η ∼ 0.8, which is roughly the position of the Feshbach resonance (see also, ref. [22]).  For a trapped interacting Fermi gas, we may calculate the contact by using the LDA density distribution, where we have summed the local contact density, I(r)/∆V = [I/(N k 2 F )] (r) × n(r)k 2 F (r), over the whole space. It is easy to see that, the two-body contact is not affected by the density average, so that I T,2b = I 2b .
However, the many-body part may be significantly affected. In Fig. 6, we present the result for the contact of a trapped 2D Fermi gas at the BEC-BCS crossover, with the many-body contact shown by a red dashed line. It has roughly the same shape as the many-body contact of a homogeneous gas, with a peak appearing at η ∼ 0.8. However, the peak is about half as high due to the average over the density distribution.
B. Spin-antiparallel static structure factor Tan's contact for a 3D strongly interacting Fermi gas has been measured in a number of ways. One appealing method is to measure the spin-antiparallel static structure factor by using Bragg spectroscopy at large momentum [34], which has a 1/q tail with a prefactor given by Tan's contact [35]. In 2D, we may make a similar prediction. It has been shown that the 2D pair correlation function n (2) (r) ∝ I ln 2 (ρ/a 2D ) [32]. The ln 2 (ρ/a 2D ) dependence can be qualitatively understood from the twobody relative wave-function ψ rel (r) ∼ ln(r/a 2D ), since n (2) (r) ∝ |ψ rel (r)| 2 . The spin-antiparallel static structure factor is simply the Fourier transform of the pair correlation function. Thus, we find that, whereq ≡ q/k F and γ = 0.577216. Compared with the 3D case, the spin-antiparallel static structure factor in 2D decays faster with increasing momentum q (q −2 compared to q −1 ). In Fig. 7 we plot the theoretical prediction for the spinantiparallel static structure factor of a trapped 2D Fermi gas in the BEC-BCS crossover for momenta q = 3k F and q = 5k F . We split the structure factor into the two-body and many-body parts, in accord with the previous classification for the contact. Near the Feshbach resonance, the static structure factor at q = 3k F is about 0.02, whose magnitude is accessible within current experimental resolution [34].

VI. CONCLUSIONS
To conclude, we have predicted theoretically and measured experimentally the density distribution and cloud size of a low-temperature two-dimensional harmonically trapped Fermi gas in the BEC-BCS crossover. The theoretical calculations have been carried out within a local density approximation, based on the ab-initio zerotemperature equation of state obtained from the fixednode diffusion Monte Carlo simulations [22]. The experimental measurements were performed using a single twodimensional Fermi cloud of 6 Li atoms near a Feshbach resonance. We have found good qualitative agreement between theory and experiment.
We have also calculated an important many-body parameter, Tan's contact, and have proposed that it can be straightforwardly measured using Bragg spectroscopy for the spin-antiparallel static structure factor as in three dimensions. In future studies, it will be interesting to study both experimentally and theoretically the density distributions at finite temperatures, which may elucidate the fermionic Berezinskii-Kosterlitz-Thouless transition in two dimensions.