Bogoliubov-de Gennes study of trapped spin-imbalanced unitary Fermi gases

It is quite common that several different phases exist simultaneously in a system of trapped quantum gases of ultra-cold atoms. One such example is the strongly-interacting Fermi gas with two imbalanced spin species, which has received a great amount of attention due to the possible presence of exotic superfluid phases. By employing novel numerical techniques and algorithms, we self-consistently solve the Bogoliubov de-Gennes equations, which describe Fermi superfluids in the mean-field framework. From this study, we investigate the novel phases of spin-imbalanced Fermi gases and examine the validity of the local density approximation (LDA), which is often invoked in the extraction of bulk properties from experimental measurements within trapped systems. We show how the validity of the LDA is affected by the trapping geometry, number of atoms and spin imbalance.


Introduction
Interest in spin-imbalanced Fermi superfluids dates back over a half century to the appearance of the Bardeen-Cooper-Schrieffer (BCS) theory. Several pioneers theoretically considered the fate of a superconductor in the presence of magnetic impurities which could disrupt the spin balance [1,2,3,4] prescribed by BCS. In the ensuing decades the BCS picture has emerged as an important paradigm in many branches of physics [5], such as nuclear physics, quantum chromodynamics, and ultracold atoms. The ultra-cold atoms, in particular, provide a highly controllable clean system to investigate the effect of spin imbalance on the nature of fermionic pairing. Over the past few years, experiments on polarized Fermi gases at Rice University [6,7,8], MIT [9,10,11,12] and ENS [13] have stimulated a flurry of theoretical activity. These efforts represent an important area within the general goal of emulating complicated many-body systems using cold atoms.
All cold atom experiments are necessarily performed in the presence of trapping potentials that hold the atoms together in an inhomogeneous environment. In order to extract the bulk properties of the system (e.g, the equations of state) from measurement on a trapped sample, a local density approximation (LDA) is employed to account for the effects of the trapping potential. The LDA states that the system can be treated locally as a part of an infinite system. In other words, on a length scale which is short in comparison to the spatial variation of the external potential, the external potential V ext (r) essentially acts as an offset to the chemical potential. In practice, one defines a local chemical potential µ(r) = µ − V ext (r), with µ being the global chemical potential. In the context of cold atoms, the LDA is often an accurate approximation when there are no phase boundaries present. However, because of the large change in the density which occurs from the center to the edge of the trap, there is often more than one phase present and care needs to be taken in application of the LDA. This is especially true at phase boundaries where the LDA usually fails. Even in the circumstances where there is only one phase present in the trap, corrections to the LDA become crucial when the number of particles in the sample is small and finite-size effects are significant [14].
The trapped polarized Fermi gas represents just such a system -in all experimental observations, phase separation between a region with vanishing spin population and one with finite spin population has been observed. The purpose of this paper is to examine the validity of the LDA as the particle number and aspect ratio of the trap is changed. To this end we employ the Bogoliubov-de Gennes (BdG) equations [15], which is a powerful mean-field tool particularly suitable for inhomogeneous Fermi superfluids, and has been recently adopted by many to study trapped ultra-cold Fermi gases [16,17,18,19,20].
In a recent paper [21] we analyzed discrepancies observed in pioneering experiments on polarized fermionic superfluids [9, 10, 6, 7] which appear to center on physics arising when the containing trap is highly elongated. In our analysis we found that as the trap becomes increasingly elongated, the solutions of the BdG show a tendency towards metastable behavior which could lead to the observation of states which are not necessarily the ground state. However we also observed that one class of solutions with structure similar to the LDA solution was consistently the lowest in energy within our analysis. Our conclusions have since been confirmed by similar calculations [22] using a density functional formulation [23] which is more sophisticated than the BdG and accounts for quantum fluctuations and interactions within the normal fluid. In this paper we focus on this class of solutions and examine how well such solutions match those obtained from the LDA approximation as a function of particle number and trap aspect ratio.
The content of the paper is organized as follows. In Sec. 2, we present the BdG formulation and describe the numerical techniques used to solve the BdG equations. In Sec. 3 we discuss our results for a relatively small number of particles and for different trapping geometry. In Sec. 4, we focus on a very elongated cigar-like trap but vary the number of atoms. Finally, a concluding remark is presented in Sec. 5.

Formulation
We consider a gas of spin-polarized fermionic atoms interacting through a contact potential (g i<j δ 3 ( r i − r j )) and confined to a harmonic trap defined in cylindrical coordinates (ρ, φ, z) by V ext (ρ, z) = m 2 (ω 2 ⊥ ρ 2 + ω 2 z z 2 ) with axial and radial frequencies denoted by (ω z , ω ⊥ ). We work at unitarity (a s → ∞) and within a cigar-shaped trap with aspect ratio α = ω ⊥ /ω z . This system of N = N ↑ + N ↓ atoms interacting through contact interaction is described by a HamiltonianĤ = d R (H 0 + H I ) with non-interacting H 0 and interacting H I portions given by: where ψ σ ( R) represents the fermionic field operators, m the mass and µ σ the chemical potential of atomic species with spin σ. The coupling constant is defined as g = 4π 2 a m . Henceforth, we work in trap units for which: m = ω z = = k B = 1. This implies that energies will be measured in units of ω z , lengths in units of l 0 = mωz and temperature (T ) in units of ω z /k B .
The Hamiltonian (1) will be treated within the mean-field BdG approximation for which there are many excellent references [15,24,25]. Here we simply state the BdG equations for the pair wave functions u j ( R) and v j ( R) which decoupleĤ: In the above coupled set of equations, u j ( R) and v j ( R) are two components of the quasiparticle wavefunction associated with energy E j . The single particle Hamiltonian H s σ is defined in our trap units by: and includes the trapping potential, the chemical potential µ σ and the Hartree meanfield potential is given by density ρ σ ( R) = Ψ † σ ( R)Ψ σ ( R) . In accordance with fermionic commutation relations [15], the quasi-particle amplitudes are normalized as: and are related to the spin densities through : where f (E) = 1/(e E/k B T + 1) is the Fermi-Dirac distribution function. In the unitary limit the Hartree terms gρσ on the diagonal of Eq. (2) do not actually diverge but are unitarity limited [26]. How to incorporating the Hartree term in the unitarity limit is beyond the mean-field BdG formalism. For this reason we ignore this term in our calculations. The paring field or gap paramter is give by Since the Hartree terms are ignored in our analysis, Eqs. (2), (4) and (6) constitute a closed set of nonlinear equations which we solve self-consitently. However as presented above, our formulation has one problem which arises from a nasty side effect of the contact interaction. The contact interaction assumes wrongly that all states are scattered in the same way regardless of their incoming energy and consequently sums in contributions from collisions at arbitrarily high energy which creates an ultra-violet divergence. Hence the gap equation 6 needs to be properly regularized as we now discuss.

Regularizing the BdG equations
Due to the assumption of contact interaction, the gap ∆ is a function of the center-ofmass coordinate of the pair, R. To discuss the regularization, it is more convenient to re-introduce back the relative coordinate r, with which the gap is defined as which diverges as ∆( R) 2πr when r → 0 [24]. To regularize Eq. (6), one simply subtracts off the 1/r divergence to obtain the regularized equation [24]: In practice, the convergene of the sum above is quite slow and we discuss here a numerically efficient way of evaluating ∆( R) to sufficient accuracy without undue effort. First an energy cutoff E c is used to break the sum of Eq. (6) into two pieces as ff.
The second term ∆ c ( R) is an approximation to the sum above the cutoff using the LDA result for the pairing field [25] which can also be written as : where k c is the momentum cutoff related to E c . This leads to a computationaly efficient form of the gap equation: Here we have employed the identity: to subsume the LDA approximation of the gap (∆ c ) into an effective interaction defined by: Below this cutoff E c , the quasiparticle states are calculated exactly by solving Eqs. (2), (10) and (12) self-consistently along with the normalization conditions: which conserve total particle number N = N ↑ + N ↓ and overall polarization P = (N ↑ − N ↓ )/N. The iterative solution of these equations is achieved using a modified Broyden's approach [27] which is a nonlinear mixing scheme allowing the formation of polarized regions even if they were not present in the initial condition. In this scheme convergence was achieved when the root mean squared difference between ∆ at different iterations was below some tolerance i.e., where j is the position index and i represents the iteration number. We should point out that at unitarity this method is analogous to a descent technique where the step is optimized through the residuals stored from a few previous steps. As mentioned previouly we choose as our initial condition ∆ LDA , the LDA solution to the BdG equations. This is an important point becasue we found in previous work [21] that at larger particle numbers the solution of the BdG can be quite sensitive to the initial conditions.

Special features
We discretize using a linear triangular finite element mesh in the ρ-z plane which anticipates that our results will retain the cylindrical symmetry of the confining potential V ext . The accuracy of these calculations are controlled by the density of the trianglular mesh and the cut-off E c used in the hybrid scheme. Both of these are changed in successive solutions until the free-energy or relevant observable converges to a sufficient accuracy. Experience has taught us that this simple renormalization scheme typically converges when the cutoff is of the order 6E F (where E F is the Fermi energy) which implies that the number of quasiparticle states to be directly calculated by Eq. (2) is about 6N. Note that this puts a constraint on the density of the discretizing mesh. Thus, for moderate system sizes, we are still presented with a very large problem. For example, for N ∼ 10 3 particles, one essentially needs to calculate ∼ 10 5 quasiparticle states at each iteration.
One important consequence of our finite element discretization is that it yields sparse matrices which are suitable to massively parallel matrix computations. This is of key importance given that the slow convergence of the sum in Eq. (10) condemns us to calculate a very large number of quasi-particle states. This is true inspite of our efficient hybrid scheme, without which calculation would be prohibitive. It is immediately obvious that these difficulties will increase with the number of particles N, and will make the problem impractical for even moderate particle numbers without very careful formulation. In our case these difficulties are inescapable since the issues to be addressed occur in the presence of finite size effects and confinement. Hence it was crucial to develop the ability to perform calculations with realistic particle numbers because it is not a priori obvious how physical properties will scale with system size. At each iteration, we need to find a large number of eigenvalues and eigenfunctions for large matrices. To this end, we use a novel shift-and-invert scheme which we developed independently but is very similar to the one described in Ref. [28].
Briefly, the scheme involves partitioning the sought spectrum amongst groups of processors working independently. The size of the group is determined by the minimum number of processors with enough total memory to store the inverted matrix which is required for building the local Krylov basis. The main challenge here is bookkeeping to prevent over-counting of states and balancing the load amongst the processor groups. It is conceivable that this method could have issues in cases where the equations support huge degenerate subspaces. In our particular formulation we exploited cylindrical symmetry and parity along the long trap axis to reduce the problem. Consequently we only had to contend with accidental degeneracies. A good analysis of these issues may be found in Ref. [28] but a more thorough description of the numerical details will be presented elsewhere. However we note here that this parallelization scheme is very efficient on distriubted computing systems and scales easily to thousands of CPU's which is as high as we have tested. Potentially it can be used to study much larger systems than we have reported here or in [21].

Validity of the BdG
We should devote a few lines here to comment on the validity of the BdG theory at unitarity becasue it is expected that quantum fluctuations and other effects due to the strong interactions could be significant in this regime. The main drawbacks of the BdG is that it fails to account for phase fluctations. At unitarity it has an additional disadvantage in that it also fails to account for interactions within the normal fluid which is unitarily limited [26]. However the BdG is widely expected to yield qualitatively reliable answers for two main reasons. First because of the finite size of these experiments, the trapped gas enjoys protection from fluctuations of arbitrarily low energy or of very long wavelengths. Secondly due to experimental evidence for superfluidity at unitarity, it is quite clear that interactions within the normal fluid are not so great that the order parameter cannot form or will be destroyed. Thus the failure to account for these effects is not expected to change the topology of the phase diagram but at most would slightly shift the phase boundaries. Since our purpose to examine the suitability of the LDA is qualitative, we are confident that the BdG can account for the essential physics. Nevertheless, due to the limitations within the BdG formalism, in particular the neglect of interaction in the normal fluid, our calculation fails to quantitatively locate position of the Clogston limit. In addition, it cannot be applied to study a system with extremely large polarization (i.e., P ≈ 1) where the polaron physics will dominate [13,29].

Results for N = 200
In this section, we focus on a relatively small particle number of N = 200. As we will show, the system is rather sensitive to the trap geometry. In the following, we first briefly discuss the case of a spherical trap with α = 1 and then concentrate on elongated cigar-like trapping potentials with α > 1 and then discuss them in detail.

Spherical trap
Liu et al. studied a spin-imbalanced Fermi gas confined in a spherical harmonic trap in Ref. [16]. To benchmark our work, we first did a series of calculations for this geometry and found our results in perfect agreement with those reported in Ref. [16]. In this case, even though we anticipate only cylindrical geometry, the density profiles always obey the spherical symmetry. Note that the authors of Ref. [16] solved the one-dimensional (1D) radial equation, hence the spherical symmetry of the cloud is automatically imposed. We refer the readers to Ref. [16] for details; here we give just a brief description of the key features. The density profiles indicate a phase-separation scenario: a fully paired BCS superfluid core at the trap center surrounded by a fully polarized shell composed of excess majority spins. A thin layer of partially polarized gas forms the interface between the superfluid core and the normal shell. In this intermediate regime, the minority density and the order parameter sharply drop to zero. Here and in other cases, we always found that the profile of the order parameter closely follows that of the density of the minority spin component. Furthermore, in this case, the LDA gives very good agreement with the full BdG calculation even for particle numbers as small as a few hundred.

Cigar trap
So far, all the experiments on spin-imbalanced Fermi gases have been performed in cigar-like traps with aspect ratio α > 1. For a given particle number, the 1D regime is eventually encountered in this geometry by increasing α and creates the possibility to study the 3D-1D dimensional crossover. Figure 1 illustrates several examples of the density profiles for N = 200 atoms confined in a moderately elongated trap with α = 5 (this trap aspect ratio is close to what has been used in the MIT experiments). We find it convenient to express our results in terms of the Fermi energy E F = (3N) 1/3 α 2/3 , central number density (2E F ) 3/2 /(6π 2 ), and the Thomas-Fermi radius along the z-axis Z F = 2E F for a single species ideal Fermi gas of N/2 particles in a trap with identical parameters. The upper row of Fig. 1 shows the density profiles of a system with a relatively small polarization P = 0.2. Here the axial spin density n 1D (z) exhibits a double-horn structure and vanishes near z = 0. This is a clear violation of the LDA which predicts that n 1D should be flat topped [30]. Fig. 1 can be examined in tandem with Fig. 2 where for a closer inspection, we plot the densities and the order parameter along the axial and radial axis for two different polarizations. Fig. 2(a) displays results for P = 0.2. The density profiles along the z-axis show clearly a phase separated threeregion structure -moving from the center to the edge of the trap, we encounter a fully paired superfluid core, a partially paired intermediate region and a fully polarized normal gas, just like in the previous case of spherical trap. In stark contrast, the density profiles for the two components along the ρ-axis are completely overlapped. In fact, this matching of the radial profiles occur for |z| ≤ 0.1. As a consequence, the axial spin density vanishes near z = 0 as shown in the upper row of Fig. 1.
That the majority and minority densities overlap along the radial direction can be understood from an argument invoking the surface energy. When induced phase separation occurs, there is an accompanying surface energy associated with the interface between the two phases. The system will then try to minimize the interface in order to reduce the associated energy. For a cigar-like trap as we study here, the superfluidnormal gas interface area can be efficiently reduced if the two spin components match their densities radially. The authors of Ref. [31,32] devised phenomenological theories to include the surface term variationally to explain the breakdown of the LDA observed in the Rice experiment [6,7]. In our calculation, the surface energy is automatically included from the self-consistent BdG formulation [33].
As polarization increases, eventually it becomes energetically unfavorable to have this radial overlap. This is illustrated in Fig. 2(b) for P = 0.6. Consequently, the axial spin density no longer vanishes near z = 0 and the LDA becomes more accurate (see the middle and bottom rows of Fig. 1). In addition, it is quite noticeable that, particularly for large P , the minority component density has a steeper down turn along the axial axis than along the radial axis. Moreover, in the partially polarized intermediate region, the order parameter has a small oscillation along the axial axis, but not along the radial axis. Similar order parameter oscillations were also found in the spherical trap case [16]. This is a consequence of the proximity effect which, in the context of superconductor, occurs when a superconductor is in contact with a normal metal, the Cooper pairs from the superconductor diffuse into the normal component.
Next, we keep N fixed at 200 but increase the trap aspect ratio to α = 50 which represents a much more elongated cigar trap and close to what is used in the Rice experiment. A similar display of the column and axial spin density profiles for different polarization as in Fig. 1 is shown in Fig. 3. In this very elongated trap, the majority and minority components have their densities matched along the radial axis up to the highest polarization we have calculated which is P = 0.7, and the minority component has a boxy-looking density profile. This further confirms that the system is able to greatly reduce the effective surface area between the normal state and the superfluid state in anisotropic cigar-like traps. Another marked feature for such an elongated trap is the prominent oscillations of the order parameter along the z-axis.
As demonstrated in Fig. 4, these oscillations are quite generic features in such a trap with finite P . As P increases, both the amplitude and the spatial extension of the oscillations increase. As shown in Fig. 4(b), at large polarizations, the axial length of the partially polarized intermediate region becomes comparable to or even larger than that of the BCS core. Accompanied by the oscillation in the order parameter, the density profiles (in particular, the minority density) also exhibit strong oscillations. Such oscillations are reminiscent of the FFLO pairing state predicted by Fulde, Ferrel, Larkin and Ovchinnikov [3,4] in which the Cooper pairs possess finite momentum and the order parameter in the bulk develops sinusoidal oscillations that break the spatial translation symmetry.
Our calculations also show that these axial oscillations are aligned along the radial axis, as shown in Fig. 4(c). We even intentionally started from an initial ansatz of ∆ where the axial oscillations are present but with the nodes mis-aligned in the radial direction, the BdG iterations eventually converge to a state where the nodes are perfectly aligned radially. This radial alignment has important impact in detecting the oscillations in column density profiles where the densities are integrated along one radial axis: Due to the radial alignment, the oscillations are not washed out and can be easily observed, for example, in the doubly integrated axial spin density n 1D , as illustrated in Fig. 4(d).
It is interesting to compare our result with the recent work by Bulgac and Forbes [23] who, using a density functional theory, argued that the FFLO pairing phase occupies a larger phase space region than people previously thought for a 3D homogeneous system. The FFLO state found in Ref. [23] is also associated with large-amplitude density oscillations, particularly in the minority component. Another perspective on the order parameter oscillations and its potential connection with the FFLO phase is dimensionality. It is well known that the partially polarized phase with FFLO-like oscillations is prominently featured in the phase space of 1D systems [25,34,35]. That the reduced dimensionality favors such a state can be understood from the argument of Fermi surface nesting or alternatively from the cost of creating domain walls. The cigarlike traps used in our calculation mimic a quasi-1D system and may be the reason that we see pronounced oscillations in our calculation. If this latter explanation is correct, i.e., the partially polarized region featuring FFLO-like oscillations is due to the effective reduction of the spatial dimension, we then expect to see these oscillations diminish as N is increased while the trap aspect ratio is fixed, which makes the system more 3D-like. To confirm this, we now turn to the next section where we keep α = 50 but vary the total particle number N.

Results for large particle numbers at α = 50
We now consider a trapped system in a very elongated cigar trap with α = 50. Fig. 5 shows the density and order parameter along the axial and radial axis at a fixed polarization P = 0.3 but different values of total particle number N. As one can clearly see, the oscillations in both the order parameter and the density profiles diminish as N is increased and the LDA approximation becomes more and more accurate, which indicates that the FFLO-like region observed above for small N does not represent a bulk 3D phase. Rather, it is a finite-size effect due to the effective reduction of the spatial dimension.
Nevertheless, we have discovered that as N increases, the system exhibits a tendency towards metastability [21]. Numerically, by starting from different initial anstze for the order parameter ∆( r), the BdG solution may converge to different final states. Our calculations show that among these different states, the one that closely resembles the LDA solution always has the lowest energy as long as N is sufficiently large (N 10 4 ), but that there may exist several metastable states with energies just slightly larger that violate the LDA. The observed LDA-violating states at Rice are most likely these metastable states. Experimentally, whether the ground state or a metastable state will be realized may depend upon how the evaporative cooling procedure is implemented [36]. This has been confirmed very recently in a new experiment by Hulet group [37].

Conclusion
In conclusion, we have carried out a systematic study of a trapped spin-imbalanced Fermi gas in the unitary limit up to a total number N ∼ 10 5 atoms. We study a class of solutions which has recently been identified as having the lowest energy [21,22] through a self-consistently solution of the Bogoliubov-de Gennes equations using stateof-the-art numerical techniques. For a given set of trapping parameters, the LDA will eventually become accurate for sufficiently large N. However, the validity of the LDA is also sensitive to the trap geometry: Traps with small anisotropy favor the LDA. Our calculations show that for a relatively small number of atoms in a very elongated trapping potential, the system contains three phases: an unpolarized BCS phase, a partially polarized FFLO-like phase and a normal phase. That the FFLO region exists may be understood from the view point of reduced effective spatial dimension. As N is increased while all other parameters remain fixed, the FFLO-like region eventually disappears. A detailed analysis of dimensional crossover will be reported in the future.