View source for Attractively interacting Fermi gases - Pairing Instability

=== BCS superfluidity ===
Superfluidity of boson was first discovered in <math>^4 He</math> system at a critical temperature of <math>T_C \sim 2.2K</math>. This was connected to the formation of <math>^4 He</math> condensates. Superfluidity of fermions, the electrons, was first discovered in Mecury at a transition temperature <math>T_C \sim 4.2K</math>, which is known as the `superconductivity' of metals. 

In the early age, there are two major confusions about the fermionic superfluidity
* what is the mechanism for superfluidity of fermions (electrons)?
**It is intuitive to suggest that two electrons could form tightly bounded pairs (Schafroth pairs) and then form condensates. However, there was no known interaction which is strong enough to overcome the ''Coulomb'' repulsion.
* why does it happen at such low temperature compared with <math>T_F </math> (typically <math>\sim 10^4 K</math> in metal)?
**For bosonic case in <math>^4 He</math>, we can estimate the transition temperature <math>T_{BEC} </math> (assuming phase space density 1 and typical Helium density) to be <math>T_{BEC} \sim 3K</math> which is consistent with the experimental findings. However, the fermi temperature in a fermionic system in Mercury is much higher (10^4) than the observed superfluidity transition temperature. 

The two puzzles remain unresolved until 1956 when Bardeen, Cooper and Schrieffer proposed the BCS theory. In short:
* It is correct to think of fermion (electron) pairs. However, instead of the tightly bound pairs, the pair here is the loosely bound BCS pair of electrons formed due to the effective attractive interaction mediated by the hosting lattice.
* The temperature scale is the ''Debye temperature'' because of the involvement of the hosting lattice in the pairing mechanism. This temperature is further modified by the pairing energy and the density of states on the Fermi sea.
=== Cooper instability in a Fermi gas with attractive interactions ===
In contrast to bosons, the non-interacting Fermi gas does not show any phase transition down to zero temperature. One might assume that this qualitative fact should not change as interactions are introduced, at least as long as they are weak. This is essentially true in the case of repulsive interactions. For attractive interactions, the situation is, however, dramatically different. Even for very weak attraction, the fermions form pairs and become superfluid, due to a generalized form of pair condensation.

The idea of pairing might be natural, as tightly bound pairs of fermions can be regarded as point-like bosons, which should form a Bose-Einstein condensate. However, for weak attractive interaction -- as is the case for the residual, phonon-induced electron-electron interaction in metals -- it is not evident that a paired state exists. Indeed, we will see in the following that in three dimensions there is no bound state for two isolated particles and arbitrarily weak interaction. However, by discussing exact solutions in 1D and 2D, where bound states exist for weak interactions, we gain insight into how a modified density of states will lead to bound states even in 3D --  this is the famous Cooper instability.

==== Two-body bound states in 1D, 2D and 3D ====
Localizing a quantum-mechanical particle of mass <math>\mu = m/2</math> to a certain range <math>R</math> leads to an increased momentum uncertainty of <math>p \sim \hbar/R</math> at a kinetic energy cost of about <math>E_{R} = p^2/m = \hbar^2 / mR^2</math>. Clearly, a shallow potential well of size <math>R</math> and depth <math>V</math> with <math>V/E_R \equiv \epsilon \ll 1</math> cannot confine the particle within its borders. But we can search for a bound state at energy <math>|E_B| \ll E_R</math> of much larger size <math>r_B = 1/\kappa \equiv\sqrt{\hbar^2/m |E_B|}  \gg R</math>.
* In 1D and 2D, there is always a bound state for an arbitrarily weak attractive interaction.
* In 3D: For a spherically symmetric well, the Schr\"odinger equation for the wave function transforms into an effective one-dimensional problem for the wave function <math>u = r \psi</math> c). We might now be tempted to think that there must always be a bound state in 3D, asthe  case in 1D. However, the boundary condition on <math>u(r)</math> is now to vanish linearly at <math>r=0</math>, in order for <math>\psi(0)</math> to be finite. Outside of the potential well, we still have <math>u \propto e^{-\kappa r}</math> for a bound state. Inside the well the wave function must fall off to zero at $r=0$ and necessarily has to change its slope from <math>-\kappa</math> outside to <math>\sim 1/R</math> inside the well over a distance <math>R</math>. This costs the large kinetic energy <math>\sim\hbar^2 u''/2m u \approx \hbar^2 /m R^2 = E_R</math>. If the well depth is smaller than a critical depth <math>V_c</math> on the order of <math>E_R</math>, the particle cannot be bound. 

The results hold for quite general shapes <math>V(r)</math> of the (purely attractive) potential well.

Applying these results to the equivalent problem of two interacting particles colliding in their center-of-mass frame, we see that in 1D and 2D, two isolated particles can bind for an arbitrarily weak purely attractive interaction. Hence in 1D and 2D, pairing of fermions can be understood already at the two-particle level. Indeed, one can show that the existence of a two-body bound state for isolated particles in 2D is a necessary and sufficient condition for the instability of the many-body Fermi sea (Cooper instability, see below)~\cite{rand89bound}. In 3D, however, there is a threshold interaction below which two isolated particles are unbound. We conclude that if pairing and condensation occur for arbitrarily weak interactions in 3D, then this must entirely be due to many-body effects.

==== Density of states ====
We show that the "density of states" in the different dimensions as the decisive factor for the existence of bound states.

Searching for a shallow bound state of energy <math>E = -\frac{\hbar^2\kappa^2}{m}</math> (<math>m/2</math> is the reduced mass), we start by writing the Schrodinger equation for the relative wave function momentum space (<math>n</math>-dimensional):
:<math>
    \psi_\kappa(\vec{q}) = - \frac{m}{\hbar^2}\frac{1}{q^2 + \kappa^2} \int \frac{d^n q'}{(2\pi)^n} V(\vec{q}-\vec{q}') \psi_\kappa(\vec{q}')
\,.
</math>
For a short-range potential, <math>V(\vec{q})</math> is practically constant <math>V_0</math> for all relevant <math>q</math>, and falls off to zero on a large <math>q</math>-scale of </amth>\approx 1/R</math>.  Thus
:<math>
    \psi_\kappa(\vec{q}) \approx  - \frac{mV_0}{\hbar^2}\frac{1}{q^2 + \kappa^2} \int_{q'\lesssim\frac{1}{R}} \frac{d^n q'}{(2\pi)^n} \psi_\kappa(\vec{q}')
</math>
We integrate once more over <math>\vec{q}</math>, applying the same cut-off <math>1/R</math>, and then divide by the common factor <math>\int_{q\lesssim\frac{1}{R}}\frac{d^n q}{(2\pi)^n} \psi_\kappa(\vec{q})</math>. We obtain the following equation for the bound state energy <math>E</math>:
:<math>
   - \frac{1}{V_0} \;=\; \frac{m}{\hbar^2}\int_{q\lesssim\frac{1}{R}} \frac{d^n q}{(2\pi)^n} \frac{1}{q^2 + \kappa^2}\; =\; \frac{1}{\Omega}\int_{\epsilon<E_R} d\epsilon \frac{\rho_n(\epsilon)}{2\epsilon+\left|E\right|}
</math>
with the density of states in <math>\rho_n(\epsilon)</math>,the energy cut-off <math>E_R = \hbar^2/m R^2</math> and the volume <math>\Omega</math> of the system. The question on the existence of bound states for arbitrarily weak interaction has now been reformulated: As <math>|V_0| \rightarrow 0</math>, the left hand side diverges. This equation has a solution for small <math>|V_0|</math> only if the right hand side also diverges for vanishing bound state energy <math>|E| \rightarrow 0</math>, and this involves an integral over the density of states.  In 1D and 2D, there is always bound states. However, in 3D the integral is finite for vanishing <math>|E|</math>, and there is a threshold for the interaction potential to bind the two particles.

These results give us an idea why there might be a paired state for two fermions immersed in a 3D medium, even for arbitrarily weak interactions: It could be that the density of available states to the two fermions is altered due to the presence of the other atoms.

==== Pairing of fermions -- The Cooper problem ====
Consider now two weakly interacting spin 1/2 fermions not in vacuum, but on top of a (non-interacting) filled Fermi sea, the Cooper problem.  For weak
interactions, the particles' momenta are essentially confined to a narrow shell above the Fermi surface. The density of states at the Fermi surface is <math>\rho_{\rm 3D}(E_F)</math>, which is a constant just like in two dimensions. We should thus find a "bound state" for the two-particle system for
"arbitrarily weak attractive interaction".

In principle, the two fermions could form a pair at any finite momentum. However, considering the discussion in the previous section, the largest binding energy can be expected for the pairs with the largest density of scattering states. For zero-momentum pairs, the entire Fermi surface is available for
scattering. If the pairs have finite center-of-mass momentum <math>\vec{q}</math>, the number of contributing states is strongly reduced, as they are confined to a circle. Consequently, pairs at rest experience the strongest binding. In the following, we will calculate this energy.

We can write the Schr\"odinger equation for the two interacting particles as before, but now we need to search for a small binding
energy <math>E_B = E-2E_F<0</math> on top of the large Fermi energy of the two particles. The equation for <math>E_B</math> is
:<math>
    -\frac{1}{V_0} = \frac{1}{\Omega}\int_{E_F<\epsilon<E_F+E_R} d\epsilon \frac{\rho_{\rm 3D}(\epsilon)}{2(\epsilon-E_F)+\left|E_B\right|}
</math>
The effect of Pauli blocking of momentum states below the Fermi surface is explicitly included by only integrating over energies
<math>\epsilon > E_F</math>.

In conventional superconductors, the natural cut-off energy <math>E_R</math> is given by the Debye frequency <math>E_R = \hbar \omega_D</math>, corresponding to the highest frequency at which ions in the crystal lattice can respond to a bypassing electron. Since we have
<math>\hbar \omega_D \ll E_F</math>, we can approximate <math>\rho_{\rm 3D}(\epsilon) \approx \rho_{\rm 3D}(E_F)</math> and find:
:<math>
    E_B = - 2 \hbar \omega_D e^{-2 \Omega/\rho_{\rm 3D}(E_F) \left|V_0\right|}\\
</math>

In the case of an atomic Fermi gas, we should replace $1/V_0$ by the physically relevant scattering length $a < 0$. The equation for the
bound state becomes
:<math>
    -\frac{m}{4\pi\hbar^2 a} = \frac{1}{\Omega}\int_{E_F}^{E_F+E_R} d\epsilon \frac{\rho_{\rm 3D}(\epsilon)}{2(\epsilon-E_F)+\left|E_B\right|}
- \frac{1}{\Omega}\int_0^{E_F+E_R} d\epsilon \frac{\rho_{\rm
3D}(\epsilon)}{2\epsilon}
</math>
The right hand expression is now finite as we let the cut-off <math>E_R \rightarrow \infty</math>, the result being (one assumes <math>\left|E_B\right|\ll E_F</math>)
:<math>
    -\frac{m}{4\pi\hbar^2 a} = \frac{\rho_{\rm 3D}(E_F)}{2\Omega}
\left(-\log\left(\frac{\left|E_B\right|}{8E_F}\right) - 2\right)
</math>
Inserting <math>\rho_{\rm 3D}(E_F) = \frac{\Omega m k_F}{2\pi^2\hbar^2}</math> with the Fermi wave vector <math>k_F = \sqrt{2mE_F/\hbar^2}</math>, one arrives at
:<math>
    E_B = - \frac{8}{e^2} E_F\, e^{-\pi/k_F \left|a\right|}
\,.
</math>
The binding energies can be compared with the result for the bound state of two particles in 2D. The role of the constant density of states <math>\rho_{\rm 2D}</math> is here played by the 3D density of states at the Fermi surface, <math>\rho_{\rm 3D}(E_F)</math>.

We, therefore, see that two weakly interacting fermions on top of a Fermi sea form a bound state due to Pauli blocking. However, in this artificial problem we neglected the interactions between particles "inside" the Fermi sea.  As we ``switch on'' the interactions for all particles from top to the bottom of the Fermi sea, the preceding discussion indicates that the gas will reorder itself into a completely new, paired state. The Fermi sea is thus unstable towards pairing (Cooper instability). The full many-body description of such a paired state, including the necessary anti-symmetrization of the
full wave function, was achieved by Bardeen, Cooper and Schrieffer (BCS) in 1957. As we will see in the next section, the self-consistent inclusion of all fermion pairs leads to more available momentum space for pairing. The effective density of states is then twice as large, giving a superfluid gap <math>\Delta</math> that differs from above result by a factor of 2 in the exponent.

It should be noted that the crucial difference to the situation of two particles in vacuum in 3D is the "constant density" of states at the Fermi energy (and not the 2D character of the Fermi surface). Therefore, if we were to consider the Cooper problem in higher dimensions and have two
fermions scatter on the <math>(n-1)</math> dimensional Fermi surface, the result would be similar to the 2D case (due to the constant density of states), and not to the case of <math>(n-1)</math> dimensions.

The conclusion of this section is that Cooper pairing is a many-body phenomenon, but the binding of two fermions can still be understood by two-body quantum mechanics, as it is similar to two isolated particles in two dimensions.  To first order, the many-body physics is not the modification of interactions, but rather the modification of the density of states due to Pauli blocking.

=== Many-body Hamiltonian with Pairing ===
It turns out that the transition from tightly bound molecules to cooper pairs can be described by the same wavefunction and there is ''no'' phase transition in between. Therefore, it is called Crossover. In this section, we discussed two approaches which describe the crossover physics.
*we start from an Ansatz wavefunction and minimize the free energy
*we start with a mean field theory similar to the weakly interaction boson case.
==== Crossover wave function ====
In 3D, two fermions in isolation can form a molecule for strong enough attractive interaction. When the size of the molecule is smaller than the interatomic distance of its constituents in the fermi sea <math>\sim 1/k_F</math>, the molecules act as bosons and the ground state of the system should be a Bose-Einstein condensate of these tightly bound pairs. This is equivalent to say that the binding energy <math>E_B > E_F</math>.

For too weak an attraction there is no bound state for two isolated fermions, but Cooper pairs can form in the medium as discussed above. The ground state of the system turns out to be a condensate of Cooper pairs as described by BCS theory. In contrast to the physics of molecular condensates, however, <math>E_B\ll E_F</math> and therefore Pauli pressure plays a major role.

The crossover from the BCS- to the BEC-regime is smooth. This is somewhat surprising since the two-body physics shows a threshold behavior at a critical interaction strength, below which there is no bound state for two particles. In the presence of the Fermi sea, however, we simply cross over from a regime of tightly bound molecules to a regime where the pairs are of much larger size than the interparticle spacing. 

Rather than the interaction strength, we take the scattering length as the parameter that ``tunes'' the interaction.  For positive $a>0$, there is a two-body bound state available at 
:<math>E_B
= -\hbar^2/m a^2
</math>
while small and negative $a<0$ corresponds to weak attraction where Cooper pairs can form in the medium. In either case, for <math>s</math>-wave interactions, two spins form a singlet and the orbital part of the pair wave function <math>\varphi(\vec{r}_1,\vec{r}_2)</math> will be symmetric under exchange,and in a uniform system, will only depend on their distance <math>\left|\vec{r}_1-\vec{r}_2\right|</math>. We will explore the many-body wave function
:<math>
    \Psi\left(\vec{r}_1,\dots,\vec{r}_N\right) = \mathcal{A}\left\{\varphi(\left|\vec{r}_1-\vec{r}_2\right|)\chi_{12}\dots \varphi(\left|\vec{r}_{N-1}-\vec{r}_N\right|)\chi_{N-1,N}\right\}
</math>
that describes a condensate of such fermion pairs, with the operator <math>\mathcal{A}</math> denoting the correct antisymmetrization of
all fermion coordinates, and the spin singlet <math>\chi_{ij} = \frac{1}{\sqrt{2}}(\left|\uparrow\right>_i \left|\downarrow\right>_j - \left|\downarrow\right>_i \left|\uparrow\right>_j)</math>. In the experiment, ``spin up'' and ``spin down'' will correspond to two atomic hyperfine states.

In second quantization notation, we write
:<math>
    \left|\Psi\right>_N = \int \prod_i d^3 r_i \,\varphi(\vec{r}_1 - \vec{r}_2) \Psi_\uparrow^\dagger(\vec{r}_1) \Psi_\downarrow^\dagger(\vec{r}_2) \dots \varphi(\vec{r}_{N-1} - \vec{r}_N) \Psi_\uparrow^\dagger(\vec{r}_{N-1}) \Psi_\downarrow^\dagger(\vec{r}_N) \left|0\right>
</math>
where the fields <math>\Psi_\sigma^\dagger(\vec{r}) = \sum_k c_{k\sigma}^\dagger \frac{e^{-i \vec{k} \cdot \vec{r}}}{\sqrt{\Omega}}</math>. With the
Fourier transform <math>\varphi(\vec{r}_1-\vec{r}_2) = \sum_k \varphi_k \frac{e^{i \vec{k} \cdot (\vec{r}_1- \vec{r}_2)}}{\sqrt{\Omega}}</math> we can introduce the pair creation operator
:<math>
    b^\dagger = \sum_k \varphi_k c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger
</math> and find
:<math>
    \left|\Psi\right>_N = {b^\dagger}^{N/2} \left|0\right>\,.
</math>
This expression for  <math>\left|\Psi\right>_N</math> is formally identical to the Gross-Pitaevskii ground state of a condensate of bosonic particles. However, the operators <math>b^\dagger</math>  obey bosonic commutation relations only in the limit of tightly bound pairs.  For the commutators, we obtain
:<math>
\begin{eqnarray}
    \left[b^\dagger,b^\dagger\right]_- &= \sum_{k k'} \varphi_k \varphi_{k'} \left[ c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger,c_{k'\uparrow}^\dagger c_{-k'\downarrow}^\dagger \right]_- &= 0  \\
    \left[b,b\right]_- &= \sum_{k k'} \varphi^*_k \varphi^*_{k'} \left[c_{-k\downarrow} c_{k\uparrow},c_{-k'\downarrow} c_{k'\uparrow}\right]_- &= 0 \nonumber\\
\left[b,b^\dagger\right]_- &= \sum_{k k'} \varphi^*_k \varphi_{k'}
\left[c_{-k\downarrow} c_{k\uparrow},c_{k'\uparrow}^\dagger
c_{-k'\downarrow}^\dagger \right]_- &= \sum_k |\varphi_k|^2 (1 -
n_{k\uparrow} - n_{k\downarrow}) \nonumber
\end{eqnarray}
</math>
The third commutator is equal to one only in the limit where the pairs are tightly bound and occupy a wide region in momentum space.  In this case, the occupation numbers <math>n_k</math> of any momentum state <math>k</math> are very small and <math>\left[b,b^\dagger\right]_- \approx \sum_k |\varphi_k|^2 = \int {\rm d}^3 r_1\int {\rm d}^3 r_2 \,|\varphi(\vec{r}_1,\vec{r}_2)|^2 = 1</math>. (or is it better to say when <math>|\varphi_k|</math> spread out to a rigion much wider than <math>|k_F|</math> in the momentum space?)

It is easier to use the grand canonical formalism, not fixing the number of atoms but the chemical potential <math>\mu</math>. A separate, crucial step is to define a many-body state which is a superposition of states with different atom numbers.  In the BEC limit, this is analogous to the use of coherent states (vs. Fock states) in quantum optics. Let <math>N_p = N/2</math> be the number of pairs. Then,
:<math>
\begin{eqnarray}
    \label{e:coherentstate}
    \mathcal{N}\left|\Psi\right> &= \sum_{J_{\rm even}} \frac{N_p^{J/4}}{(J/2)!} \left|\Psi\right>_J &= \sum_M \frac{1}{M!} {N_p^{M/2}\; b^\dagger}^M \left|0\right> = e^{\sqrt{N_p} \;b^\dagger} \left|0\right>\nonumber \\
&= \prod_k e^{\sqrt{N_p}\;  \varphi_k c_{k\uparrow}^\dagger
c_{-k\downarrow}^\dagger} \left|0\right> &= \prod_k (1 + \sqrt{N_p}\;
\varphi_k\, c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger)
\left|0\right>
\end{eqnarray}
</math>
The second to last equation follows since the operators $c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger$ commute for different $k$, and the last equation follows from $c_k^{\dagger 2} = 0$. After normalization we have,
:<math>
\begin{equation}
    \left|\Psi_{\rm BCS}\right> = \prod_k (u_k + v_k c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger) \left|0\right>
    \label{e:BCSstate}
\end{equation}
</math>
with <math>v_k = \sqrt{N_p}\,\varphi_k u_k</math> and <math>|u_k|^2 + |v_k|^2 = 1</math>. 

This is the famous BCS wave function which is the exact solution of the simplified Hamiltonian. It is a product of wave functions referring to the occupation of pairs of single-particle momentum states, <math>(\vec{k},\uparrow,-\vec{k},\downarrow)</math>.

As a special case, it describes a non-interacting Fermi sea, with all momentum pairs occupied up to the Fermi momentum (<math>u_k=0,v_k=1</math> for <math>k<k_F</math> and <math>u_k=1, v_k=0</math> for <math>k>k_F</math>). In general, arbitrary <math>v_k, u_k</math> describe a ``molten'' Fermi sea, modified by the coherent scattering of pairs with zero total momentum. Pairs of momentum states are seen to be in a superposition of being fully empty and fully occupied.

It is worth noticing that we didn't make any further assumption on the component <math>\varphi_k</math>, therefore this wave function encompasses the entire regime of pairing, from point bosons (small molecules) to weakly and non-interacting fermions.

==== Gap and number equation ====

The variational parameters <math>u_k, v_k</math> are derived in the standard way by minimizing the free energy <math>E - \mu N =
\left<\hat{H} - \mu \hat{N}\right></math>. The many-body Hamiltonian for the system is
:<math>
    \hat{H} = \sum_{k,\sigma} \epsilon_k c_{k\sigma}^\dagger c_{k\sigma} + \frac{V_0}{\Omega}\sum_{k,k',q} c_{k+\frac{q}{2} \uparrow}^\dagger c_{-k+\frac{q}{2}\downarrow}^\dagger c_{k'+\frac{q}{2}\downarrow} c_{-k'+\frac{q}{2}\uparrow}
</math>

The dominant role in superfluidity is played by fermion pairs with zero total momentum (Cooper pairs with zero momentum have the largest binding energy.) Therefore, we simplify the mathematical description by neglecting interactions between pairs at finite momentum. This is a very drastic simplification, as hereby density fluctuations are eliminated. For neutral superfluids, sound waves  (the Bogoliubov-Anderson mode, ) are eliminated by this approximation.
The approximate Hamiltonian (``BCS Hamiltonian``) reads
:<math>
    \hat{H} = \sum_{k,\sigma} \epsilon_k c_{k\sigma}^\dagger c_{k\sigma} + \frac{V_0}{\Omega}\sum_{k,k'} c_{k \uparrow}^\dagger c_{-k\downarrow}^\dagger c_{k'\downarrow} c_{-k'\uparrow}
</math>
Minimizing <math>E-\mu N</math> leads to
:<math>
\begin{eqnarray}
    v_k^2 &=& \frac{1}{2}\left(1 - \frac{\xi_k}{E_k}\right) \nonumber \\
    u_k^2 &=& \frac{1}{2}\left(1 + \frac{\xi_k}{E_k}\right) \nonumber \\
    \mbox{with }\;E_k &=& \sqrt{\xi_k^2 + \Delta^2}
\end{eqnarray}
</math>
where <math>\Delta</math> is given by the "gap equation"
:<math>
\Delta \equiv \frac{V_0}{\Omega} \sum_k \left<c_{k\uparrow}
c_{-k\downarrow}\right> = - \frac{V_0}{\Omega} \sum_k u_k v_k = -
\frac{V_0}{\Omega} \sum_k \frac{\Delta}{2 E_k}
</math> or
:<math>
    -\frac{1}{V_0} = \int \frac{d^3 k}{\left(2\pi\right)^3} \;\frac{1}{2E_{k}}
</math>
Note the similarity to the bound state equation in free space and in the simplified Cooperproblem, Eq.~\ref{e:Cooper}. An additional constraint is given by the "number equation" for the total particle density <math>n = N / \Omega</math>
:<math>
n = 2 \int \frac{d^3 k}{\left(2\pi\right)^3} \; v_k^2
</math>
Gap and number equations have to be solved simultaneously to yield the two unknowns <math>\mu, \Delta</math>. 
With <math>V_0</math> replaced by scattering length, and using  <math>E_F, k_F</math> as the relevant scales, the gap and number equation give
:<math>
\begin{eqnarray}
    -\frac{1}{k_F a} &=& \frac{2}{\pi} \left(\frac{2}{3 I_2\left(\frac{\mu}{\Delta}\right)}\right)^{1/3} I_1\left(\frac{\mu}{\Delta}\right)\\
    \frac{\Delta}{E_F} &=& \left(\frac{2}{3 I_2\left(\frac{\mu}{\Delta}\right)}\right)^{2/3}
\end{eqnarray}
</math>
The result for <math>\mu</math> and <math>\Delta</math> as a function of <math>1/k_F a</math> is shown in Fig.~\ref{f:Deltamu}. It is possible to obtain analytic expressions for the solutions in terms of complete elliptic integrals

==== Mean field theory and Bogoliubov transformation ====
In this derivation, we have combined the simplified Hamiltonian with the BCS variational Ansatz.  Alternatively one can
apply a decoupling (mean field) approximation to the Hamiltonian similar to the mean-field theory for weakly interacting bosons. Expecting that there will be some form of pair condensate, we assume that the pair creation and annihilation operator only weakly fluctuates around its non-zero expectation value
:<math>
\begin{eqnarray}
    C_k=  \left<c_{k\uparrow} c_{-k\downarrow}\right>= -\left< c_{k\uparrow}^\dagger
    c_{-k\downarrow}^\dagger\right>
\end{eqnarray}
</math>
chosen to be real (since the relative phase of states which differ in particle number by two can be arbitrarily chosen).
That is, we write
:<math>
\begin{equation}
c_{k\uparrow} c_{-k\downarrow} = C_k + (c_{k\uparrow} c_{-k\downarrow} - C_k)
\end{equation}
</math>
with the operator in parentheses giving rise to fluctuations that are small on the scale of <math>C_k</math>.
The gap parameter <math>\Delta</math> is now defined as
:<math>
\Delta=\frac{V_0}{\Omega} \sum_{k} C_k
</math>
We only include terms in the interaction part of the Hamiltonian which
involve the $C_k$'s at least once. That is, we neglect the
correlation of fluctuations of the pair creation and annihilation
operators. One obtains
:<math>
\begin{equation}
    \hat{H} = \sum_{k} \epsilon_k (c_{k\uparrow}^\dagger c_{k\uparrow}+ c_{k\downarrow}^\dagger c_{k\downarrow}) -\Delta \sum_k\left( c_{k \uparrow}^\dagger c_{-k\downarrow}^\dagger + c_{k \downarrow}
    c_{-k\uparrow} + \sum_{k'} C_{k'}\right)
\end{equation}
</math>

This Hamiltonian  is bilinear in the creation and annihilation operators and can easily be solved by a Bogoliubov transformation  to new quasi-particle operators 
<math>\gamma_{k\downarrow}</math>:
:<math>
\begin{eqnarray}
    \gamma_{k\uparrow} &=& u_k c_{k\uparrow} - v_k c_{-k\downarrow}^\dagger \\
    \gamma_{-k\downarrow}^\dagger &=& u_k c_{-k\downarrow}^\dagger + v_k c_{k\uparrow} \nonumber
\end{eqnarray}
</math>
The <math>u_k</math> and <math>v_k</math> are determined from the requirements that the new operators fulfill fermionic commutation relations and that the transformed Hamiltonian is diagonal with respect to the quasiparticle operators. This solution is identical to the one obtained before for the <math>u_k</math> and <math>v_k</math>, and the transformed Hamiltonian becomes
:<math>
\begin{equation}
\hat{H} - \mu \hat{N} = - \frac{\Delta^2}{V_0/\Omega} + \sum_k (\xi_k - E_k) + \sum_k E_k
(\gamma_{k\uparrow}^\dagger \gamma_{k\uparrow} +
\gamma_{k\downarrow}^\dagger \gamma_{k\downarrow})
\end{equation}
</math>
The first two terms give the free energy of the pair condensate when the correct <math>u_k</math> and <math>v_k</math> are inserted. The third term represents the energy of excited quasi-particles, and we identify <math>E_k</math> as excitation energy of a quasi-particle.  The superfluid ground state is the quasi-particle vacuum.

This approach via the pairing field is analogous to the Bogoliubov treatment of an interacting Bose-Einstein condensate: There, the creation and annihilation operators for atoms with zero momentum are replaced by <math>\sqrt{N_0}</math>, the square root of the number of condensed atoms (i.e.~a coherent field). In the interaction term of the Hamiltonian, the Hamiltonian is solved by keeping only certain pair interactions, either by using a variational pairing wave function, or by introducing a mean pairing field. It should be noted that these approximations are not even necessary, as the BCS wave function can be shown to be the {\it exact} solution to the
reduced Hamiltonian.

=== Discussion of the three regimes -- BCS, BEC and crossover ===
==== BCS limit ====
In the BCS-limit of weak attractive interaction, <math>k_F a \rightarrow 0_-</math>, we have
:<math>
\begin{eqnarray}
    \mu &\approx& E_F \\
    \Delta &\approx& \frac{8}{e^2} e^{-\pi/2k_F\left|a\right|}
\end{eqnarray}
</math>
*The first equation tells us that adding a spin up and spin down particle to the system costs a Fermi energy per particle (with the implicit assumption that both a spin up and a spin down particle are added, raising the total energy by <math>2 \mu</math>).  In the weakly interacting BCS limit Pauli blocking still dominates over interactions, and hence the particles can only be added at the Fermi surface. 
*The second equation is the classic result of BCS theory for the superfluid gap. Compared to the bound state energy for a single Cooper pair on top of a non-interacting Fermi sea, the gap is larger (the negative exponent is smaller by a factor of two), as the entire collection of particles now takes part in the pairing. However, the gap is still exponentially small compared to the Fermi energy: Cooper pairing is fragile.

The ground state energy of the BCS state can be calculated from mean-field and is
:<math>
\begin{equation}
    E_{\rm G,\, BCS} = \frac{3}{5} N E_F - \frac{1}{2}\,\rho(E_F)\, \Delta^2
\end{equation}
</math>
The first term is the energy of the non-interacting normal state. The second term is the energy due to condensation, negative as it should be, indicating that the BCS state is energetically favorable compared to the normal state. Though the kinetic energy of the system increases (populating states above the fermi sea), the total energy drops due to the attractive interaction.

The second term can be interpreted in two ways. One way refers to the number of pairs $N/2$ times the energy per pair <math>- \frac{3}{4} \Delta^2/E_F</math>. The other interpretation refers to pairing on the surface only but with a pairing energy <math>\Delta</math>. The number of pairs is the number of states within <math>E_F\pm\Delta</math> which is <math>\sim \rho(E_F) \, \Delta\sim N \Delta/E_F</math>.  The second interpretation justifies the picture of a Cooper pair condensate. In the solution of the Cooper problem, the pair wave function has a peak occupation per momentum state of <math>\sim 1/ \rho(E_F) \Delta</math>. Therefore, one can stack up  <math>\sim  \rho(E_F) \Delta</math> pairs with zero total momentum and construct a Bose-Einstein condensate of Cooper pairs. However, the Cooper pairs are not bosons, as shown by commutator. However, if there were only a few Cooper pairs, much less than <math>\rho(E_F) \Delta</math>, the occupation of momentum states <math>n_k</math> would still be very small compared to 1 and these pairs would be to a good approximation bosons.

It depends on the experiment whether it reveals a pairing energy of <math>\frac{1}{2} \Delta^2/E_F</math> or of <math>\Delta</math>.  In RF spectroscopy, all momentum states can be excited, and the spectrum shows the previous one. Tunnelling experiments in superconductors probe the region close to the Fermi surface, and show a pairing gap of <math>\Delta</math>.

To give a sense of scale, Fermi energies in dilute atomic gases are on the order of a <math>\mu\rm K</math>, corresponding to <math>1/k_F \sim 4\,000\, a_0</math>. In the absence of scattering resonances, a typical scattering length would be about <math>50-100\; a_0</math> (on the order of the van der Waals-range). Therefore, the realization of superfluidity in Fermi gases requires scattering or Feshbach resonances to increase the scattering length, bringing the gas into the strongly interacting regime where <math>k_F \left|a\right| > 1</math>. In this case, the above mean field theory predicts <math>\Delta > 0.22 \;E_F</math>  and this is the regime where current experiments are operating.

==== BEC limit ====
In the BEC limit of tightly bound pairs, for <math>k_F a \rightarrow 0_+</math>, one finds
:<math>
\begin{eqnarray}
    \mu &=& -\frac{\hbar^2}{2 m a^2} + \frac{\pi \hbar^2 a n}{m}\\
    \Delta &\approx& \sqrt{\frac{16}{3\pi}} \frac{E_F}{\sqrt{k_F a}}
\end{eqnarray}
</math>
*The first term for the chemical potential is the binding energy per fermion in a tightly bound molecule. This reflects again the implicit assumption (made by using the wave function in that we always add two fermions of opposite spin at the same time to the system.
*The second term for <math>\mu</math> is a mean field contribution describing the repulsive interaction ''between'' molecules in the gas, indicating a scattering length of <math>2a</math>. Exact calculation for the interaction between four fermions shows <math>a_M = 0.6\, a</math>. 

Here <math>\Delta</math> signifies neither the binding energy of molecules nor does it correspond to a gap in the excitation spectrum. Indeed, in the BEC-regime, as soon as $\mu <0$, there is no longer a gap at non-zero <math>k</math> in the single-fermion excitation spectrum . Instead, we have for the quasi-particle energies <math>E_k = \sqrt{(\epsilon_k - \mu)^2 + \Delta^2} \approx |\mu| + \epsilon_k + \frac{\Delta^2}{2|\mu|}</math>. So only the combination <math>\Delta^2/|\mu|</math> is important. 
:<math>    
\frac{\Delta^2}{2|\mu|} = \frac{8}{3\pi} \frac{E_F^2}{k_F a} \frac{2 m a^2}{\hbar^2} = \frac{4}{3\pi}\frac{\hbar^2}{m} k_F^3 a = \frac{4\pi\hbar^2}{m} n\; a
</math>
which is two times the molecular mean field. In fact, it can be interpreted here as the mean field energy experienced by a single fermion in a gas of molecules.
It might surprise that the simplified Hamiltonian contains interactions between two molecules or between a molecule and a single fermion at all. In fact, a crucial part of the simplification has been to explicitly neglect such three- and four-body interactions. This effective repulsion is due to Pauli blocking where a molecule or a fermion cannot enter the region already taken by another molecule.

We see that the only way interactions between pairs, or between a pair and a single fermion, enter in the simplified description of the BEC-BCS crossover is via the anti-symmetry of the many-body wave function.

==== Evolution from BCS to BEC ====
Our variational approach smoothly interpolates between the two known regimes of a BCS-type superfluid and a BEC of molecules. It is a crossover, which occurs approximately between <math>1/k_F a = -1</math> and +1 and is fully continuous. The occupation of momentum states <math>n_k = v_k^2</math> evolves smoothly from the step-function of a degenerate Fermi gas, broadened over a width <math>\Delta \ll E_F</math> due to pairing, to that of <math>N_p</math> molecules, namely the number of molecules <math>N_p</math> times the probability <math>|\varphi_k|^2</math> to find a molecule with momentum <math>k</math>. 

It is also interesting to follow the evolution of the ``Cooper pair'' wave function, which would be <math>v_k/u_k\sqrt{N_p}</math>. The definition given here is the two-point correlation function. Both definitions for the Cooper pair wave function show a sharp feature, either a peak or an edge at the Fermi surface, of width <math>\sim \delta k</math>, thus giving similar behavior for the real space wave function. both in <math>k</math>-space, where it is given by
<math>\left<\Psi_{BCS}\right|c^\dagger_{k\uparrow} c^\dagger_{-k\downarrow}\left|\Psi_{BCS}\right> = u_k v_k</math>, and in real space, where it is
:<math>
\begin{eqnarray}
\psi(\vec{r}_1,\vec{r}_2) &=&
\left<\Psi_{BCS}\right|\Psi_{\uparrow}^\dagger(\vec{r}_1)\Psi_{\downarrow}^\dagger(\vec{r}_2)\left|\Psi_{BCS}\right>
= \frac{1}{\Omega}\sum_k u_k v_k e^{-i \vec{k} \cdot (\vec{r}_1-\vec{r}_2)} \nonumber \\
&=&\frac{1}{\Omega}\sum_k \frac{\Delta}{2E_k}\, e^{-i \vec{k} \cdot
(\vec{r}_1-\vec{r}_2)} \label{e:cooperpairwavefunction}
\end{eqnarray}
</math>
In the BCS limit, the pairing occurs near the Fermi surface <math>k= k_F</math>. Therefore, the spatial wave function of Cooper pairs has a strong modulation at the inverse wave vector <math>1/k_F</math>, and an overall extent of the inverse width of the pairing region, <math>\sim1/\delta k \sim \frac{\hbar v_F}{\Delta} \gg 1/k_F</math>.  The characteristic size of the Cooper pair, or the ''two-particle correlation length'' <math>\xi_0</math>, can be
defined as <math>\xi_0^2=\frac{\left<\psi(\vec{r})\right|r^2\left|\psi(\vec{r})\right>}{\left<\psi(\vec{r})|\psi(\vec{r})\right>}</math>,
and this gives indeed <math>\xi_0 \sim 1/\delta k</math>,
:<math>
    \xi_0 \approx \xi_{BCS} \equiv \frac{\hbar v_F}{\pi \Delta} \gg 1/k_F \quad \mbox{in the BCS-limit}
</math>

In the BEC limit, <math>u_k v_k \propto \frac{1}{1+(k a)^2}</math>, and so
:<math>
    \psi(\vec{r}_1,\vec{r}_2) \sim \frac{e^{- \left|\vec{r}_1 - \vec{r}_2\right|/a}}{\left|\vec{r}_1 - \vec{r}_2\right|}
</math>
which is simply the wave function of a molecule of size <math>\sim a</math>. The two-particle correlation length  <math>\xi_{phase}</math> that is associated with spatial fluctuations of the order parameter. The two length scales coincide in the BCS-limit, but differ in the BEC-limit, where <math>\xi_{phase}</math> is given by the healing length <math>\propto\frac{1}{\sqrt{n a}}</math>, is thus <math>\xi_0 \sim a</math>.
Figs.~\ref{f:pairwavefunction} and~\ref{f:pairsize} summarize the
evolution of the pair wave function and pair size throughout the
crossover.
 Back to: [[Quantum gases]]