BCS pairing in a trapped dipolar Fermi gas

We present a detailed study of the BCS pairing transition in a trapped polarized dipolar Fermi gas. In the case of a shallow nearly spherical trap, we find a decrease in the transition temperature as a function of the trap aspect ratio and predict the existence of the optimal trap geometry. The latter corresponds to the highest critical temperature of the BCS transition for a given number of particles. We also derive the phase diagram for an ultracold trapped dipolar Fermi gas in the situation where the trap frequencies can be of the order of the critical temperature of the BCS transition in the homogeneous case, and determine the critical value of the dipole–dipole interaction energy below which the BCS transition ceases to exist. The critical dipole strength is obtained as a function of the trap aspect ratio. Alternatively, for a given dipole strength, there is a critical value of the trap anisotropy for the BCS state to appear. The order parameter calculated at criticality exhibits novel non-monotonic behaviour resulting from the combined effect of the confining potential and the anisotropic character of the interparticle dipole–dipole interaction.


Introduction
One of the most challenging goals of modern atomic, molecular and optical physics is to observe the superfluid (BCS) transition [1] in a trapped Fermi gas. The possibility of such a transition for gases with attractive short-range interactions has been predicted in [2] and, in connection with trapped gases, in [3], and has been a subject of very intensive experimental investigations since then (see [4] for the latest experimental results). In typical experiments, evaporative cooling is used to cool fermions. However, since the Pauli principle forbids the s-wave scattering for fermions in the same internal state, Fermi-Fermi [5] or Fermi-Bose [6] mixtures have to be used to ensure collisional thermalization of the gas. Such a combination of evaporation and sympathetic cooling allows temperatures of T 0.1T F to be reached, where T F is the Fermi temperature at which the gas exhibits quantum degeneracy. Unfortunately, critical temperatures for the BCS transition, T c , are predicted to be much smaller than T F . Nowadays, the standard method employs a Feshbach resonance in order to increase the atomic scattering length a s to larger negative values. Such a 'resonance superfluidity' should lead to superfluid transition at T c ≈ 0.1T F [7]. In the most promising scenario, one starts with a molecular condensate formed for a s 0 and modifies a s towards the negative values [4,10]. Another way to achieve the BCS regime is to use a cooling scheme that can overcome the Pauli blocking, such as appropriately designed laser cooling [8]. Yet another promising pathway is to place the Fermi gas in an optical lattice and enter the 'high T c ' regime [9].
The temperature of BCS transition in a two-component Fermi gas depends dramatically on the difference in concentration of the two components, which presents another experimental obstacle [10]. This problem, however, is not relevant for a polarized Fermi gas with long-range interactions, such as dipole-dipole ones. The possibility of the Cooper pairing has been predicted in [11] and the critical temperature (including many-body corrections), as well as the order parameter, have been obtained for a homogeneous gas in [12]. Possible realizations of dipolar gases include ultracold heteronuclear molecular gases [13], atomic gases in a strong dc electric field [14] and atomic gases with laser-induced dipoles [15] or with magnetic dipoles [16]. For dipolar moments d of the order of one Debye and densities n of 10 12 cm −3 , T c should be in the 100 nK range, i.e. experimentally feasible.
Dipole-dipole interaction is not only long-range, but also anisotropic, i.e. partially attractive and partially repulsive. Thus, the nature of the interaction for trapped gases may be controlled by the geometry of the trap. For a dipolar Bose gas in a cylindrical trap with axial (radial) frequency denoted by ω z (ω ρ ), there exists a critical aspect ratio λ = (ω z /ω ρ ) 1/2 above which the Bosecondensed gas collapses if the atom number is too large [17] and below which the condensate exhibits roton-maxon instability [18]. The trap geometry is also expected to control the physics of trapped dipolar Fermi gases. So far, however, only partial results have been obtained [19]: analytic corrections to T c in 'loose' traps and solution of the case of an infinite 'slab' with ω ρ = 0 and with ω z finite. In the latter case, there exists a critical frequency above which the superfluid phase does not exist. Very recently, we have reported results for the case of a general trap [20], and predicted the critical dipole strength below which the BCS pairing transition ceases to exist.
In this paper, we present a detailed derivation of the results of [19,20] and study the effect of trap geometry on the BCS transition in trapped dipolar Fermi gases. We first consider the case of a shallow, nearly spherical trap with trap frequencies ω z , ω ρ much smaller than the critical temperature T c of a spatially homogeneous gas with the density equal to the maximal density of a trapped gas sample, ω z , ω ρ T c . In this case, the presence of a confining potential results in a decrease of the critical temperature when compared with the spatially homogeneous gas. In the case of a strong confining potential where ω z , ω ρ ∼ T c , we calculate the phase diagram in the plane − λ −1 , where = 36nd 2 /πµ is the strength of the dipole-dipole interaction relative to the chemical potential µ. Below the critical value of the interaction, < c , the BCS transition does not take place. Similarly, for a given dipole interaction strength, there is a critical value of λ −1 above which the BCS state appears. We determine the dependence of c on λ −1 and calculate the order parameter at the criticality. The order parameter exhibits a novel non-monotonic behaviour in strongly elongated cylindrical traps.

BCS pairing in a dipolar Fermi gas
We consider a dipolar Fermi gas polarized along the z-direction in a cylindrical trap. The corresponding Hamiltonian reads where ψ † (r) and ψ(r) are the canonical fermionic creation and annihilation operators, respectively, of particles with mass m and dipolar moment d; The BCS pairing is characterized by the order parameter (r 1 , , which attains nonzero values below the critical temperature T c . Just below T c , the order parameter is a non-trivial solution of the BCS gap equation (see, e.g., [1]) and is an extremum of the functional The kernel K in the above expression is where is the Matsubara Green function of the Fermi gas in the normal phase. In this formula, ω n = πT(2n + 1), n = 0, ±1, ±2, . . . , ψ ν (r) and ε ν are the eigenfunctions and the corresponding eigenvalues of the free particle Hamiltonian in the trap [−h 2 ∇ 2 /2m + V trap (r)]ψ ν (r) = ε ν ψ ν (r).
The solution of the general problem given by equations (2)-(4) is very difficult even for numerical methods due to the large number of variables, relatively low symmetry of the system and a long-range character of the interparticle interaction. Remarkably, in the case of a shallow, nearly spherical trap, the solution can be found analytically, whereas the general case can be treated by using a variational approach.

Shallow nearly spherical trap
We begin with the case of a weakly deformed spherical trap with frequencies much lower than the critical temperature: ω ρ , ω z T c . In the new variables the second term in equation (3) reads where The kernel K depends on variables R and R only due to the presence of the trapping potential; but, as a function of r and r , the kernel decays rapidly for r, r > ξ 0 , where 2mµ being the Fermi momentum) determines the characteristic scale for pairing correlations. Under the condition ω ρ , ω z T c , the gap (r 1 , r 2 ) = (R 12 , r 12 ) is a slowly varying function (on the scale ξ 0 ) of R 12 = (r 1 + r 2 )/2 (see the end of this section). At the same time, the Fourier transform of with respect to r 12 = r 1 − r 2 , varies on a scale of the order of the Fermi momentum, p ∼ p F ∼hn 1/3 , see [12]. It is therefore convenient to write equation (5) in the form where R c = (R + R )/2 = (R 12 + R 34 )/2 and r c = R − R and, keeping in mind that the pairing takes place in the central region of the gas sample, where U trap (R) µ, we can expand the order parameter in powers of (r + r )/4.
The leading term of this expansion is with G ω n (R, P) = r G ω n (R, r) exp(−iPr).
In the case ω ρ , ω z T c , the Green function G ω n (R, P) can be approximated as the integration over r c in equation (9) gives The term containing U trap (∇ q /2) can be neglected because ∇ q ∼ 1/p F , and, therefore, equation (9) can finally be written in the form This expression corresponds to the local density approximation with µ → µ(R c ) = µ − U trap (R c ) expanded in powers of U trap (R c )/µ up to the first order. The next term of the expansion of equation (8) in powers of (r + r )/4 is quadratic (the linear in (r + r )/4 contribution vanishes due to the symmetry of the problem) and has the form where ← − ∇ i and − → ∇ i are the ith component of the gradient ∇ R c acting on the left (on * ) and on the right (on ), respectively. After neglecting the q dependence of , the integrations over r c and q are straightforward and equation (14) can written as One can show with the help of the explicit form of the Green functions, equation (10), that the main contribution from the integrals over r and r is with n being the unit vector in the direction of P, n = P/P, and ξ P (R c ) = P 2 /2m − µ(R c ). This expression decays exponentially for |ξ P (R c )| > T c , and, therefore, can be approximated by the simpler one in integrals over P with a slow varying functions of P where ζ(z) is the Riemann zeta-function. With this simplification, equation (15) becomes where ν F (µ(R c )) = mp F (R c )/2π 2 is the density of states on the local Fermi surface.
After combining together equations (13) and (16) and performing the variation with respect to * (R c , P), we obtain the gap equation in the form The first line of equation (17) coincides with the gap equation in the spatially homogeneous case, see [12], with R c being a parameter. Therefore, following the method developed in [12], we find that the order parameter on the local Fermi surface has the form (R c , [12] for more details) and the function (R c ) obeys the equation where C = 0.5772 is the Euler constant, λ 0 (R) = (mp F (R)/2π 2 ) d , with d being the dipoledipole scattering amplitude ( d = −8d 2 /π in the Born approximation). With an account of the expressions for the critical temperature T (0) c and explicit form of ϕ 0 (n) in the spatially homogeneous case, see [12], equation (18) for (R c ) takes the final form where Note that, in obtaining this equation from equation (19), we also expand λ 0 (R c ) up to the first order in U trap (R c )/µ, similar to equation (13). Equation (19) is equivalent to the Schrödinger equation for a three-dimensional (3D) anisotropic harmonic oscillator. As a result, the shift in the critical temperature due to the presence of the trapping potential is given by the lowest eigenvalue and equals In the considered case ω i /T c 1, the critical temperature in the trap is only slightly lower than that in the spatially homogeneous case.
Just below T c , the order parameter has the Gaussian form Thomas-Fermi radius of the gas sample in the ith direction. This justifies our assumption that pairing takes place only in the central part of the sample. On the other hand, we have l i ξ 0 and, therefore, the gradient expansion of the order parameter in powers of r + r is legitimate.
We see that there exists the optimal value of the trap aspect ratio λ * = 0.81, which corresponds to the highest transition temperature in the trap. The existence of the optimal value for the trap aspect ratio is a result of the competition between the anisotropic character of the dipole-dipole interparticle interaction and the confinement of the gas sample in the radial direction. The former becomes predominantly attractive for smaller values of λ (cylindrical traps) and, therefore, favours the BCS pairing. The latter, in contrast, acts destructively on the pairing due to the size effect and, therefore, is less pronounced for larger values of λ.

Critical aspect ratio
We have seen in section 3 that the trap geometry has a more pronounced influence on the BCS pairing in a dipolar Fermi gas compared with a two-component Fermi gas with a short-range interaction. This is due to the fact that the states, which form Cooper pairs in a trapped dipolar Fermi gas, have different quantum numbers n z . Therefore, their energies are different, at least by the amount of ω z . When this difference becomes of the order of T c , pairing is obviously impossible. As a result, superfluid transition in a trapped dipolar Fermi gas can take place only in traps with ω z < ω zc , where the critical frequency ω zc is found to be ω zc = 1.8 T c [19]. For ω z < ω zc , as can be seen from equation (21), confinement to the radial direction decreases the critical temperature as well. Therefore, one would expect the existence of a critical aspect ratio λ c such that pairing is possible only in traps with λ < λ c [20].
In this section, we study the BCS pairing in the case of a cylindrical trap where the trap frequencies can be of the order of the critical temperature in the spatially homogeneous gas, ω ρ , ω z ∼ T (0) c , but still much lower than the chemical potential, ω ρ , ω z µ. In this case, the BCS gap equation (2) can hardly be tractable even numerically without additional simplifications.
As it was shown in [12], the BCS pairing is dominated by the p-wave scattering with zero projection of the angular momentum on the z-axis, l z = 0, in which the interaction is attractive. Contributions of higher angular momenta, although present due to the long-range character of the dipole-dipole interaction, are numerically small (see also [11]). In the p-wave channels with l z = ±1, the interaction is repulsive and, therefore, leads only to renormalizations of a Fermiliquid type, and will be neglected here. Therefore, for the considered pairing problem, we can model the dipole-dipole interaction by the following (off-shell) scattering amplitude where p is the incoming momentum, p the outgoing one and γ 1 (E) some function of the energy E. The amplitude 1 describes anisotropic scattering only in the p-wave channel with the projection of the angular momentum l z = 0. The function γ 1 (E) obeys the equation that follows from the Lipmann-Schwinger equation for the off-shell scattering amplitude [21]; γ 1 (E) is assumed to be negative in order to guarantee the BCS pairing. The cut-off parameter ensures convergence of the integral and, in fact, can be expressed in terms of the observable scattering data corresponding to the on-shell scattering amplitude with p = p and E = p 2 /m. It follows from equation (23) that γ 1 (E) is inversely proportional to E, γ 1 (E) = γ 1 (2mE) −1 , with some coefficient γ 1 . Therefore, the on-shell amplitude is energy-independent, as it should be for low-energy scattering on the dipole-dipole potential (see [22]).
The coefficient γ 1 determines the value of the critical temperature T c of the BCS transition in a spatially homogeneous gas and can be expressed through the dipole moment d using the results of [12]. In a homogeneous gas, the order parameter has the form (p) = p z 0 with some constant 0 , and the linearized gap equation implies where ξ p = p 2 /2m − µ, and the bare interaction is renormalized in terms of the scattering amplitude with γ 1 (µ) = γ 1 /p 2 F at the Fermi energy ε F = µ = p 2 F /2m along the lines of [2] (p F is the Fermi momentum). After integrating over p, we obtain the equation on T c : where ν F = mp F /2π 2 is the density of states at the Fermi energy. After comparing the solution of equation (25) with the result of [12] for T c , we find γ 1 = −24d 2 /π. In the ordinary space, the scattering amplitude 1 is where r and r are the relative distances between the two incoming and the two outgoing particles, respectively. Therefore, the order parameter in the trapped gas, (r 1 , r 2 ) ∼ ψ(r 1 )ψ(r 2 ) , has the form where r = r 1 − r 2 , R = (r 1 + r 2 )/2, and the corresponding equation for the critical temperature is Here ξ i = ξ(n i ), n = (n x , n z , n z ), is the harmonic oscillator quantum number, ξ(n) =h[ω z (n z + 1/2) + ω ρ (n x + n y + 1)] − µ, and the function M n 1 n 2 (R) is defined as with ϕ n (z) being the harmonic oscillator wavefunction.
The gap equation (28) is still hardly tractable numerically and, hence, we employ additional approximations. We assume a large number of particles such that the chemical potential µ is much larger than the trap frequencies, µ ω z , ω ρ . Therefore, while calculating the functions M n 1 n 2 (R), we can use the WKB approximation for the wave functions ϕ n of the most important for the BCS pairing states with energies near the Fermi energy ε F = µ. Another simplification is due to the fact that the BCS order parameter 0 (R) varies slowly on an interparticle distance scale n −1/3 ∼h/p F , where p F = √ 2mµ is now the Fermi momentum in the centre of the trap in the Thomas-Fermi approximation. As a result, the pairing correlations are pronounced only between states that are close in energy. This allows us to neglect q 2 /4m in the denominator of the second term in equation (28) together with rapidly oscillating terms in the functions M n 1 n 2 (R). We then obtain (see the appendix) and where n ≡ |n 1 − n 2 | N ≡ (n 1 + n 2 )/2, l iN = √ 2Nh/nω i = l 0i √ 2N, whereas U n (z) = sin [(n + 1) arccos z]/ sin(arccos z) and T n (x) = cos(n arccos x) are the Chebyshev polynomials. The functions M (z) n 1 n 2 (z) and M (ρ) n 1 n 2 (x) fulfil the following completeness relations: with δ 0 = 1 and δ n>0 = 2, which follow from the completeness of the Chebyshev polynomials. It is convenient to rewrite equation (28) in the following way: because the sum ξ 1 + ξ 2 does not depend on n 1 − n 2 and, therefore, with the help of formulae (33) and (34), the summation over n can easily be performed in kernels K 2 and K 3 . On the other hand, the kernel K 1 is determined entirely by the states near the chemical potential µ.
The calculation of the sums in the kernel K 3 (R, R ) gives where µ(R) = µ − i mω 2 i R 2 i /2, and we have replaced the discrete sums over N i with the integrals over continuous variables p i = √ 2mω i N i . As a result, the kernel K 3 (R, R ) is Following equation (23), we see that the kernel K 3 results in the replacement µ → µ(R) in the scattering amplitude γ 1 (µ) on the left-hand side of the gap equation (35). We are interested in the critical value λ c of the aspect ratio, above which the BCS pairing does not take place for a given strength of the dipole interaction. Therefore, this value corresponds to vanishing critical temperature. We therefore calculate the kernels K 1 and K 2 in the limit T ω i . The summation over ξ(N) = (ξ 1 ) + ξ 1 )/2 in the kernel K 2 is then within the limits −µ ξ(N) −ω z /2, where −ω z /2 is the upper limit due to the fact that the function M (z) n 1 n 2 is non-zero only when |n 1z − n 2z | 1. These sums can again be replaced by integrals with the following result: where p F (R) = √ 2mµ(R) is the local Fermi momentum in the Thomas-Fermi approximation. In order to calculate the kernel K 1 in the limt T ω i , we note that non-zero contributions to the sums over n 1 and n 2 originate only from the region ω z /2 |ξ(N)| ε(n)/2, where n i = |n 1i − n 2i | and ε(n) =h i ω i n 1 . As a result, the kernel K 1 (R, R ) can be written as where the functions V are defined as An extremum of this form at T c corresponds to the eigenvector of the matrix M m z ,m ρ ,n z ,n ρ with a zero eigenvalue.
The condition that the matrix M has a zero eigenvalue imposes a constraint on the interaction parameter and the trap frequencies ω i and, therefore, allows to determine the critical aspect ratio. Indeed, the parameter can be written as = 36n(0)d 2 /πµ, where n(0) = (2mµ) 3/2 /6π 2h3 is the gas density in the centre of the trap and, hence, for a given dipole moment d, depends only on the chemical potential µ. On the other hand, the chemical potential µ and the total number of particles in the trap N are also related, 3N = µ 3 /ω z ω 2 ρ . As a result, for a given and N, the product of the trap frequencies ω z ω 2 ρ is completely determined, and the only free parameter is the trap aspect ratio λ = (ω z /ω ρ ) 1/2 . We may thus determine its critical value λ c from the condition that the lowest eigenvalue is zero.
The calculation of the matrix elements M m z ,m ρ ,n z ,n ρ is naturally divided into two parts (see equation (39)). The contribution with the local kernel L(r) is a 3D integral that can easily be computed using, for instance, the Monte Carlo integration routine, such as the VEGAS algorithm from the GSL library [23]. The part with the non-local kernel K(r, r ) is a triple sum, whose elements are 8D integrals. The straightforward approach based on the same numerical method fails in this case because it is too time-consuming. To overcome this problem, we calculate the functions V n z ,n ρ n (α, β, ζ) in the following way. We integrate numerically over r for fixed α, β and n values (the value of ζ is then fixed by the δ-function) and use these data for a 2D spline interpolation of the integrand for the last integrations over α and β. In this way, the time required to compute the matrix elements M m z ,m ρ ,n z ,n ρ reduces to about 72 h on a standard workstation. This procedure gives convergent results with respect to the highest powers of the polynomials M ζ and M ρ in our ansatz, equation (42). We also checked that our solutions do not depend on the number of shots in the Monte Carlo algorithm and on the number of points chosen for interpolation. They are also a proper asymptotic behaviour for r → 0.
The results of our calculations are presented in two figures. Figure 2 shows the desired relation between the interaction strength and the inverse critical aspect ratio λ −1 c . The three curves correspond to three different numbers of particles. As it could be expected, for a larger number of particles, the critical aspect ratio λ c is larger because the interaction is stronger. For a pancake trap, λ −1 < 1, the interaction is predominantly repulsive, and higher values of for a fixed λ are required to achieve the BCS transition. On the other hand, for a cigar trap, λ −1 > 1 , the interaction is predominantly attractive and the BCS transition takes place at smaller values of . The existence of the critical interaction strength (for a given value of λ) is a result of the discreteness of the spectrum in the trap and of the specific structure of the order parameter (∼p z ). The latter allows pairing only between particles in the states, where quantum numbers n z differ by an odd number (intershell pairing). Therefore the pairing correlations have to be strong enough to overcome the corresponding energy difference. The solid line shows 0 (z, ρ = 0) and the broken line corresponds to 0 (z = 0, ρ). Figure 3 shows the order parameter 0 (r) for the cigar trap with λ −1 = 2.2. The order parameter is a non-monotonic function of the distance from the trap centre, in contrast to the case of the BCS order parameter in a two-component Fermi gas with short-range interactions [24] (see, however, [25], in which an oscillating and highly non-monotonic order parameter was obtained for the case of an intershell pairing). This effect persists, despite being less pronounced, for the case of a pancake geometry. In the axial direction, the order parameter (z, ρ = 0) develops a minimum at ρ < 1, whereas, in the radial direction, (z = 0, ρ) becomes negative in the outer part of the cloud. This behaviour, originating from the anisotropy of the interparticle interaction, can have profound consequences for the form and spectrum of the elementary excitations.

Conclusions
We have presented a detailed theory and analysed the influence of a trapping potential on the BCS pairing in a dipolar single-component Fermi gas. We have determined the phase diagram for trapped dipolar Fermi gases in the − λ −1 plane, where measures the dipole strength and λ is the trap aspect ratio. BCS transition at a finite temperature T is possible iff > c (λ). We have calculated the critical line c (λ) and the order parameter at criticality.
where we neglect the rapidly oscillating contribution cos [ N+n/2 (x) + N−n/2 (x) − π]. In the case n N, the phase difference N+n/2 (x) − N−n/2 (x) can be simplified as Note also that one obtains the same result if, instead of differentiating the WKB wavefunction equation (A.1), one uses the well-known relations between the wave functions of a harmonic oscillator and its derivatives.