BEC-BCS Crossover
BCS superfluidity
Superfluidity of boson was first discovered in system at a critical temperature of . This was connected to the formation of condensates. Superfluidity of fermions, the electrons, was first discovered in Mecury at a transition temperature , 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 (typically in metal)?
- For bosonic case in , we can estimate the transition temperature (assuming phase space density 1 and typical Helium density) to be 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.
Pairing on the Fermi surface
\section{Theory of the BEC-BCS crossover} \label{c:BECBCStheory}
\subsection{Elastic collisions} \label{s:elasticcollisions}
Due to their diluteness, most properties of systems of ultracold atoms are related to two-body collisions. If we neglect the weak magnetic dipole interaction between the spins, the interatomic interaction is described by a central potential $V(r)$. At large distances from each other, atoms interact with the van der Waals-potential $-C_6/r^6$ as they experience each other's fluctuating electric dipole\footnote{For distances on the order of or larger than the characteristic wavelength of radiation of the atom, $\lambda \gg r_0$, retardation effects change the potential to a $-1/r^7$ law.}. At short distances on the order of a few Bohr radii $a_0$, the two electron clouds strongly repel each other, leading to ``hard-core repulsion. If the spins of the two valence electrons (we are considering alkali atoms) are in a triplet configuration, there is an additional repulsion due to Pauli's exclusion principle. Hence, the triplet potential $V_T(r)$ is shallower than the singlet one $V_S(r)$. The exact inclusion of the interatomic potential in the description of the gas would be extremely complicated. However, the gases we are dealing with are ultracold and ultradilute, which implies that both the de Broglie wavelength $\lambda_{dB}$ and the interparticle distance $n^{-1/3} \sim 5\,000-10\,000\, a_0$ are much larger than the range of the interatomic potential $r_0$ (on the order of the van der Waals length $r_0 \sim \left(\mu C_6 / \hbar^2\right) \sim 50\, a_0$ for \li). As a result, scattering processes never explore the fine details of the short-range scattering potential. The entire collision process can thus be described by a single quantity, the {\it scattering length}.
Since the description of Feshbach resonances and of the BCS-BEC crossover require the concept of the effective range and renormalization of the scattering length, we quickly summarize some important results of scattering theory.
The Schr\"odinger equation for the reduced one-particle problem in the center-of-mass frame of the colliding atoms (with reduced mass $m/2$, distance vector $r$, and initial relative wave vector $\vect{k}$) is \begin{equation}
(\nabla^2 + k^2)\Psi_{\vect{k}}(\vect{r}) = v(r)\Psi_{\vect{k}}(\vect{r}) \quad\mbox{with } k^2 = \frac{m E}{\hbar^2} \quad \mbox{and } v(r) = \frac{m V(r)}{\hbar^2}
\label{e:schrodinger} \end{equation} Far away from the scattering potential, the wave function $\Psi_{\vect{k}}(\vect{r})$ is given by the sum of the incident plane wave $e^{i \vect{k} \cdot \vect{r}}$ and an outgoing scattered wave: \begin{equation}
\Psi_{\vect{k}}(\vect{r}) \sim e^{i \vect{k} \cdot \vect{r}} + f(\vect{k}',\vect{k}) \frac{e^{i k r}}{r}
\label{e:psiasymptotic}
\end{equation} $f(\vect{k}',\vect{k})$ is the scattering amplitude for scattering an incident plane wave with wave vector $\vect{k}$ into the direction $\vect{k}' = k\, \vect{r}/r$ (energy conservation implies $k' = k$).
Since we assume a central potential, the scattered wave must be axially symmetric with respect to the incident wave vector $\vect{k}$, and we can perform the usual expansion into partial waves with angular momentum $l$~\cite{land77qm}. For ultracold collisions, we are interested in describing the scattering process at {\it low momenta} $k \ll 1/r_0$, where $r_0$ is the range of the interatomic potential. In the absence of resonance phenomena for $l \ne 0$, {\it $s$-wave scattering} $\,l=0$ is dominant over all other partial waves (if allowed by the Pauli principle): \begin{equation} f \approx f_s = \frac{1}{2ik}(e^{2i\delta_s}-1) = \frac{1}{k \cot \delta_s - i k} \label{e:scattamp} \end{equation} where $f_s$ and $\delta_s$ are the $s$-wave scattering amplitude and phase shift, resp.~\cite{land77qm}. Time-reversal symmetry implies that $k\cot\delta_s$ is an even function of $k$. For low momenta $k \ll 1/r_0$, we may expand it to order $k^2$: \begin{equation} k \cot \delta_s \approx -\frac{1}{a} + r_{\rm eff} \frac{k^2}{2} \end{equation} which defines the {\it scattering length} \begin{equation} a = -\lim_{k \ll 1/r_0} \frac{\tan \delta_s}{k}, \end{equation} and the effective range $r_{\rm eff}$ of the scattering potential. For example, for a spherical well potential of depth $V \equiv \hbar^2 K^2/m$ and radius $R$, $r_{\rm eff} = R - \frac{1}{K^2 a} - \frac{1}{3} \frac{R^3}{a^2}$, which deviates from the potential range $R$ only for $|a| \lesssim R$ or very shallow wells. For van der Waals potentials, $r_{\rm eff}$ is of order $r_0$~\cite{flam99scatt}. With the help of $a$ and $r_{\rm eff}$, $f$ is written as~\cite{land77qm} \begin{equation} f(k) = \frac{1}{-\frac{1}{a} + r_{\rm eff} \frac{k^2}{2} - ik} \label{e:scattamplitude} \end{equation} In the limit $k|a| \ll 1$ and $|r_{\rm eff}| \lesssim 1/k$, $f$ becomes independent on momentum and equals $-a$. For $k|a| \gg 1$ and $r_{\rm eff} \ll 1/k$, the scattering amplitude is $f = \frac{i}{k}$ and the cross section for atom-atom collisions is $\sigma = \frac{4\pi}{k^2}$. This is the so-called unitarity limit. Such a divergence of $a$ occurs whenever a new bound state is supported by the potential (see section~\ref{s:squarewell}).
\subsection{Pseudo-potentials} \label{s:renormalization} If the de Broglie wavelength $\frac{2\pi}{k}$ of the colliding particles is much larger than the fine details of the interatomic potential, $1/k \gg r_0$, we can create a simpler description by modifying the potential in such a way that it is much easier to manipulate in the calculations, but still reproduces the correct $s$-wave scattering. An obvious candidate for such a ``pseudo-potential is a delta-potential $\delta(\vect{r})$.
However, there is a subtlety involved which we will address in the following. The goal is to find an expression for the scattering amplitude $f$ in terms of the potential $V(r) = \frac{\hbar^2 v(r)}{m}$, so that we can try out different pseudo-potentials, always ensuring that $f \rightarrow -a$ in the $s$-wave limit. For this, let us go back to the Schr\"odinger equation Eq.~\ref{e:schrodinger}. If we knew the solution to the following equation: \begin{equation}
(\nabla^2 + k^2)G_k(\vect{r}) = \delta(\vect{r})
\label{e:defGreen}
\end{equation} we could write an integral equation for the wave function $\Psi_{\vect{k}}(\vect{r})$ as follows: \begin{equation}
\Psi_{\vect{k}}(\vect{r}) = e^{i \vect{k}\cdot\vect{r}} + \int d^3 r' G_k(\vect{r}-\vect{r}')v(\vect{r}')\Psi_{\vect{k}}(\vect{r'})
\label{e:integralequation} \end{equation} This can be simply checked by inserting this implicit solution for $\Psi_{\vect{k}}$ into Eq.~\ref{e:schrodinger}. $G_k(\vect{r})$ can be easily obtained from the Fourier transform of Eq.~\ref{e:defGreen}, defining $G_k(\vect{p}) = \int d^3 r e^{-i \vect{p} \cdot \vect{r}} G_k(\vect{r})$: \begin{equation} (-p^2 + k^2)G_k(\vect{p}) = 1 \end{equation} The solution for $G_k(\vect{r})$ is \begin{equation} \label{e:Green} G_{k,+}(\vect{r}) = \int \frac{d^3 p}{(2\pi)^3} \frac{e^{i \vect{p} \cdot \vect{r}}}{k^2 - p^2 + i \eta} = -\frac{1}{4\pi}\frac{e^{ikr}}{r} \end{equation} where we have chosen (by adding the infinitesimal constant $i \eta$, with $\eta>0$ in the denominator) the solution that corresponds to an outgoing spherical wave. $G_{k,+}(\vect{r})$ is the {\it Green's function} of the scattering problem. Far away from the origin, $|\vect{r}-\vect{r'}| \sim r - \vect{r'}\cdot \vect{u}$, with the unit vector $\vect{u} = \vect{r}/r$, and \begin{equation} \Psi_{\vect{k}}(\vect{r}) \sim e^{i \vect{k} \cdot \vect{r}} - \frac{e^{i k r}}{4\pi r} \int d^3 r' e^{-i \vect{k}'\cdot \vect{r}'} v(\vect{r}')\Psi_{\vect{k}}(\vect{r'}) \end{equation} where $\vect{k}' = k \vect{u}$. With Eq.~\ref{e:psiasymptotic}, this invites the definition of the scattering amplitude via \begin{equation} f(\vect{k}',\vect{k}) = -\frac{1}{4\pi} \int d^3 r\, e^{-i \vect{k}'\cdot \vect{r}} v(\vect{r})\Psi_{\vect{k}}(\vect{r}) \end{equation} Inserting the exact formula for $\Psi_{\vect{k}}(\vect{r})$, Eq.~\ref{e:integralequation}, combined with Eq.~\ref{e:Green}, leads to an integral equation for the scattering amplitude \begin{eqnarray} f(\vect{k}',\vect{k}) = -\frac{v(\vect{k}'-\vect{k})}{4\pi} +\int \frac{d^3 q}{(2\pi)^3} \, \frac{v(\vect{k}'-\vect{q})f(\vect{q},\vect{k})}{k^2-q^2+i\eta} \label{e:lippmannschwinger} \end{eqnarray} where $v(\vect{k})$ is the Fourier transform of the potential $v(\vect{r})$ (which we suppose to exist). This is the Lippmann-Schwinger equation, an exact integral equation for $f$ in terms of the potential $v$, useful to perform a perturbation expansion. Note that it requires knowledge of $f(\vect{q},\vect{k})$ for $q^2 \ne k^2$ (``off the energy shell). However, the dominant contributions to the integral do come from wave vectors $\vect{q}$ such that $q^2 = k^2$. For low-energy $s$-wave scattering, $f(\vect{q},\vect{k}) \rightarrow f(k)$ then only depends on the magnitude of the wave vector $\vect{k}$. With this approximation, we can take $f(k)$ outside the integral. Taking the limit $k \ll 1/r_0$, dividing by $f(k)$ and by $v_0 \equiv v(\vect{0})$, we arrive at \begin{eqnarray} \frac{1}{f(k)} \approx -\frac{4\pi}{v_0} + \frac{4\pi}{v_0} \int \frac{d^3 q}{(2\pi)^3}\, \frac{v(-\vect{q})}{k^2-q^2+i\eta} \label{e:scattampintegral} \end{eqnarray} If we only keep the first order in $v$, we obtain the scattering length in {\it Born approximation}, $a = \frac{v_0}{4\pi}$. For a delta-potential $V(\vect{r}) = V_0\, \delta(\vect{r})$, we obtain to first order in $V_0$ \begin{equation} V_0 = \frac{4\pi \hbar^2 a}{m} \end{equation} However, already the second order term in the expansion of Eq.~\ref{e:scattampintegral} would not converge, as it involves the divergent integral $\int \frac{d^3 q}{(2\pi)^3} \frac{1}{q^2}$. The reason is that the Fourier transform of the $\delta$-potential does not fall off at large momenta. Any physical potential {\it does} fall off at some large momentum, so this is not a ``real problem. For example, the van-der-Waals potential varies on a characteristic length scale $r_0$ and will thus have a natural momentum cut-off $\hbar/r_0$. A proper regularization of contact interactions employs the pseudo-potential~\cite{huan87} $V(\vect{r})\psi(\vect{r}) = V_0 \delta(\vect{r})\frac{\partial}{\partial r} (r \psi(\vect{r}))$. It leads exactly to a scattering amplitude $f(k) = -a/(1+ i k a)$ if $V_0 = \frac{4\pi\hbar^2 a}{m}$.
Here we will work with a Fourier transform that is equal to a constant $V_0$ at all relevant momenta in the problem, but that falls off at very large momenta, to make the second order term converge. The exact form is not important. If we are to calculate physical quantities, we will replace $V_0$ in favor of the observable quantity $a$ using the formal prescription \begin{equation} \frac{1}{V_0} = \frac{m}{4\pi\hbar^2 a} - \frac{m}{\hbar^2}\int \frac{d^3 q}{(2\pi)^3} \frac{1}{q^2} \label{e:renormalize} \end{equation} We will always find that the diverging term is exactly balanced by another diverging integral in the final expressions, so this is a well-defined procedure~\cite{melo93,haus99}.
Alternatively, one can introduce a ``brute force energy cut-off $E_R = \hbar^2/m R^2$ (momentum cut-off $\hbar/R$), taken to be much larger than typical scattering energies. Eq.~\ref{e:scattampintegral} then gives \begin{eqnarray} \frac{1}{f(k)} \approx -\frac{4\pi}{v_0} - \frac{2}{\pi} \frac{1}{R} + \frac{2 R}{\pi} k^2 - i k \label{e:scattampcutoff} \end{eqnarray} This is now exactly of the form Eq.~\ref{e:scattamplitude} with the scattering length \begin{eqnarray} a = \frac{\pi}{2}\frac{R}{1+\frac{2\pi^2 R}{v_0}} \label{e:acutoff} \end{eqnarray} For any physical, given scattering length $a$ we can thus find the correct strength $v_0$ that reproduces the same $a$ (provided that we choose $R \ll a$ for positive $a$). This approach implies an effective range $r_{\rm eff} = \frac{4}{\pi}R$ that should be chosen much smaller than all relevant distances. Note that as a function of $v_0$, only one pole of $a$ and therefore only one bound state is obtained, at $v_0 = -2\pi^2 R$.
This prompts us to discuss the relation between Eq.~\ref{e:renormalize} and Eq.~\ref{e:lippmannschwinger}: The Lippmann-Schwinger equation is an exact reformulation of Schr\"odinger's equation for the scattering problem. One can, for example, exactly solve for the scattering amplitude in the case of a spherical well potential~\cite{bray71}. In particular, all bound states supported by the potential are recovered. However, to arrive at Eq.~\ref{e:renormalize}, one ignores the oscillatory behavior of both $v(\vect{q})$ and $f(\vect{q},\vect{k})$ and replaces them by $\vect{q}$-independent constants. As a result, Eq.~\ref{e:renormalize}, with a cut-off for the diverging integral at a wave vector $1/R$, only allows for {\it one} bound state to appear as the potential strength is increased (see Eq.~\ref{e:acutoff}).
We will analyze this approximation for a spherical well of depth $V$ and radius $R$. The true scattering length for a spherical well is given by~\cite{land77qm} \begin{equation}
\frac{a}{R} = 1 - \frac{\tan(K R)}{K R}
\end{equation} with $K^2 = m V/\hbar^2$. which one can write as \begin{eqnarray}
\frac{a}{R} &=& 1 - \frac{\prod_{n=1}^\infty (1 - \frac{K^2 R^2}{n^2 \pi^2})}{\prod_{n=1}^\infty(1 - \frac{4K^2 R^2}{(2n-1)^2\pi^2})} \quad \left.%
\begin{array}{ll}
\leftarrow \mbox{Zeros of }$a-R$ &\\
\leftarrow \mbox{Resonances of }a &\\
\end{array}% \right. \end{eqnarray} In contrast, Eq.~\ref{e:renormalize} with $V_0 = - \frac{4\pi}{3} V R^3$ and the ``brute force cut-off at $1/R$ gives \begin{equation}
\frac{a}{R} = \frac{K^2 R^2}{\frac{2}{\pi}K^2 R^2 - 3}
\end{equation} The sudden cut-off strips the scattering length of all but one zero (at $V = 0$) and of all but one resonance. For a shallow well that does not support a bound state, the scattering length still behaves correctly as $a = -\frac{1}{3} \frac{V}{E_R} R$. However, the sudden cut-off $v(\vect{q}) \approx {\rm const.}$ for $q \le \frac{1}{R}$ and 0 beyond results in a shifted critical well depth to accommodate the first bound state, $V = \frac{3\pi}{2} E_R$, differing from the exact result $V = \frac{\pi^2}{4} E_R$. This could be cured by adjusting the cut-off. But for increasing well depth, no new bound state is found and $a$ saturates at $\sim R$, contrary to the exact result.
At first, such an approximation might be unsettling, as the van-der-Waals potentials of the atoms we deal with contain many bound states. However, the gas is in the ultracold regime, where the de Broglie-wavelength is much larger than the range $r_0$ of the potential. The short-range physics, and whether the wave function has one or many nodes within $r_0$ (i.e. whether the potential supports one or many bound states), is not important. All that matters is the phase shift $\delta_s$ {\it modulo $2\pi$} that the atomic wave packets receive during a collision. We have seen that with a Fourier transform of the potential that is constant up to a momentum cut-off $\hbar/R$, we can reproduce any low-energy scattering behavior, which is described by the scattering length $a$. We can even realize a wide range of combinations of $a$ and the effective range $r_{\rm eff}$ to capture scattering at finite values of $k$. An exception is the situation where $0 < a \lesssim r_{\rm eff}$ or potentials that have a negative effective range. This can be cured by more sophisticated models (see the model for Feshbach resonances in chapter~\ref{c:feshbach}).
\subsection{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~\footnote{Repulsive interactions still allow for the possibility of induced $p$-wave superfluidity (Kohn and Luttinger~\cite{kohn65}, also see~\cite{bara98}) however at very low temperatures $T_C \approx E_F \exp[-13(\pi/2k_F|a|)^2]$.}. 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 from 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.
\subsubsection{Two-body bound states in 1D, 2D and 3D} \label{s:boundstates}
\begin{figure}
% Requires \usepackage{graphicx}
\includegraphics[width=5.5in]{figs_BECBCSCrossover/boundstates.eps}\\
\caption[Bound state wave functions in 1D, 2D and 3D]{Bound
state wave functions in 1D, 2D and 3D for a potential well of size $R$ and depth $V$. In 1D and 2D, bound states exist for arbitrarily shallow wells. In terms of the small parameter $\epsilon = V/E_R$ with $E_R = \hbar^2/m R^2$, the size of the bound state in 1D is $R/\epsilon$. In 2D, the bound state is exponentially large, of size $R e^{-1/\epsilon}$. In 3D, due to the steep slope in $u(r) = r \psi(r)$, bound states can only exist for well depths $V_{\rm 3D}$ larger than a certain threshold $V_c \approx E_R$. The size of the bound state diverges as $R E_R / (V_{\rm 3D}-V_c)$ for $V_{\rm 3D}>V_c$. }\label{f:squarewell} \end{figure}
Localizing a quantum-mechanical particle of mass $\mu = m/2$ to a certain range $R$ leads to an increased momentum uncertainty of $p \sim \hbar/R$ at a kinetic energy cost of about $E_{R} = p^2/m = \hbar^2 / m R^2$. Clearly, a shallow potential well of size $R$ and depth $V$ with $V/E_R \equiv \epsilon \ll 1$ cannot confine the particle within its borders. But we can search for a bound state at energy $|E_B| \ll E_R$ of much larger size $r_B = 1/\kappa \equiv \sqrt{\hbar^2/m |E_B|} \gg R$.
\begin{itemize} \item {\bf 1D}: The bound state wave function far away from the well necessarily behaves like $e^{\pm \kappa x}$ for negative (positive) $x$ (see Fig.~\ref{f:squarewell}a). As we traverse the well, the wave function has to change its slope by $2\kappa$ over a range $R$. This costs kinetic energy $\approx \hbar^2 \kappa/m R$ that has to be provided by the potential energy $-V$. We deduce that $\kappa \approx m R V / \hbar^2 = \epsilon/R$, where $\epsilon=V/E_R$ is a small number for a weak potential. The size of the bound state $r_B \approx R / \epsilon$ is indeed much larger than the size of the well, and the bound state energy $E_B \approx - E_R\, \epsilon^2/2$ depends quadratically on the weak attraction $-V$. Importantly, we can {\it always} find a bound state even for arbitrarily weak (purely) attractive potentials.
\item {\bf 2D}: For a spherically symmetric well, the Schr\"odinger equation for the radial wave function $\psi(r)$ {\it outside} the well reads $\frac{1}{r}\partial_r(r\partial_r \psi) = \kappa^2 \psi$. The solution is the modified Bessel function which vanishes like $e^{-\kappa r}$ as $r \gg 1/\kappa$ (see Fig.~\ref{f:squarewell}b). For $R\ll r \ll 1/\kappa$, we can neglect the small bound state energy $E_B \propto -\kappa^2$ compared to the
kinetic energy and have $\partial_r(r\psi') = 0$ or $\psi(r) \approx \log (\kappa r)/\log(\kappa R)$, where $1/\kappa$ is the natural scale of evolution for
$\psi(r)$ and we have normalized $\psi$ to be of order 1 at $R$. Note that in 2D, it is not the change in the slope $\psi'$ of the wave function
which costs kinetic energy, but the change in $r \psi'$. {\it Inside} the well, we can assume $\psi(r)$ to be practically constant as $V \ll E_R$.
Thus, $r \psi'$ changes from $\approx 1/\log \kappa R$ (outside) to $\approx 0$ (inside) over a distance $R$. The corresponding kinetic energy cost
is $\frac{\hbar^2}{m r}\partial_r(r\psi')/\psi \approx \hbar^2/m R^2 \log (\kappa R) = E_R /\log (\kappa R)$, which has to be provided by the potential
energy $-V$. We deduce $\kappa \approx \frac{1}{R}\, e^{-c E_R/V}$ and $E_B \approx -E_R\, e^{-2c E_R / V}$ with $c$ on the
order of 1. The particle is extremely weakly bound, with its bound state energy depending exponentially on the shallow potential $-V$.
Accordingly, the size of the bound state is exponentially large, $r_B \approx R\, e^{c E_R/V}$. Nevertheless, we can {\it always} find
this weakly bound state, for arbitrarily small attraction.
\item {\bf 3D}: For a spherically symmetric well, the Schr\"odinger equation for the wave function $\psi$ transforms into an effective one-dimensional problem for the
wave function $u = r \psi$ (see Fig.~\ref{f:squarewell}c). We might now be tempted to think that there must always be a bound state in 3D, as we already found this to be the
case in 1D. However, the boundary condition on $u(r)$ is now to vanish linearly at $r=0$, in order for $\psi(0)$ to be finite. Outside the potential well, we still have
$u \propto e^{-\kappa r}$ 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 $-\kappa$ outside to
$\sim 1/R$ inside the well over a distance $R$. This costs the large kinetic energy $\sim\hbar^2 u/2m u \approx \hbar^2 /m R^2 = E_R$. If the well depth $V$ is smaller than a {\it critical depth} $V_c$ on the order of $E_R$, the particle cannot be bound. At $V=V_c$, the first
bound state enters at $E=0$. As $\kappa=0$, $u$ is then constant outside the well. If the potential depth is further increased by a small amount $\Delta V \ll V_c$, $u$ again
falls off like $e^{-\kappa r}$ for $r > R$. This requires an additional change in slope by $\kappa$ over the distance $R$, provided by $\Delta V$. So we find analogously to the
1D case $\kappa \sim m R \Delta V / \hbar^2$. Hence, the bound state energy $E_B \approx - \Delta V^2 / E_R$ is quadratic in the ``detuning $\Delta V = (V-V_C)$, and the size
of the bound state diverges as $r_B \approx R E_R / (V - V_C)$. We will find exactly the same behavior for a weakly bound state when discussing Feshbach resonances
in chapter~\ref{c:feshbach}.
\end{itemize}
\begin{table} \centering \begin{tabular}{c|c|c|c}
% after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
& 1D & 2D & 3D \\ \hline
$V$ & $\ll E_R$ & $\ll E_R$ & $> V_c \approx E_R$ \\[2pt] \hline
$\psi(r>R)$ & $e^{-r/r_B}$ & $K_0(\frac{r}{r_B}) = \left\{%
\begin{array}{ll}
-\log r/r_B, & \hbox{$R \ll r \ll r_B$} \\
e^{-r/r_B}, & \hbox{$r \gg r_B$} \\
\end{array}% \right.$ & $\frac{e^{- r / {r_B}}}{r}$ \\[12pt] \hline
$r_B$ & $R \frac{E_R}{V}$ & $R\, e^{c E_R / V}$ & $R \frac{E_R}{V-V_c}$ \\[4pt]\hline
$E_B = -\frac{\hbar^2}{m r_B^2}$ & $-V^2 / E_R$ & $-E_R e^{-2 c E_R / V}$ & $-(V-V_c)^2 / E_R$ \\
\end{tabular} \caption{Bound states in 1D, 2D and 3D for a potential well of size $R$ and depth $V$. $\psi(r>R)$ is the wave function outside the well, $r_B$ is the size of the bound state, and $E_B$ its energy ($E_R = \hbar^2/m R^2$).}\label{t:boundstate} \label{t:boundstates} \end{table}
The analysis holds for quite general shapes $V(r)$ of the
(purely attractive) potential well (in the equations, we only need to replace $V$ by
its average over the well - if it exists -,
$\frac{1}{R}\int_{-\infty}^\infty V(x) dx$ in 1D,
$\frac{1}{R^2}\int_0^\infty r V(r) dr$ in 2D etc.).
Table~\ref{t:boundstates} summarizes the different cases.
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.
\subsubsection{Density of states}
\begin{table} \centering \begin{tabular}{c|c|c|c}
% after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
& 1D & 2D & 3D \\ \hline
$\frac{\hbar^2}{m\Omega}\rho(\epsilon)$ & $\frac{1}{\pi}\sqrt{\frac{\hbar^2}{2m\epsilon}}$ & $\frac{1}{2\pi}$ & $\frac{1}{2\pi^2}\sqrt{\frac{2m\epsilon}{\hbar^2}}$ \\[2pt] \hline
$\frac{1}{\left|V_0\right|} = \frac{1}{\Omega}\int_{\epsilon<E_R} d\epsilon \frac{\rho_n(\epsilon)}{2\epsilon+\left|E\right|}$ & $\sqrt{\frac{m}{4\hbar^2 \left|E\right|}}$ & $\frac{m}{4\pi\hbar^2}\log\frac{2E_R + |E|}{|E|}$ & $\frac{1}{2\pi^2}\frac{m^{3/2}}{\hbar^3}(\sqrt{2E_R} - \frac{\pi}{2} \sqrt{|E|})$ \\[4pt]\hline
$E = -\frac{\hbar^2 \kappa^2}{m}$ & $-\frac{m}{4\hbar^2} V_0^2$ & $-2 E_R\, e^{-\frac{4\pi\hbar^2}{m\left|V_0\right|}}$ & $-\frac{8}{\pi^2}E_R\,\frac{(\left|V_0\right| - V_{0c})^2}{\left|V_0\right|^2} = -\hbar^2/m a^2$ \\ \end{tabular} \caption{Link between the density of states and the existence of a bound state for arbitrarily weak interaction. The table shows the density of states, $\rho(\epsilon)$, the equation relating the bound state energy $E$ to $V_0$, and the result for $E$. It is assumed that $E_R \gg |E|$. To compare with table~\ref{t:boundstates} note that $|V_0| \sim V R^n$. $V_{0c} = \sqrt{2}\,\pi^2 E_R R^3$ is the threshold interaction strength for the 3D case. The formula for the 3D bound state energy follows from the renormalization procedure outlined in section~\ref{s:renormalization}, when expressing $V_0$ in terms of the scattering length $a$ using Eq.~\ref{e:renormalize}.} \label{t:momentumboundstates} \end{table}
What physical quantity decides whether there are bound states or
not? To answer this question, we formulate the problem of two
interacting particles of mass $m$ in momentum space. This allows a
particularly transparent treatment for all three cases (1D, 2D,
3D) at once, and identifies the {\it 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 $E = -\frac{\hbar^2 \kappa^2}{m}$ ($m/2$ is the reduced mass), we start by writing the Schr\"odinger equation for the relative wave function $\frac{\hbar^2}{m}(\nabla^2-\kappa^2)\psi = V \psi$ in ($n$-dimensional) momentum space: \begin{equation}
\psi_\kappa(\vect{q}) = - \frac{m}{\hbar^2}\frac{1}{q^2 + \kappa^2} \int \frac{d^n q'}{(2\pi)^n} V(\vect{q}-\vect{q}') \psi_\kappa(\vect{q}')
\end{equation} For a short-range potential of range $R \ll 1/\kappa$, $V(\vect{q})$ is practically constant, $V(\vect{q}) \approx V_0$, for all relevant $q$, and falls off to zero on a large $q$-scale of $\approx 1/R$. For example, for a potential well of depth $V$ and size $R$, we have $V_0 \sim - V R^n$. Thus, \begin{equation}
\psi_\kappa(\vect{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(\vect{q}')
\end{equation} We integrate once more over $\vect{q}$, applying the same cut-off $1/R$, and then divide by the common factor $\int_{q\lesssim\frac{1}{R}} \frac{d^n q}{(2\pi)^n} \psi_\kappa(\vect{q})$. We obtain the following equation for the bound state energy $E$: \begin{equation}
- \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|}
\label{e:densityboundstates} \end{equation} with the density of states in $n$ dimensions $\rho_n(\epsilon)$, the energy cut-off $E_R = \hbar^2/m R^2$ and the volume $\Omega$ of the system (note that $V_0$ has units of energy times volume). The question on the existence of bound states for arbitrarily weak interaction has now been reformulated: As $|V_0| \rightarrow 0$, the left hand side of Eq.~\ref{e:densityboundstates} diverges. This equation has a solution for small $|V_0|$ only if the right hand side also diverges for vanishing bound state energy $|E| \rightarrow 0$, and this involves an integral over the density of states. Table~\ref{t:momentumboundstates} presents the different cases in 1D, 2D, 3D. In 1D, the integral diverges as $1/\sqrt{|E|}$, so one can always find a bound state solution. The binding energy depends quadratically on the interaction, as we had found before. In 2D, where the density of states $\rho_{\rm 2D}$ is {\it constant}, the integral still diverges logarithmically as $|E|\rightarrow 0$, so that again there is a solution $|E|$ for any small $|V_0|$. The binding energy now depends exponentially on the interaction and $\rho_{\rm 2D}$: \begin{equation}
E_{\rm 2D} = - 2 E_R \, e^{-\frac{2\Omega}{\rho_{\rm 2D} \left|V_0\right|}}
\label{e:boundstate2D}
\end{equation} However, in 3D the integral is finite for vanishing $|E|$, 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 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. This is exactly what happens, as will be discussed in the next section.
\subsubsection{Pairing of fermions -- The Cooper problem} \label{s:cooperproblem}
\begin{figure}
% Requires \usepackage{graphicx}
\centering
\includegraphics[width=4in]{figs_BECBCSCrossover/cooperproblem.eps}\\
\caption[Cooper problem: Two particles scattering on top of a Fermi sea]{Cooper problem: Two particles scattering on top of a Fermi sea. a) Weakly interacting particles with equal and opposite momenta can scatter into final states in a narrow shell (blue-shaded) on top of the Fermi sea (gray shaded), which blocks possible final momentum states. b) For non-zero total momentum $2\vect{q}$, particles can scatter only in a narrow band around a circle with radius $\sqrt{k_F^2 - q^2}$.}\label{f:cooperproblem}
\end{figure}
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~\cite{coop56}. Momentum states below the Fermi surface are not available to the two scattering particles due to Pauli blocking (Fig.~\ref{f:cooperproblem}a). 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 $\rho_{\rm 3D}(E_F)$, which is a constant just like in two dimensions. We should thus find a {\it bound state} for the two-particle system {\it 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, as we can see from Fig.~\ref{f:cooperproblem}a. If the pairs have finite center-of-mass momentum $\vect{q}$, the number of contributing states is strongly reduced, as they are confined to a circle (see Fig.~\ref{f:cooperproblem}b). 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 $E_B = E-2E_F<0$ on top of the large Fermi energy $2E_F$ of the two particles. The equation for $E_B$ is \begin{equation}
-\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|}
\label{e:Cooper} \end{equation} The effect of Pauli blocking of momentum states below the Fermi surface is explicitly included by only integrating over energies $\epsilon > E_F$.
In conventional superconductors, the natural cut-off energy $E_R$ is given by the Debye frequency $\omega_D$, $E_R = \hbar \omega_D$, corresponding to the highest frequency at which ions in the crystal lattice can respond to a bypassing electron. Since we have $\hbar \omega_D \ll E_F$, we can approximate $\rho_{\rm 3D}(\epsilon) \approx \rho_{\rm 3D}(E_F)$ and find: \begin{equation}
E_B = - 2 \hbar \omega_D e^{-2 \Omega/\rho_{\rm 3D}(E_F) \left|V_0\right|}\\
\label{e:coopersuperconductor}
\end{equation}
In the case of an atomic Fermi gas, we should replace $1/V_0$ by the physically relevant scattering length $a < 0$ using the prescription in Eq.~\ref{e:renormalize}. The equation for the bound state becomes \begin{equation}
-\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} \label{e:cooperrenorm} \end{equation} The right hand expression is now finite as we let the cut-off $E_R \rightarrow \infty$, the result being (one assumes $\left|E_B\right|\ll E_F$) \begin{equation}
-\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) \end{equation} Inserting $\rho_{\rm 3D}(E_F) = \frac{\Omega m k_F}{2\pi^2 \hbar^2}$ with the Fermi wave vector $k_F = \sqrt{2mE_F/\hbar^2}$, one arrives at \begin{equation}
E_B = - \frac{8}{e^2} E_F\, e^{-\pi/k_F \left|a\right|}
\label{e:cooperproblem}
\end{equation} The binding energies Eqs.~\ref{e:coopersuperconductor} and \ref{e:cooperproblem} can be compared with the result for the bound state of two particles in 2D, Eq.~\ref{e:boundstate2D}. The role of the constant density of states $\rho_{\rm 2D}$ is here played by the 3D density of states at the Fermi surface, $\rho_{\rm 3D}(E_F)$.
The result is remarkable: 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 {\it in} 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~\cite{bard57}. 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 $\Delta$ that differs from $|E_B|$ (Eq.~\ref{e:cooperproblem}) by a factor of 2 in the exponent: \begin{equation}
\Delta = \frac{8}{e^2} E_F\, e^{-\pi/2 k_F \left|a\right|}
\end{equation}
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 $n$ and have two fermions scatter on the $(n-1)$ 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 $(n-1)$ 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.
\subsection{Crossover wave function} \label{s:crossoverwavefunction} From section~\ref{s:boundstates} we know that in 3D, two fermions in isolation can form a molecule for strong enough attractive interaction. The ground state of the system should be a Bose-Einstein condensate of these tightly bound pairs. However, if we increase the density of particles in the system, we will ultimately reach the point where the Pauli pressure of the fermionic constituents becomes important and modifies the properties of the system. When the Fermi energy of the constituents exceeds the binding energy of the molecules, we expect that the equation of state will be fermionic, i.e. the chemical potential will be proportional to the density to the power of 2/3. Only when the size of the molecules is much smaller than the interparticle spacing, i.e. when the binding energy largely exceeds the Fermi energy, is the fermionic nature of the constituents irrelevant -- tightly bound fermions are spread out widely in momentum space and do not run into the Pauli limitation of unity occupation per momentum state.
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, the binding energy of these pairs is much less than the Fermi energy and therefore Pauli pressure plays a major role.
It was realized by Leggett~\cite{legg80}, building upon work by Popov~\cite{popo66}, Keldysh~\cite{keld68} and Eagles~\cite{eagl69}, that 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. Closely following Leggett's work~\cite{legg80}, and its extension to finite temperatures by Nozi\`eres and Schmitt-Rink~\cite{nozi85}, we will describe the BEC-BCS crossover in a simple ``one-channel model of a potential well. Rather than the interaction strength $V_0$ as in section~\ref{s:boundstates}, we will take the scattering length $a$ as the parameter that ``tunes the interaction. The relation between $V_0$ and $a$ is given by Eq.~\ref{e:renormalize} and its explicit form Eq.~\ref{e:acutoff}. For positive $a>0$, there is a two-body bound state available at $E_B = -\hbar^2/m a^2$ (see table~\ref{t:momentumboundstates}), while small and negative $a<0$ corresponds to weak attraction where Cooper pairs can form in the medium. In either case, for $s$-wave interactions, the orbital part of the pair wave function $\varphi(\vect{r}_1,\vect{r}_2)$ will be symmetric under exchange of the paired particles' coordinates and, in a uniform system, will only depend on their distance $\left|\vect{r}_1-\vect{r}_2\right|$. We will explore the many-body wave function \begin{equation}
\Psi\left(\vect{r}_1,\dots,\vect{r}_N\right) = \mathcal{A}\left\{\varphi(\left|\vect{r}_1-\vect{r}_2\right|)\chi_{12}\dots \varphi(\left|\vect{r}_{N-1}-\vect{r}_N\right|)\chi_{N-1,N}\right\}
\label{e:fermicondensatepsi} \end{equation} that describes a condensate of such fermion pairs, with the operator $\mathcal{A}$ denoting the correct antisymmetrization of all fermion coordinates, and the spin singlet $\chi_{ij} = \frac{1}{\sqrt{2}}(\left|\uparrow\right>_i \left|\downarrow\right>_j - \left|\downarrow\right>_i \left|\uparrow\right>_j)$. In the experiment, ``spin up and ``spin down will correspond to two atomic hyperfine states.
In second quantization notation we write \begin{equation}
\left|\Psi\right>_N = \int \prod_i d^3 r_i \,\varphi(\vect{r}_1 - \vect{r}_2) \Psi_\uparrow^\dagger(\vect{r}_1) \Psi_\downarrow^\dagger(\vect{r}_2) \dots \varphi(\vect{r}_{N-1} - \vect{r}_N) \Psi_\uparrow^\dagger(\vect{r}_{N-1}) \Psi_\downarrow^\dagger(\vect{r}_N) \left|0\right>
\end{equation} where the fields $\Psi_\sigma^\dagger(\vect{r}) = \sum_k c_{k\sigma}^\dagger \frac{e^{-i \vect{k} \cdot \vect{r}}}{\sqrt{\Omega}}$. With the Fourier transform $\varphi(\vect{r}_1-\vect{r}_2) = \sum_k \varphi_k \frac{e^{i \vect{k} \cdot (\vect{r}_1- \vect{r}_2)}}{\sqrt{\Omega}}$ we can introduce the pair creation operator \begin{equation}
b^\dagger = \sum_k \varphi_k c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger
\end{equation} and find \begin{equation}
\left|\Psi\right>_N = {b^\dagger}^{N/2} \left|0\right>
\end{equation} This expression for $\left|\Psi\right>_N$ is formally identical to the Gross-Pitaevskii ground state of a condensate of bosonic particles. However, the operators $b^\dagger$ obey bosonic commutation relations only in the limit of tightly bound pairs. For the commutators, we obtain \begin{eqnarray} \label{e:commutators}
\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}
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 $n_k$ of any momentum state $k$ are very small (see section~\ref{s:evolution} below), and $\left[b,b^\dagger\right]_- \approx \sum_k |\varphi_k|^2 = \int {\rm d}^3 r_1\int {\rm d}^3 r_2 \,|\varphi(\vect{r}_1,\vect{r}_2)|^2 = 1$.
Working with the $N$-particle state $\left|\Psi\right>_N$ is inconvenient, as one would face a complicated combinatoric problem in manipulating the sum over all the $c_k^\dagger$'s (as one chooses a certain $k$ for the first fermion, the choices for the second depend on this $k$, etc.). It is preferable to use the grand canonical formalism, not fixing the number of atoms but the chemical potential $\mu$. 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 $N_p = N/2$ be the number of pairs. Then, \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} 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$. If we choose the constant $\mathcal{N} = \prod_k \frac{1}{u_k} = \prod_k \sqrt{1 + N_p |\varphi_k|^2}$, then $\left|\Psi\right>$ becomes a properly normalized state \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} with $v_k = \sqrt{N_p}\,\varphi_k u_k$ and $|u_k|^2 + |v_k|^2 = 1$. This is the famous BCS wave function, first introduced as a variational Ansatz, later shown to be the exact solution of the simplified Hamiltonian Eq.~\ref{e:Hsimplified} (below). It is a product of wave functions referring to the occupation of pairs of single-particle momentum states, $(\vect{k},\uparrow,-\vect{k},\downarrow)$. As a special case, it describes a non-interacting Fermi sea, with all momentum pairs occupied up to the Fermi momentum ($u_k=0, v_k=1$ for $k<k_F$ and $u_k=1, v_k=0$ for $k>k_F$). In general, for a suitable choice of the $v_k$'s and $u_k$'s, it describes 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. The above derivation makes it clear that this wave function encompasses the entire regime of pairing, from point bosons (small molecules) to weakly and non-interacting fermions.
\subsection{Gap and number equation}
The variational parameters $u_k$ and $v_k$ are derived in the standard way by minimizing the free energy $E - \mu N = \left<\hat{H} - \mu \hat{N}\right>$. The many-body Hamiltonian for the system is \begin{equation}
\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}
\label{e:Hamiltonian} \end{equation}
The dominant role in superfluidity is played by fermion pairs with zero total momentum. Indeed, as we have seen in section~\ref{s:cooperproblem}, Cooper pairs with zero momentum have the largest binding energy. Therefore, we simplify the mathematical description by neglecting interactions between pairs at finite momentum, i.e. we only keep the terms for $\vect{q} = 0$. This is a very drastic simplification, as hereby density fluctuations are eliminated. It is less critical for charged superfluids, where density fluctuations are suppressed by Coulomb interactions. However, for neutral superfluids, sound waves (the Bogoliubov-Anderson mode, see section~\ref{s:collexcitations}) are eliminated by this approximation. The approximate Hamiltonian (``BCS Hamiltonian) reads \begin{equation}
\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}
\label{e:Hsimplified} \end{equation} The free energy becomes \begin{eqnarray}
\label{e:free-energy}
\mathcal{F} = \left<\hat{H} - \mu \hat{N}\right> &=& \sum_k 2 \xi_k v_k^2 + \frac{V_0}{\Omega}\sum_{k,k'} u_k v_k u_{k'} v_{k'}\\
\mbox{with }\;\xi_k &=& \epsilon_k - \mu \nonumber
\end{eqnarray} Minimizing $E-\mu N$ leads to \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}
\label{e:ukvk}
\end{eqnarray} where $\Delta$ is given by the {\it gap equation} $\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}$ or \begin{equation}
-\frac{1}{V_0} = \int \frac{d^3 k}{\left(2\pi\right)^3} \;\frac{1}{2E_{k}}
\end{equation} Note the similarity to the bound state equation in free space, Eq.~\ref{e:densityboundstates}, and in the simplified Cooper problem, Eq.~\ref{e:Cooper}. An additional constraint is given by the {\it number equation} for the total particle density $n = N / \Omega$ \begin{equation} n = 2 \int \frac{d^3 k}{\left(2\pi\right)^3} \; v_k^2 \end{equation} Gap and number equations have to be solved simultaneously to yield the two unknowns $\mu$ and $\Delta$. We will once more replace $V_0$ by the scattering length $a$ using prescription Eq.~\ref{e:renormalize}, so that the gap equation becomes (compare Eq.~\ref{e:cooperrenorm}) \begin{equation} -\frac{m}{4\pi\hbar^2 a} = \int \frac{d^3 k}{\left(2\pi\right)^3} \left(\frac{1}{2 E_k} - \frac{1}{2 \epsilon_k}\right) \end{equation} where the integral is now well-defined. The equations can be rewritten in dimensionless form with the Fermi energy $E_F = \hbar^2 k_F^2 / 2m$ and wave vector $k_F = (3 \pi^2 n)^{1/3}$~\cite{orti05bcs} \begin{eqnarray} -\frac{1}{k_F a} &= &\frac{2}{\pi} \sqrt{\frac{\Delta}{E_F}}\; I_1\left(\frac{\mu}{\Delta}\right) \\ 1 &=& \frac{3}{2}\left(\frac{\Delta}{E_F}\right)^{3/2} I_2\left(\frac{\mu}{\Delta}\right)\\ \mbox{with }\; I_1(z) &=& \int_0^\infty dx \;x^2 \left(\frac{1}{\sqrt{\left(x^2 - z\right)^2 + 1}} - \frac{1}{x^2}\right)\\ \mbox{and }\; I_2(z) &=& \int_0^\infty dx \;x^2 \left(1 - \frac{x^2 - z}{\sqrt{\left(x^2-z\right)^2 + 1}}\right) \end{eqnarray} This gives \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} The first equation can be inverted to obtain $\mu / \Delta$ as a function of the {\it interaction parameter} $1/k_F a$, which can then be inserted into the second equation to yield the gap $\Delta$. The result for $\mu$ and $\Delta$ as a function of $1/k_F a$ is shown in Fig.~\ref{f:Deltamu}. It is possible to obtain analytic expressions for the solutions in terms of complete elliptic integrals~\cite{mari98becbcs}.
\begin{figure} \centering
\includegraphics[width=3in]{figs_BECBCSCrossover/deltamu.eps}\\
\caption[Chemical potential and gap in the BEC-BCS crossover]{Chemical potential (dotted line) and gap (straight line, red) in the BEC-BCS crossover as a function of the interaction parameter $1/k_Fa$. The BCS-limit of negative $1/k_F a$ is to the right on the graph. The resonance where $1/k_F a = 0$ is indicated by the dashed line.}\label{f:Deltamu}
\end{figure}
In this derivation, we have combined the simplified Hamiltonian, Eq.~\ref{e:Hsimplified} with the BCS variational Ansatz. Alternatively one can apply a decoupling (mean field) approximation to the Hamiltonian~\cite{peth02bec}. 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 \begin{eqnarray}
C_k= \left<c_{k\uparrow} c_{-k\downarrow}\right>= -\left< c_{k\uparrow}^\dagger
c_{-k\downarrow}^\dagger\right>
\end{eqnarray} 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 \begin{equation} c_{k\uparrow} c_{-k\downarrow} = C_k + (c_{k\uparrow} c_{-k\downarrow} - C_k) \end{equation} with the operator in parentheses giving rise to fluctuations that are small on the scale of $C_k$. The gap parameter $\Delta$ is now defined as \begin{equation} \Delta=\frac{V_0}{\Omega} \sum_{k} C_k \end{equation} 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 \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}
This Hamiltonian is bilinear in the creation and annihilation operators and can easily be solved by a Bogoliubov transformation~\cite{bogo58,vala58,peth02bec} from the particle operators $c_{k\downarrow}$ and $c_{k\uparrow}$ to new quasi-particle operators $\gamma_{k\uparrow}$ and $\gamma_{k\downarrow}$: \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
\label{e:quasiparticles}
\end{eqnarray} The $u_k$ and $v_k$ 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 $u_k$ and $v_k$, and the transformed Hamiltonian becomes \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}) \label{e:BogoliubovH} \end{equation} The first two terms give the free energy $E - \mu N$ of the pair condensate, identical to Eq.~\ref{e:free-energy} when the correct $u_k$ and $v_k$ are inserted. The third term represents the energy of excited quasi-particles, and we identify $E_k$ as excitation energy of a quasi-particle. The superfluid ground state is the quasi-particle vacuum: $\gamma_{k\uparrow} \left|\Psi\right> = 0 = \gamma_{k\downarrow} \left|\Psi\right>$.
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 $\sqrt{N_0}$, the square root of the number of condensed atoms (i.e.~a coherent field). In the interaction term of the Hamiltonian all terms are dropped that contain less than two factors of $\sqrt{N_0}$. In other words, the Hamiltonian (Eq.~\ref{e:Hsimplified}) 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 $Eq.~\ref{e:Hsimplified}$~\cite{duke04review}.
\subsection{Discussion of the three regimes -- BCS, BEC and crossover} \subsubsection{BCS limit} \label{s:BCSlimit} In the BCS-limit of weak attractive interaction, $k_F a \rightarrow 0_-$, we have\footnote{This follows by substituting $\xi = x^2 - z$ in the integrals and taking the limit $z\rightarrow \infty$. One has $I_1(z) \approx \sqrt{z}\left(\log(8z) - 2\right)$ and $I_2(z) = \frac{2}{3} z^{3/2}$.} \begin{eqnarray}
\mu &\approx& E_F \\
\Delta &\approx& \frac{8}{e^2} e^{-\pi/2k_F\left|a\right|}
\label{e:BCSLimit}
\end{eqnarray}
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 $2 \mu$). 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\footnote{The present mean field treatment does not include density fluctuations, which modify the prefactor in the expression for the gap $\Delta$ ~\cite{gork61,peth02bec}.}. Compared to the bound state energy for a single Cooper pair on top of a non-interacting Fermi sea, Eq.~\ref{e:cooperproblem}, 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\footnote{In the self-consistent BCS solution, not only the momentum states above the Fermi surface contribute to pairing, but also those {\it below} it, in a symmetric shell around the Fermi momentum. In the Cooper problem the states below the Fermi surface were excluded, reducing the effective density of states by a factor of two.}. 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 Eq.~\ref{e:free-energy} and is \begin{equation}
E_{\rm G,\, BCS} = \frac{3}{5} N E_F - \frac{1}{2}\,\rho(E_F)\, \Delta^2
\end{equation}
The first term is the energy of the non-interacting normal state, where $\frac{3}{5}E_F$ is the average kinetic energy per fermion in the Fermi sea. 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.
Although the total kinetic energy of the Fermi gas has been increased (by populating momentum states above $E_F$), the total energy is lower due to the gain in potential energy. This is valid for any kind of pairing (i.e. proton and electron forming a hydrogen atom), since the localization of the pair wave function costs kinetic energy.
The energy of the BCS state, $- \frac{1}{2}\,\rho(E_F)\, \Delta^2$ can be interpreted in two ways. One way refers to the wave function Eq.~\ref{e:fermicondensatepsi}, which consists of $N/2$ identical fermion pairs. The energy per pair is then $- \frac{3}{4} \Delta^2/E_F$. The other interpretation refers to the BCS wave function Eq.~\ref{e:BCSstate}. It is essentially a product of a ``frozen Fermi sea (as $v_k \approx 1$, $u_k \approx 0$ for low values of $k$) with a paired component consisting of $\sim \rho(E_F) \, \Delta \sim N \Delta/E_F$ pairs, located in an energy shell of width $\Delta$ around the Fermi energy. They each contribute a pairing energy on the order of $\Delta$. The second interpretation justifies the picture of a Cooper pair condensate. In the solution of the Cooper problem (section~\ref{s:cooperproblem}), the pair wave function has a peak occupation per momentum state of $\sim 1/ \rho(E_F) \Delta$. Therefore, one can stack up $\sim \rho(E_F) \Delta$ pairs with zero total momentum without getting into serious trouble with the Pauli exclusion principle and construct a Bose-Einstein condensate consisting of $\sim \rho(E_F) \Delta$ Cooper pairs~\footnote{Similarly to the fermion pairs described by the operator $b^\dagger$, the Cooper pairs from section~\ref{s:cooperproblem} are not bosons, as shown by the equivalent of Eq.~\ref{e:commutators}. However, if there were only a few Cooper pairs, much less than $\rho(E_F) \Delta$, the occupation of momentum states $n_k$ 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 $\frac{1}{2} \Delta^2/E_F$ or of $\Delta$. In RF spectroscopy, all momentum states can be excited (see section~\ref{e:chap2RFspectroscopy}), and the spectrum shows a gap of $\frac{1}{2} \Delta^2/E_F$ (see section~\ref{s:RFspectrum}). Tunnelling experiments in superconductors probe the region close to the Fermi surface, and show a pairing gap of $\Delta$.
The two interpretations for the BCS energy carry along two possible choices of the pairing wave function (see section~\ref{s:evolution}). The first one is $\varphi_k = u_k/v_k\sqrt{N_p}$, which can be shown to extend throughout the whole Fermi sea from zero to slightly above $k_F$, whereas the second one, $\psi(k) = u_k v_k$, is concentrated around the Fermi surface (see Fig.~\ref{f:excitation}).
To give a sense of scale, Fermi energies in dilute atomic gases are on the order of a $\mu\rm K$, corresponding to $1/k_F \sim 4\,000\, a_0$. In the absence of scattering resonances, a typical scattering length would be about $50-100\; a_0$ (on the order of the van der Waals-range). Even if $a < 0$, this would result in a vanishingly small gap $\Delta/k_B \approx 10^{-30}\dots 10^{-60}\, \rm K$. 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 $k_F \left|a\right| > 1$ (see chapter~\ref{c:feshbach}). In this case, the above mean field theory predicts $\Delta > 0.22 \;E_F$ or $\Delta/k_B > 200\, \rm nK$ for $k_F |a| > 1$, and this is the regime where current experiments are operating.
\subsubsection{BEC limit} \label{s:BEClimit}
In the BEC limit of tightly bound pairs, for $k_F a \rightarrow 0_+$, one finds\footnote{This result follows from the expansion of the integrals for $z<0$ and $|z|\rightarrow \infty$. One finds $I_1(z) = -\frac{\pi}{2}\sqrt{|z|} -\frac{\pi}{32}\frac{1}{|z|^{3/2}}$ and $I_2(z) = \frac{\pi}{8}\frac{1}{\sqrt{|z|}}$.} \begin{eqnarray}
\label{e:BEClimit}
\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}}
\label{e:BEClimitgap} \end{eqnarray}
The first term in the expression for the chemical potential is the binding energy per fermion in a tightly bound molecule (see table~\ref{t:momentumboundstates}). This reflects again the implicit assumption (made by using the wave function in Eq.~\ref{e:fermicondensatepsi}) that we always add {\it two} fermions of opposite spin at the same time to the system.
The second term is a mean field contribution describing the repulsive interaction between molecules in the gas. Indeed, a condensate of molecules of mass $m_M = 2m$, density $n_M = n/2$ and a molecule-molecule scattering length $a_M$ will have a chemical potential $\mu_M = \frac{4\pi \hbar^2 a_M n_M}{m_M}$. Since $\mu_M$ is twice the chemical potential for each fermion, we obtain from the above expression the molecule-molecule scattering length $a_M = 2 a$. However, this result is not exact. Petrov, Shlyapnikov and Salomon~\cite{petr04dimers} have performed an exact calculation for the interaction between four fermions and shown that $a_M = 0.6\, a$. The present mean field approach neglects correlations between different pairs, or between one fermion and a pair. If those are included, the correct few-body physics is recovered~\cite{pier00becbcs,holl04bosefermi,hu06becbcs}.
The expression for the quantity $\Delta$ 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 $k$ in the single-fermion excitation spectrum (see Fig.~\ref{f:pairwavefunction} below). Instead, we have for the quasi-particle energies $E_k = \sqrt{(\epsilon_k - \mu)^2 + \Delta^2} \approx |\mu| + \epsilon_k + \frac{\Delta^2}{2|\mu|}$. So $\Delta$ itself does not play a role in the BEC-regime, but only the combination $\Delta^2/|\mu|$ is important. As we see from Eq.~\ref{e:BEClimit}, \begin{equation}
\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
\label{e:Deltasquared} \end{equation} which is two times the molecular mean field. In fact, we will show in section~\ref{s:excitations} that 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 Eq.~\ref{e:Hsimplified} 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 {\it neglect} such three- and four-body interactions. The solution to this puzzle lies in the Pauli principle, which acts as an effective repulsive interaction: In a molecule, each constituent fermion is confined to a region of size $\sim a$ around the molecule's center of mass (see next section). The probability to find another like fermion in that region is strongly reduced due to Pauli blocking. Thus, effectively, the motion of molecules is constrained to a reduced volume $\Omega' = \Omega - c N_M a_M^3$, with the number of molecules $N_M$ and $c$ on the order of 1. This is the same effect one has for a gas of hard-sphere bosons of size $a_M$, and generally for a Bose gas with scattering length $a_M$. An analogous argument leads to the effective interaction between a single fermion and a 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.
\subsubsection{Evolution from BCS to BEC} \label{s:evolution} \begin{figure} \centering
\includegraphics[width=3in]{figs_BECBCSCrossover/nk.eps}\\
\caption[Occupation $n_k$ of momentum states $k$ in the BEC-BCS crossover]{Occupation $n_k$ of momentum states $k$ in the BEC-BCS crossover. The numbers give the interaction parameter $1/k_F a$. After~\cite{nozi85}.}\label{f:momoccupation}
\end{figure}
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 $1/k_F a = -1$ and +1 and is fully continuous. The occupation of momentum states $n_k = v_k^2$ evolves smoothly from the step-function $\Theta(k_F - k)$ of a degenerate Fermi gas, broadened over a width $\Delta \ll E_F$ due to pairing, to that of $N_p$ molecules, namely the number of molecules $N_p$ times the probability $|\varphi_k|^2$ to find a molecule with momentum $k$ (we have $\varphi_k = \frac{(2 \pi a)^{3/2}}{\sqrt{\Omega}}\frac{1}{\pi}\frac{1}{1+k^2 a^2}$) (see Fig.~\ref{f:momoccupation}). It is also interesting to follow the evolution of the ``Cooper pair wave function\footnote{Note that this definition is not equal to the Fourier transform of the pair wave function $\varphi(\vect{r})$ introduced in Eq.~\ref{e:fermicondensatepsi}, which would be $v_k/u_k\sqrt{N_p}$. 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 $\sim \delta k$, thus giving similar behavior for the real space wave function.} both in $k$-space, where it is given by $\left<\Psi_{BCS}\right|c^\dagger_{k\uparrow} c^\dagger_{-k\downarrow}\left|\Psi_{BCS}\right> = u_k v_k$, and in real space, where it is \begin{eqnarray} \psi(\vect{r}_1,\vect{r}_2) &=& \left<\Psi_{BCS}\right|\Psi_{\uparrow}^\dagger(\vect{r}_1)\Psi_{\downarrow}^\dagger(\vect{r}_2)\left|\Psi_{BCS}\right> = \frac{1}{\Omega}\sum_k u_k v_k e^{-i \vect{k} \cdot (\vect{r}_1-\vect{r}_2)} \nonumber \\ &=&\frac{1}{\Omega}\sum_k \frac{\Delta}{2E_k}\, e^{-i \vect{k} \cdot (\vect{r}_1-\vect{r}_2)} \label{e:cooperpairwavefunction} \end{eqnarray}
\begin{figure}
\centering
\includegraphics[width=5.5in]{figs_BECBCSCrossover/pairwavefunction.eps}\\
\caption[Evolution of the spatial pair wave function $\psi(r)$ in the BEC-BCS crossover]{Evolution of the spatial pair wave function $\psi(r)$ in the BEC-BCS crossover. The inset shows the Fourier transform $\psi(k)$, showing clearly that in the BCS-limit, momentum states around the Fermi surface make the dominant contribution to the wave function. In the crossover, the entire Fermi sphere takes part in the pairing. In the BEC-limit, $\psi(k)$ broadens as the pairs become more and more tightly bound. $\psi(r)$ was obtained via numerical integration of $\int_{-\mu}^\infty d\xi \frac{\sin(r \sqrt{\xi + \mu})}{\sqrt{\xi^2 + \Delta^2}}$ (here, $\hbar = 1 = m$), an expression that follows from Eq.~\ref{e:cooperpairwavefunction}.}\label{f:pairwavefunction}
\end{figure}
\begin{figure} \centering
\includegraphics[width=3in]{figs_BECBCSCrossover/pairsize.eps}\\
\caption[From tightly bound molecules to long-range Cooper pairs]{From tightly bound molecules to long-range Cooper pairs. Evolution of the pair size $\xi_0 =\sqrt{\frac{\left<\psi(\vect{r})\right|r^2\left|\psi(\vect{r})\right>}{\left<\psi(\vect{r})|\psi(\vect{r})\right>}}$ as a function of the interaction parameter $1/k_F a$. On resonance (dashed line), the pair size is on the order of the inverse wave vector, $\xi_0(0) \sim \frac{1}{k_F}$, about a third of the interparticle spacing.}\label{f:pairsize}
\end{figure}
In the BCS limit, the pairing occurs near the Fermi surface $k= k_F$, in a region of width $\delta k \sim \frac{\partial k}{\partial \epsilon} \delta \epsilon \approx \frac{\Delta}{\hbar v_F}$, where $v_F$ is the velocity of fermions at the Fermi surface. Therefore, the spatial wave function of Cooper pairs has a strong modulation at the inverse wave vector $1/k_F$, and an overall extent of the inverse width of the pairing region, $\sim 1/\delta k \sim \frac{\hbar v_F}{\Delta} \gg 1/k_F$. More quantitatively, Eq.~\ref{e:cooperpairwavefunction} gives (setting $r = \left|\vect{r}_1-\vect{r}_2\right|$)~\cite{bard57} \begin{equation}
\psi(r) = \frac{k_F}{\pi^2 r}\frac{\Delta}{\hbar v_F} \sin(k_F r)\, K_0\left(\frac{r}{\pi \xi_{BCS}}\right)\; \stackrel{r\rightarrow \infty}{\sim}\;\sin\left(k_F r\right)\, e^{-r/(\pi \xi_{BCS})}
\end{equation} where $K_0(k r)$ is the modified Bessel function that falls off as $e^{-k r}$ at infinity. We have encountered a similar exponential envelope function for a two-body bound state (see table~\ref{t:boundstate}). The characteristic size of the Cooper pair, or the {\it two-particle correlation length} $\xi_0$, can be defined as $\xi_0^2 =\frac{\left<\psi(\vect{r})\right|r^2\left|\psi(\vect{r})\right>}{\left<\psi(\vect{r})|\psi(\vect{r})\right>}$, and this gives indeed $\xi_0 \sim 1/\delta k$, \begin{equation}
\xi_0 \approx \xi_{BCS} \equiv \frac{\hbar v_F}{\pi \Delta} \gg 1/k_F \quad \mbox{in the BCS-limit}
\end{equation}
In the BEC limit, $u_k v_k \propto \frac{1}{1+(k a)^2}$, and so \begin{equation}
\psi(\vect{r}_1,\vect{r}_2) \sim \frac{e^{- \left|\vect{r}_1 - \vect{r}_2\right|/a}}{\left|\vect{r}_1 - \vect{r}_2\right|}
\end{equation} which is simply the wave function of a molecule of size $\sim a$ (see table~\ref{t:boundstate}). The two-particle correlation length\footnote{This length scale should be distinguished from the {\it coherence length} $\xi_{phase}$ 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 $\xi_{phase}$ is given by the healing length $\propto \frac{1}{\sqrt{n a}}$. See~\cite{pist94xi} for a detailed discussion.} is thus $\xi_0 \sim a$. 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
