Coexistence of pairing gaps in three-component Fermi gases

We study a three-component superfluid Fermi gas in a spherically symmetric harmonic trap using the Bogoliubov-deGennes method. We predict a coexistence phase in which two pairing field order parameters are simultaneously nonzero, in stark contrast to studies performed for trapped gases using local density approximation. We also discuss the role of atom number conservation in the context of a homogeneous system.

Ultracold Fermi gases have opened up a way to explore multi-component gases experimentally. Recently, a degenerate three-component gas was successfully created [4,5]. The stability of three-component gases is hindered by the three-body recombination, reducing the pairing in the gas and the lifetime of the sample [6][7][8]. However, there are ways to stabilize the gas against such losses using, for example, optical lattices. Optical lattices are interesting also due to the rich phase diagram: theoretical investigations have found that color superconductivity competes with normal phase and formation of trions [9][10][11][12][13][14][15].
The hyperfine energy spacing in alkaline atoms stabilizes the atom numbers, making the atom number of each component separately conserved. This is the case for most two-component Fermi gas experiments and a natural extension of these studies is a non-SU(3) symmetric case in which the atom numbers in all three components are fixed. Here we study such a threecomponent system with fixed atom numbers in all three components in a spherically symmetric harmonic trap using the Bogoliubov-deGennes (BdG) equations. We study the coexistence of the pairing gaps in these systems and discuss the scaling of the system size up to the thermodynamical limit.
In section 2, we give an overview of the three-component system and the corresponding mean-field theory. In section 3, we consider the BCS-type mean-field theory in homogeneous space and describe the effect of boundary conditions on the stability of different phases. In section 4, we consider the effects of trapping potential using the BdG method. In section 5, we show the main results obtained from the BdG method, especially regarding the coexistence of two pairing gaps. We conclude with a discussion in section 6.

The system setup
The general mean-field Hamiltonian for a three-component system in the contact interaction potential approximation is (up to a constant) where the first term of the Hamiltonian includes contributions from the kinetic energy, the external trapping potential V σ (r) (that can depend on the component |σ ), and the chemical potentials µ σ , respectively. Interactions are described by the two-body scattering T-matrix. In the contact interaction potential approximation, it can be written as 3 where a σ σ is the scattering length between atoms in hyperfine states |σ and |σ , and m r = 2m σ m σ /(m σ + m σ ) is twice the reduced mass. The Hartree fields are denoted by W σ (r) = σ =σ U σ σ (r) n σ (r) and the densities are n σ (r) = † σ (r) σ (r) . The pairing (mean-)field σ σ (r) =Ũ σ σ (r) σ (r) σ (r) includes a renormalized interactionŨ σ σ (r) that is used to remove the ultraviolet divergence following the standard procedure (see below). In our model, we neglect the possibility of three-body bound states and other three-body effects that can affect the lifetime of the gas [7].
A three-component system has three possible pairing fields corresponding to the three interaction channels U 12 , U 13 , U 23 , and these can be combined into a total pairing field vector = [ 23 , − 13 , 12 ] T . Identical properties make the system SU(3)-symmetric and the pairing vector can be reduced to a single gap by a simple unitary transformation. The orientation of the pairing vector corresponds to a choice of the global gauge [16][17][18][19], [26], and the simplest choice is the one where only one of the pairing gaps, say 12 , is non-zero. This, of course, makes atoms in component |3 effectively non-interacting. Indeed, in an SU(3)-symmetric case, there always exists a gapless branch describing unpaired atoms.
In our study, one interaction is always suppressed (we choose U 13 = 0). This is the case, for example, in 6 Li where Feshbach resonances between the three lowest hyperfine states |1 − |2 (at B = 834 G) and |2 − |3 (at B = 811 G) lie close to each other, while the resonance |1 − |3 (the nearest one at B = 690 G) is sufficiently far away [28]. Similar behavior occurs in 40 K [29], where the richness of the hyperfine level structure allows even more freedom in choosing the suitable interaction strengths. Moreover, mixtures of 6 Li and 40 K offer interesting possibilities [30,31]. Also, we do not consider here interactions between |1 and |3 induced by the component |2 [32]. Thus, neglecting the |1 − |3 interaction channel altogether, the symmetry is broken at least to SU(2) × SU(1). Analogous to the SU(3) symmetric case, the total pairing field is now a two-dimensional (2D) vector = [ 12 , 23 ] T , and in the symmetric case, where components |1 and |3 are identical, it is preserved under spin rotations of the hyperfine states |1 and |3 . This symmetry implies that the ground state is degenerate with respect to the orientation of the total pairing field vector. However, the degeneracy is lifted by changing masses, chemical potentials or interaction strengths and, as we will soon show, also by imposing boundary conditions, such as fixing the number of atoms in different components.
Boundary conditions, such as fixed particle numbers or fixed chemical potentials, manifest themselves in different ways in atomic gases. While the total particle numbers are, in practice, fixed, the local densities are not, as the particles are allowed to move around in the trap. Hence, from the local density approximation point of view, locally the relevant boundary condition appears to be a fixed chemical potential. However, globally the relevant boundary condition is the fixed particle number, and in the BdG method we indeed fix the mean particle number. Below we will also discuss how the two pictures merge in the limit of large system size N → ∞.

Homogeneous system
For a homogeneous system, the densities and gaps lose their spatial dependence. This corresponds to the usual BCS approximation in which pairing can only occur between atoms with opposite momenta |k, σ and | − k, σ . The mean-field Hamiltonian can be written in matrix form as where C is a constant and the single-particle dispersion is ξ σ k =¯h 2 k 2 2m σ − µ σ . We have here neglected the Hartree fields since at the level of our approximation they provide only a constant energy shift in a homogeneous system. The standard approach calls for diagonalizing this using the Bogoliubov transformation and iteratively solving for the pairing fields 12 and 23 . To satisfy the fixed mean atom number boundary condition, the iteration must adjust the chemical potentials in a self-consistent manner as well. Note that the inherent atom number fluctuations implied by the mean-field theory play no role here as long as the typical fluctuations (scaling as √ N ) are much smaller than the atom numbers in different species (scaling as N ).
As discussed above, in the symmetric case (µ 1 = µ 3 , m 1 = m 3 and U 12 = U 23 ), the Hamiltonian has SU(2)×SU(1) symmetry and all of the pairing fields with 2 12 + 2 23 constant yield the same total energy. The ground state is thus degenerate. However, different orientations of the pairing field vector yield different atom numbers in components |1 and |3 . Thus, fixing the numbers of atoms N 1 and N 3 breaks the degeneracy and a well-defined energy minimum is found. Figure 1 shows typical energy landscapes H MF + σ µ σ N σ as a function of the pairing fields 12 and 23 for equal interaction strengths U 12 = U 23 . The Fermi momentum k F here and throughout this work is defined as the Fermi momentum of component |2 , F is the Fermi energy of component |σ . The energies have been calculated for fixed atom numbers: the chemical potentials µ σ are solved for every point ( 12 , 23 ) so that the atom number constraints are satisfied. The figures show clearly how the ground state becomes non-degenerate and realizes itself in a particular combination of pairing fields. In the special case where there is an equal number of atoms in all three components, N 1 = N 2 = N 3 , the ground state is still degenerate. However, if the number of atoms in component |2 is changed (keeping N 1 = N 3 but N 2 = N 1 ), the degeneracy is broken and a non-degenerate energy minimum appears. This is in stark contrast to the case where the chemical potentials are kept constant and the atom numbers are allowed to vary. In such a case the ground state remains degenerate as long as the chemical potentials for components |1 and |3 are equal, µ 1 = µ 3 .
In the case of a number mismatch or difference in the interaction strengths of components |1 and |3 , the energy minimum will be shifted from the equal pairing case. Depending on the number of atoms in component |2 , the minimum appears either at the edge of the energy landscape (yielding either of the two pairing fields 12 or 23 zero) or somewhere in between. This too is an important difference from the case of fixed chemical potentials where the breaking of the symmetry (by either changing chemical potentials or interaction strengths) always results in either of the two pairing fields dominating and the other becoming zero. Thus, for fixed chemical potentials, one does not observe coexistence of the two pairing fields 12 and 23 except possibly in the symmetric, or degenerate, case, whereas for fixed atom numbers the coexistence phase (described by two non-vanishing pairing field order parameters 12 and 23 ) is very real.
The pairing scheme is revealed by the momentum distribution of each state. Figures 2 and 3 show the momentum distributions of the three components for equal pairing gaps ( 12 = 23 ) and for projected pairing gaps (obtained by setting 23 = 0, which is always an allowed solution, and minimizing the energy by varying only 12 ). In the first case, equal pairing gaps imply that there are equal numbers of 12 and 23 pairs. Since only zero-momentum Cooper pairs are considered here, one can filter the paired atoms from the momentum distributions and determine the momentum distribution of unpaired atoms by calculating the difference n σ k − n 2k /2 for σ = 1, 3. The distribution of unpaired atoms is seen to form a clear Fermi sphere, but with a maximal occupation probability of 0.5. In the case of a projected pairing gap in figure 3, the pairing atoms |1 and |2 form a breach due to a number mismatch between the two components.
The boundary condition of fixed atom numbers for each component separately is natural for atomic gas experiments. However, the results for a homogeneous gas must be approached with caution since experiments are always conducted in non-uniform trapping potentials. Using these homogeneous system results in conjunction with local density approximation means locally fixing the chemical potentials instead of the atom numbers. This discrepancy concerning which boundary condition to use can be solved by treating the trapping effects explicitly using the BdG method.

Harmonic trap-the Bogoliubov-deGennes (BdG) method
In order to consider trapped systems, we use the BdG method, which allows the inclusion of trap effects exactly. The mean-field BdG method is not expected to be able to capture all relevant  physics in the strongly interacting regime. However, in an unbalanced two-component system, it has been shown [33] that for small polarizations and symmetric trap geometries, there is good agreement between the mean-field BdG approach and the real-space dynamical meanfield theory. We solve the three-component mean-field system in a spherically harmonic trap 7 V σ (r) = 1 2 m σ ω 2 σ r 2 using the eigenbasis of the 3D harmonic oscillator, where Y lm are the spherical harmonics and the radial wavefunctions are given by Here, L l+1/2 n (r 2 σ ) is the associated Laguerre polynomial andr σ ≡ r √ m σ ω σ /h. The mean-field Hamiltonian separates for different l-quantum numbers H = l H l + C, such that [H , H l ] = 0 and C is constant. Introducing a finite cutoff energy E c and keeping only single-particle states with energy less than the cutoff allows the writing of each H l using block matrices, where J l 2 = J l 21 + J l 23 . The block matrices are defined as with the Hartree shift, and the pairing field, F l σ σ nn = ∞ 0 drr 2 R l σ n (r ) σ σ (r )R l σ n (r ).
The connection between the interaction strength U H σ σ used in the Hartree shift and the bare interaction strength U σ σ will be discussed below. The energy matrix l σ is diagonal with elements σ nl =hω σ (2n + l + 3/2) − µ σ and the operator vectors are c σ l = [c σ 0 l0 · · · c σ N l l0 ] T . We denote the number of single-particle states with fixed l, whose energy is below the cutoff, Similarly to free space, the H l matrices can be diagonalized using the Bogoliubov transformation, which is provided by unitary 3N l c × 3N l c -matrices W l . By inserting the identity operator (W l ) † W l between the matrix and the operator vectors, we have the quasiparticle basis as (γ 1 l γ † 2 l γ 3 l ) T = W l (c 1 l c † 2 l c 3 l ) T . The rotation matrix W l is chosen such that the matrix in (7) is diagonalized.
The equations for the pairing fields and the densities are where n F is the Fermi distribution andσ = σ − 1, σ ∈ {1, 2, 3}. The total number of particles in each component |σ is obtained by integration N σ = ∞ 0 dr r 2 n σ (r ). As in usual BCS theory, the gap equation is ultraviolet divergent, hence the energy cutoff is E c . In order to make the model cutoff independent, we follow a standard approach [34] and use a renormalized interactionŨ σ σ but now generalized to a three-component system, 1 where Here the momentum cutoff k c,σ σ and the local Fermi momentum k F,σ σ are defined as where the local average chemical potential is We solve these gap and number equations self-consistently using fixed point iteration. For every iteration step in the gap equation, we solve the chemical potentials µ σ to keep the particle numbers constant. The iteration is terminated when the subsequent gap profiles in the iteration differ by at most 5 × 10 −5h ω. The cutoff energy is chosen to be 2.5 × max{E σ F } (with a higher cutoff, the results do not qualitatively change). If the convergence is slow, we try different initial values to ensure that the final result is correct. We use a small finite temperature (T = 10 −3 T F ) to smooth the Fermi distribution and to help solve the number equations in the presence of a discrete energy spectrum. However, we have checked that the results are unchanged even for zero temperature; for example, in figures 6 and 7. We do not consider higher temperatures in this work. However, we have checked that our results are sufficiently robust to survive low but experimentally relevant temperature T = 0.05 T F . An example of the effect of the temperature is shown in figure 4, where exactly the same parameters as these in figure 6 are used except that T = 0.05 T F . The results, especially the coexistence region, are practically identical; only the minor features at the edge of the gas have been smoothed.
The Hartree fields become infinite with a diverging scattering length a σ σ → ∞. This unphysical effect is caused by improper treatment of two-body scattering effects, and in practice these energy shifts are limited by the Fermi energy. Monte Carlo results on a two-component Fermi gas suggest that the Hartree fields at unitarity do not exceed |W | ≈ 0.5 E F [35,36]. We limit the Hartree field interaction to be smaller than this value by imposing a hard cutoff on the Hartree interaction strength. That is, instead of the bare interaction U σ σ we use U H σ σ which is limited from above by Note that component |2 experiences two Hartree fields due to the two components, |1 and |3 , and thus the total Hartree shift experienced by this component can be up to double the above cutoff. Since the Hartree shifts induce a mismatch between the Fermi surfaces, the pairing amplitude is reduced. Note that this is in contrast to the balanced two-component case where both components experience the same Hartree potential and the densities remain thus equal. In three-component systems, the inclusion of the third interaction U 13 would reduce the mismatch, but because it is usually weaker, the mismatch does not totally disappear. However, the mismatch can be countered by careful choice of interaction strengths and atom numbers, so that local density imbalances are reduced. The Hartree shift can be seen to produce interesting shell structures for certain parameters. Such exotic shell structures created by the Hartree shift will be considered elsewhere, requiring a more complete treatment of the Hartree effect to confirm their validity. In the present work, we have chosen to focus on a parameter range in which these peculiar features are not present, and in all of the results shown in this work except in figure 5, we neglect the Hartree effects; that is, we use U H σ σ = 0 for all σ, σ . Our choice does not significantly limit the parameter range, because these effects appeared always only in tiny islands in the parameter space, typically at the edges of the trap. Indeed, we have checked that the inclusion of the Hartree shift in the way described above does not qualitatively change the results presented here. An example is shown in figure 5, which presents the case of figure 6 but with Hartree fields included. The qualitative behavior is the same, although the numerical values of the order parameters are smaller. Importantly, the Hartree fields do not affect the coexistence of the two-order parameters.  inside 23 is relatively constant when increasing the system size (increasing atom numbers N ), the characteristic length scale given by the oscillator length r osc = √h /m 2 ω 2 . Locally, the dominating pairing channel is the one for which the atom densities are least mismatched, the strength of the interaction being only a secondary factor. This local nature of pairing allows interesting shell structures [37], as shown in figure 7. However, we do not pursue these issues here but rather concentrate on more general features. Note that the number of atoms in component |3 has been chosen to be slightly smaller than that in components |1 and |2 , but the features shown here are very general.    (18) Figure 8(c) shows this parameter as a function of interaction strength and position, revealing a large coexistence region in the somewhat narrow interaction strength window −0.51 < (k F a 23 ) −1 < −0.47, but also a coexistence region close to the edge of the trap across a wide range of interactions. In figure 8(d), we show how the coexistence parameter at the center of the trap scales with increasing system size N (the atom numbers N 1 and N 3 are scaled correspondingly so that the relative polarizations are fixed). The coexistence area is suppressed as N grows large, implying that the coexistence may vanish in the thermodynamic N → ∞ limit. However, with sufficiently accurate choice of interaction strengths, the coexistence region should be experimentally accessible with reasonable sized atom gases. We have not studied the scaling of the coexistence regions at the edge of the trap, but since the penetration length in figure 7 is given by the oscillator length, we expect the N → ∞ limit to yield a phase separation into a 23 core and a surrounding 12 shell.

Results
To better understand the nature of the pairing scheme in the coexistence areas, figure 9(a) shows the occupation numbers of different n-quantum number states for a symmetric case scheme, as in the homogeneous system (see figure 2). The unpaired atoms are distributed among the two components |1 and |3 and form a step at the Fermi surface.

Conclusions
We have studied the pairing of a three-component Fermi gas when two of the interspecies interaction channels are dominant. In a homogeneous system, we showed that different boundary conditions, namely fixing the chemical potential or fixing the atom numbers, produce qualitatively different results and phases. We have not considered the possibility of a phase separation in which the two pairing fields would be spatially separated in an otherwise uniform system.
For trapped systems, our BdG study reveals an interesting coexistence region where both pairing channels 12 and 23 are present. This is a mesoscopic effect and likely to vanish in the limit of a large system N → ∞, resulting in phase separation into shells of different pairing fields. However, the coexistence region is present at atom numbers relevant for atom gas experiments, making the observation of this intriguing double-gap prediction feasible.
There is already a wide range of standard experimental techniques for detecting such pairing correlations. For example, one could use radio-frequency spectroscopy [38][39][40] to drive atoms from the hyperfine states |1 and |3 separately into some fourth non-interacting state |e . Other possibilities include transforming pairing correlations into molecular 12 and/or 23 pairs through magnetic field sweeps [41] or optical molecular spectroscopy [42].