Vortex formation in a rotating two-component Fermi gas

# Vortex formation in a rotating two-component Fermi gas

Harmen J. Warringa    Armen Sedrakian Institut für Theoretische Physik, Goethe-Universität Frankfurt am Main, Max-von-Laue-Straße 1, 60438 Frankfurt am Main, Germany
July 14, 2019
###### Abstract

A two-component Fermi gas with attractive -wave interactions forms a superfluid at low temperatures. When this gas is confined in a rotating trap, fermions can unpair at the edges of the gas and vortices can arise beyond certain critical rotation frequencies. We compute these critical rotation frequencies and construct the phase diagram in the plane of scattering length and rotation frequency for different total number of particles. We work at zero temperature and consider a cylindrically symmetric harmonic trapping potential. The calculations are performed in the Hartree-Fock-Bogoliubov approximation which implies that our results are quantitatively reliable for weak interactions.

## I Introduction

A characteristic feature of superfluids is the appearance of vortices when they are rotated. This fact has been used to demonstrate that a two-component atomic Fermi gas becomes a superfluid at sufficiently low temperatures Zwierlein05 (); Zwierlein06 (); Zwierlein06b (); Schunck07 (). A superfluid state can be created in such a gas by trapping fermionic atoms in two distinct hyperfine states. The interaction strength between the two components can be controlled by an external magnetic field. When the interactions are tuned to be attractive, and the atoms are cooled to sufficiently low temperatures, the components will form pairs via the Cooper instability. Due to this pair formation the Fermi gas becomes a Bardeen-Cooper-Schrieffer (BCS) superfluid as was envisaged in Refs. Stoof96 (); Houbiers97 ().

The response of the superfluid to rotation can be investigated by rotating the trapping potential with a certain frequency . For a non-rotating trap the entire gas will form a superfluid without vortices. Let us now imagine increasing the rotation frequency at zero temperature. Up to a certain critical rotation frequency, the superfluid will stay in the vortex-free state carrying zero angular momentum. For low temperatures the angular momentum will be quenched, as has also been observed experimentally AngularMomentumQuenching (). Above the critical frequency, angular momentum will be inserted in the gas by either unpairing the fermions near the edges of the gas, by formation of vortices, or by the combination of both effects. The goal of this paper is to compute the rotation frequencies at which these transitions take place.

Besides the experimental investigations Zwierlein05 (); Zwierlein06 (); Zwierlein06b (); Schunck07 (); AngularMomentumQuenching (), various theoretical studies of rotating two-component Fermi gases have been performed (for a review see e.g. Ref. Giorgini08 ()). The profile of a single vortex was analyzed in this context for the first time in Ref. Rodriguez01 () using a Ginzburg-Landau approach and in Ref. Nygaard03 () by solving the Bogoliubov-de Gennes equation. In Ref. Bulgac03 () it was concluded that a single vortex can induce a sizable density depletion at its core. This was also found to be the case for a vortex lattice Feder04 (). Density depletion is important, since it allows the vortices to be detected experimentally Zwierlein05 (). The vortex profile was investigated in a population imbalanced gas in Ref. Takahashi06 () and in situation in which the two components have unequal mass in Ref. Iskin08 (). Real-time dynamics of vortices has been studied in Ref. Bulgac11 ().

Vortex lattices in two-component Fermi gases have been examined in several ways in Refs. Feder04 (); Botelho06 (); Tonini08 (). At high rotation frequencies a completely unpaired phase is preferred over a vortex lattice. The critical rotation frequency corresponding to this transition was computed in Refs. Zhai06 (); Veillette06 (). A vortex lattice can also be destroyed by heating the gas. The corresponding critical temperature was computed in Refs. Feder04 (); Moller07 (). When the number of trapped components is unequal, a vortex lattice might be formed within the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase. Its melting temperature was investigated in Ref. Shim06 () and its upper critical rotation frequency in Ref. Kulic09 ().

The first analysis of the critical frequency for vortex formation in a two-component Fermi gas was performed in Ref. Bruun01 (). To obtain the Helmholtz free energy difference between a vortex with unit angular momentum located at the center of the trap and the vortex-free superfluid was estimated at . For a cylindrically symmetric infinite well as the trapping potential, was obtained by solving the Bogolibuov-de Gennes equation in Refs. Nygaard03 (); Nygaard04 (). The free energy difference arises from the loss of condensation energy at the vortex core, the kinetic energy of fermions circulating around the vortex core, and the energy needed to expand the cloud to accommodate the excess particles removed from the vortex core (the latter effect was not considered in Ref. Bruun01 (), but was taken into account in Refs. Nygaard03 (); Nygaard04 ()). When rotating the trap, the free energy decreases with , where is the angular momentum contained in the gas. For the vortex-free superfluid , while for the vortex with unit angular momentum , with the number of trapped particles. These estimates are correct if rotating the trap does not cause unpairing near the edges of the gas. The critical rotation frequency in this case is . This frequency is similar to the lower critical magnetic field in type-II superconductors Bardeen69 ().

Fermions confined in a rotating trap can unpair near the edges Bausmerth08 (); Urban08 (); Iskin09 (). This effect was not considered in Refs. Bruun01 (); Nygaard03 (); Nygaard04 (). In this paper we will take into account this possibility in order to obtain a more reliable value of . We will make detailed study of how varies with the interaction strength and the number of trapped particles. Furthermore we will compute the critical rotation frequency for unpairing. We consider a cylindrically symmetric harmonic trapping potential. Our calculations are carried out in the Hartree-Fock-Bogoliubov approximation. Therefore, we expect that our results are reliable only in the weak coupling limit.

Let us finally remark that vortices have also been observed in rotating Bose-Einstein condensates (BEC) of bosonic atoms BECvortices () and of bosonic dimers composed of fermionic atoms Zwierlein05 () (for a review see e.g. Ref. Fetter09 ()). The behavior of vortices through the BEC-BCS crossover was investigated experimentally in Ref. Zwierlein05 (). In this article we only discuss the BCS regime. A computation of the critical rotation frequency for vortex formation in a BEC is discussed in Ref. Stringari99 (). Theoretical studies of behavior of vortices through the BEC-BCS crossover have been performed in Refs. Kawaguchi04 (); Chien06 (); Sensarma06 ().

This article is organized as follows. In Sec. II we will introduce the action from which one can derive the properties of the Fermi gas. To achieve this in practice we will adopt the two-particle irreducible (2PI) effective action, which we explain in Sec. III. From the 2PI effective action one obtains the Dyson-Schwinger equation which is the main equation we have to solve. This can be achieved by finding the solution of the Bogoliubov-de Gennes equation, which is explained in Sec. IV. The numerical methods by which we have solved the Bogoliubov-de-Gennes and the Dyson-Schwinger equation are discussed in Secs. V and VI respectively. The reader who is not interested in the details of the calculation can immediately go to Sec. VII where we present the results. The phase diagrams presented in Figs. 14 and 15 are our main results. We draw our conclusions in Sec. VIII. Several details are relegated to the appendices. In Appendix A we review the 2PI effective action. A derivation of the Bogoliubov-de Gennes equation is presented in Appendix B. To solve the Bogoliubov-de Gennes equation numerically, we will use a basis based on Maxwell polynomials. We discuss the computation of the quadrature weights and nodes of these polynomials in Appendix C. In Appendix D we derive the representation of the single particle Hamiltonian in the basis we will employ.

## Ii Setup

Let us consider a two-component Fermi gas in which -wave interactions are dominant and label its components by . Typically in experimental realizations the higher partial waves can be neglected and the superfluidity is driven by attractive -wave interactions. The interactions among the same components can be neglected since the Pauli principle admits -wave interactions only between the different species. We will denote the -wave interaction potential as and specify it below. Under these assumptions the interacting Fermi gas is described by the following action (see e.g. Ref. Stoof ()), where

 Skin=∑α=↑,↓∫dXψ∗α(X)[ℏ∂∂τ+H(Ω)−μα]ψα(X), (1) Sint=12∑α=↑,↓∫dX∫dYψ∗α(X+)ψ∗−α(Y+)×V(x−y)δ(τx−τy)ψ−α(Y)ψα(X). (2)

Here is the (path-integral) quantum field corresponding to the component, and . We write integration over spatial coordinates and imaginary time as with the inverse temperature . Here and in the rest of the article and is an infinitesimal small positive number. We have inserted and in Eq. (2) in order to maintain the correct ordering of the fields in the path-integral. We will achieve this in a different way for the kinetic term and explain this at the end of Appendix B. Furthermore is equivalent to , denotes the chemical potential, and is the single-particle Hamiltonian. We will assume that the particles are trapped in a potential that is rotating in the - plane with angular frequency . It is then convenient to perform the calculation in the rotating frame. The single-particle Hamiltonian in the rotating frame reads (see e.g. Ref. Stoof ())

 H(Ω)=p22M+U(x)−ΩLz, (3)

where is the fermion mass and is the -component of the angular momentum. The trapping potential realized in experiments is typically harmonic. In this paper we will study a cylindrically shaped trap given by the potential

 U(x)=12Mω2(x2+y2), (4)

which implies confinement in the - plane, and infinite extension in the -direction. In an experiment this regime can be reached by choosing the trapping frequency in the -direction much smaller than in the --direction. In cylindrical coordinates, , the single-particle Hamiltonian reads

 H(Ω)=ℏ22M(−d2dρ2−1ρddρ+L2zℏ2ρ2−d2dz2)+12Mω2ρ2−ΩLz, (5)

where . The normalized eigenfunctions of this Hamiltonian can be written as a product of three functions,

 ψ0nmpz(x)=1√LRnm(ρ)fm(ϕ)eipzz, (6)

with the radial quantum number , the angular momentum quantum number , and the momentum in the -direction with . The constant denotes the length of the system in the -direction. In this article we will consider the limit . Then . The function is a normalized eigenfunction of and is given explicitly by

 fm(ϕ)=1√2πeimϕ. (7)

The radial eigenfunctions are given by

 Rnm(ρ)=cnmL|m|n(¯ρ2)¯ρ|m|e−¯ρ2/2, (8)

where denotes the generalized Laguerre polynomial which has degree . Furthermore with the harmonic oscillator length . Normalization gives . The energy spectrum of is given by

 ϵ0nmpz=ℏω(1+2n+|m|)−ℏΩm+p2z2M. (9)

At low enough temperatures and densities, the typical wavelength of the particles will be much longer than the range of the interaction potential. In that case the detailed structure of the potential is unimportant and the only relevant interaction parameter is the scattering length . To perform a calculation in this situation, one can just choose the most convenient potential that has scattering length . Following Ref. Bruun99 () we will use the Huang-Yang potential Huang57 ()

 V(r)=gδ(r)∂∂rr, (10)

where the coupling constant . Pairing between fermions requires attractive interactions, that is and hence . Note that the Huang-Yang potential is not equivalent to an ordinary -function potential , because the derivative operator also acts on the fields in Eq. (2). The advantage of the Huang-Yang potential is that all relevant physical quantities become automatically convergent Bruun99 (). Furthermore, the coupling constant does not have to be renormalized so that the foregoing relation between the scattering length and the coupling constant always holds Bruun99 (). In the case of an ordinary -function potential one will encounter divergences. This will require a regularization prescription and a renormalization of the coupling constant.

In order to perform calculations it is convenient to rewrite the action in the Nambu-Gor’kov basis. For that purpose we introduce the Nambu-Gor’kov fields

 Ψ(X)=(ψ↑(X)ψ∗↓(X)), (11)

and write the kinetic part of the action, Eq. (1), as

 Skin=−ℏ∫dX∫dX′Ψ†(X)G−10(X,X′)Ψ(X′), (12)

where the bare inverse Nambu-Gor’kov propagator reads

 G−10(X,X′)=−1ℏ⎛⎝ℏ∂∂τ+H(Ω)−μ↑00ℏ∂∂τ−H(Ω)∗+μ↓⎞⎠×δ(X−X′). (13)

Here we used the fact that . In the Nambu-Gor’kov basis the interaction part of the action becomes

 Sint=−∑α=±∫dXdYdX′dY′Ψ†(X′)σαΨ(X)×Ψ†(Y′)σ−αΨ(Y)Vα(X,Y;X′,Y′), (14)

where and . The potential is given by

 (15)

Because the Huang-Yang potential contains a derivative operator, we had to introduce several -functions in order to separate the potential operator from the quantum fields.

## Iii 2PI effective action

To study the interacting Fermi gas we will compute the resummed propagator and the grand potential by using the two-particle irreducible (2PI) effective action 2PI (). Readers who are not interested in the details of this formalism can immediately continue with Sec. IV where we discuss the Bogoliubov-de Gennes equation which follows from the 2PI effective action.

The 2PI effective action is also known as the Cornwall-Jackiw-Tomboulis formalism and is equivalent to the Luttinger-Ward functional approach LW (). It has been applied previously to investigate pairing in atomic gases Haussmann07 () and in quark matter CSC2PI ().

The main advantage of the 2PI effective action is that it is a functional method which generates the resummed Nambu-Gor’kov propagator and the corresponding grand potential. Here, we will truncate the 2PI effective action at order , which is equivalent to the Hartree-Fock-Bogoliubov approximation. This leads to the well known Bogoliubov-de Gennes equation. Another advantage of the 2PI method is that any truncation can be systematically improved straightforwardly by taking into account higher order diagrams. This is necessary when extending our results to a strongly coupled Fermi gas.

The 2PI effective action reads (see Appendix A for details) 2PI ()

 Γ[G]=−TrlogG−1−Tr(G−10G−1)+Γ2[G], (16)

where is the sum of all 2PI diagrams generated from with propagators . The interaction vertex can be directly read off from Eq. (14). We have displayed the Feynman rules for computing the 2PI effective action in Fig. 1. All diagrams contributing to up to order are displayed in Fig. 2.

By minimizing with respect to one obtains the Dyson-Schwinger equation

 G−1=G−10−Σ[G], (17)

where the 1PI self-energy is . The grand potential is equal to the minimum value of and reads

 ΦG=−1βTrlogG−1−1βTr(Σ[G]G)+1βΓ2[G], (18)

where is now the solution of the Dyson-Schwinger equation, Eq. (17).

In order to perform a practical calculation, one has to truncate at some order. In this paper we will take into account the contributions of order to , which are represented by the first two diagrams in Fig. 2. Such truncation leads to the Hartree-Fock-Bogoliubov approximation. This approximation is expected to give accurate results for small coupling . In general the Hartree-Fock-Bogoliubov approximation can also be applied to strongly correlated systems if the interaction kernel is renormalized appropriately. Such renormalization entails a resummation of ladder diagrams in the Lippman-Schwinger equation. Since we are interested in the weakly interacting BCS limit there is no need to do this.

Applying the Feynman rules to the first two diagrams in Fig. 2 we find that to order , is given by

 Γ2[G]=−1ℏ∑α=±∫dXdYdX′dY′Vα(X,Y;X′,Y′)×tr[G(X,X′)σα]tr[G(Y,Y′)σ−α]+1ℏ∑α=±∫dXdYdX′dY′Vα(X,Y;X′,Y′)×tr[G(X,Y′)σαG(Y,X′)σ−α]. (19)

The relative minus sign arises because the first diagram in Fig. 2 contains two closed loops whereas the second diagram contains only one closed loop when following the arrows.

The self-energy can be obtained by differentiating Eq. (19) with respect to . This yields

 Σ[G](X,X′)=−2ℏ∑α=±σa∫dYdY′tr[G(Y,Y′)σ−α]×Vα(X′,Y;X,Y′)+2ℏ∑α=±∫dYdY′σαG(Y,Y′)σ−α×Vα(X′,Y;Y′,X). (20)

The first term is the Hartree self-energy, the second one is the pairing contribution.

The last two equations can be simplified by computing the traces and inserting the explicit expressions for and . One has then to act the Huang-Yang potential on the propagators. It can be shown Bruun99 () that the diagonal components of are finite in the limit where . For these components the Huang-Yang potential acts as an ordinary -function potential. We will use this fact to simplify the expressions involving the diagonal components of . On the other hand the off-diagonal components of will have in general a singularity proportional to when Bruun99 (). In such situations the full form of the Huang-Yang potential needs to be taken into account.

For later convenience we will now define the pairing field as

 Δ(x)≡∫d3x′V(x−x′)G↑↓(x,τ;x′,τ). (21)

Following Ref. Bruun99 () we will split the off-diagonal component of the Nambu-Gor’kov propagator into a singular and regular part

 limx′→xG↑↓(x,τ;x′,τ)=Cr+Greg↑↓(X,X), (22)

where is some constant and the superscript reg indicates the regular part of . By inserting the Huang-Yang potential in Eq. (21) one can see that the pairing field does not contain any singularity Bruun99 ()

 Δ(x)=gGreg↑↓(X,X). (23)

The diagonal components of the propagator can be expressed in terms of number densities. To find the relation between the propagator and the densities we use the definition of the propagator in terms of field operators

 Gij(X,X′)=−⟨Tτ^Ψi(X)^Ψ†j(X′)⟩=θ(τ′−τ)⟨^Ψ†j(X′)^Ψi(X)⟩−θ(τ−τ′)⟨^Ψi(X)^Ψ†j(X′)⟩, (24)

here denotes time ordering in imaginary time. Hence the number densities of the two species are related to as

 n↑(x) = ⟨^ψ†↑(x,τ)^ψ↑(x,τ)⟩=G↑↑(X,X+), (25) n↓(x) = ⟨^ψ†↓(x,τ)^ψ↓(x,τ)⟩=−G↓↓(X,X−). (26)

From Eq. (24) it follows that and . This can be used to show that which implies that

 (27)

Here we used the fact that the singular part of is a function of , so that and could be interchanged without changing the result.

We can now use the above definitions to simplify Eqs. (19) and (20). After computing the traces we find that to order , is given by

 1βΓ2[G]=g∫d3xn↑(x)n↓(x)+∫d3xG↑↓(X,X)∗Δ(x). (28)

The 1PI self-energy to order is

 Σ[G](X,X′)=1ℏ(gn↓(x)Δ(x)Δ∗(x)−gn↑(x))δ(X−X′). (29)

It follows that so that the grand potential , Eq. (18), to order becomes

 ΦG=−1βTrlogG−1−g∫d3xn↑(x)n↓(x)−∫d3xG↑↓(X,X)∗Δ(x). (30)

Since the particle number in the trap is fixed we solve the equations for the chemical potentials to obtain the desired number of fermions in each hyperfine state. The appropriate thermodynamic potential is then the Helmholtz free energy, which reads

 F=ΦG+μ↑N↑+μ↓N↓, (31)

where denote the total number of particles of a particular species. Since we consider a cylindrically shaped trap, and are proportional to the length of the trap in the -direction . Since is taken to be infinite is convenient to consider instead the free energy and particle number per unit of the harmonic oscillator length in the -direction. For this reason we define

 F=FL/λ,N↑,↓=N↑,↓L/λ. (32)

Furthermore we will define to be the total number of particles per unit of length in the -direction, .

## Iv Bogoliubov-de Gennes equation

To proceed, we will insert the explicit expression for the 1PI self-energy, Eq. (29), into Eq. (17). This yields the Dyson-Schwinger equation for the Nambu-Gor’kov propagator,

 G−1(X,X′)=−1ℏ(ℏ∂∂τ+H)δ(X−X′), (33)

with

 H=(H(Ω)−μ↑+gn↓(x)Δ(x)Δ∗(x)−H(Ω)∗+μ↓−gn↑(x)). (34)

To solve the Dyson-Schwinger equation, one first inverts both the left- and right-hand sides of Eq. (33). As explained in detail in Appendix B this can be achieved by solving the Bogoliubov-de Gennes equation Gennes64 ()

 H(ui(x)vi(x))=Ei(ui(x)vi(x)). (35)

The functions and have to be normalized as . Using the explicit expression of , Eq. (100), and Eqs. (25), (26) one can now read off the densities

 n↑(x) = ∑if(Ei)|ui(x)|2, (36) n↓(x) = ∑if(−Ei)|vi(x)|2, (37)

where is the Fermi-Dirac distribution function.

As follows from Eqs. (23) and (100) to obtain we need to extract the regular part of the propagator

 G↑↓(x,τ;x′,τ)=∑if(Ei)ui(x)v∗i(x′), (38)

in the limit . To do so we will use the method proposed in Ref. Bruun99 () with the improvements suggested in Refs. Bulgac02 (); Grasso03 ().

The sum over all modes in Eq. (38) is logarithmically divergent for . The singular part of arises from the modes in the integrand with large negative energy. To obtain the regular part in the limit we will first subtract this large energy contribution. For this purpose we define

 νc(x)=∑|Ei|

where denotes an energy cutoff introduced to regulate the logarithmic divergence in . The part of dominated by the modes with large negative energies can be approximated as Bruun99 (); Bulgac02 (); Grasso03 () with

 K(x,x′;Ec)=∑Es<ϵi

where here and below and denote the eigenvectors and eigenvalues of the Hartree-Fock Hamiltonian

 HHF=H(Ω=0)−μ+gn(x). (41)

Here and . The low-energy cut-off in Eq. (40) is arbitrary, except that it should be chosen positive in order to avoid singularities arising from the Fermi surface. As an alternative to introducing a low-energy cutoff, one can add a small imaginary part to as done in Refs. Bruun99 (); Bulgac02 (); Grasso03 (). In that case the integrand of has a peak near the Fermi surface, which makes the numerical integration over difficult. One can reduce this peak by increasing the magnitude of the imaginary part. However, that will worsen the large negative energy approximation of . The low-energy cutoff which we apply here does not suffer from these problems.

Let us next define . Because contains the logarithmic divergent part of , the difference is finite, and hence converges for large enough . There is some freedom in choosing . For example, we could have left out the term in the Hartree-Fock Hamiltonian as in Refs. Bruun99 (); Bulgac02 (). However, by including this term, approximates much better, which implies that a much smaller value of is sufficient to compute accurately Grasso03 ().

Summarizing the discussion above, we found that in the limit , we can write . Following Refs. Bulgac02 (); Grasso03 () the singular part can now be obtained by making use of the Thomas-Fermi approximation. In the limit one finds that

 K(x,x′;∞)=K(x,x;Ec′)−12∫k

Here is an infinitesimal small positive number and is a second energy cutoff that can be chosen different from . For large enough the sum of first two terms in Eq. (42) is convergent. The inhomogeneous wavevector cut-off can be found from

 ℏ2k2c′(ρ)2M+12Mω2ρ2−μ+gn(ρ)=Ec′. (43)

The last term of Eq. (42) contains the singularity. One can now perform the integration over analytically, which in the limit gives Bulgac02 (); Grasso03 ()

 G↑↓(x,τ;x′,τ)=−MΔ(x)4πℏ21|x−x′|+Greg↑↓(x,τ;x,τ), (44)

where

 Greg↑↓(x,τ;x,τ)=νs(x)−Δ(x)K(x,x;Ec′)+Δ(x)M2π2ℏ2[kc′(ρ)−12kF(ρ)log(kc′(ρ)+kF(ρ)kc′(ρ)−kF(ρ))]. (45)

Here we have introduced the length of the Fermi wavevector which can be found from the equation

 ℏ2k2F(ρ)2M=μ−gn(ρ)−12Mω2ρ2−iγ. (46)

To obtain we have to multiply Eq. (45) by as follows from Eq. (23).

Inserting Eq. (104) into Eqs. (30) and (31) gives the Helmholtz free energy

 F=−∑i[|Ei|2+1βlog(1+e−β|Ei|)]+∑iϵi−∫d3xG↑↓(x,τ;x,τ)∗Δ(x)−g∫d3xn↑(x)n↓(x)+μ↑N↑+μ↓N↓. (47)

Although some of the individual terms in the last equation are ultraviolet divergent, their sum, and hence the Helmholtz free energy, is ultraviolet finite. The divergence present in the sum over the eigenvalues of the Bogoliubov-de Gennes matrix is canceled by the sum over the eigenvalues of the Hartree-Fock Hamiltonian and by the logarithmic divergence originating from . This can be made clearer by expressing in terms of and . One finds

 F=−∑|Ei|

In the absence of a trapping potential is known analytically; then it can be seen that Eq. (48) is ultraviolet finite. As we will see in the next section, is also finite in the general case.

## V Solving the Bogoliubov-de Gennes-equation

In the previous section we have reduced the Dyson-Schwinger equation to a nonlinear equation of the form

 (Δ(x),n↑,↓(x),N↑,↓)=F(Δ(x),n↑,↓(x),μ↑,↓). (49)

We will now discuss how to compute the function in practice. In the next section we will solve the Dyson-Schwinger equation.

First we will use the symmetries of our problem to simplify the analysis. For zero rotation frequency the superfluid is in a vortex-free phase. Then, the pairing field will be a function of the radial coordinate only, i.e. , where .

We will assume that the first vortex that appears when increasing the rotation frequency carries one unit of angular momentum and is located at the center of the trap. The pairing field for such vortex has the following form: . After the single vortex has appeared, a vortex lattice can be formed by further increasing the rotation frequency.

For these reasons we will make the following ansatz for the pairing field

 Δ(x)=~Δ(ρ)exp(ikϕ), (50)

where is the winding number (unit of angular momentum) of the vortex at the center. Hence the case corresponds to the vortex-free phase. To determine the onset of the vortex phase, we have to compare the Helmholtz free energies with and . Because of the cylindrical symmetry of the trap and the fact that a possible vortex is located at the origin, the number densities are a function of only, i.e. .

In this case the solutions of the Bogoliubov-de Gennes equation, Eq. (35), are of the following form

 ui(x) = 1√Leipzz1√2πeimϕ1√λρunmpz(ρ), (51) vi(x) = 1√Leipzz1√2πei(m−k)ϕ1√λρvnmpz(ρ), (52)

which can be verified by inserting these expressions into Eq. (35). This also yields the Bogoliubov-de Gennes equation for and ,

 Hmpz(unmpz(ρ)vnmpz(ρ))=Enmpz(unmpz(ρ)vnmpz(ρ)), (53)

where

 Hmpz=(Hm(Ω)−μ↑+gn↓(ρ)~Δ(ρ)~Δ(ρ)−Hk−m(Ω)+μ↓−gn↑(ρ)), (54)

and

 Hm(Ω)=ℏ22M(−d2dρ2−14ρ2+m2ρ2+λ2η2)+12Mω2ρ2+p2z2M−mℏΩ. (55)

Due to the factor in Eqs. (51) and (52) is a bit different from the single particle Hamiltonian defined in Eq. (5). We have inserted this factor for later convenience. Furthermore, by inserting the term we have modified the centrifugal potential in such a way that it becomes regular at . The original centrifugal potential is reproduced for . As we explain in Appendix D, this slight modification is necessary for the numerical computation of the wavefunctions and energies. In our computations we have taken and checked that the results are completely stable if is varied by orders of magnitude around this value.

From the normalization condition on and it follows that and have to be normalized as

 1λ∫∞0dρ[unmpz(ρ)2+vnmpz(ρ)2]=1. (56)

To solve the Bogoliubov-de Gennes equation numerically, we have to discretrize Eq. (53). One can do this by expanding the wavefunctions and in a certain basis. For practical purposes, this basis has to be truncated, i.e., we will represent and by a finite number of basis functions . More specifically we will write

 unmpz(ρ)=N∑i=1ciℓi(ρ)vnmpz(ρ)=N∑i=1diℓi(ρ), (57)

here and are the expansion coefficients, which depend on the quantum numbers , and . We will require the basis functions to be orthonormal in the following way

 1λ∫∞0dρℓi(ρ)ℓj(ρ)=δij. (58)

From Eq. (56) it then follows that the expansion coefficients have to be normalized as

 N∑i=1(c2i+d2i)=1. (59)

Equation (53) can now be transformed into an ordinary eigenvalue equation for a matrix which reads

 (¯Hm−¯μ↑¯Δ¯Δ−¯Hk−m+¯μ↓)(cd)=Enmpz(cd). (60)

Here , and are real symmetric matrices which are given by

 (¯Hm)ij = 1λ∫∞0dρℓi(ρ)Hm(Ω)ℓj(ρ), (61) (¯μα)ij = 1λ∫∞0dρℓi(ρ)[μα−gn−α(ρ)]ℓj(ρ), (62) (¯Δ)ij = 1λ∫∞0dρℓi(ρ)~Δ(ρ)ℓj(ρ). (63)

Similarly we define as the representation of the Hartree-Fock Hamiltonian defined in Eq. (41). Once these matrices have been computed explicitly, Eq. (60) can be solved numerically using standard linear algebra routines. From the solutions the wavefunctions and the expressions for and can be constructed. We will give the explicit expressions at the end of this section.

Since truncating the basis is an approximation, we have to choose the basis carefully in order to make sure that the basis functions can describe the exact solution to good accuracy. We can always improve the accuracy of the truncation by taking a larger value of , but the drawback is that this increases the computational cost as well.

A good basis has to be able to describe the wave-functions in the case , which are the radial single particle wave functions given in Eq. (8) times . For that reason we will choose a basis of the following form

 ℓi(ρ)=√w(ρ/λ)li(ρ/λ), (64)

where and is a set of linear independent polynomials of degree . In this basis all single particle wave functions can be represented exactly by a finite number of basis functions. For the calculation of vortices this is a desirable feature, because in that case mixing between states with different angular momentum quantum numbers occurs.

From Eq. (58) it follows that the polynomials have to be chosen orthogonal with respect to weight function , i.e.

 ∫∞0dxw(x)li(x)lj(x)=δij. (65)

The set of polynomials of increasing degree that are orthonormal to each other with the weight function on the interval are called Maxwell polynomials Shizgal81 (). We will write the Maxwell polynomials for as , where denotes the degree.

One could have chosen as basis functions the Maxwell polynomials directly, for example . However, the computation of the Bogoliubov-de Gennes matrix, Eq. (60), becomes much easier if we use Lagrange interpolating functions as basis functions. These Lagrange interpolating functions are a particular linear combination of Maxwell polynomials and will be specified below. This approach is generally known as the Discrete Variable Representation (DVR) method, or alternatively as the Lagrange mesh discretization Light85 (); Baye86 (); Szalay93 (). By applying this method one can obtain highly accurate values for the energies and wavefunctions Baye02 (). The DVR method can be applied to any set of orthogonal polynomials Baye99 (), and is not restricted to Maxwell polynomials. To our knowledge, this is the first time that the DVR method is used with Maxwell polynomials.

The DVR method is based on the Gaussian quadrature formula

 ∫∞0dxw(x)f(x)≈N∑n=1wnf(xn). (66)

Here and in the following the nodes are the roots of the Maxwell polynomial and the corresponding quadrature weights. The integration formula Eq. (66) is exact for all polynomials of degree less than . From the properties of the orthogonal polynomials, one can show that all roots are real and all weights are positive. The nodes and the weights are the only non-trivial properties of the Maxwell polynomials one needs to know in order to apply the DVR method. Since the Maxwell polynomials are non-standard polynomials, we review the computation of its nodes and weights in Appendix C.

In the DVR method one chooses the functions to be the Lagrange interpolating functions through the nodes . Explicitly, these functions read

 li(x)=1√wiN∏n=1,n≠ix−xnxi−xn. (67)

It follows directly that the polynomials satisfy the following useful property, . Because the combination is a polynomial of degree , the Gaussian quadrature is exact for this combination. Hence we can use the Gaussian quadrature to show that satisfies the required orthonormality condition given in Eq. (65),

 ∫∞0dxw(x)li(x)lj(x)=N∑k=1wk√wiwjδikδjk=δij. (68)

We will now compute the matrices that appear in the Bogoliubov-de Gennes equation, Eqs. (61)-(63) explicitly. In the DVR method one uses the Gaussian quadrature to compute the integrals. In this way we find that

 (¯μα)ij ≈ μαδij−gn−α(λxi)δij, (69) (¯Δ)ij