Construction of a giant vortex state in a trapped Fermi system

A superfluid atomic Fermi system may support a giant vortex if the trapping potential is anharmonic. In such a potential, the single-particle spectrum has a positive curvature as a function of angular momentum. A tractable model is presented in which the lowest and next lowest Landau levels are occupied. Different parameter regimes are identified and characterized. Due to the dependence of the interaction on angular momentum quantum number, the Cooper pairing is at its strongest not only close to the Fermi level, but also close to the energy minimum. It is shown that the gas is superfluid in the interior of the toroidal density distribution and normal in the outer regions. Furthermore, the pairing may give rise to a localized density depression in configuration space.


2
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Introduction
Among the features that distinguish a superfluid from a normal fluid is the quantization of fluid circulation. This implies that a superfluid responds to external rotation by forming one or more vortices, each of which carries an integer number of circulation quanta. In an infinite homogeneous fluid, it is well known that singly quantized vortices are energetically favourable for a given rotation rate, since the angular momentum is proportional to the quantum number and the energy is proportional to its square [1,2]. As a result, a rotating homogeneous superfluid does normally not form giant vortices, i.e., vortices with quantum number larger than one, but instead a lattice of singly quantized vortices will form. Such lattices are also observed in bosonic [3,4] and fermionic [5] superfluid gases contained in harmonic trapping potentials. The parameters in these experiments are such that the system is much larger than the core of a single vortex, so that the energy minimization argument sketched above goes through.
The situation is different if the length scale of a vortex is comparable to the system radius. In finite superconductors, it has been found numerically that giant vortices form in small samples at strong fields [6]. For a Bose condensate, the quantum number of a vortex depends on the shape of the potential used to confine the atoms as well as on parameters associated with rotation and inter-particle interactions. When the potential in the radial direction is harmonic, an array of singly quantized vortices is the energetically favourable state for all parameter values [7,8], but in an anharmonic potential, a Bose condensate will develop a giant vortex if the rotation speed is high enough and the interaction energy is weak enough [9]- [18].
Because a condensed Bose gas at zero temperature obeys a nonlinear Schrödinger equation that in the non-interacting limit reduces to the single-particle Schrödinger equation [2,8], it is possible to put up an exact argument proving the energetic stability of giant vortices in anharmonic traps and their instability in harmonic traps [9]. In a Fermi gas, no such criteria have yet been established. Therefore, as a first step towards understanding the rotational properties in superfluid Fermi gases, this paper shows that it is possible, within a BCS ansatz, to construct a situation where a giant vortex is the energetically favourable state. In order to make the problem tractable, the study is confined to the socalled lowest Landau level (LLL), i.e., the space of states without radial excitations [19]; the extension to the next lowest Landau level (NLLL) will also be studied. In section 2, the physical setting is described and definitions are introduced. In section 3, the solution of the equations is carried out. Section 4 discusses the physics of the solution in different parameter regimes. In section 5, the effects of extending the space of occupied states and approaching more realistic parameter values are discussed. Finally, section 6 provides a conclusion and outlook.

Lowest Landau level
Consider a system of spin-1/2 fermions confined in an anharmonic, i.e. quadratic plus quartic, trapping potential. The fermions are subject to a rotational force with a rotation frequency . (Alternatively, can be seen as a Lagrange multiplier applied to ensure a finite angular momentum L z .) The number of particles in each spin state is N, so that the total number is 2N. The confinement in the direction along the rotation axis is assumed to be so tight that the motion in that direction is frozen out, so that the system can be assumed to be two-dimensional. The Hamiltonian is where σ runs over the two spin indices ↑, ↓, and the interaction strength g is assumed negative, g = −|g|. Choosing units so that the particle mass and Planck's constanth are unity, the singleparticle Hamiltonian is The potential is so that ω is the trap frequency associated with the harmonic part and a is the degree of anharmonicity. The passage to dimensionless units is completed by setting ω = 1. The singleparticle spectrum associated with the Hamiltonian (2) can be solved numerically, resulting in a positive curvature for the energy as a function of angular quantum number m for a fixed radial quantum number n r , as seen in figure 1. In order to make the problem tractable one needs to assume that no states are occupied except those in the LLL, i.e. the levels with no radial nodes, n r = 0. The restrictions on physical parameters needed for this assumption to hold will shortly be explored. In general, the single-particle spectrum in the LLL around the minimum can be 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT approximated by a quadratic expression, where κ = m − m 0 is a shifted azimuthal quantum number; the offset m 0 is the quantum number at which the energy has its minimum; and the curvature α has to, in general, be computed numerically. When the anharmonicity a is very small, the energy spectrum can be solved perturbatively using harmonic-oscillator eigenfunctions, resulting in Thus, in this limit the minimum-energy quantum number is given by the curvature is α = a, and the ground-state energy 0 = 1 − m 2 0 + a. For not so small a, or high enough quantum numbers m, the corrections to the perturbative expressions become important, but they only result in quantitative shifts. In the following, we shall use the effective curvature parameter α, which may be slightly different from a, in the calculations, but the more physical discussions and order-of-magnitude estimates will be stated in terms of the actual anharmonicity a.
The number of particles is determined by the chemical potential µ, which in the weakly interacting limit is equal to the Fermi energy. For convenience, a Fermi quantum number κ F is defined through the relation so that the actual number of particles per spin state is N = 2κ F in the weakly interacting limit, but one must bear in mind that it may differ from that when interactions are finite. The parameters are assumed to be such that µ < 0, implying m 0 − κ F > 0. That way, the occupied part of the spectrum is bounded from below at some finite angular momentum quantum number m > 0, and one can treat the spectrum as if it is extended indefinitely in both the positive and negative κ-direction.
The coupling constant g introduced in (1) is a pseudopotential whose relation to the actual scattering length a s is in general nontrivial [20], but in the weak-coupling limit it reduces to the simple form g = 4πa s . The conventional dimensionless coupling parameter in the study of trapped Fermi gases is the product ζ = k F a s , the Fermi wavenumber times the s-wave scattering length. In terms of our Fermi quantum number κ F and coupling constant g, we obtain ζ = α 1/2 κ F g/4π, and we shall confine the study to small values of ζ, i.e. weak coupling, in order not to violate the conditions for the model to be valid. We now discuss the conditions for the LLL assumption to hold. Firstly, the NLLL, with n r = 1, must lie higher than the chemical potential. The spacing between the Landau levels is 2 (in the present units where ω = 1), so the condition is µ < 0 + 2 or in terms of physical quantities,

5
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Furthermore the condition µ < 0 stated above needs to be improved. Interactions broaden the occupation of fermions and smear out the Fermi level on a scale |g|, the interaction strength.
The distance between the lower Fermi level and the edge of the spectrum at m = 0 thus needs to be larger than |g|, that is, Note that according to the first inequality (8), a must be of the order N −2 or smaller, which means that the second restriction (9) is automatically fulfilled as long as − 1 is positive and not smaller than |g|N −2 , which is an extremely small quantity. The first inequality does, however, put quite severe limitations on the experimental implementation of the state, since it forces the anharmonicity to be less than or of the order of the inverse square of the number of particles. These matters will be discussed further in the concluding section.

Bogoliubov-de Gennes (BdG) analysis
Through a number of assumptions on the physical parameters, we have arrived at a description of a Fermi gas with a quadratic single-particle spectrum. This is the textbook case, and in principle, one may simply apply the usual textbook BCS expressions known for superconductors [1]. Two things remain to be sorted out: firstly, the gap must be assumed to be larger than the level spacing, which permits us to treat the levels as a continuum (for paired Fermi systems with a discrete spectrum, see [21]). Secondly, one must consider the variation of the coupling with quantum number κ. This will have an effect on the density profile and quasi-particle energy spectrum, as we shall see.

The BdG equation in two dimensions reads
where we have anticipated that the solutions u κ , v κ can be labelled by the shifted angular momentum quantum number κ. When the single-particle energy spectrum is symmetric around m = m 0 (or κ = 0), the gap function can be written This amounts to assuming a cylindrically symmetric gap function and corresponds to the usual assumption of a homogeneous gap function in an infinite system. In order to match the θ dependence, the amplitudes take on the form where φ m are single-particle eigenfunctions in the anharmonic trap and u κ and v κ are spaceindependent amplitudes. The BdG equation (10) is now multiplied through by eigenfunctions φ m and integrated, resulting in and The solution for the energies and amplitudes is similar to the traditional BCS expressions: Now turn to the self-consistency equation for the gap function, Inserting the calculated Bogoliubov amplitudes into the self-consistency equation, one obtains Multiply by φ m 0 +κ (r)φ m 0 −κ (r), integrate and rename the subscripts. The result is where the coupling matrix element was defined as In order to be able to produce a definite expression, let us make use of the results for a harmonic trap. Since (8) confines the parameters to lie in the regime of a very weak anharmonicity, the anharmonic trap can at worst bring in small quantitative corrections. The matrix element for the contact potential in the space of harmonic-oscillator eigenfunctions is [8] m, n|k, l = 1 2π Note that the potential is separable: 7 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT and furthermore, the potential is in the limit of large m 0 approximated by A separable potential presents a considerable simplification. The solution for the gap is where the constant is found by solving the self-consistency equation in its final form, This sum has, in general, to be computed numerically. Nevertheless, there is plenty of information to be extracted from its qualitative form, as we shall explore in the next section.

Parameter regimes
Consider the self-consistency equation (25). If the coupling is weak, is expected to be small as well. In this case, the largest contribution to the self-consistency sum comes from the states close to the Fermi surface, and the usual textbook BCS solution is recovered [1], where ω c is the cut-off frequency for the sum and N(0) is the density of states at the Fermi level. However, that solution is only valid when the interaction strength at the Fermi surface is larger than the separation of the energy levels; if it is smaller, one cannot neglect discreteness effects [21]. This criterion translates to On the other hand, if |g| is too large, pairing will happen away from the Fermi surface, and as will shortly be shown, this occurs when |g| > ακ 2 F . Conversely, the textbook case of only on-Fermi surface pairing needs the inequality |g| < ακ 2 F to be fulfilled. Combining the two inequalities, one finds that a criterion for avoiding both discreteness and off-Fermi surface pairing is which cannot be fulfilled for any value of κ F . We conclude that the standard textbook solution is not applicable in the present system. We now derive the criterion already used above, for states away from the Fermi surface to contribute to the self-consistency sum. When the sum is exhausted by terms far from the Fermi surface we can write 8 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Setting = 0 yields an inequality, The sum is estimated by inserting the maximum value of the summand and multiplying it by the effective range of the potential V κ , where √ 2m 0 is the width of the function V κ and V 2 0 = (4π 3 m 0 ) −1/2 is its value in the region where the summand is at its largest. The factors involving m 0 cancel out, so the critical coupling strength for an off-Fermi surface contribution reduces to where g c ∝ |ξ 0 | or where K is an unknown numerical constant which has to be included because of the order-ofmagnitude character of the estimation leading up to the inequality (31). In figure 2, the numerically computed value of is plotted in various slices of parameter space. When (31) is fulfilled and pairing happens off the Fermi surface, the magnitude of the gap is expected to become of order unity or larger in the present units. Indeed, it is seen how the = 0.01 level curve can be reasonably well fitted to the curve g = g c , with K ≈ 5.5, in all three plots. Beyond this curve, the discreteness of the spectrum is expected to be important and the BCS solution cannot be trusted to be quantitatively correct. The quasi-particle spectrum differs qualitatively between the different parameter regimes. This is illustrated in figure 3. When |g| is less than or of order g c , as in (a), one recognizes the textbook situation: pairing takes place mostly at the Fermi surface and the spectrum is equal to the absolute value of the single-particle spectrum except close to the Fermi surface at κ = ±60, where a gap opens. In figure 3(a), the full line representing the quasi-particle energy appears to have a cusp at the Fermi level, because the gap is too small to be visible on this scale; on a finer scale, the cusp is rounded off. The plot also includes the κ-dependent gap function κ and the summand in the self-consistency equation (25), S κ = V 2 κ /E κ . It is seen that the latter term, represented by the dotted line, has peaks at the Fermi level, which means that an appreciable contribution to the sum comes from terms close to the Fermi surface. If |g| were even smaller, the sum would be entirely exhausted by the terms close to the Fermi surface, but as was found above, the BCS solution would be altered by effects associated with the discreteness of the spectrum; therefore the dotted line also has a maximum near κ = 0. For stronger coupling, |g| > g c , the pairing happens predominantly close to κ = 0, where the interaction is at its maximum, as seen in figures 3(c) and (d). The quasi-particle spectrum is affected in this region, but at larger κ it resembles the usual BCS spectrum. The spatial profile of the gap function is completely insensitive to the profile in κ-space of the summand as long as the population is confined to the LLL. It is given by the self-consistency equation in configuration space (18), Thus, the gap function is a constant times the 2m 0 th radial harmonic-oscillator eigenfunction. At the same time, the density is much broader as a function of r, and it does depend on the coupling.  The density is with |v κ | 2 = (1 − ξ κ /E κ )/2. It is a peculiar feature of the LLL states that each eigenfunction φ m 0 −κ (r) is confined to a very narrow, shell-like region in configuration space, centred at r = (m 0 − κ) 1/2 . As a result, the distribution of particles in κ-space is mirrored by the spatial density profile, giving direct observational access to the distribution of particles in angular momentum space. For weak coupling, all states below the Fermi surface contribute, just as in ordinary BCS theory. The spatial profile as a function of the radial coordinate is a flattened, almost square profile, nicely mirroring the filled Fermi sea. As a result, only the interior of the toroidal density distribution is superfluid, while the outer regions are normal. For stronger coupling, the Bogoliubov amplitude v κ is depleted around κ = 0 and the density profile has a dip in the region where the gap function is at its largest. Mathematically, it is easy to see why the density is depressed at the trap bottom when the coupling is strong: the Bogoliubov amplitude v κ is equal to half its peak value when E κ = 2 κ + ξ κ ξ κ . Physically, it is the strong interactions that alter the occupation of the single-particle levels. As a result, one obtains a spatial region of strong Cooper pairing that shows up as a density depression. In the next section, we shall examine the conditions for this effect to persist when higher Landau levels are included.

Next lowest Landau level
In order to approach slightly more realistic territory in parameter space, let us now study what happens when the NLLL is brought into the calculation. Inspection of (8) reveals that this means allowing the product aN 2 /16 to be slightly larger than unity. This system is still not easily realizable in experiment, but it is a first correction to the even more difficult LLL case studied above. The new Landau level introduces a new parameter into the system, namely, the energy difference between Landau levels, which is equal to 2ω = 2 in dimensionless units. Thus, the balance between three main parameters determines the physics of the system as well as the validity of the truncation of the basis. The three are the chemical potential µ, the coupling strength g, and the energy difference 2ω. In addition, the balance between the Fermi level κ F and the width of the coupling, √ m 0 , plays a role in determining the physics.
The single-particle eigenfunctions in the NLLL are The perturbative expression for the energies in the quartic trap are (see figure 1) The position of the energy minimum of the NLLL coincides with that of the LLL, m 0 , to within 3/2 quanta. That shift can safely be neglected here. On the other hand, when many Landau levels are brought into the calculation, the resulting asymmetry of the single-particle energy spectrum is likely to break the cylindrical symmetry such that the order parameter becomes a superposition of several angular-momentum states. The resulting profile is expected to be an annulus pierced by a ring of singly quantized vortices [14]. Now consider the coupling matrix elements. They are 12 The other matrix elements can be obtained from these by a sign change, e.g. V 0010 κ,q = V 0001 κ,−q . It is seen that these new matrix elements are not separable in the way that simplified the calculations within the LLL. Moreover, the matrix elements V 0001 and V 1101 are smaller than the others by a factor m 1/2 0 ; however, that does not mean that the inter-level coupling cannot be neglected, since the matrix elements V 0011 and V 0101 are still as large as the intra-level matrix elements.
Employing the cylindrically symmetric assumption for (r) (keeping in mind that this may change when more Landau levels are brought into the problem), the ansatz for the Bogoliubov amplitudes has to be labelled by an angular momentum quantum number κ and a branch index j = 0, 1 (enumerating the two solutions that exist for every κ), The Bogoliubov equation for the (k, j)th mode is We have defined ξ κn = κn − µ and the matrix elements The self-consistency equation is in general and integrating one obtains The coupled equations are solved numerically. Figure 4 displays a few illustrative cases. It has been checked in trial calculations that the inclusion of more Landau levels does not alter the density profiles shown in this paper, and only marginally distorts the quasi-particle energy curves. Basically, there are four parameter regimes of interest. When µ 1 and |g| 1, the NLLL is not occupied at all. This case was studied in the previous section. When µ ∼ 1 and |g| 1, the Fermi surface cuts through both energy bands, but the coupling is not too strong, as in figures 4(a) and (b). There is a gap at the bottom of each of the two bands. The NLLL has a larger gap because the coupling is stronger for smaller κ. The associated density profile is expected to be bimodal, with an enhanced density in the interior of the toroidal cloud where two Landau levels are populated (cf [22]). Because the Fermi level has been limited to κ F = 60 in the numerical calculation, this 'wedding cake' structure does not come out as clearly in figure 4(b) as it would in a bigger system. Figures 4(c) and (d) illustrates the case where the chemical potential is small, but the coupling is comparable to the separation between the energy bands; µ 1 and |g| ∼ 1. The Fermi level lies far below the NLLL, but because of the strong coupling, both levels are deformed. The resulting density profile is still depleted in the region of strong Cooper pairing, because the contribution from the NLLL is modest. For stronger interactions, more Landau levels are coupled and the NLLL approximation would fail. Note that the effect of the strong coupling on the quasi-particle spectrum is much more pronounced than the effect on the density. Finally, in figures 4(e) and (f), both µ and |g| are comparable to the inter-level gap. The energy levels are heavily distorted by the strong coupling. The density near r = m 1/2 0 is no longer depleted, because the NLLL states now occupy that regime. The gap function takes on a slightly nontrivial shape in configuration space; this is again an effect of the relatively few single-particle levels involved in the calculation. Trial calculations indicate that when the number of states is increased, the gap function acquires a smoother shape, but the contribution from the radially excited states broadens its profile compared to the LLL case.
Extrapolating from the results presented here, one may observe that the density is affected by two competing effects: as the Fermi surface includes more and more Landau levels, a 'wedding cake' structure is expected in the density for weak interactions, but at the same time, strong interactions will deplete the interior of the toroidal density distribution.

Conclusion and outlook
This paper has shown that a Fermi superfluid trapped in an anharmonic potential can develop a giant vortex if it is rotated at high enough angular velocity, since the single-particle energy spectrum will then have its global minimum at a finite angular momentum. A giant vortex state has been constructed in the basis of the lowest radially excited states, the LLL. Although attaining the LLL regime experimentally puts severe constraints on precision, the theoretical construction of such a state demonstrates that there exist parameter regimes in anharmonic traps for which a Fermi gas develops a giant vortex. As long as only the LLL is occupied, the density and gap function are cylindrically symmetric and form a toroidal profile. The superfluid region is confined to the interior of the torus, while the edges are normal. Since Cooper pairing is strongest in the centre of the toroidal profile, the density is depleted in that region. The effect on the density appears, however, to be exclusive to the LLL. As more Landau levels are occupied, singly quantized vortices may form within the torus. Considering the general similarity between superfluid Fermi and Bose gases, it is not far-fetched to suspect that in the more general case, when all Landau levels contribute to the gap and density, one will find the same type of phase diagram that have been predicted to exist in Bose gases in anharmonic traps, containing vortex lattices, vortex lattices with a giant vortex in the centre, and a pure giant vortex [9,14]. At the present, the features of this vortex phase diagram can only be subject to speculation.
In order to have only the LLL populated in the quadratic plus quartic trap, the anharmonicity and rotation frequency have to be extremely fine-tuned. In addition, any mechanical stirring method [3] will bring in corrections to the anharmonic trap that may be hard to control. It is therefore probably better to control the angular momentum of the gas by first stirring it in a harmonic potential and subsequently applying the anharmonicity. As an alternative to an extremely small quartic trapping term, an optical hard-wall potential at a suitably large distance from the trap centre may be a more practical choice. It is possible that such a setup, with extensive amounts of fine-tuning, may be a more practical route towards accessing the LLL.
Nevertheless, this study should first and foremost be seen as a demonstration by construction that trapped Fermi superfluids can sustain giant vortices. Several important issues remain to be studied, notably the size and shape of the allowed parameter regime for giant vortices; the physics in the BEC-BCS crossover regime; and whether giant vortices can ever be the favourable rotational configuration in purely harmonic traps. We hope to address these issues in a future study.