Super-shell structures and pairing in ultracold trapped Fermi gases

TL;DR

This paper investigates super-shell structures and pairing phenomena in ultracold trapped Fermi gases, combining analytical and numerical methods to understand their energy levels and pairing gaps under various trapping potentials.

Contribution

It provides an analytical approach to super-shell structures and pairing gaps in ultracold Fermi gases, extending previous Hartree-Fock results with periodic orbit and WKB calculations.

Findings

Super-shell structures are analytically characterized.

Level densities influence pairing gaps.

Super-shell effects are predicted for attractive interactions.

Abstract

We calculate level densities and pairing gaps for an ultracold dilute gas of fermionic atoms in harmonic traps under the influence of mean field and anharmonic quartic trap potentials. Super-shell structures, which were found in Hartree-Fock calculations, are calculated analytically within periodic orbit theory as well as from WKB calculations. For attractive interactions, the underlying level densities are crucial for pairing and super-shell structures in gaps are predicted.

Figures

Equations

$$H$$

$$U(r)=\frac{2\pi {\hslash }^{2}a}{m}\rho (r).$$

$$E_F=\frac{{\hslash }^2{k}_F^2(r)}{2m}+\frac{1}{2}m{\omega }^2{r}^2+U(r).$$

$$\rho (r)$$

$$U(r)\simeq \epsilon r^4,$$

$$g(E)=\frac{E^2}{{\left(\hslash \omega \right)}^3}\left(1+Re\sum _{k=-\infty }^{\infty }{\left(-1\right)}^k\mathcal{M}{e}^{i2\pi kE/\hslash \omega }\right).$$

$$\mathcal{M}=\frac{\hslash }{k\sigma }\left(e^{-i2k\sigma /\hslash -i\pi /2}+e^{-i3k\sigma /\hslash +i\pi /2}\right),$$

$$g(E)=\frac{E^2}{(\hslash \omega {)}^3}+\frac{4}{\pi \epsilon }\sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k}{k}\times$$

$$\cos \left(\frac{k}{\hslash }\left(\frac{2\pi E}{\omega }-\frac{5}{2}\sigma \right)\right)\sin \left(\frac{k\sigma }{2\hslash }\right).$$

$$n_F=E/\hslash \omega =\sqrt{2S\hslash \omega /|\epsilon |}.$$

$${\mathcal{R}}_{n_{F}l}(r)\simeq \frac{2}{\sqrt{\pi }}\frac{\sin (k_l(r)r+\theta )}{k_l^{1/2}(r)r},$$

$$k_l^2(r)=2{\tilde{n}}_F-r^2-l(l+1)/r^2.$$

$$E_{n_F,l}$$

$$=$$

$$D\equiv |E_{n_F,l\equiv 0}-E_{n_F,l=n_F}|=\epsilon n_F^2/2$$

$$E_{n_F,l}-{\tilde{n}}_F\hslash \omega$$

$$I={\int }_{-1}^1{\left(1-x\sqrt{1-l(l+1)/{\tilde{n}}_F^2}\right)}^{3/2}\frac{dx}{\sqrt{1-x^2}}$$

$$D=2/\left(3{\pi }^2\right)a{n}_F^{3/2}\hslash \omega \left(8\sqrt{2}/3-\pi \right).$$

$$\left(8\sqrt{2}/3-\pi \right)2/\left(3{\pi }^2\right)a{n}_F^{3/2}=S.$$

$$I$$

$$E_{n_F,l}-\tilde{n}\hslash \omega$$

$$1=\frac{G}{n_F^2}{\int }_0^{\sim 2n_F}\frac{g(E)dE}{\sqrt{(E-\mu {)}^2+\Delta (\mu {)}^2}}.$$

Full text

Ultracold atomic gases have recently been used to create novel quantum many-body systems such as strongly interacting high temperature superfluids of fermions, Bose-Einstein condensates, Mott insulators in optical lattices, etc. These lab phenomena have a strong overlap with condensed matter BCS , nuclear BM and neutron star physics Pines . Finite fermion systems such as atoms in traps, nuclei, helium and metal clusters, semiconductor quantum dots, superconducting grains, etc., have additional interesting quantum structures such as level spectra, densities and pairing. These will be observable as temperatures are further lowered in atomic trap experiments. The high degree of control over physical parameters, including interaction strength and density, makes the atomic traps marvelous model systems for general quantum phenomena.

The purpose here is to calculate the level spectra, densities and pairing for zero-temperature Fermi gases in harmonic oscillator (HO) traps with anharmonic and mean field perturbations, and to show that novel super-shell structures appear in both level densities and pairing. In calculating level spectra by analytical periodic orbit theory and WKB as well as numerical Hartree-Fock, we also relate these different theoretical approaches to one another.

We treat a gas of $N$ fermionic atoms of mass $m$ in a HO potential at zero temperature, interacting via a two-body interaction with s-wave scattering length $a$. We shall mainly discuss a spherically symmetric trap and a dilute gas (i.e. where the density $\rho$ obeys the condition $\rho|a|^{3}\ll 1$) of particles with two spin states of equal population. The Hamiltonian is then given by

[TABLE]

We will consider both external anharmonic potentials of the form $U=\varepsilon r^{4}$ and particle interactions: $U({\bf r}_{i})=(2\pi\hbar^{2}a/m)\sum_{j\neq i}\delta^{3}({\bf r}_{i}-{\bf r}_{j})$. When interactions are weak, the latter can be approximated by the mean field potential

[TABLE]

For a large number of particles and $U=0$ the Fermi energy is $E_{F}=\tilde{n}_{F}\hbar\omega$ where $n_{F}=\tilde{n}_{F}-3/2\simeq(3N)^{1/3}$ is the HO quantum number at the Fermi surface. The HO shells are highly degenerate with states having angular momenta $l=n_{F},n_{F}-2,...,mod(n_{F},2)$, due to the $U(3)$ symmetry of the 3D spherically symmetric HO potential. However, interactions split this degeneracy. In the Thomas-Fermi (TF) approximation (see, e.g., PS ) the Fermi energy is

[TABLE]

The density $\rho(r)=k_{F}^{3}(r)/3\pi^{2}$ is

[TABLE]

inside the cloud $r\leq R_{TF}=a_{osc}\sqrt{2\tilde{n}_{F}}$, where $\rho_{0}=(2\tilde{n}_{F})^{3/2}/3\pi^{2}a_{osc}^{3}$ is the central density TF3 . For convenience we set the oscillator length $a_{osc}=\sqrt{\hbar/m\omega}=1$ in the following.

Taylor expanding the density and thereby also the mean field of Eq. (2) around the center gives

[TABLE]

the first term will simply incorporate a constant shift in energies whereas the term quadratic in radius renormalizes the HO frequency as $\omega_{{\it eff}}=\omega\sqrt{1-6\pi a\rho_{0}/R_{TF}^{2}}$. The third term is quartic in radius and is therefore also of the same form as the external potential

[TABLE]

with $\varepsilon=(3\pi\hbar^{2}a/4m)\rho_{0}/R_{TF}^{4}$. Both the pure quartic potential and the mean field potential of Eq. (2) are anharmonic and change the level density by splitting the $l$ degeneracy of the HO shell $n_{F}$ at the Fermi surface.

We will now calculate analytically the level spectra from perturbative periodic orbit theory for the quartic potential and subsequently within semiclassical WKB wavefunctions for both the quartic and the mean field potential of Eq. (2). We will start with repulsive interactions where pairing is not present.

In periodic orbit theory BB , the level density can be written (to leading order in $\hbar^{-1}$) in terms of a perturbative HO trace formula Brack ; Creagh

[TABLE]

For the unperturbed HO ($U=0$) the modulation factor is ${\cal M}=1$. For a quartic perturbed potential, as in Eq. (6), the modulation factor was calculated in Brack

[TABLE]

with $\sigma=\varepsilon\pi E^{2}/\hbar^{2}\omega^{3}$, being a small classical action. The two terms arise from the change in actions for the circle and diameter orbits respectively due to the quartic potential Brack . The resulting level density can be written in the factorised form Ogren

[TABLE]

[TABLE]

Here, the first term is the average level density, the cosine factor gives the rapid HO shell oscillations (modified by the perturbation) which, however, are slowly modulated by the sine factor resulting in a beating pattern. Moreover, the non-perturbed HO limit, equivalent to ${\cal M}=1$ in Eq. (7), is recovered in the limit of $|\varepsilon|\rightarrow 0$, where the $U\left(3\right)$ symmetry is restored. The $k=1$ term in Eq. (9) gives the major oscillations in the level density and is shown in Fig. 1 (a). The beating pattern or super-shells is clearly observed. The shell oscillations vanish when the argument of the sine in Eq. (9) is an integer $S=1,2,3,...$ times $\pi$, i.e. $|\varepsilon|E^{2}/2(\hbar\omega)^{3}=S$. This gives the supernode condition

[TABLE]

We now turn to an alternative calculation of the level density with WKB. The splitting of the HO shells degenerate levels $l=n_{F},n_{F}-2,...,mod\left(n_{F},2\right)$ in the shell $n_{F}$ by the mean-field potential can be calculated perturbatively in the dilute limit. An excellent approximation for the radial HO wave function with angular momentum $l$ and $(n_{F}-l)/2$ radial nodes in the HO shell when $n_{F}\gg 1$ is the WKB one HM ; Heiselberg :

[TABLE]

between turning points $r^{2}_{\pm}=\tilde{n}_{F}\pm\sqrt{\tilde{n}_{F}^{2}-l(l+1)}$. Here, $\tilde{n}_{F}=n_{F}+3/2$ and the WKB wave number $k_{l}(r)$ is

[TABLE]

When $n_{F}\gg 1$ the wave function has many nodes $1\ll l\ll n_{F}$ and the oscillations in ${\cal R}_{nl}^{2}(r)$ can be averaged $\langle{\sin^{2}(k_{l}(r)r)}\rangle=1/2$ HM . The phase $\theta$ is then unimportant. The single-particle energies for the anharmonic potential of Eq. (6) are simply

[TABLE]

It is special for the quartic perturbation that the level energies are linear in $l(l+1)$. The resulting level spacing increases as $(2l+1)$ just as the level degeneracy for $SO(3)$ symmetry. Therefore the level density is constant within the bandwidth

[TABLE]

on energy scales larger than $2D/n_{F}$ but smaller than $D$. The level density vanishes between the bandwidths of two neighbouring $n$ shells and therefore it generally has a strong oscillatory behavior as shown in Fig. 1 (a). Its amplitude is largest when $D\sim\hbar\omega/2$. However, when $D\simeq\hbar\omega$ the level density is constant and the oscillatory behavior vanishes. This phenomenon repeats when $D=S\hbar\omega$ since the level spectra then overlap $S$ times. With the bandwidth of Eq. (15) under this condition, we obtain exactly the same supernode condition as for periodic orbit theory, Eq. (10). We conclude that Craig’s perturbative periodic orbit theory Creagh is in exact agreement with perturbative WKB for a quartically perturbed spherical symmetric HO in three dimensions.

We now turn to the slightly more complicated mean field potential of Eq. (2). Its level spectrum can also be calculated from the WKB wave functions of Eq. (11). Inserting them in Eq. (13), we obtain

[TABLE]

Here, the integral $I$ is

[TABLE]

where $x=(r^{2}-\tilde{n}_{F})/\sqrt{\tilde{n}_{F}^{2}-l(l+1)}$. This integral is $I=\pi$ for $l\simeq n_{F}$ and $I=8\sqrt{2}/3$ for $l=0$. The bandwidth is therefore

[TABLE]

Inserting this bandwidth in the supernode condition $D=S\hbar\omega$ gives

[TABLE]

For example in the case $2\pi a=1$ the supernodes $S=1,2,3,..$ should occur when $n_{F}\sim 28,44,58$, etc. The Hartree-Fock (HF) calculations of the oscillating part of the total energy, which is proportional to the level density at the Fermi level BB , result in slightly higher supernodes, as in Fig 1 (b). The differences arise because the WKB calculations are perturbative in the interaction strength, whereas in the HF calculation the MF potential $U$ includes a large scattering length which, e.g., leads to corrections for the effective oscillator frequency. Also for the purely quartic term the perturbative approach underestimates the exact supernodes (see Fig. 3 of Brack ). For weaker interactions $2\pi a=0.1$, the first supernode $S=1$ should occur at $n_{F}=130$ according to the condition of Eq. (19), in closer agreement with the HF result of Fig. 1 (c).

For comparison, the Taylor expansion of the mean field potential leads to the supernode condition of Eq. (10) with $\varepsilon=(3\pi\hbar^{2}a/4m)\rho_{0}/R_{TF}^{4}$. It differs from Eq. (19) by the prefactor, which is $\sim 34$% smaller. It is a better approximation to expand e.g. around $r=R_{TF}/2\sqrt{2}=\sqrt{n_{F}}/2$, where the corresponding prefactor is only $\sim 8$% smaller, such that the supernode in Fig 1 (c) is predicted to $n_{F}=137$. Now expanding $I$ of Eq. (17) for small $l\ll n_{F}$, one finds

[TABLE]

resulting in the level spectrum HM

[TABLE]

This level density is constant at low $l$ as for the potential in Eq. (14). However, near $l\sim n_{F}$ the density of levels is slightly smaller as can be seen from the bandwidth corresponding to Eq. (21), which is $\sim 12$% larger, for a given $n_{F}$, than the bandwidth of Eq. (18). That the level density is not completely constant within the bandwidth has the effect that a small periodicity remains even at the super-shell condition $D=S\hbar\omega$. Therefore the shell oscillations do not disappear completely at the supernodes, as can be seen in Fig. 1 (b,c), whereas for the purely quartic case (a) the oscillations disappear completely at the supernodes.

Most atomic traps are not spherical but cigar shaped (prolate) with \omega_{z}\raisebox{-2.15277pt}{\stackrel{{\scriptstyle<}}{{\sim}}}\omega_{\perp}. The unperturbed HO energies $E=n_{z}\hbar\omega_{z}+n_{\perp}\hbar\omega_{\perp}$ will generally lead to a constant level density for energy scales larger than $\hbar\omega_{z}$ but smaller than $\hbar\omega_{\perp}$. When the oscillator frequency ratio $\omega_{\perp}/\omega_{z}$ is a rational number, level degeneracies and larger oscillations will occur on the scale $\hbar\omega_{z}$. Interactions will, however, smear this level density. In any case, super-shell structure is not expected as in the spherical symmetric case. In very oblate traps $\omega_{z}\gg\omega_{\perp}$ the mean field potential is effectively two-dimensional and quadratic, i.e. it does not split the HO shells HM ; Zyl . Thus we may expect strong oscillations in the level density on the scale $\hbar\omega_{\perp}$, but again no super-shell structure.

Attractive interactions lead to pairing by an amount that is exponentially sensitive to the underlying level density near the Fermi surface HM ; BM ; BH ; Heiselberg . The level density is the same for repulsive and attractive interactions except that the levels are reversed when the sign of $\varepsilon$ (Eqs. (9)) and (13)) or $a$ is changed (Eq. (16)). Therefore we can use the level densities and bandwidths calculated above for pairing calculations. Pairing in finite systems is described by the Bogoliubov-de Gennes (BdG) equations deGennes and take place between time-reversed states. As shown in BH these states can be approximated by HO wave functions in dilute HO traps as long as the gap does not exceed the oscillator energy, \Delta\raisebox{-2.15277pt}{\stackrel{{\scriptstyle<}}{{\sim}}}\hbar\omega. Solving BdG for such finite systems is numerically complicated and we shall therefore apply further simplifying approximations, namely that the pairing gap $\Delta_{nl}$ and the wavefunction overlap matrix elements vary slowly with level $l$ in a shell $n$. Both approximations are fair for the trapped atoms as argued in Heiselberg and deviations can be understood. As result we arrive at a much simplified gap equation

[TABLE]

Here, the supergap $G=32\sqrt{2n_{F}}|a|\hbar\omega/15\pi^{2}$ was calculated in HM as the pairing gap when all states in a shell can pair; this is the case for a region of interaction strengths and particle number where the gap is large as compared to the level splitting, yet small compared to the shell splitting $\hbar\omega$. $\Delta(\mu)=\Delta_{nl}$ is the gap at the Fermi surface. $g(E)=n_{F}^{2}/D$ is the level density within each bandgap $D$ around every shell n=0,1,...,\sim$$2n_{F} but vanishes between the bandgaps. The gap equation thus reduces to $1=(G/D)\sum_{n=0}^{\sim 2n_{F}}\int_{0}^{D}dE/\sqrt{(E+n\hbar\omega-\mu)^{2}+\Delta^{2}}$. The chemical potential $\mu$ can be determined from the level spectrum; as we gradually fill particles into the shell $n_{F}$ at the Fermi surface, $\mu$ increases from $n_{F}\hbar\omega$ to $n_{F}\hbar\omega+D$. The cut-off n\raisebox{-2.15277pt}{\stackrel{{\scriptstyle<}}{{\sim}}}2n_{F} in the sum of the gap equation models as a first approximation the more rigorous regularization procedure described in Ref. BruunBCS that is required for a delta-function pseudo-potential.

By solving this simplified gap equation of Eq. (22), we find that it still contains and displays the essential interplay between the variation in level density and pairing. To illustrate the super-shell structure in pairing, we take the strongly anharmonic trap potential used for the level spectra in Fig. 1 (a), and calculate the pairing arising from a weak attractive scattering length $a<0$. For sufficiently weak interactions such that pairing only takes place in the shell at the Fermi surface, we obtain the expected result from the gap equation: $\Delta=G$ when $D\ll\Delta$, whereas for $D\gg\Delta$ we get $\Delta=D\exp(-D/2G)$ midshell ($\mu=n_{F}\hbar\omega+D/2$) and $\Delta=2D\exp(-D/G)$ endshell ($\mu=n_{F}\hbar\omega$ or $\mu=n_{F}\hbar\omega+D$). Pairing is thus stronger at midshell than at endshell, where there are fewer states to pair Heiselberg , and strong shell oscillations follow as shown in Fig. 2. For stronger interactions, pairing also takes place between states in shells around the Fermi shell and Eq. (22) gives: $\Delta=G/\left(1-2G\ln\left(n_{F}\right)/\hbar\omega\right)$ for small bandwidth BH . In Fig. 2 this curve is compared with the finite bandwidth result, which has strong oscillations except at the supernodes where the level density is continuous. At a supernode $D=\hbar\omega$ and the gap equation (22) leads to a gap $\Delta=2n_{F}\hbar\omega\exp(-\hbar\omega/2G)$ Heiselberg .

In summary, level densities, shell-oscillations and super-shell structures in anharmonic traps calculated from numerical Hartree-Fock and analytical periodic orbit theory as well as WKB were found to match to leading order. Analogous super-shell structures were found in pairing from an approximated BdG calculation. The mean field in atomic nuclei also have a large anharmonic potential and the HO shells start to overlap (the first supernode) already for heavy nuclei with $n_{F}\sim 5-6$. The interplay of level spectra and multishell pairing is, however, difficult to disentangle in nuclear pairing due to strong spin-orbit effect and small particle number. Ultracold atomic traps, however, provide ideal systems for observing the rich quantum structures such as level densities and pairing.

Discussions with Matthias Brack on periodic orbit theory, Ben Mottelson on (nuclear) shell theory and pairing, and proof reading by Joel Corney, are gratefully acknowledged.