Critical Temperature in the BCS-BEC Crossover with Spin-Orbit Coupling

Abstract

We review the study of the superfluid phase transition in a system of fermions whose interaction can be tuned continuously along the crossover from Bardeen–Cooper–Schrieffer (BCS) superconducting phase to a Bose–Einstein condensate (BEC), also in the presence of a spin–orbit coupling. Below a critical temperature the system is characterized by an order parameter. Generally a mean field approximation cannot reproduce the correct behavior of the critical temperature $T_c$ over the whole crossover. We analyze the crucial role of quantum fluctuations beyond the mean-field approach useful to find $T_c$ along the crossover in the presence of a spin–orbit coupling, within a path integral approach. A formal and detailed derivation for the set of equations useful to derive $T_c$ is performed in the presence of Rashba, Dresselhaus and Zeeman couplings. In particular in the case of only Rashba coupling, for which the spin–orbit effects are more relevant, the two-body bound state exists for any value of the interaction, namely in the full crossover. As a result the effective masses of the emerging bosonic excitations are finite also in the BCS regime.

1. Introduction

Experimental developments in confining, cooling and controlling the strength of the interaction of alkali atoms brought a lot of attention to the physics of the crossover between two fundamental and paradigmatic many-body systems, the Bose–Einstein condensation (BEC) and the Bardeen–Cooper–Schrieffer (BCS) superconductivity.

A system of weakly attracting fermions can be treated in the context of the well-known BCS theory [1], originally formulated to describe superconductivity in some materials where an effective attractive interaction between electrons can arise from electron-phonon interaction. This theory predicts a phase transition and the formation of the so-called Cooper pairs [2] with the appearance of a gap in the single particle spectrum. The fermions composing these weakly bounded pairs are spatially separated since their correlation length is larger than the mean inter-particle distance; for this reason the pairs cannot be considered as bosons. By increasing the strength of the attraction among the fermions, the system could be seen as made of more localized Cooper pairs, or almost bosonic molecules, which may undergo into a true Bose–Einstein condensation. Since there is no evidence of a broken symmetry by continuously varying the width of the couples, this process realizes a crossover between a BCS and a BEC state. The idea of this crossover preceded by far its experimental realization [3] and became a subject of active investigations by several authors. Those studies were first motivated by the purpose of understanding superconductivity in metals at very low electron densities, since in this diluted regime Cooper pairs are smaller than the mean inter-particle distance and, therefore, can be treated like bosons. In a particularly relevant pioneering work [4] a system of interacting fermions has been investigated to obtain the evolution of the critical temperature of the system along the crossover by changing the strength of the interaction between the BCS and the BEC limits. [5] The discovery of the high-temperature superconductors [6], where the coupling between electrons is supposed to be stronger than that in conventional superconductors, renewed the interest in the BCS-BEC crossover with the aim to explore the regime of stronger interaction. Few years later, an important work [7] provided for the first time a path integral formulation of the problem, which we will review here, describing the crossover at finite temperature with tunable strength of the coupling between pairs. The scientific interest in this field was first purely theoretical because the interaction strength in real solid state materials cannot be tuned easily. The situation changed drastically after the first realizations of magnetically confined ultracold alkali gases. A typical alkali gas obtained in laboratory is made by ${10}^5$ to ${10}^9$ atoms, with a density n inside the bulk of a realistic trap in the range from ${10}^{12}$ to ${10}^{15}$ cm${}^{-3}$. An important property is their diluteness: these densities are such that the mean inter-particle distance is larger than the length scale associated to the atom-atom potential and implying that only two-body physics is relevant. Below those densities, the thermalization of the gas after a perturbation becomes too slow, while at higher densities three body collisions become relevant. At low temperatures, namely at low kinetic energies, the scattering properties of the atoms can be characterized by a single parameter, the scattering length. These gases can be furthermore cooled at such temperatures to enter in a quantum degenerate regime characterized by a thermal de Broglie wavelength of the same order of the interatomic distance, i.e., ${\lambda }_{th}^2=\frac{2\pi {\hslash }^2}{m{k}_{B}T}\sim n^{-2/3}$. The degenerate temperatures, for the densities mentioned before, are smaller than a $\mu K$. These are the conditions achieved in the first experimental realization of a Bose–Einstein condensation [8]. The importance of such experimental setups does not rely only on these extreme conditions, but also on the possibility of fine tuning the parameters which describe the gas. In particular the mechanism of the so-called Feshbach resonance allows one to control the strength of the interaction between two atoms by simply a magnetic field, for example, making it sufficiently strong to support a new bound state. This phenomenon was first predicted and then observed experimentally [9,10] in bosonic systems. Soon after these techniques have been extended to create and manipulate ultracold gases for fermionic systems [11,12,13,14] employing atomic alkali gases with two different components, a mixture of atoms in two different spin or pseudospin states, where the quantum degeneracy is reached when lowering the temperature the energy of the system ceases to depend on the temperature. In this system by Feshbach resonance the interaction among the atoms can be made sufficiently strong to allow the formation of a two-body bound state, called Feshbach molecule, with a lifetime compatible to the time needed to cool the system down and to observe a Bose–Einstein condensation. Later the experimental observation of fermionic pairs was extended to the whole region of the crossover [15,16], studying also the thermodynamical properties [17,18]. A renewed interest in the study of the BCS-BEC crossover arose in 2011 after the experimental realization, in the context of ultracold gases, of a synthetic spin–orbit (SO) interaction [19,20,21,22]. In contrast to the spin–orbit interaction in solid state systems, in atomic gases this coupling is produced by means of laser beams and therefore is perfectly tunable. The realization of this artificial spin–orbit coupling is a bright example of the versatility and controllability of the ultracold techniques, showing once more that ultracold gases are the perfect experimental playgrounds for testing physical models hardly realizable by means of other solid state setups. These experimental achievements stimulated intense theoretical effort to understand the spin–orbit effects with the so-called Rashba [23] and Dresselhaus [24] terms in the evolution from BCS to BEC superfluidity, considering many SO configurations, in two and three dimensions, at zero and finite temperature, exploring eventually novel phases of matter or topological properties [25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67]. In particular SO coupling induces triplet pairing [43] which can produce nodes in the quasiparticle excitation spectrum. This allows for detecting bulk topological phase transitions of the Lifshitz type in polarized systems [41] or, more recently, in heterostructures [68]. Moreover the experimental realization of 2D degenerate Fermi gases for ultracold atoms in a highly anisotropic disk-shaped potential [69] is one of the reasons of the growing interest for fermions in reduced dimensionality, where exotic gapless superfluid states and Fulde–Ferrell–Larkin–Ovchinnikov (FFLO) pairing states can be observed [55].

The aim of this paper is to review some basic concepts on the BCS-BEC crossover, whose properties will require taking into account quantum fluctuations in order to get a better theoretical predictions when comparing with experimental results [70]. We will reproduce many results beyond the mean-field level, and extend the study in the presence of a spin–orbit interaction, in a detailed and comprensive way, by means of a path integral approach, with a particular emphasis to the derivation of the critical temperature $T_c$. We will show how $T_c$ can be strongly enhanced by the SO couplings in the BCS side, already at the mean-field level, while its increase is softened by the effects of the quantum fluctuations specially in the intermediate region of the crossover and in the strong coupling regime.

1.1. Basic Concepts in Ultracold Atomic Physics

To understand the structure of an atom it is necessary to consider the forces among its constituents. For alkali atoms (Li, Na, K, Rb...), the relevant force is due to Coulomb interaction that acts between the off-shell electron and the effective charge $Ze$ of the nucleus. This central force produces the following Bohr energy spectrum

$$E_n=-\frac{m_r{c}^2}{2}\frac{Z^2{\alpha }^2}{n^2}$$

(1)

with $n=1,2,\dots$, where $\alpha =\frac{e^2}{4\pi {ϵ}_0\hslash c}\approx \frac{1}{137}$ is the so-called fine structure constant and $m_r\approx m_e$ the reduced mass, approximately given by the mass of the electron. This spectrum is degenerate in the quantum numbers L and $L_z$, associated to the angular momentum operators ${\widehat{L}}^2$ and ${\widehat{L}}_z$.

The Coulomb interaction, due to its symmetry, preserves the angular momentum and the spin separately. However, in principle, they are coupled via a spin–orbit term, ${\widehat{H}}_{so}=\frac{{\alpha }_{so}}{{\hslash }^2}\widehat{L}·\widehat{S}$, therefore, only the total angular momentum $\widehat{J}=\widehat{L}+\widehat{S}$ is conserved. This additional term causes the known correction to the spectrum, called fine structure, shifting the degeneracy in the energy spectrum by $\Delta E_J=\frac{{\alpha }_{so}}{2}\left[J(J+1)-L(L+1)-S(S+1)\right]$. However, for ultracold alkali gases, at very low temperatures, only the fundamental state is occupied, so that the off-shell electron is in the s orbital ($L=0$). In this state $J=S$ and the spin–orbit fine structure can be neglected, $\Delta E_J=0$. However, also the coupling between the spins of the electrons $\widehat{S}$ and of the nucleus $\widehat{I}$ can remove the degeneracies and it is generally approximated by a dipole–dipole interaction

$${\widehat{H}}_{hf}=\frac{{\alpha }_h}{{\hslash }^2}\widehat{I}·\widehat{S}$$

(2)

This term, again, does not lead to the conservation of both $\widehat{I}$ and $\widehat{S}$ separately but only of the total spin of the atom $\widehat{F}=\widehat{I}+\widehat{S}$, since $L=0$. Then it is possible to label the states of an alkali atom by the quantum numbers F and $F_z$, related to the operators ${\widehat{F}}^2$ and ${\widehat{F}}_z$. The effect of this contribution is a small correction to the Bohr spectrum, usually called hyperfine structure

$$\Delta E_F=\frac{{\alpha }_h}{2}\left[F(F+1)-I(I+1)-S(S+1)\right].$$

(3)

Since $S=\frac{1}{2}$ and because of the properties of the sums of angular momentums, the quantum number F is allowed only to have the values $F=I\pm \frac{1}{2}$, so the hyperfine correction is given by

$$\Delta E_{F=I\pm \frac{1}{2}}=\pm \frac{{\alpha }_h}{2}\left(I+\frac{1}{2}\mp \frac{1}{2}\right)$$

(4)

In the presence of an homogeneous magnetic field, this coupling between the latter and the electronic and nuclear spins gives rise to a shift of the energy levels of the atoms, called the Zeeman effect. An alkali atom at very low energy and in a magnetic field can be described, therefore, by the following Hamiltonian

$$\widehat{H}={\widehat{H}}_{hf}-{\gamma }_{N}B{\widehat{I}}_z+{\gamma }_{E}B{\widehat{S}}_z$$

(5)

where the Zeeman couplings are inversely proportional to the mass of the particles, therefore ${\gamma }_N\ll {\gamma }_E$, since the nucleus is heavier than the electrons. As long as the effect of the magnetic field is small, F and $F_z$ can be considered almost good quantum numbers, while in general the eigenstates are labeled by $I_z$ e $S_z$, associated to ${\widehat{I}}_z$ and ${\widehat{S}}_z$. Diagonalizing the Hamiltonian in Equation (5) we can find the corrections to the Bohr spectrum. For example, in the case of an hydrogen atom $I=\frac{1}{2}$ and we get

$$\Delta E_{\pm }=-\frac{{\alpha }_h}{4}\pm \frac{{\alpha }_h}{2}\sqrt{1+{\left(\frac{\hslash ({\gamma }_N+{\gamma }_E)B}{{\alpha }_h}\right)}^2}.$$

(6)

The dependence on the magnetic field is quadratic for low intensities and linear for high field. An important experimental feature is that varying the magnetic field it is possible to vary the energy levels of the single atom, a property that can be exploited to induce a confining potential for atoms, realizing an atomic trap.

There are two main kinds of atomic traps available, the magnetic and the optical dipole traps, which are briefly discussed below.

• Magnetic Traps

This technique is based on the concrete possibility of creating position dependent magnetic fields which, as a consequence, make every hyperfine energies also position dependent. In these conditions an alkali atom has a kinetic energy plus an effective potential, that, properly taylored, can be seen as a confining potential for the gas,

$$E_{\mathrm{atom}}\left(\mathbf{r}\right)=\frac{{\hslash }^2{p}^2}{2m}+V\left(B\left(\mathbf{r}\right)\right)$$

(7)

with a confining force $\mathbf{F}\left(\mathbf{r}\right)=-\nabla V\left(\mathbf{r}\right)$. Trapping a gas in a magnetic field is the key to achieve degeneracy temperatures; in these conditions the gas is isolated from the external world getting rid of the main source of heating which is due to the walls of the container. Typically the magnetic trap is modeled to obtain an harmonic potential, however there are several experimental possibilities depending only on the ability to create non homogeneous magnetic fields. In particular these potentials highly depend on the hyperfine states of the atom and this characteristic can be exploited to select particular states in the trap. Another interesting possibility is to create highly anisotropic potentials by strong confinement in some dimensions, realizing one- or two-dimensional systems. This technique is compatible with the cooling methods of atomic gases but not with the mechanism of the magnetic Feshbach resonance. For the latter reason an alternative trapping method is required.

• Optical Traps

The optical traps are based upon the interaction between an atom and light of a laser beams. The key features of this method are the great stability of the trapped gas and the possibility of making a lot of different geometries. Under proper conditions, called far-detuning, the confinement is independent from the internal state of the atom. When an atom is put in a laser beam, the electric field $\mathbf{E}\left(\mathbf{r}\right)={\mathbf{E}}_0\left(\mathbf{r}\right)e^{i\omega t}$ induces an electric dipole moment $\mathbf{p}\left(\mathbf{r}\right)={\alpha }_{\omega }{\mathbf{E}}_0\left(\mathbf{r}\right)e^{i\omega t}$, where ${\alpha }_{\omega }$ is the polarizability. The effective potential given by the iteration between the dipole moment and the electric field is given by a temporal mean value

$$V\left(\mathbf{r}\right)=-\frac{1}{2}\overline{\mathbf{p}\left(\mathbf{r}\right)·\mathbf{E}\left(\mathbf{r}\right)}$$

(8)

which is simply directly proportional to the intensity $\mathrm{I}\left(\mathbf{r}\right)$ of the beam. It is then possible to produce an artificial potential, as before, considering a position dependent focalization of the laser so that there is a force acting on the single atom given by $\mathbf{F}\left(\mathbf{r}\right)=-\nabla V\left(\mathbf{r}\right)\propto \nabla I\left(\mathbf{r}\right)$. As mentioned, this technique makes possible to easily create different geometries. Indeed, employing two counter-propagating laser in more dimensions for example, lattice potentials can be realized, particularly interesting since they can simulate solid state systems in a fully controllable environment. For a more complete review on this confining technique one can refer to Refs. [71,72].

The confined gas has to be cooled down in order to reach the degeneracy temperatures. This can be done in two steps, by laser cooling and evaporative cooling [73].

• Laser or Doppler Cooling

The term laser cooling refers to several techniques employed to cool atomic gases, taking advantage of the interaction between light and matter. The most relevant one takes advantage also from the Doppler effect. This technique is applied to slow down the mean velocity of the atoms realizing a situation in which an atom absorbs more radiation in a direction and less in the opposite one. This can be done slightly detuning the frequency of the laser below the resonance of the exited state of the atom. Because of the Doppler effect, the atom will see a frequency closer to the resonant one and so it will have an higher absorption probability. In the opposite direction the effect will be reversed. The atom, which is now in an exited state, can spontaneously emit a photon. The overall result of the absorption and emission process is to reduce the momentum of the atom, therefore its speed. If the absorption and emission are repeated many times, the average speed, and therefore the kinetic energy of the atoms, will be reduced, and so the temperature. By this technique one can reach low temperatures of the order of some $\mu K$, still not enough to get a degenerate gas.

• Evaporative Cooling

The degenerate temperatures can be reached by the so-called forced evaporative cooling method. A fundamental feature of every trap discussed before is the possibility of varying its depth, defined as the energy difference between the bottom and the scattering threshold of the potential. Lowering the depth of the trap it is possible to let the most energetic atoms escape; this correspond to extracting energy from the gas at the cost of loosing a considerable quantity of matter. After the thermalization of the gas, the temperature is lowered and this procedure can be repeated. The lower limit in temperature now is related to the rate between inelastic collisions, that fix the lifetime of the sample, and the elastic collisions, that assure thermalization and fix the time needed for the cooling process. Loading more atoms in a trap lowers the limit of the lowest temperature achievable but also makes the gas unstable. With the correct compromise this scheme allows to go well below the degenerate temperatures. For alkali atoms the limits in temperatures achievable by evaporative cooling are at the order of T∼pK [74]. Once a degenerate gas is obtained, letting the gas expand in free flight by turning off the trap, from the absorption spectrum one can access to the temperature of the gas or to the density profile in position space that can be compared with theoretical results.

1.2. Artificial Spin-Orbit Interaction in Ultracold Gases

Laser techniques can be employed also to build a tunable synthetic spin–orbit coupling between two pseudo-spin states of an atom. In a previous section it was shown that spin–orbit coupling (fine structure) for atoms in the fundamental state ($n=0$, $L=0$) vanishes. In order to induce it in laboratory it is necessary to select a gas with two components. In the first realization of this system [19] two states ($|F,F_z\rangle$) of ${}^{87}$Rb were chosen with pseudo-spins $|\uparrow \rangle =|1,0\rangle$ and $|\downarrow \rangle =|1,-1\rangle$. In order to understand schematically how it works, let us suppose that the atom is in an initial state $|\downarrow \rangle$, then, two lasers can be tuned in such a way that one induces a transition to an intermediate exited state and the other induces a stimulated emission from the latter state to $|\uparrow \rangle$. The inverse process is equally allowed. Due to the Doppler effect, this coupling will be necessarily momentum-sensitive and, as shown in Refs. [19,20], the effective generated term is equivalent to a spin–orbit coupling of a spinful particle moving in a static electric field $\mathbf{E}$, let us say, along the z axis. An effective additional term in the Hamiltonian is $-{\mu }_s·{\mathbf{B}}_{SO}$, where ${\mathbf{B}}_{SO}=\mathbf{v}\times \mathbf{E}$ is the magnetic field seen by the moving particle and ${\mu }_s$ the electron magnetic moment, which can be written as

$$H_R=v_R({\sigma }_x{k}_y-{\sigma }_y{k}_x)$$

(9)

This is called Rashba coupling, in analogy to the same coupling occurring in crystals. With the same technique other interactions can be produced in the laboratory such as the Dresselhaus term

$$H_D=v_D({\sigma }_x{k}_y+{\sigma }_y{k}_x)$$

(10)

or a Weyl term $H_W=\lambda \sigma ·\mathbf{k}$. The two lasers are detuned by a frequency $2\delta$ from the Raman resonance and couple the states $|\uparrow \rangle$ and $|\downarrow \rangle$ with the Raman strength $2\mathsf{\Omega }$. This generates additional terms, $\mathsf{\Omega }{\sigma }_z$ and $\delta {\sigma }_y$. We will neglect the detunig term and will consider the Raman coupling as an effective artificial Zeeman term

$$H_Z=h{\sigma }_z.$$

(11)

By this procedure one can study the effects of a spin–orbit interaction in ultracold gases, which has the great advantage of being fully under control, allowing for fine tuning experiments where the couplings can be easily manipulated, not possible otherwise in solid state systems.

2. Two-Body Scattering Problem

In the previous section we said that only two-body physics is relevant in a dilute ultracold gas, for that reason we will briefly review the two body problem in a generic central potential. We are going to consider only scattering properties after two-body collisions neglecting the three-body ones because of the low probability of finding three particles within the range of interaction. Generally interatomic potentials are not known analytically and their approximations are not so easy to use in theoretical and numerical evaluations. The density of the gas can fix the interparticle distances and the range of the potential. This is what typically happens in the experiments with ultracold dilute gases. In this situation an approximate potential, or pseudo-potential, can be employed, which should exhibit some general properties of the full scattering problem. A detailed analysis of the scattering problems, specifically for ultracold gases, can be found in Refs. [75,76,77]. Scattering theory provides the tools useful to characterize the interaction between particles in terms of few parameters that can be used to create suitable approximate potentials. The Schrödinger equation is employed to determine the wavefunctions and the scattering amplitudes, then we require that a simple pseudo-potential should reproduce the same results in a low energy limit, namely in the ultracold limit. We will assume that the two-body interaction is described by a short-range and spherically symmetric potential.

It is known that Schrödinger equation for the wavefunction of two colliding particles can be decomposed in one equation that describes the center-of-mass motion, a free particle equation with mass $M=m_1+m_2$, and another that describes the relative motion. With the definition $m=\frac{m_1{m}_2}{m_1+m_2}$ for the reduced mass, in the relative frame the Schödinger equation reads

$$\left[-\frac{{\hslash }^2{\nabla }^2}{m}+V\left(\mathbf{r}\right)\right]\psi \left(\mathbf{r}\right)=\left[{\widehat{H}}_0+V\left(\mathbf{r}\right)\right]\psi \left(\mathbf{r}\right)=E\psi \left(\mathbf{r}\right).$$

(12)

Let us suppose that the potential is short-ranged, so that there exists a characteristic length $r^{*}$ such that the potential can be neglected outside this radius,

$$V\left(\mathbf{r}\right)\approx 0,\forall |\mathbf{r}|>r^{*}$$

(13)

In this case Equation (12) becomes a free Schrödinger equation whose solution is given by the composition of an incoming plane wave with momentum $\mathbf{p}$ and a scattering state with momentum ${\mathbf{p}}^{\prime }$. It is relevant to look for elastic collisions, that can be studied fixing the scattering energy to $E=2{ϵ}_{\mathbf{p}}=\frac{p^2}{m}$, equal to the incoming wave energy. The solution is given by

$$\left[2{ϵ}_{\mathbf{p}}-{\widehat{H}}_0\right]|{\psi }_{\mathbf{p}}\rangle =\widehat{V}|{\psi }_{\mathbf{p}}\rangle .$$

(14)

The homogeneous solution ($\widehat{V}=0$) is simply $|{\psi }_{\mathbf{p}}\rangle =|\mathbf{p}\rangle$, therefore, it is possible to write an implicit solution, known as the Lippmann–Schwinger equation

$$|{\psi }_{\mathbf{p}}\rangle =|\mathbf{p}\rangle +\frac{1}{2{ϵ}_{\mathbf{p}}-{\widehat{H}}_0+iϵ}\widehat{V}|{\psi }_{\mathbf{p}}\rangle$$

(15)

with $ϵ$ an infinitesimally small positive real number. In the coordinate representation, it reads

$$\begin{array}{ccc}\hfill {\psi }_{\mathbf{p}}\left(\mathbf{r}\right)& =& \langle \mathbf{r}|\mathbf{p}\rangle +\langle \mathbf{r}|\frac{1}{2{ϵ}_{\mathbf{p}}-{\widehat{H}}_0+iϵ}\widehat{V}|{\psi }_{\mathbf{p}}\rangle \hfill \\ & =& \frac{e^{\frac{i}{\hslash }\mathbf{r}·\mathbf{p}}}{{\left(2\pi \hslash \right)}^{3/2}}-\frac{m}{4\pi {\hslash }^2}\int d^3{r}^{\prime }\frac{e^{\frac{i}{\hslash }p|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}\langle {\mathbf{r}}^{\prime }|\widehat{V}|{\psi }_{\mathbf{p}}\rangle \hfill \end{array}$$

(16)

At long distances, (see Appendix A for details) Equation (16) can be written as

$${\psi }_{\mathbf{p}}\left(\mathbf{r}\right)=\frac{1}{{\left(2\pi \hslash \right)}^{3/2}}\left\{e^{\frac{i}{\hslash }\mathbf{r}·\mathbf{p}}+f(\mathbf{q},\mathbf{p})\frac{e^{\frac{i}{\hslash }pr}}{r}\right\}+O\left(r^{-2}\right)$$

(17)

after defining $\mathbf{q}\equiv p\widehat{\mathbf{r}}=\frac{p}{r}\mathbf{r}$, and the scattering amplitude

$$f(\mathbf{q},\mathbf{p})=-2{\pi }^2\hslash m\langle \mathbf{q}|\widehat{V}|{\psi }_{\mathbf{p}}\rangle .$$

(18)

The scattering solution at large distances is then composed by an incoming plane wave and an outgoing spherical wave, weighted by the scattering amplitude. In terms of the latter it is also possible to define the differential cross section of the process, at large distances,

$$\frac{d\sigma }{d\mathsf{\Omega }}={|f(\mathbf{q},\mathbf{p})|}^2.$$

(19)

Due to the spherically symmetric potential, the scattering amplitude depends only on the modulus p and the angle $\theta =arccos\left(\frac{\mathbf{q}·\mathbf{p}}{p^2}\right)=arccos(\widehat{\mathbf{r}}·\widehat{\mathbf{p}})$. It is convenient, therefore, to perform a spherical wave expansion

$$f(\mathbf{q},\mathbf{p})=f(p,\theta )=\sum _{\ell =0}^{\infty }(2\ell +1)f_{\ell }\left(p\right)P_{\ell }(cos\theta )$$

(20)

with $P_{\ell }$ the Legendre polynomials. In order to preserve the normalization of the wavefunction, one has to impose the following condition for the coefficients $f_{\ell }\left(p\right)$ (see Appendix B)

$$|1+2ip{f}_{\ell }\left(p\right)|=1\Rightarrow 1+2ip{f}_{\ell }\left(p\right)\equiv e^{2i{\delta }_{\ell }\left(p\right)}$$

(21)

which defines the phase shift ${\delta }_{\ell }\left(p\right)$ and where we rescaled, here and in what follows, $p\to \hslash p$ for simplicity, to get rid of ℏ. This relation imposes that when the relative momentum vanishes, also the phase shift should be zero, ${\delta }_{\ell }\left(0\right)=0$. The scattering amplitude, at fixed angular momentum, can be rewritten as

$$f_{\ell }\left(p\right)=\frac{1}{pcot\left[{\delta }_{\ell }\left(p\right)\right]-ip}$$

(22)

In the radial Schrödinger equation (see Appendix C), for $\ell \ne 0$, besides the potential $V\left(r\right)$ there is also the so-called centrifugal barrier, that depends on the angular momentum at which the scattering happens and on the relative distance. In ultracold atoms, at sufficiently low temperatures, the atoms may not have enough energy to overcome this centrifugal barrier, therefore the only relevant contribution to the scattering amplitude will come from s-wave scattering ($\ell =0$). To clarify this point one can estimate the angular momentum as the product of the range of the interaction $r^{*}$ times the momentum given by the inverse of the thermal de Broglie length $\ell \approx \frac{r^{*}}{{\Lambda }_{th}}$. For ultracold gases the temperatures are typically $T\sim nK$ and, since the thermal wavelength ${\Lambda }_{th}\sim T^{-1/2}$, it is allowed to consider $\ell \approx 0$, namely, the scattering amplitude is dominated by the s-wave contribution

$$f(\mathbf{q},\mathbf{p})\approx f_0\left(p\right)=\frac{1}{pcot\left[{\delta }_0\left(p\right)\right]-ip}$$

(23)

It is expected that scattering processes modified by the presence of $V\left(r\right)$ for $\ell \ne 0$ are greatly suppressed at sufficiently low energy. Moreover, also the $\ell =0$ contribution is affected by the low energy hypothesis. As shown in Ref. [76], it is possible to define the scattering length at fixed ℓ considering the low energy limit

$$f_{\ell }\left(p\right)\underset{p\to 0}{\sim }-a_{\ell }{p}^{2\ell }$$

(24)

which comes from the low energy limit of the phase shift

$${\delta }_{\ell }\left(p\right)\underset{p\to 0}{\sim }{p}^{2\ell +1}$$

(25)

This can be derived by solving the Schrödinger equation as shown in Appendix C. It is remarkable that ultracold gases can be cooled down to the point where only one partial wave in the two-body problem becomes dominant, the s-wave contribution. At low energies ($p\sim 1/\Lambda$) it is possible to expand the phase ${\delta }_0\left(p\right)=-a_{s}p+O\left(p^2\right)$, so defining the scattering length as

$$a_s\equiv -\underset{p\to 0}{lim}\frac{{\delta }_0\left(p\right)}{p}$$

(26)

Actually only $a_s=a_{\ell =0}$ has the dimension of a length, contrary to the other terms $a_{\ell }$ with $\ell >0$. Considering the relation in Equation (22) and the following expansion

$$pcot\left({\delta }_0\left(p\right)\right)\approx -\frac{1}{a_s}+\frac{r_{\mathrm{eff}}}{2}{p}^2$$

(27)

the scattering amplitude is then given, at low energy, by the expression

$$f_0\left(p\right)\approx -\frac{a_s}{1-a_s\frac{r_{\mathrm{eff}}}{2}{p}^2+i{a}_{s}p}$$

(28)

which shows explicitly that for ultracold gases, under the hypothesis of short-range central potential, at low energy, the two-body scattering depends only on two parameters: $a_s$, the s-wave scattering length, and $r_{\mathrm{eff}}$, an effective range proportional to $r^{*}$. We will not treat $r_{\mathrm{eff}}$, the effective range of the potential, more deeply here, because, for what follows, only the first order expansion of the phase shift is needed. Finally also the cross section admits a low energy limit. Writing $\sigma \left(p\right)={\sum }_{\ell }{\sigma }_{\ell }\left(p\right)$, we have ${\sigma }_{\ell \ne 0}\left(p\right)=\frac{8\pi }{p^2}(2\ell +1){sin}^2\left({\delta }_{\ell }\left(p\right)\right)\underset{p\to 0}{\sim }8\pi (2\ell +1)p^{4\ell }$ and ${\sigma }_{\ell =0}\left(p\right)\underset{p\to 0}{\sim }8\pi a_s^2$. In particular, from Equations (A15) and (A17), one gets $\sigma \left(p\right)\underset{p\to 0}{\sim }8\pi a_s^2$ for spinless bosons, $\sigma \left(p\right)\underset{p\to 0}{\sim }{p}^4$ for spinless (polarized) fermions, ${\sigma }_S\left(p\right)\underset{p\to 0}{\sim }8\pi a_s^2$ for fermions in a singlet state, ${\sigma }_{T\left(p\right)}\underset{p\to 0}{\sim }{p}^4$ for fermions in a triplet state. These results show that the interaction between polarized fermions is suppressed, while interaction in the singlet channel is dominant in the low energy limit, consistently with the Pauli exclusion principle.

2.1. Two-Body Scattering Matrix

It is possible to mimic a realistic potential using the pseudo-potential with the same scattering properties a low energies, namely the same s-wave scattering length. Let us introduce the operator ${\widehat{T}}_{2B}$, called T scattering matrix,

$$\widehat{V}|{\psi }_{\mathbf{p}}\rangle ={\widehat{T}}_{2B}|\mathbf{p}\rangle$$

(29)

By this definition it is possible to rewrite the scattering amplitude as a function of ${\widehat{T}}_{2B}$

$$f(\mathbf{q},\mathbf{p})=-2{\pi }^{2}m\langle \mathbf{q}|\widehat{V}|{\psi }_{\mathbf{p}}\rangle =-2{\pi }^{2}m\langle \mathbf{q}|{\widehat{T}}_{2B}|\mathbf{p}\rangle .$$

(30)

In the same way also the Lippmann–Schwinger equation, introduced earlier, can be written as

$${\widehat{T}}_{2B}\left(z\right)=\widehat{V}+\widehat{V}\frac{1}{z-{\widehat{H}}_0+iϵ}{\widehat{T}}_{2B}\left(z\right)$$

(31)

whose solution, in terms of ${\widehat{T}}_{2B}$, is the following

$${\widehat{T}}_{2B}\left(z\right)=\widehat{V}+\widehat{V}\frac{1}{z-\widehat{H}+iϵ}\widehat{V}.$$

(32)

This solution can be expressed as an expansion in the potential $\widehat{V}$

$${\widehat{T}}_{2B}\left(z\right)=\widehat{V}\sum _{n=0}^{\infty }{\left[\frac{1}{z-{\widehat{H}}_0+iϵ}\widehat{V}\right]}^n.$$

(33)

Given a complete set of eigenstates of $\widehat{H}$, inserting a completeness relation, one can write ${\widehat{T}}_{2B}\left(z\right)=\widehat{V}+{\sum }_{\alpha }\widehat{V}\frac{|{\psi }_{\alpha }\rangle \langle {\psi }_{\alpha }|}{z-{ϵ}_{\alpha }}\widehat{V}+\frac{\mathsf{\Omega }}{{\left(2\pi \right)}^3}\int d^{3}p\widehat{V}\frac{|{\psi }_{\mathbf{p}}\rangle \langle {\psi }_{\mathbf{p}}|}{z-2{ϵ}_{\mathbf{p}}}\widehat{V}$, where ${ϵ}_{\alpha }<0$ are the energies of the bound states, $|{\psi }_{\mathbf{p}}\rangle$ the scattering states and $\mathsf{\Omega }$ is a volume. It is then clear that the two-body scattering matrix $T_{2B}\left(z\right)$ has simple poles associated to bound states and a branch cut on the real axis caused by the continuum of the scattering states. Indeed the transition amplitude given by Equation (28) can be also rewritten in terms of the scattering energy $E=\frac{{\hslash }^2{p}^2}{m}$ so that

$$T_{2B}\left(E\right)\approx \frac{1}{2{\pi }^{2}m}\frac{a_s}{1-i{a}_s\sqrt{\frac{mE}{{\hslash }^2}}-\frac{m{a}_s{r}_{\mathrm{eff}}}{2{\hslash }^2}E}$$

(34)

At low energies ($E\to 0^{-}$) it is found that, for positive scattering length, $a_s>0$, this quantity has a simple pole at

$$E_b=-\frac{{\hslash }^2}{m{a}_s^2}.$$

(35)

This pole indicates the presence of a bound state close to the scattering threshold ($E=0$) and gives an additional meaning to the scattering length $a_s$.

As a final remark, by the Fano–Feshbach resonance method one can control by a magnetic field the value of $a_s$. We refer to Refs. [11,12,13,14] for a detailed description of this technique.

2.2. Renormalization of the Contact Potential

Let us consider the contact potential in the position representation

$${\widehat{V}}_{\mathrm{eff}}\left(\mathbf{r}\right)=-g{\delta }^{\left(3\right)}\left(\mathbf{r}\right)$$

(36)

The action of this operator on the eigenstates of the momentum operator is

$$\langle \mathbf{q}|{\widehat{V}}_{\mathrm{eff}}|\mathbf{p}\rangle =\langle \mathbf{q}|\int d^{3}r{\widehat{V}}_{\mathrm{eff}}|\mathbf{r}\rangle \langle \mathbf{r}|\mathbf{p}\rangle =-g\langle \mathbf{q}|0\rangle \langle 0|\mathbf{p}\rangle =-\frac{g}{{\left(2\pi \right)}^3}$$

(37)

We will want this potential to reproduce the scattering amplitude of a realistic potential. Exploiting the Lippmann–Schwinger equation for $T_{2B}$, written in Equation (31), we have

$$\langle \mathbf{q}|{\widehat{T}}_{2B}|\mathbf{p}\rangle =\langle \mathbf{q}|{\widehat{V}}_{\mathrm{eff}}|\mathbf{p}\rangle +\int d^{3}k\langle \mathbf{q}|{\widehat{V}}_{\mathrm{eff}}\frac{1}{E-{\widehat{H}}_0+iϵ}|\mathbf{k}\rangle \langle \mathbf{k}|{\widehat{T}}_{2B}|\mathbf{p}\rangle$$

(38)

We now impose that, for the contact potential, the zero energy scattering amplitude should give the same result of a realistic potential. From Equation (34) we have

$$T_{2B}(E=0)=\frac{a_s}{2{\pi }^{2}m}$$

(39)

while from Equation (38), using Equation (37), we get

$$\frac{1}{g}=-\frac{m}{4\pi a_s}+\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\frac{1}{2{ϵ}_{\mathbf{k}}}$$

(40)

This relation may seem inconsistent at first sight due to the divergence of the last term. However, as we will see in what follows, this result will turn out to be very helpful since this divergence cancels out another divergence appearing in the so-called gap equation, for the order parameter, resulting as an elegant renormalization procedure for the contact potential.

2.3. With Spin-Orbit Coupling

As shown in Refs. [44,45], in the presence of a spin–orbit coupling, the equation useful to cure the ultraviolet divergences is formally the same as Equation (40), (see Equation (A82)), with the only difference that $a_s$ is replaced by the scattering length $a_r$ which contains also a dependence on the spin–orbit parameters (see Appendix D for more details).

3. Fermi Gas with Attractive Potential

In the previous introductory sections we showed that the extreme regimes of ultracold gases, namely the diluteness and the low temperatures, allow us to consider only the two-body physics and simplify the real interaction replacing it with a contact potential that is able to reproduce the same scattering properties. For these reasons this kind of systems can be described by a set of fermions in three dimensions with a BCS-like Hamiltonian

$$H=H_0+H_I=\int d^{3}r\left\{\sum _{\sigma }{\psi }_{\sigma }^{†}\left(\overrightarrow{r}\right)\left(-\frac{{\nabla }^2}{2m}\right){\psi }_{\sigma }\left(\overrightarrow{r}\right)-g{\psi }_{\uparrow }^{†}\left(\overrightarrow{r}\right){\psi }_{\downarrow }^{†}\left(\overrightarrow{r}\right){\psi }_{\downarrow }\left(\overrightarrow{r}\right){\psi }_{\uparrow }\left(\overrightarrow{r}\right)\right\}$$

(41)

The interaction coupling is positive ($g>0$) for an attractive potential. In the following section we will review the study of this model by a path integral approach, focusing on the evaluation of the critical temperature as a function of the coupling g, along the crossover from weak coupling (BCS regime) to strong coupling (BEC regime). The Hamiltonian in Equation (41) is associated to the grand canonical partition function in the path integral formalism which is

$$\mathcal{Z}=\mathrm{Tr}\left(e^{-\beta (H-\mu N)}\right)=\int D{\psi }^{†}D\psi e^{-S[{\psi }^{†},\psi ]}$$

(42)

where the action is give by

$$S[{\psi }^{†},\psi ]=\int d^{3}r{\int }_0^{\beta }d\tau \left\{\sum _{\sigma }{\psi }_{\sigma }^{†}\left(\overrightarrow{r}\right)\left({\partial }_{\tau }-\frac{{\nabla }^2}{2m}-\mu \right){\psi }_{\sigma }\left(\overrightarrow{r}\right)-g{\psi }_{\uparrow }^{†}\left(\overrightarrow{r}\right){\psi }_{\downarrow }^{†}\left(\overrightarrow{r}\right){\psi }_{\downarrow }\left(\overrightarrow{r}\right){\psi }_{\uparrow }\left(\overrightarrow{r}\right)\right\}$$

(43)

It is possible to take advantage of the properties of Gaussian integrals (see for instance Ref. [77]) to handle the interacting term introducing the complex auxiliary field $\Delta (\mathbf{r},\tau )$

$$e^{g\int d^{3}r{\int }_0^{\beta }d\tau {\psi }_{\uparrow }^{†}{\psi }_{\downarrow }^{†}{\psi }_{\downarrow }{\psi }_{\uparrow }}\propto \int D{\Delta }^{*}D\Delta e^{\int d^{3}r{\int }_0^{\beta }d\tau \left[-\frac{{|\Delta |}^2}{g}+{\Delta }^{*}{\psi }_{\downarrow }{\psi }_{\uparrow }+\Delta {\psi }_{\uparrow }^{†}{\psi }_{\downarrow }^{†}\right]}$$

(44)

The grand canonical partition function can be rewritten, therefore, as it follows

$$\mathcal{Z}=\int D{\psi }^{†}D\psi \int D{\Delta }^{*}D\Delta e^{-S[{\psi }^{†},\psi ,\Delta ]}$$

(45)

with the use of a new action, dependent now also on an auxiliary field

$$S[{\psi }^{†},\psi ,\Delta ]=\int d^{3}r{\int }_0^{\beta }d\tau \left\{\sum _{\sigma }{\psi }_{\sigma }^{†}\left({\partial }_{\tau }-\frac{{\nabla }^2}{2m}-\mu \right){\psi }_{\sigma }+\frac{{|\Delta |}^2}{g}-{\Delta }^{*}{\psi }_{\downarrow }{\psi }_{\uparrow }-\Delta {\psi }_{\uparrow }^{†}{\psi }_{\downarrow }^{†}\right\}$$

(46)

The physical meaning of $\Delta$ will be clear from the saddle point equation. It is more convenient to introduce the Nambu spinors

$$\mathsf{\Psi }=\left(\begin{array}{c}{\psi }_{\uparrow }\\ {\psi }_{\downarrow }^{†}\end{array}\right)\overline{\mathsf{\Psi }}=\left({\psi }_{\uparrow }^{†},{\psi }_{\downarrow }\right)$$

(47)

so that the action of the system can be written as

$$S[\overline{\mathsf{\Psi }},\mathsf{\Psi },\Delta ]=\int d^{3}r{\int }_0^{\beta }d\tau \left[\frac{{|\Delta |}^2}{g}+\overline{\mathsf{\Psi }}{\mathcal{G}}^{-1}\mathsf{\Psi }\right]+\beta \sum _{\mathbf{k}}{\xi }_{\mathbf{k}}$$

(48)

with ${\mathcal{G}}^{-1}$, the inverse of the interacting Green function, defined as

$${\mathcal{G}}^{-1}=\left(\begin{array}{cc}-{\partial }_{\tau }+\frac{{\nabla }^2}{2m}+\mu & \Delta \\ {\Delta }^{*}& -{\partial }_{\tau }-\frac{{\nabla }^2}{2m}-\mu \end{array}\right)$$

(49)

Actually in order to put $S[{\psi }^{†},\psi ,\Delta ]$ in a quadratic form in the spinor representation we exploited the anticommutation relations of the Grassmann fields, getting also the last term in Equation (48). Defining $z=e^{-\beta {\sum }_{\mathbf{k}}{\xi }_{\mathbf{k}}}$, this term can be put outside the functional integral. The grand canonical partition function is the following

$$\mathcal{Z}=z\int D{\Delta }^{*}D\Delta e^{-\int d^{3}r{\int }_0^{\beta }d\tau \frac{{\Delta }^2}{g}}\int D\overline{\mathsf{\Psi }}D\mathsf{\Psi }{e}^{-\int d^{3}r{\int }_0^{\beta }d\tau \overline{\mathsf{\Psi }}{\mathcal{G}}^{-1}\mathsf{\Psi }}$$

(50)

Performing the standard Gaussian integral over the fermionic degrees of freedom one gets

$$\mathcal{Z}=z\int D{\Delta }^{*}D\Delta e^{-\int d^{3}r{\int }_0^{\beta }d\tau \frac{{|\Delta |}^2}{g}+\mathrm{Tr}[ln\left({\mathcal{G}}^{-1}\right)]}\equiv \int D{\Delta }^{*}D\Delta e^{-S_e[\Delta ]}$$

(51)

which defines the effective action depending only on the auxiliary field

$$S_e[\Delta ]=\int d^{3}r{\int }_0^{\beta }d\tau \frac{{\left|\Delta (\overrightarrow{r},\tau )\right|}^2}{g}-\mathrm{Tr}[ln\left({\mathcal{G}}^{-1}\right)]+\beta \sum _{\mathbf{k}}{\xi }_{\mathbf{k}}.$$

(52)

3.1. Gap Equation

Starting from this action it is possible to write down the saddle point equation and look for an homogeneous solution. With $\Delta (\overrightarrow{r},\tau )={\Delta }_0$ the effective action is given by

$$S_e\left[{\Delta }_0\right]=\frac{\beta \mathcal{V}}{g}{|{\Delta }_0|}^2-\mathrm{Tr}[ln\left({\mathcal{G}}_0^{-1}\right)]+\beta \sum _{\mathbf{k}}{\xi }_{\mathbf{k}}$$

(53)

where $\mathcal{V}$ is the volume. In the momentum space and in the Matsubara frequencies the matrix ${\mathcal{G}}^{-1}$ defined as in Equation (49) is

$${\mathcal{G}}_0^{-1}(\mathbf{k},i{\omega }_n)=\left(\begin{array}{cc}i{\omega }_n-\frac{{\mathbf{k}}^2}{2m}+\mu & {\Delta }_0\\ {\Delta }_0^{*}& i{\omega }_n+\frac{{\mathbf{k}}^2}{2m}-\mu \end{array}\right)=\left(\begin{array}{cc}i{\omega }_n-{\xi }_{\mathbf{k}}& {\Delta }_0\\ {\Delta }_0^{*}& i{\omega }_n+{\xi }_{\mathbf{k}}\end{array}\right)$$

(54)

which has poles in real frequencies at $\pm E_{\mathbf{k}}$ where

$$E_{\mathbf{k}}=\sqrt{{\xi }_{\mathbf{k}}^2+{|{\Delta }_0|}^2}$$

(55)

The minimum of the spectrum is at ${ϵ}_{\mathbf{k}}=\mu$ with the minimum energy $|{\Delta }_0|$, characterized by the existence of a gap compared to the minimum of the spectrum of a free particle, which is zero. This implies that there is a minimum energetic cost for the creation of an elementary excitation. The physical interpretation is that the minimum energy requested for breaking a pair is $2|{\Delta }_0|$. The value of this minimum energy can be obtained imposing $\frac{\delta S_e}{\delta {\Delta }_0}=0$, getting

$$\frac{1}{g}=\frac{1}{\beta \mathcal{V}}\sum _{\mathbf{k},i{\omega }_n}\frac{1}{{\omega }_n^2+E_{\mathbf{k}}^2}$$

(56)

Performing the standard sum over the Matsubara frequencies and going to the continuum limit ${\displaystyle \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\to \frac{1}{{\left(2\pi \right)}^3}\int d^{3}k}$, one gets the so called gap equation

$$\frac{1}{g}=\frac{1}{{\left(2\pi \right)}^3}\int d^{3}k\frac{tanh(\beta E_{\mathbf{k}}/2)}{2E_{\mathbf{k}}}$$

(57)

The contact potential is a great simplification in the model, however it gives rise to ultraviolet divergences. Actually, Equation (57) diverges linearly with the ultraviolet cut-off.

BCS Superconductors

The Hamiltonian and the obtained gap equation are the same as those of the BCS theory for conventional superconductors. However in the BCS theory the ultraviolet behavior is regularized by the existence of a natural cut-off at the Debye frequency ${\omega }_D$ of the underlying lattice. Typically, in classical superconductors $\hslash {\omega }_D\ll {ϵ}_{\mathrm{F}}=\mu$ and a finite gap can exists only for particles near the Fermi energy. Exploiting the approximation of constant density of states (DOS), $\frac{1}{{\left(2\pi \right)}^3}\int d^{3}k\approx {\nu }_0\int d\xi$, where ${\nu }_0$ is the DOS at the Fermi energy, the gap equation becomes

$$\frac{1}{g}\approx {\nu }_0\underset{0}{\overset{\hslash {\omega }_D}{\int }}d\xi \frac{tanh\left(\frac{({\xi }^2+{|\Delta |}^2{)}^{1/2}}{2T_c}\right)}{({\xi }^2+{|\Delta |}^2{)}^{1/2}}$$

(58)

Solving the gap equation Equation (58) one gets, at $T=0$,

$$|{\Delta }_0|=2\hslash {\omega }_D{e}^{-\frac{1}{{\nu }_{0}g}}.$$

(59)

Supposing that at a critical temperature, $T=T_c$, the thermal fluctuations spoil the superconductivity, ${\Delta }_0\left(T_c\right)=0$, one gets

$$T_c=2\frac{e^{\gamma }}{\pi }\hslash {\omega }_D{e}^{-\frac{1}{{\nu }_{0}g}},$$

(60)

with $\gamma \approx 0.5772$, the Euler-Mascheroni number. Finally, for $T<T_c$ the gap is

$$|{\Delta }_0\left(T\right)|\propto \sqrt{T_c(T_c-T)}.$$

(61)

3.2. Number Equation

The second equation that we should considered describes the mean number of particles, usually called number equation. Calculations with the full effective action (52) are difficult, so one has to resort to an approximated scheme expanding the action close to the classical solution of the gap equation. Considering $\Delta \left(x\right)={\Delta }_0+\delta \Delta \left(x\right)$, the action can be expanded in the fluctuations of the auxiliary field at the desired order. The zeroth order is equivalent to the mean field approximation. We will then include the fluctuations at the Gaussian level so that

$$\begin{array}{c}{S}_e[\Delta ]=S_e\left[{\Delta }_0\right]+\underset{=0}{\underbrace{\int dx{\left.\frac{\delta S_e[\Delta ]}{\delta \Delta \left(x\right)}\right|}_{{\Delta }_0}}}\delta \Delta \left(x\right)+\hfill \\ \hfill +\frac{1}{2}\int d{x}_{1}d{x}_2{\left.\frac{{\delta }^2{S}_e[\Delta ]}{\delta \Delta \left(x_1\right)\delta \Delta \left(x_2\right)}\right|}_{{\Delta }_0}\delta \Delta \left(x_1\right)\delta \Delta \left(x_2\right)+o\left({\Delta }^3\right)\end{array}$$

(62)

3.2.1. Mean Field

At the mean field level the action is proportional to the grand canonical potential

$$\mathsf{\Omega }[{\Delta }_0,]=\beta S_e\left[{\Delta }_0\right]$$

(63)

from which is possible to derive the equation for the mean number of particles

$$N=-\frac{\partial \mathsf{\Omega }}{\partial \mu }=\sum _{\mathbf{k}}\left[1+\frac{1}{\beta }\sum _{i{\omega }_n}\mathrm{Tr}\left({\mathcal{G}}_0\frac{\partial {\mathcal{G}}_0^{-1}}{\partial \mu }\right)\right]=\sum _{\mathbf{k}}\left[1-\frac{1}{\beta }\sum _{i{\omega }_n}\left(\frac{2{\xi }_{\mathbf{k}}}{{\omega }_n^2+E_{\mathbf{k}}^2}\right)\right]$$

(64)

Summing over the Matsubara frequencies one gets

$$N=\sum _{\mathbf{k}}\left[1-\frac{{\xi }_{\mathbf{k}}}{E_{\mathbf{k}}}tanh\left(\frac{E_{\mathbf{k}}}{2T}\right)\right]$$

(65)

Let us consider the cases at $T=0$ and $T=T_c$.

• At $T=0$.

In order to study the behaviors of the gap ${\Delta }_0$ and of the chemical potential $\mu$ as functions of the coupling, at $T=0$, at the mean field level, the equations to solve are

$$\begin{array}{ccc}\hfill \frac{1}{g}& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{1}{2E_{\mathbf{k}}}\hfill \end{array}$$

(66)

$$\begin{array}{ccc}\hfill n& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left(1-\frac{{\xi }_{\mathbf{k}}}{E_{\mathbf{k}}}\right)\hfill \end{array}$$

(67)

where $n=N/\mathcal{V}$, the particle density.

• At $T=T_c$.

On the other hand if we want to calculate the critical temperature $T_c$, putting by definition ${\Delta }_0\left(T_c\right)=0$, at the mean field level the equation to solve in terms of $T_c$ are

$$\begin{array}{ccc}\hfill \frac{1}{g}& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}}{2T_c}\right)}{2{\xi }_{\mathbf{k}}}\hfill \end{array}$$

(68)

$$\begin{array}{ccc}\hfill n& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[1-tanh\left(\frac{{\xi }_{\mathbf{k}}}{2T_c}\right)\right].\hfill \end{array}$$

(69)

3.2.2. Inclusion of the Gaussian Fluctuations

Let us now include quantum fluctuations beyond the mean field approach at the Gaussian level. The fermionic propagator appearing in Equation (52), in momentum space, can be written as

$${\mathcal{G}}^{-1}=\left(\begin{array}{cc}i{\omega }_n-{\xi }_{\mathbf{k}}& \Delta \\ {\Delta }^{*}& i{\omega }_n+{\xi }_{\mathbf{k}}\end{array}\right)=\left(\begin{array}{cc}i{\omega }_n-{\xi }_{\mathbf{k}}& 0\\ 0& i{\omega }_n+{\xi }_{\mathbf{k}}\end{array}\right)+\left(\begin{array}{cc}0& \Delta \\ {\Delta }^{*}& 0\end{array}\right)\equiv {\widehat{G}}^{-1}+\widehat{\Delta }$$

(70)

Expanding the logarithm for small $\Delta$ one gets

$$ln\left({\mathcal{G}}^{-1}\right)=ln\left({\widehat{G}}^{-1}(\mathbb{I}+\widehat{G}\widehat{\Delta })\right)=ln\left({\widehat{G}}^{-1}\right)+\widehat{G}\widehat{\Delta }-\frac{1}{2}\widehat{G}\widehat{\Delta }\widehat{G}\widehat{\Delta }+O\left({\Delta }^3\right)$$

(71)

The linear term in $\widehat{\Delta }$ has a vanishing trace, therefore, up to second order the action is given by

$$S_e[\Delta ]=S_e^{\left(0\right)}+\int d^{3}r{\int }_0^{\beta }d\tau \frac{|\Delta (\overrightarrow{r},\tau ){|}^2}{g}+\frac{1}{2}\mathrm{Tr}\left(\widehat{G}\widehat{\Delta }\widehat{G}\widehat{\Delta }\right)+o\left({\Delta }^3\right)$$

(72)

Making explicit the quadratic term in $\Delta$ we have

$$\widehat{G}\widehat{\Delta }\widehat{G}\widehat{\Delta }=\left(\begin{array}{cc}0& G_{11}\Delta \\ G_{22}{\Delta }^{*}& 0\end{array}\right)\left(\begin{array}{cc}0& G_{11}\Delta \\ G_{22}{\Delta }^{*}& 0\end{array}\right)=\left(\begin{array}{cc}{G}_{11}\Delta G_{22}{\Delta }^{*}& 0\\ 0& G_{22}{\Delta }^{*}{G}_{11}\Delta \end{array}\right)$$

All the calculation can be done in momentum space

$$\frac{1}{2}\mathrm{Tr}\left(\widehat{G}\widehat{\Delta }\widehat{G}\widehat{\Delta }\right)=\mathrm{Tr}(G_{22}{\Delta }^{*}{G}_{11}\Delta )=\frac{1}{\beta \mathcal{V}}\sum _{\mathsf{q}}{\Delta }^{*}\left(\mathsf{q}\right)\left[\sum _{\mathsf{k}}{G}_{11}\left(\mathsf{k}\right)G_{22}(\mathsf{k}-\mathsf{q})\right]\Delta \left(\mathsf{q}\right)$$

(73)

where $\mathsf{k}=(\mathbf{k},i{\omega }_n)$ and $\mathsf{q}=(\mathbf{q},i{\nu }_m)$ are the four-momenta and $\widehat{G}\left(\mathsf{k}\right)$ is given by

$$\widehat{G}\left(\mathsf{k}\right)=\widehat{G}(\mathbf{k},i{\omega }_n)=\left(\begin{array}{cc}{G}_{11}& 0\\ 0& G_{22}\end{array}\right)=\left(\begin{array}{cc}\frac{1}{i{\omega }_n-{\xi }_{\mathbf{k}}}& 0\\ 0& \frac{1}{i{\omega }_n+{\xi }_{\mathbf{k}}}\end{array}\right)$$

(74)

whose components fulfill the property $G_{11}\left(\mathsf{k}\right)=-G_{22}(-\mathsf{k})$ that can be exploited to rewrite the quadratic term. Calling for simplicity $G_{11}\left(\mathsf{k}\right)=G\left(\mathsf{k}\right)$, we have

$$\frac{1}{2}\mathrm{Tr}\left(\widehat{G}\widehat{\Delta }\widehat{G}\widehat{\Delta }\right)=\sum _{\mathsf{q}}{\Delta }^{*}\left(\mathsf{q}\right)\left[-\frac{1}{\beta \mathcal{V}}\sum _{k}G\left(\mathsf{k}\right)G(\mathsf{q}-\mathsf{k})\right]\Delta \left(\mathsf{q}\right)\equiv \sum _{\mathsf{q}}{\Delta }^{*}\left(\mathsf{q}\right)\chi \left(\mathsf{q}\right)\Delta \left(\mathsf{q}\right)$$

(75)

where we defined the function $\chi \left(\mathsf{q}\right)$ which, after performing the summation over the Matsubara frequencies $i{\omega }_n$, reads

$$\chi \left(\mathsf{q}\right)=\frac{1}{\beta \mathcal{V}}\sum _{\mathbf{k},i{\omega }_n}\frac{1}{i{\omega }_n-{\xi }_{\mathbf{k}}}·\frac{1}{i{\omega }_n-i{\nu }_m+{\xi }_{\mathbf{q}-\mathbf{k}}}=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{q}-\mathbf{k}}\right)}{i{\nu }_m-{\xi }_{\mathbf{k}}-{\xi }_{\mathbf{q}-\mathbf{k}}}$$

(76)

where $n_{\mathrm{F}}$ is the Fermi distribution. The action at second order in $\Delta$, corresponding to the inclusion of the Gaussian fluctuation around mean field solution, is

$$S_e[\Delta ]=S_e^{\left(0\right)}+\sum _{\mathsf{q}}{\Delta }^{*}\left(\mathsf{q}\right){\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)\Delta \left(\mathsf{q}\right)+o\left({\Delta }^3\right)$$

(77)

where we introduced the function

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)=\frac{1}{g}+\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}}\right)}{i{\nu }_m-{\xi }_{\mathbf{k}}-{\xi }_{\mathbf{k}-\mathbf{q}}}$$

(78)

The partition function at the Gaussian level is, therefore,

$${\mathcal{Z}}_G={\mathcal{Z}}_0\int D{\Delta }^{*}D\Delta e^{-{\sum }_{\mathsf{q}}{\Delta }^{*}\left(\mathsf{q}\right){\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)\Delta \left(\mathsf{q}\right)}={\mathcal{Z}}_{0}det\left(\mathsf{\Gamma }\right)$$

(79)

with the definition ${\mathcal{Z}}_0=e^{-S_e^{\left(0\right)}}$, which can be written as

$${\mathcal{Z}}_{0}det\left(\mathsf{\Gamma }\right)={\mathcal{Z}}_0{e}^{-ln[det\left({\mathsf{\Gamma }}^{-1}\right)]}={\mathcal{Z}}_0{e}^{-\mathrm{Tr}[ln\left({\mathsf{\Gamma }}^{-1}\right)]}=e^{-S_e^{\left(0\right)}-\mathrm{Tr}[ln\left({\mathsf{\Gamma }}^{-1}\right)]}$$

(80)

the grand canonical potential is shown to be

$${\mathsf{\Omega }}_G=-\frac{1}{\beta }ln\left({\mathcal{Z}}_G\right)=\frac{1}{\beta }{S}_e^{\left(0\right)}+\frac{1}{\beta }\mathrm{Tr}[ln\left({\mathsf{\Gamma }}^{-1}\right)]$$

(81)

The number equation, at the Gaussian level, is given by

$$n=-\frac{1}{\mathcal{V}}\frac{\partial {\mathsf{\Omega }}_G}{\partial \mu }=n^{\left(0\right)}+n^{\left(2\right)}$$

(82)

where

$$n^{\left(2\right)}=-\frac{1}{\beta \mathcal{V}}\frac{\partial }{\partial \mu }\mathrm{Tr}[ln\left({\mathsf{\Gamma }}^{-1}\right)]=-\frac{1}{\beta \mathcal{V}}\frac{\partial }{\partial \mu }\sum _{\mathbf{q},i{\nu }_m}ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},i{\nu }_m)\right)$$

(83)

The first term in Equation (82) is simply the density of particles at the mean field level while the second one is due to the quantum corrections. In general if $ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},z)\right)$, has poles in the complex plane at $z=z_j$ ($j=0,1,2\cdots$) and a branch cut on the real axis, the sum over Matsubara frequencies can be rewritten as

$$\begin{array}{c}\frac{1}{\beta }\sum _{i{\nu }_n}ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},i{\nu }_n)\right)=\sum _j{n}_{\mathrm{B}}\left(z_j\right)\mathrm{Res}\left[ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},z_j)\right)\right]-\hfill \\ \hfill -\frac{1}{2\pi i}{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\left[ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega +iϵ)\right)-ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega -iϵ)\right)\right]\end{array}$$

(84)

In our case the integrand has poles only on the cut, so that this formula reduces to

$$\frac{1}{\beta }\sum _{i{\nu }_n}ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},i{\nu }_n)\right)=\frac{1}{2\pi i}{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\left[ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega -iϵ)\right)-ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega +iϵ)\right)\right]$$

(85)

where $n_{\mathrm{B}}$ is the Bose distribution. Since $\mathsf{\Gamma }(\mathbf{q},\omega )$ has to be real on the real axis, it can be written as a modulus times a phase. Defining a phase shift $\delta \left(\mathsf{q}\right)$ as in Ref. [4], we can write

$$\mathsf{\Gamma }{(\mathbf{q},\omega \pm iϵ)}^{-1}=\left|\mathsf{\Gamma }{(\mathbf{q},\omega \pm iϵ)}^{-1}\right|e^{\pm i\delta (\mathbf{q},\omega +iϵ)}$$

(86)

therefore

$$\frac{1}{\beta }\sum _{i{\nu }_n}ln\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},i{\nu }_n)\right)==-\frac{1}{\pi }{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\delta (\mathbf{q},\omega )$$

(87)

At the Gaussian level, therefore, the density of particles can written in the following form

$$n=n^{\left(0\right)}+\frac{1}{\pi \mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }.$$

(88)

The number and the gap equations allow us to find the critical temperature and the chemical potential as functions of the interaction strength. In the following section we will solve those equations, first at the mean field level and then with Gaussian fluctuations. However, before we proceed with the calculation, due to the presence of a contact potential, we notice that the right-hand-side of Equation (57) has an ultraviolet divergence that has to be cured.

4. BCS-BEC Crossover

In this section we will derive the critical temperature $T_c$ and the chemical potential at $T_c$ as functions of the interaction coupling, discussing the limits of strong and weak attraction, namely the BEC and the BCS limits, respectively, solving the gap equation and the number equation introduced previously. In Section 2.2 we showed that an effective contact potential, $V\left(\mathbf{r}\right)=-g\delta \left(\mathbf{r}\right)$, can reproduce the same scattering properties of a realistic potential, in the low energy limit, if g and the scattering length $a_s$ satisfy the relation in Equation (40), reported here for convenience

$$\frac{m}{4\pi a_s}=-\frac{1}{g}+\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{1}{2{ϵ}_{\mathbf{k}}}$$

(89)

Let us substitute the coupling g from Equation (89) in the gap equation Equation (57), linking in this way the toy model discussed until now with a realistic ultracold gas, through the experimentally accessible scattering length $a_s$. In this way the gap equation reads

$$-\frac{m}{4\pi a_s}=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left(\frac{tanh({\xi }_{\mathbf{k}}/2T_c)}{2{\xi }_{\mathbf{k}}}-\frac{1}{2{ϵ}_{\mathbf{k}}}\right)$$

(90)

The divergence in the gap equation is exactly canceled out by the divergence in Equation (89). The limits of strong ($g\to \infty$) and weak coupling ($g\to 0$) have to be reviewed in terms of $a_s$ which is now the tuning parameter along the crossover. For this purpose it is useful to introduce an ultraviolet cut-off $\Lambda$ so that Equation (89) becomes

$$\frac{m}{4\pi a_s}=-\frac{1}{g}+\frac{1}{\mathcal{V}}\sum _{|\mathbf{k}|<\Lambda }\frac{1}{2{ϵ}_{\mathbf{k}}}$$

(91)

In the limit $g\to 0$, namely in the BCS limit, the cut-off $\Lambda$ does not play any role and eventually can be sent to infinity, therefore the scattering length is negative and $\frac{1}{a_s}\to -\infty$.

In the limit $g\to \infty$, namely in the BEC limit, we get, instead,

$$\frac{m}{4\pi a_s}=\frac{1}{\mathcal{V}}\sum _{|k|<\Lambda }\frac{1}{2{ϵ}_{\mathbf{k}}}=\frac{1}{2{\pi }^2}{\int }_0^{\Lambda }dk\frac{k^2}{2{ϵ}_{\mathbf{k}}}=\frac{m}{2{\pi }^2}{\int }_0^{\Lambda }dk=\frac{m}{2{\pi }^2}\Lambda$$

(92)

Sending the cut-off to infinity the limit of strong interaction can be captured for $\frac{1}{a_s}\to +\infty$. Since the scattering length comes from the low energy limit of the interaction, it is an intrinsic parameter and not an adjustable property of the potential. However, as already mentioned, by the mechanism of the Feshbach resonance, the scattering length can be modified on a wide range of values, therefore it is the tunable parameter along the crossover.

4.1. Mean Field Theory

At the mean field level one has to simultaneously solve the gap and the number equations found previously,

$$\begin{array}{ccc}\hfill -\frac{m}{4\pi a_s}& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left(\frac{tanh({\xi }_{\mathbf{k}}/2T_c)}{2{\xi }_{\mathbf{k}}}-\frac{1}{2{ϵ}_{\mathbf{k}}}\right)\hfill \end{array}$$

(93)

$$\begin{array}{ccc}\hfill n& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[1-tanh\left(\frac{{\xi }_{\mathbf{k}}}{2T_c}\right)\right]\hfill \end{array}$$

(94)

We will discuss in what follows the weak (BCS) and the strong (BEC) interacting limits.

4.1.1. Weak Coupling

Let us consider the number equation Equation (94), which, in the continuum limit, reminding that ${ϵ}_{\mathbf{k}}=\frac{k^2}{2m}$ and introducing the density of state, becomes

$$n=\frac{{\left(2m\right)}^{\frac{3}{2}}}{2{\pi }^2}{\int }_0^{\infty }\frac{\sqrt{ϵ}}{e^{\frac{ϵ-\mu }{T_c}}+1}dϵ$$

(95)

For $T_c$ very small one can resort to the well-known Sommerfield expansion, getting for the density of particle at low temperature

$$n={\left(2m\right)}^{\frac{3}{2}}\left[\frac{1}{3{\pi }^2}{\mu }^{\frac{3}{2}}+\frac{1}{24}\frac{T_c^2}{{\mu }^{1/2}}+o\left(T_c^4\right)\right]$$

(96)

At fixed density $n=\frac{{\left(2m{ϵ}_{\mathrm{F}}\right)}^{3/2}}{3{\pi }^2}$, which is the typical experimental situation, this equation of state implies ${ϵ}_{\mathrm{F}}^{3/2}={\mu }^{\frac{3}{2}}+\frac{{\pi }^2}{8}\frac{T_c^2}{{\mu }^{1/2}}+o\left(T_c^4\right)$, and solving for $\mu$ we find an explicit relationship between the chemical potential and the temperature

$$\mu ={ϵ}_{\mathrm{F}}-\frac{{\pi }^2}{12}\frac{T_c^2}{{ϵ}_{\mathrm{F}}}+o\left(T_c^4\right)$$

(97)

The chemical potential is, therefore, almost constant, $\mu \simeq {ϵ}_{\mathrm{F}}$. Thermal corrections occur only at second order. Exploiting this result, we can solve the gap equation Equation (93), as done in Ref. [77]. With constant chemical potential, Equation (93) fixes the behavior of the critical temperature $T_c$ as a function of the scattering length. In the continuum limit, after introducing the density of state and after a change of variables, $x=ϵ/\mu$, Equation (93) becomes

$$-\frac{\pi }{k_{\mathrm{F}}{a}_s}={\int }_0^{\infty }dx\sqrt{x}\left(\frac{tanh(\mu (x-1)/2T_c)}{x-1}-\frac{1}{x}\right)$$

(98)

At $\mu ={ϵ}_{\mathrm{F}}$ and supposing $T_c\ll {ϵ}_{\mathrm{F}}$ one gets $\frac{\pi }{k_{\mathrm{F}}{a}_s}=-2ln\left(\frac{8e^{\gamma }{ϵ}_{\mathrm{F}}}{e^2\pi T_c}\right)$, with $\gamma$ being the Eulero–Mascheroni number. We remember that $a_s$ has a small negative value in the weak coupling limit. Solving in terms of $T_c$ we obtain

$$T_c=\frac{8e^{\gamma }}{\pi e^2}{ϵ}_{\mathrm{F}}{e}^{-\frac{\pi }{2k_{\mathrm{F}}|a_s|}}$$

(99)

a critical temperature exponentially small in the weak coupling limit, in agreement with the BCS theory and consistently with the hypothesis $T_c\ll {ϵ}_{\mathrm{F}}$.

4.1.2. Strong Coupling

In the strong coupling limit, we will solve first the gap equation Equation (93) to obtain the chemical potential $\mu$ which is expected to be very large and negative, as in the case of free fermions at high temperatures, and that $|\mu \left(T_c\right)|\gg T_c$. Taking advantage of this ansatz, whose validity will be proved a posteriori, we have

$$-\frac{m}{4\pi a_s}=\frac{m^{3/2}{T}_c^{1/2}}{2{\pi }^2}{\int }_0^{\infty }dx\sqrt{x}\left(\frac{tanh(x-\tilde{\mu })}{x-\tilde{\mu }}-\frac{1}{x}\right)$$

(100)

where we defined $\tilde{\mu }=\frac{\mu \left(T_c\right)}{2T_c}$. For $|\tilde{\mu }|\gg 1$, we can approximate $\frac{tanh(x-\tilde{\mu })}{x-\tilde{\mu }}\approx \frac{1}{x-\tilde{\mu }}$, so that, rescaling $x\to -\tilde{\mu }x$, we get

$$-\frac{m}{4\pi a_s}\approx \frac{m^{3/2}{T}_c^{1/2}}{2{\pi }^2}\sqrt{-\tilde{\mu }}{\int }_0^{\infty }dx\frac{\sqrt{x}}{x(x+1)}=\frac{m^{3/2}}{2\pi }\sqrt{-\mu }$$

(101)

The relation between the chemical potential and the scattering length, in the strong coupling limit, is, therefore, given by

$$\mu \approx -\frac{1}{2m{a}_s^2}=\frac{E_b}{2}$$

(102)

where we used Equation (35) (here $\hslash =1$), reminding, from Section 2.1, that for positive scattering length the two-body problem supports the existence of a bound state with energy $E_b$. We will use Equation (102) in the number equation Equation (94) to derive the critical temperature. After introducing the density of state and after a rescaling, in the limit $\frac{1}{a_s}\to +\infty$, we have

$$\begin{array}{c}\hfill n=\frac{{\left(2m\right)}^{\frac{3}{2}}}{2{\pi }^2}{\int }_0^{\infty }\frac{\sqrt{ϵ}}{e^{\frac{ϵ-\mu }{T_c}}+1}dϵ\approx \frac{{\left(2m{T}_c\right)}^{\frac{3}{2}}}{2{\pi }^2}{e}^{\frac{\mu }{T_c}}{\int }_0^{\infty }\sqrt{x}{e}^{-x}dx=\frac{{\left(m{T}_c\right)}^{\frac{3}{2}}}{\sqrt{2}{\pi }^{\frac{3}{2}}}{e}^{-\frac{|E_b|}{2T_c}}\end{array}$$

At fixed density of particles along the crossover, $n=\frac{{\left(2m{ϵ}_{\mathrm{F}}\right)}^{3/2}}{3{\pi }^2}$, we have ${\left(\frac{{ϵ}_{\mathrm{F}}}{|E_b|}\right)}^{\frac{3}{2}}\propto {\left(\frac{T_c}{|E_b|}\right)}^{\frac{3}{2}}{e}^{-\frac{|E_b|}{2T_c}}$, therefore in the strong coupling limit the leading term for the critical temperature at the mean field level is given by

$$T_c\approx \frac{|E_b|}{3ln\left(|E_b|/{ϵ}_{\mathrm{F}}\right)},$$

(103)

namely, $T_c$ grows to infinity in the deep strong coupling regime. Actually this is the temperature at which Cooper pairs break down due to thermal fluctuations, since it was derived imposing $\Delta \left(T_c\right)=0$, namely in a system of free fermions. It is then naturally expected that when the interaction becomes strong, this approximation fails. As we have seen, at strong coupling, $T_c/{ϵ}_{\mathrm{F}}\gg 1$, the gas is out of the degenerate regime. Actually, at high temperatures the system is described by a mixture of free particles and pairs, both following almost a Maxwell–Boltzmann distribution therefore the equilibrium can be imposed by he condition ${\mu }_p=2{\mu }_f$, namely when the chemical potential of the pairs matches twice the chemical potential of the unpaired fermions. In this picture $T_c$ can be seen as the dissociation temperature for the pairs.

4.1.3. Along the Crossover

Let us consider Equations (93) and (94) in the continuum limit. Defining $y=\frac{1}{k_{\mathrm{F}}{a}_s}$ and $\tilde{\mu }=\frac{\mu }{2T_c}$, working at fixed density, we can write

$$\begin{array}{c}\hfill y{\left(\frac{{ϵ}_{\mathrm{F}}}{T_c}\right)}^{\frac{1}{2}}=-\frac{\sqrt{2}}{\pi }{\int }_0^{\infty }\sqrt{x}\left(\frac{tanh(x-\tilde{\mu })}{x-\tilde{\mu }}-\frac{1}{x}\right)dx\equiv I_g\left(\tilde{\mu }\right)\end{array}$$

(104)

$$\begin{array}{c}\hfill {\left(\frac{{ϵ}_{\mathrm{F}}}{T_c}\right)}^{\frac{3}{2}}=\frac{3\sqrt{2}}{2}{\int }_0^{\infty }\left(1-tanh(x-\tilde{\mu })\right)dx\equiv I_0\left(\tilde{\mu }\right)\end{array}$$

(105)

These integrals can be calculated numerically for different values of $\tilde{\mu }$, building the lists of values for $\left(\frac{\mu }{{ϵ}_{\mathrm{F}}},y\right)$ and $\left(\frac{T_c}{{ϵ}_{\mathrm{F}}},y\right)$. The graphs reported in Figure 1 are obtained in this way, confirming the analytical results at mean field level.

4.2. Beyond Mean Field: Gaussian Fluctuations

As seen previously, the second order expansion of the action is given by Equation (77), which is

$$S_e[\Delta ]=S_e^{\left(0\right)}+\sum _{\mathsf{q}}{\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right){|\Delta \left(\mathsf{q}\right)|}^2+o\left({\Delta }^3\right)$$

(106)

with ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$ given by Equation (78). After regularizing the contact interaction through Equation (89), and after a shift of the momenta $\mathbf{k}\to \mathbf{k}+\mathbf{q}/2$ in order, for simplicity, to get rid of the angular dependence in the denominator on the angle defined by $\mathbf{k}·\mathbf{q}=kqcos\theta$, we can write

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)=-\left\{\frac{m}{4\pi a_s}+\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}/2}\right)}{{\xi }_{\mathbf{k}+\mathbf{q}/2}+{\xi }_{\mathbf{k}-\mathbf{q}/2}-i{\nu }_n}-\frac{1}{2{ϵ}_{\mathbf{k}}}\right]\right\}$$

(107)

The sum is invariant under $\mathbf{k}\to -\mathbf{k}$, therefore the numerator $[1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}/2}\right)]$ can be written as $[1-2n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}\right)]=tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}}{2T_c}\right)$, getting simply

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)=-\left\{\frac{m}{4\pi a_s}+\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}}{2T_c}\right)}{{\xi }_{\mathbf{k}+\mathbf{q}/2}+{\xi }_{\mathbf{k}-\mathbf{q}/2}-i{\nu }_n}-\frac{1}{2{ϵ}_{\mathbf{k}}}\right]\right\}$$

(108)

In the continuum limit and writing explicitly the spectra ${\xi }_{\mathbf{k}\pm \mathbf{q}/2}=\frac{{|\mathbf{k}\pm \mathbf{q}/2|}^2}{2m}$ we have

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)=-\frac{m}{4\pi a_s}-\frac{m}{{\left(2\pi \right)}^2}{\int }_{-1}^{1}dcos\theta \int dk{k}^2\left[\frac{tanh\left(\frac{k^2+q^2/4-2m\mu +kqcos\theta }{4m{T}_c}\right)}{k^2+q^2/4-2m\mu -im{\nu }_n}-\frac{1}{k^2}\right]$$

(109)

The dependence on $cos\theta$ is only in the hyperbolic tangent and we can integrate over it, using ${\int }_{-1}^{1}dxtanh\left(a+bx\right)=\frac{1}{b}[ln(cosh(a+b))-ln(cosh(a-b))]$,

$${\int }_{-1}^{1}dxtanh\left(\frac{k^2+\frac{q^2}{4}+kqx-2m\mu }{4m{T}_c}\right)=\frac{4m{T}_c}{kq}\left[A(k,q)-A(k,-q)\right]$$

(110)

where

$$A(k,q)=ln\left[cosh\left(\frac{{(k+q/2)}^2}{4m{T}_c}-\frac{\mu }{2T_c}\right)\right]$$

(111)

We have, at the end, the following expression

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)=-\left\{\frac{m}{4\pi a_s}+\frac{m}{{\left(2\pi \right)}^2}{\int }_0^{+\infty }dk{k}^2\left[\frac{4m{T}_c}{kq}\frac{A(k,q)-A(k,-q)}{k^2+q^2/4-2m\mu -im{\nu }_n}-\frac{2}{k^2}\right]\right\}$$

(112)

Let us consider the weak and the strong coupling limits.

4.2.1. Weak Coupling

In the weak coupling limit the chemical potential is finite and the critical temperature and the scattering length go to zero, $a_s\to 0^{-}$. Without entering into the details we have

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)=-\frac{m}{4\pi }\left(\frac{1}{a_s}+f(\mathsf{q},\mu )\right)$$

(113)

where $f(\mathsf{q},\mu )$ some complex function of the four-momentum $\mathsf{q}$ (momentum and frequency) and chemical potential $\mu$. In the weak coupling limit, approximately ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$ does not depend on critical temperature since, for small $T_c$ we have $A(k,q)\sim \frac{1}{2T_c}\left|\frac{{(k+q/2)}^2}{2m}-\mu \right|$. Calculating explicitly the second order correction to the density of particles, we have

$$\begin{array}{ccc}\hfill n^{\left(2\right)}& =& \frac{1}{\beta \mathcal{V}}\frac{\partial }{\partial \mu }\sum _{\mathsf{q}}ln\left({\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)\right)\hfill \\ & =& \frac{1}{\beta \mathcal{V}}\sum _{\mathsf{q}}\frac{1}{\frac{1}{a_s}+f(\mathsf{q},\mu )}\frac{\partial f(\mathsf{q},\mu )}{\partial \mu }\underset{a_s\to 0^{-}}{⟶}\frac{a_s}{\beta \mathcal{V}}\sum _{\mathsf{q}}\frac{\partial f(\mathsf{q},\mu )}{\partial \mu }\propto a_s\hfill \end{array}$$

(114)

In this limit the correction to the density of particles is proportional to the scattering length and therefore is negligible. As a result the chemical potential will be almost constant, $\mu \approx {ϵ}_{\mathrm{F}}$, and since also the gap equation is unchanged, in the limit $\frac{1}{k_{\mathrm{F}}{a}_s}\to -\infty$, the critical temperature is expected to be almost the same as that obtained at the mean field level. In other words, the Gaussian fluctuations are not so effective at weak coupling. This expectation will be confirmed by numerical calculations.

4.2.2. Strong Coupling

As we will see, in the strong coupling limit the effective partition function is equivalent to that of a system of free bosons whose critical temperature is that of the BEC. In order to prove this result we first need to manipulate ${\mathsf{\Gamma }}^{-1}$. As we have seen, in the BEC limit, $\beta \mu \to -\infty$ and $1/k_{\mathrm{F}}{a}_s\to \infty$. For a large and negative chemical potential

$$A(k,q)=ln\left[cosh\left(\frac{\frac{{(k+q/2)}^2}{2m}-\mu }{2T_c}\right)\right]\underset{\mu \to -\infty }{⟶}\frac{\frac{{(k+q/2)}^2}{2m}-\mu }{2T_c}$$

(115)

therefore

$$A(k,q)-A(k,-q)\underset{\mu \to -\infty }{⟶}\frac{\frac{{(k+q/2)}^2}{2m}-\mu }{2T_c}-\frac{\frac{{(k-q/2)}^2}{2m}-\mu }{2T_c}=\frac{kq}{2m{T}_c}$$

(116)

Inserting this result in Equation (112), after integration, we get

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)\approx -\frac{m}{4\pi }\left(\frac{1}{a_s}-\sqrt{\frac{q^2}{4}-2m\mu -im{\nu }_n}\right)$$

(117)

Taking advantage of the definition of the bound state $E_b=-\frac{1}{m{a}_s^2}$, and defining ${\mu }_{\mathrm{B}}=2\mu -E_b$, we can write

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)\approx \frac{1}{R\left(\mathsf{q}\right)}\left[i{\nu }_n-\left(\frac{q^2}{4m}-{\mu }_{\mathrm{B}}\right)\right]$$

(118)

where $R\left(\mathsf{q}\right)=\frac{4\pi }{m^{3/2}}\left(\sqrt{-E_b}+\sqrt{\frac{q^2}{4m}-2\mu -i{\nu }_n}\right)$, which is simply $R\left(\mathsf{q}\right)\approx \frac{8\pi }{m^2{a}_s}$ for small fluctuations in the deep BEC limit. The quantity ${\mu }_{\mathrm{B}}$ can be seen as the chemical potential for free bosons, which makes sense only if it is negative. The partition function takes the form

$${\mathcal{Z}}_G={\mathcal{Z}}_0\int D{\Delta }^{*}D\Delta e^{-{\sum }_{\mathsf{q}}{\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right){|\Delta \left(\mathsf{q}\right)|}^2}\approx {\mathcal{Z}}_0\int D{\Delta }^{*}D\Delta e^{-{\sum }_{\mathsf{q}}\left[i{\nu }_n-\left(\frac{q^2}{4m}-{\mu }_{\mathrm{B}}\right)\right]\frac{{|\Delta \left(\mathsf{q}\right)|}^2}{R\left(\mathsf{q}\right)}}$$

(119)

Rescaling the fields $\varphi \left(\mathsf{q}\right)=\frac{\Delta \left(\mathsf{q}\right)}{\sqrt{R\left(\mathsf{q}\right)}}$ and defining $m_{\mathrm{B}}=2m$ the partition function at the Gaussian level is equivalent to that for free bosons with mass $m_{\mathrm{B}}$ (equal to the mass of two fermions)

$${\mathcal{Z}}_G\propto \int D{\varphi }^{*}D\varphi e^{-{\sum }_{\mathsf{q}}{\varphi }^{*}\left(\mathsf{q}\right)\left[i{\nu }_n-\left(\frac{q^2}{2m_{\mathrm{B}}}-{\mu }_{\mathrm{B}}\right)\right]\varphi \left(\mathsf{q}\right)}$$

(120)

Now we can take advantage of the well-known expression for the critical temperature for the Bose–Einstein condensation which is given by

$$T_c=\frac{2\pi }{m_{\mathrm{B}}}{\left(\frac{n_{\mathrm{B}}}{2\zeta (3/2)}\right)}^{\frac{2}{3}}$$

(121)

If in the deep BEC regime all fermions form pairs, therefore $n_{\mathrm{B}}=n/2$. At fixed fermionic density $n=\frac{{\left(2m{ϵ}_{\mathrm{F}}\right)}^{2/3}}{3{\pi }^2}$, we have

$$T_c=\frac{2\pi }{m_{\mathrm{B}}}{\left(\frac{n_{\mathrm{B}}}{2\zeta (3/2)}\right)}^{\frac{2}{3}}=\frac{\pi }{m}{\left(\frac{{\left(2m{ϵ}_{\mathrm{F}}\right)}^{3/2}}{6{\pi }^2\zeta (3/2)}\right)}^{\frac{2}{3}}={\left(\frac{\sqrt{2}}{3\sqrt{\pi }\zeta (3/2)}\right)}^{\frac{2}{3}}{ϵ}_{\mathrm{F}}\approx 0.218{ϵ}_{\mathrm{F}}$$

(122)

According to the BEC theory, at the critical point, the effective chemical potential approaches zero from negative values, therefore

$${\mu }_{\mathrm{B}}\left(T_c\right)=2\mu \left(T_c\right)-E_b\underset{\frac{1}{k_{\mathrm{F}}{a}_s}\to \infty }{⟶}0$$

(123)

which is exactly the solution of the gap equation in the strong coupling limit, Equation (102).

4.2.3. Full Crossover

As shown in Section 3.2.2, the set of equations to solve at second order is

$$\begin{array}{ccc}\hfill -\frac{m}{4\pi a_s}& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left(\frac{tanh({\xi }_{\mathbf{k}}/2T_c)}{2{\xi }_{\mathbf{k}}}-\frac{1}{2{ϵ}_{\mathbf{k}}}\right)\hfill \end{array}$$

(124)

$$\begin{array}{ccc}\hfill n=n^{\left(0\right)}+n^{\left(2\right)}& =& \frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[1-tanh\left(\frac{{\xi }_{\mathbf{k}}}{2T_c}\right)\right]+\frac{1}{\pi \mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }\hfill \end{array}$$

(125)

with the phase shift $\delta (\mathbf{q},\omega )=-\mathrm{Arg}\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega +iϵ)\right)$. We need, therefore, to find the second term of the number equation, $n^{\left(2\right)}$, possibly resorting also to some approximations.

• Bosonic Approximation

As we have seen before, the corrections due to quantum fluctuations are more important in the BEC regime while they are negligible in the BCS weak coupling regime. The inverse of the vertex function in the strong coupling limit is

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)\approx i{\nu }_n-\left(\frac{q^2}{2m_{\mathrm{B}}}-{\mu }_{\mathrm{B}}\right)$$

(126)

which describes a system of free bosons, with ${\mu }_{\mathrm{B}}=2\mu -E_b$ and $m_{\mathrm{B}}=2m$. We can use this simple result as a first approximation to include Gaussian fluctuations in order to calculate the critical temperature. This approximation consists of simply substituting the second order correction $n^{\left(2\right)}$ in the number equation Equation (125) with the density of free bosons whose inverse propagator is given by Equation (126). It is well known that the critical temperature $T_c$ of free bosons is given by Equation (121). Let us suppose that the number of bosonic excitations composed by pairs of fermions is $n_{\mathrm{B}}=2n^{\left(2\right)}$, while the rest of fermions remains unpaired. In this approximation, using Equation (121), we have, therefore,

$$n^{\left(2\right)}\approx 2\zeta (3/2){\left(\frac{m{T}_c}{\pi }\right)}^{3/2}$$

(127)

As a result, we have to solve Equation (124) together with the number equation which, in the bosonic approximation, reads

$$n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[1-tanh\left(\frac{{\xi }_{\mathbf{k}}}{2T_c}\right)\right]+2\zeta (3/2){\left(\frac{m{T}_c}{\pi }\right)}^{3/2}.$$

(128)

After defining the following quantities

$$y=\frac{1}{k_{\mathrm{F}}{a}_s},\overline{\mu }=\frac{\mu }{{ϵ}_{\mathrm{F}}},{\overline{T}}_c=\frac{T_c}{{ϵ}_{\mathrm{F}}}$$

(129)

in the continuum limit, Equations (124) and (128), at fixed density $n=\frac{k_{\mathrm{F}}^3}{3{\pi }^2}$ can be written as

$$\begin{array}{ccc}& & y=-\frac{2}{\pi }{\int }_0^{\infty }dx{x}^2\left(\frac{tanh\left(\frac{x^2-\overline{\mu }}{2{\overline{T}}_c}\right)}{x^2-\overline{\mu }}-\frac{1}{x^2}\right)\hfill \end{array}$$

(130)

$$\begin{array}{ccc}& & \frac{2}{3}={\int }_0^{\infty }dx{x}^2\left[1-tanh\left(\frac{x^2-\overline{\mu }}{2{\overline{T}}_c}\right)\right]+2\zeta (3/2)\sqrt{\frac{\pi }{2}}{\overline{T}}_c^{3/2}\hfill \end{array}$$

(131)

This set of equations can be easily solved numerically, finding the lists of values $({\overline{T}}_c,\overline{\mu },y)$ which fulfill the two equations. Some results are reported in Figure 2.

• Exact Gaussian Fluctuations

Let us work further on the function ${\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega +iϵ)$, as obtained in Equation (112). With the rescalings $\tilde{\mu }=\frac{\mu }{2T_c}$, $p=\frac{k}{\sqrt{4m{T}_c}}$ and $Q=\frac{q}{\sqrt{4m{T}_c}}$, the function $A(k,q)$ becomes

$$A(k,q)=ln\left[cosh\left(\frac{\frac{{(k+q/2)}^2}{2m}-\mu }{2T_c}\right)\right]=ln\left[cosh\left({(p+Q/2)}^2-\tilde{\mu }\right)\right]=\tilde{A}(p,Q)$$

(132)

Making the same substitutions in the vertex function, after a Wick rotation in the frequencies (${\nu }_n\to \omega +iϵ$, with $ϵ>0$ an infinitesimal quantity), we can define ${\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega +iϵ)={\tilde{\mathsf{\Gamma }}}^{-1}(\mathbf{Q},\nu +iϵ)$, with the rescaled real frequency $\nu =\frac{\omega }{4T_c}$,

$${\tilde{\mathsf{\Gamma }}}^{-1}(\mathbf{Q},\nu +iϵ)=-\left(\frac{m}{4\pi }\right)\left\{\frac{1}{a_s}+\frac{\sqrt{4m{T}_c}}{\pi Q}{\int }_0^{\infty }dp\left[\frac{p\left[\tilde{A}(p,Q)-\tilde{A}(p,-Q)\right]}{p^2+Q^2/4-\tilde{\mu }-\nu -iϵ}-2Q\right]\right\}$$

(133)

We remind that, looking at Equation (88), the interesting part of the function ${\tilde{\mathsf{\Gamma }}}^{-1}(\mathbf{Q},\nu +iϵ)$ is actually its argument, therefore, its prefactor is irrelevant. Within the range of integration the integrand has a pole in $p_0=\sqrt{\tilde{\mu }+\nu -Q^2/4}$, if $p_0$ real and positive. Defining the function

$$t(p,Q)=\frac{p\left[\tilde{A}(p,Q)-\tilde{A}(p,-Q)\right]-2Q(p^2-p_0^2)}{p+p_0}$$

(134)

we can write

$${\tilde{\mathsf{\Gamma }}}^{-1}(\mathbf{Q},\nu +iϵ)=-C\left[Qy{\left(\frac{{ϵ}_{\mathrm{F}}}{T_c}\right)}^{1/2}+\frac{\sqrt{2}}{\pi }{\int }_0^{+\infty }\frac{t(p,Q)}{p-p_0-iϵ}dp\right]$$

(135)

with the prefactor $C=-{\left(\frac{m}{2}\right)}^{3/2}\frac{\sqrt{T_c}}{\pi Q}$. The integral in Equation (135) can be evaluated knowing its Cauchy principal value

$${\int }_0^{+\infty }\frac{t(p,Q)}{p-p_0\pm iϵ}dp=\mathcal{P}{\int }_0^{+\infty }\frac{t(p,Q)}{p-p_0}dp\mp i\pi {\int }_0^{+\infty }\delta (p-p_0)t(p,Q)$$

(136)

For $(\tilde{\mu }+\nu -Q^2/4)>0$, namely for $p_0$ real and positive, ${\mathsf{\Gamma }}^{-1}$ has an imaginary part. We can write the vertex function ${\tilde{\mathsf{\Gamma }}}^{-1}$ and its argument $\tilde{\delta }$ as it follows

$${\tilde{\mathsf{\Gamma }}}^{-1}(\mathbf{Q},\nu )=-C\left(X+iY\right),\tilde{\delta }(\mathbf{Q},\nu )=arctan\left(\frac{Y}{X}\right)$$

(137)

where the real and imaginary parts of ${\tilde{\mathsf{\Gamma }}}^{-1}$ are given by

$$\begin{array}{ccc}& & X=Qy{\left(\frac{{ϵ}_{\mathrm{F}}}{T_c}\right)}^{1/2}+\frac{\sqrt{2}}{\pi }\mathcal{P}{\int }_0^{\infty }\frac{t(p,Q)}{p-p_0}dp\hfill \end{array}$$

(138)

$$\begin{array}{ccc}& & Y=\sqrt{2}t(p_0,Q)=\frac{1}{\sqrt{2}}\left[\tilde{A}(p_0,Q)-\tilde{A}(p_0,-Q)\right]\hfill \end{array}$$

(139)

We obtained an analytic expression for the imaginary part of the vertex function. We can now calculate the contribution to the density of particles due to quantum Gaussian fluctuations from Equation (88), which by rescaling all the parameters, reads

$$\begin{array}{ccc}\hfill n^{\left(2\right)}& =& \frac{1}{\pi \mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }\hfill \\ \hfill & =& \frac{{\left(4m{T}_c\right)}^{3/2}}{{\pi }^3}{\int }_0^{+\infty }dQ{Q}^2{\int }_{-\infty }^{+\infty }d\nu \frac{1}{e^{4\nu }-1}\frac{\partial \tilde{\delta }(\mathbf{Q},\nu )}{\partial \tilde{\mu }}\equiv \frac{{\left(2m{T}_c\right)}^{3/2}}{3{\pi }^2}{I}_2\left(\tilde{\mu },y\sqrt{\frac{{ϵ}_{\mathrm{F}}}{T_c}}\right)\hfill \end{array}$$

(140)

Writing the number of particles in terms of ${ϵ}_{\mathrm{F}}$, the final number equation can be rewritten as

$${\left(\frac{{ϵ}_{\mathrm{F}}}{T_c}\right)}^{3/2}=I_0\left(\tilde{\mu }\right)+I_2\left(\tilde{\mu },y\sqrt{\frac{{ϵ}_{\mathrm{F}}}{T_c}}\right)$$

(141)

We have, therefore, to solve this equation together with the gap equation Equation (104). The complete set of equations, obtained going beyond the mean field theory, with the inclusion of the full Gaussian fluctuations contained in the function $I_2$, is, therefore,

$$\begin{array}{ccc}& & y{\left(\frac{{ϵ}_{\mathrm{F}}}{T_c}\right)}^{1/2}=I_g\left(\tilde{\mu }\right)\hfill \end{array}$$

(142)

$$\begin{array}{ccc}& & {\left(\frac{{ϵ}_{\mathrm{F}}}{T_c}\right)}^{3/2}=I_0\left(\tilde{\mu }\right)+I_2\left(\tilde{\mu },I_g\left(\tilde{\mu }\right)\right)\hfill \end{array}$$

(143)

This set of equations can be solved numerically with root-finding algorithms, finding the values for $(y,\mu ,T_c)$ which satisfy simultaneously the two equations. The results are shown in Figure 3. In particular Equation (143) is actually an equation of state which put in relation $T_c$ and $\mu$. As one can see from Figure 3, the quantum fluctuations are more important in the BEC regime, namely in the strong coupling limit, as already discussed.

The failure of the mean field theory in correctly reproducing the behavior for the critical temperature along the crossover is due to the fact that, in that scheme, the number equation describes a set of non interacting fermions, disregarding, therefore, the two-body physics which becomes relevant specially when the attraction becomes strong. We showed that going towards strong coupling, the partition function of the system actually tends to that of a system of bosons which can be seen as composed by pairs of fermions.

5. Fermi Gas with Spin-Orbit Interaction, Zeeman Term and Attractive Potential

The aim of the next sections is to see how the critical temperature along the BCS-BEC crossover, described in the previous section, is modified by the presence of the spin–orbit (SO) coupling discussed in section Section 1.2. Let us consider the following Hamiltonian $H=H_0+H_I$, in the momentum space, made by a single-particle term

$$H_0=\sum _{\mathbf{k}}{\widehat{\mathsf{\Psi }}}^{†}\left(\mathbf{k}\right)\left\{\frac{k^2}{2m}\mathbb{I}+v_{\mathrm{R}}({\sigma }_x{k}_y-{\sigma }_y{k}_x)+v_{\mathrm{D}}({\sigma }_x{k}_y+{\sigma }_y{k}_x)+h{\sigma }_z\right\}\widehat{\mathsf{\Psi }}\left(\mathbf{k}\right)$$

(144)

where $\widehat{\mathsf{\Psi }}\left(\mathbf{k}\right)=\left(\begin{array}{c}{\widehat{\psi }}_{\uparrow }\left(\mathbf{k}\right)\\ {\widehat{\psi }}_{\downarrow }\left(\mathbf{k}\right)\end{array}\right)$ and an interacting term

$$H_I=-\frac{g}{\mathcal{V}}\sum _{\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}}{\widehat{\psi }}_{\uparrow }^{†}(\mathbf{k}+\mathbf{q}){\widehat{\psi }}_{\downarrow }^{†}(-\mathbf{k}){\widehat{\psi }}_{\downarrow }(\mathbf{q}-{\mathbf{k}}^{\prime }){\widehat{\psi }}_{\uparrow }(-{\mathbf{k}}^{\prime })$$

(145)

where ${\sigma }_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$ and ${\sigma }_y=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right)$ are Pauli matrices. The new ingredient with respect to the previous case is the spin–orbit interaction, contained in the two terms of the single-particle Hamiltonian, respectively the Rashba term, $v_{\mathrm{R}}({\sigma }_x{k}_y-{\sigma }_y{k}_x)$, Equation (9), and the Dresselhaus term, $v_{\mathrm{D}}({\sigma }_x{k}_y+{\sigma }_y{k}_x)$, Equation (10). We included also an effective artificial Zeeman term $h{\sigma }_z$, Equation (11), which can be tuned by a Raman laser [19]. This Hamiltonian is associated, in the path-integral formalism, by the Grassmann variables $\psi$ and ${\psi }^{†}$, to the partition function

$$\mathcal{Z}=\mathrm{Tr}\left(e^{-\beta (H-\mu N)}\right)=\int D{\psi }^{†}D\psi e^{-S[{\psi }^{†},\psi ]}$$

(146)

with the action

$$\begin{array}{c}S[{\psi }^{†},\psi ]=\sum _{\mathsf{k},s}{\psi }_s^{†}\left(\mathsf{k}\right)\left[i{\omega }_n+\frac{k^2}{2m}-\mu \right]{\psi }_s\left(\mathsf{k}\right)-\frac{g}{\beta \mathcal{V}}\sum _{\mathsf{k},{\mathsf{k}}^{\prime },\mathsf{q}}{\psi }_{\uparrow }^{†}(\mathsf{k}+\mathsf{q}){\psi }_{\downarrow }^{†}(-\mathsf{k}){\psi }_{\downarrow }(\mathsf{q}-{\mathsf{k}}^{\prime }){\psi }_{\uparrow }\left({\mathsf{k}}^{\prime }\right)\hfill \\ \hfill +\sum _{\mathsf{k}}\left({\psi }_{\uparrow }^{†}\left(\mathsf{k}\right),{\psi }_{\downarrow }^{†}\left(\mathsf{k}\right)\right)\left[v_{\mathrm{R}}({\sigma }_x{k}_y-{\sigma }_y{k}_x)+v_{\mathrm{D}}({\sigma }_x{k}_y+{\sigma }_y{k}_x)+h{\sigma }_z\right]\left(\begin{array}{c}{\psi }_{\uparrow }\left(\mathsf{k}\right)\\ {\psi }_{\downarrow }\left(\mathsf{k}\right)\end{array}\right)\end{array}$$

(147)

where ${\psi }_s\left(\mathsf{k}\right)$ and ${\psi }_s^{†}\left(\mathsf{k}\right)$ are Grassmann variables depending on the four-momentum $\mathsf{k}=({\omega }_n,\mathbf{k})$. As already seen, the term with four fermions can be decoupled by means of an Hubbard-Stratonovich transformation, introducing the auxiliary field $\Delta \left(\mathsf{q}\right)$

$$\begin{array}{ccc}\hfill & & e^{\frac{g}{\beta \mathcal{V}}\sum _{\mathsf{k},{\mathsf{k}}^{\prime },\mathsf{q}}{\psi }_{\uparrow }^{†}(\mathsf{k}+\mathsf{q}){\psi }_{\downarrow }^{†}(-\mathsf{k}){\psi }_{\downarrow }(-{\mathsf{k}}^{\prime }){\psi }_{\uparrow }({\mathsf{k}}^{\prime }+\mathsf{q})}\hfill \\ & & =\int D{\Delta }^{*}D\Delta e^{\sum _{\mathsf{q}}\left[-\frac{\beta \mathcal{V}}{g}{|\Delta |}^2+\sum _{\mathsf{k}}\left({\Delta }^{*}\left(\mathsf{q}\right){\psi }_{\downarrow }(-\mathsf{k}){\psi }_{\uparrow }(\mathsf{k}+\mathsf{q})+{\psi }_{\uparrow }^{†}(k+\mathsf{q}){\psi }_{\downarrow }^{†}(-\mathsf{k})\Delta \left(\mathsf{q}\right)\right)\right]}\hfill \end{array}$$

(148)

in this way the partition function takes the form $\mathcal{Z}=\int D{\psi }^{†}D\psi \int D{\Delta }^{*}D\Delta e^{-S[{\psi }^{†},\psi ,\Delta ]}$ where, after defining

$$\gamma \left(\mathbf{k}\right)=v_{\mathrm{R}}(k_x+i{k}_y)+v_{\mathrm{D}}(k_x-i{k}_y)$$

(149)

the new action reads

$$\begin{array}{c}S[{\psi }^{†},\psi ,\Delta ]=\sum _{\mathsf{k},s}{\psi }_s^{†}\left(\mathsf{k}\right)\left[i{\omega }_n+{\xi }_{\mathbf{k}}^{\left(s\right)}\right]{\psi }_s\left(\mathsf{k}\right)+\sum _{\mathsf{k}}{\psi }_{\downarrow }^{†}\left(\mathsf{k}\right){\gamma }^{*}\left(\mathbf{k}\right){\psi }_{\uparrow }\left(\mathsf{k}\right)+\sum _{\mathsf{k}}{\psi }_{\uparrow }^{†}\left(\mathsf{k}\right)\gamma \left(\mathbf{k}\right){\psi }_{\downarrow }\left(\mathsf{k}\right)\hfill \\ \hfill +\sum _{\mathsf{q}}\frac{\beta \mathcal{V}}{g}{|\Delta \left(\mathsf{q}\right)|}^2-\sum _{\mathsf{q},\mathsf{k}}{\Delta }^{*}\left(\mathsf{q}\right){\psi }_{\downarrow }(-\mathsf{k}){\psi }_{\uparrow }(\mathsf{k}+\mathsf{q})-\sum _{\mathsf{q},\mathsf{k}}{\psi }_{\uparrow }^{†}(\mathsf{k}+\mathsf{q}){\psi }_{\downarrow }^{†}(-\mathsf{k})\Delta \left(\mathsf{q}\right)\end{array}$$

(150)

where ${\xi }_{\mathbf{k}}^{(\uparrow )}={\xi }_{\mathbf{k}}+h$ and ${\xi }_{\mathbf{k}}^{(\downarrow )}={\xi }_{\mathbf{k}}-h$. It is now convenient to introduce the following Nambu–Jona-Lasinio multispinors

$$\mathsf{\Psi }\left(\mathsf{k}\right)=\left(\begin{array}{c}{\psi }_{\uparrow }\left(\mathsf{k}\right)\\ {\psi }_{\downarrow }^{†}(-\mathsf{k})\\ {\psi }_{\downarrow }\left(\mathsf{k}\right)\\ {\psi }_{\uparrow }^{†}(-\mathsf{k})\end{array}\right)\mathrm{and}\overline{\mathsf{\Psi }}\left(\mathsf{k}\right)=\left({\psi }_{\uparrow }^{†}\left(\mathsf{k}\right),{\psi }_{\downarrow }(-\mathsf{k}),{\psi }_{\downarrow }^{†}\left(\mathsf{k}\right),{\psi }_{\uparrow }(-\mathsf{k})\right)$$

(151)

so that Equation (150) can be written as

$$S[{\psi }^{†},\psi ,\Delta ]=\frac{\beta \mathcal{V}}{g}\sum _{\mathsf{q}}{|\Delta \left(\mathsf{q}\right)|}^2+\frac{1}{2}\sum _{\mathsf{k},\mathsf{p}}\overline{\mathsf{\Psi }}\left(\mathsf{k}\right){\mathcal{G}}^{-1}(\mathsf{k},\mathsf{p})\mathsf{\Psi }\left(\mathsf{p}\right)+\beta \sum _{\mathbf{k}}{\xi }_{\mathbf{k}}$$

(152)

where the sum $\beta {\sum }_{\mathbf{k}}{\xi }_{\mathbf{k}}$ comes from interchanging the fermionic fields and the corresponding equal-time limiting procedure in order to write the action in a matrix form [77], namely ${\displaystyle \sum _{\mathbf{k},s}}{\widehat{\psi }}_s^{†}\left(\mathbf{k}\right){\xi }_{\mathbf{k}}{\widehat{\psi }}_s\left(\mathbf{k}\right)=\frac{1}{2}{\displaystyle \sum _{\mathbf{k}}}{\xi }_{\mathbf{k}}[{\widehat{\psi }}_{\uparrow }^{†}\left(\mathbf{k}\right){\widehat{\psi }}_{\uparrow }\left(\mathbf{k}\right)-{\widehat{\psi }}_{\uparrow }(-\mathbf{k}){\widehat{\psi }}_{\uparrow }^{†}(-\mathbf{k})+{\widehat{\psi }}_{\downarrow }^{†}\left(\mathbf{k}\right){\widehat{\psi }}_{\downarrow }\left(\mathbf{k}\right)-{\widehat{\psi }}_{\downarrow }(-\mathbf{k}){\widehat{\psi }}_{\downarrow }^{†}(-\mathbf{k})]$$+{\displaystyle \sum _{\mathbf{k}}}{\xi }_{\mathbf{k}}$. In Equation (152) we also introduced the following inverse Green function

$${\mathcal{G}}^{-1}(\mathsf{k},\mathsf{p})=\left(\begin{array}{cccc}(i{\omega }_n+{\xi }_{\mathbf{k}}+h){\delta }_{\mathsf{k},\mathsf{p}}& -\Delta (\mathsf{p}-\mathsf{k})& \gamma \left(\mathbf{k}\right){\delta }_{\mathsf{k},\mathsf{p}}& 0\\ -{\Delta }^{*}(\mathsf{k}-\mathsf{p})& (i{\omega }_n-{\xi }_{\mathbf{k}}+h){\delta }_{\mathsf{k},\mathsf{p}}& 0& -\gamma (-\mathbf{k}){\delta }_{\mathsf{k},\mathsf{p}}\\ {\gamma }^{*}\left(\mathbf{k}\right){\delta }_{\mathsf{k},\mathsf{p}}& 0& (i{\omega }_n+{\xi }_{\mathbf{k}}-h){\delta }_{\mathsf{k},\mathsf{p}}& \Delta (\mathsf{p}-\mathsf{k})\\ 0& -{\gamma }^{*}(-\mathbf{k}){\delta }_{\mathsf{k},\mathsf{p}}& {\Delta }^{*}(\mathsf{k}-\mathsf{p})& (i{\omega }_n-{\xi }_{\mathbf{k}}-h){\delta }_{\mathsf{k},\mathsf{p}}\end{array}\right)$$

(153)

where $\gamma (-\mathbf{k})=-\gamma \left(\mathbf{k}\right)$. Being the action quadratic, we can perform a Gaussian integral in the Grassmann variables getting an effective theory which depends only on the auxiliary complex field

$$\mathcal{Z}=\int D{\psi }^{†}D\psi D{\Delta }^{*}D\Delta e^{-S[{\psi }^{†},\psi ,\Delta ]}\propto \int D{\Delta }^{*}D\Delta e^{S_e[\Delta ]}$$

(154)

obtaining an effective action depending only on the pairing function

$$S_e[\Delta ]=\frac{\beta \mathcal{V}}{g}\sum _{\mathsf{q}}{|\Delta \left(\mathsf{q}\right)|}^2-\frac{1}{2}ln\left[det\left({\mathcal{G}}^{-1}\right)\right]+\beta \sum _{\mathbf{k}}{\xi }_{\mathbf{k}}.$$

(155)

5.1. Gap Equation

As done in Section 3, we first search for the homogeneous pairing which minimizes the action, selecting only the contribution at $\mathsf{q}=0$ in momentum space. Calling $\Delta (\mathsf{q}=0)={\Delta }_0$ the determinant appearing in Equation (155) is simply

$$det\left({\mathcal{G}}^{-1}\right)=\prod _{\mathbf{k},i{\omega }_n,j}(i{\omega }_n-E_{j\mathbf{k}})$$

(156)

where the energies are

$$\begin{array}{ccc}& & E_{1\mathbf{k}}=\sqrt{{\xi }_{\mathbf{k}}^2+{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2-2\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\hfill \end{array}$$

(157)

$$\begin{array}{ccc}& & E_{2\mathbf{k}}=\sqrt{{\xi }_{\mathbf{k}}^2+{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2+2\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\hfill \end{array}$$

(158)

and where

$${\gamma }_{h\mathbf{k}}=\sqrt{{|\gamma \left(\mathbf{k}\right)|}^2+h^2}=\sqrt{{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2{k}_y^2+{(v_{\mathrm{R}}-v_{\mathrm{D}})}^2{k}_x^2+h^2}$$

(159)

with $\gamma \left(\mathbf{k}\right)=v_{\mathrm{R}}(k_y+i{k}_x)+v_{\mathrm{D}}(k_y-i{k}_x)$, Equation (149), and ${\xi }_{\mathbf{k}}=k^2/2m-\mu$. $v_{\mathrm{R}}$ and $v_{\mathrm{D}}$ are the Rashba and Dresselhaus velocities and h is the Zeeman field. The effective action, therefore, takes the form

$$S_e\left[{\Delta }_0\right]=\frac{\beta \mathcal{V}}{g}{|{\Delta }_0|}^2-\frac{1}{2}\sum _{\mathbf{k},i{\omega }_n,j}ln\left(i{\omega }_n-E_{j\mathbf{k}}\right)+\beta \sum _{\mathbf{k}}{\xi }_{\mathbf{k}}$$

(160)

The saddle point equation sets the mean field solution and is given by the condition $\frac{\partial S_e}{\partial {\Delta }_0}=0$

$$\frac{{\Delta }_0^{*}}{g}=-\frac{1}{2\beta \mathcal{V}}\sum _{\mathbf{k},i{\omega }_n,j}\frac{1}{i{\omega }_n-E_{j\mathbf{k}}}\frac{\partial E_{j\mathbf{k}}}{\partial {\Delta }_0}$$

(161)

After summing over the Matsubara frequencies the gap equation, in the presence of spin–orbit interaction and Zeeman term, reads

$$\begin{array}{c}\hfill \frac{1}{g}=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{tanh(\beta E_{1\mathbf{k}}/2)}{E_{1\mathbf{k}}}\left(1-\frac{h^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)\right.\\ \hfill \left.+\frac{tanh(\beta E_{2\mathbf{k}}/2)}{E_{2\mathbf{k}}}\left(1+\frac{h^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)\right]\end{array}$$

(162)

5.2. Number Equation

As seen in the previous case, from the full action we can not calculate the partition function exactly. The approximation scheme that we will employ is based again on the expansion of the action in the gap field, as in Equation (62). The zeroth order is equivalent to the mean field approximation. We will then include the fluctuations at the Gaussian level.

5.2.1. Mean Field

The effective action, at mean field level, is the action in Equation (160), evaluated at the homogeneous solution of the gap equation, corresponding to the saddle point approximation. The associated grand canonical potential is simply $\mathsf{\Omega }[{\Delta }_0,]=\beta S_e\left[{\Delta }_0\right]$. The mean number of particles, defined by $N=-\frac{\partial \mathsf{\Omega }}{\partial \mu }$ is, therefore, given by

$$N=\frac{1}{2\beta }\frac{\partial }{\partial \mu }\sum _{\mathbf{k},i{\omega }_n,j}ln\left(i{\omega }_n-E_{j\mathbf{k}}\right)-\frac{\partial }{\partial \mu }\sum _{\mathbf{k}}{\xi }_{\mathbf{k}}=\sum _{\mathbf{k}}\left(1-\frac{1}{2\beta }\sum _{i{\omega }_n,j}\frac{1}{i{\omega }_n-E_{j\mathbf{k}}}\frac{\partial E_{j\mathbf{k}}}{\partial \mu }\right)$$

(163)

Summing over the Matsubara frequencies, the number equation, at the mean-field level, reads

$$\begin{array}{c}\hfill N=\sum _{\mathbf{k}}\left[1-\frac{tanh(\beta E_{1\mathbf{k}}/2)}{2E_{1\mathbf{k}}}\left({\xi }_{\mathbf{k}}-\frac{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)\right.\\ \hfill \left.-\frac{tanh(\beta E_{2\mathbf{k}}/2)}{2E_{2\mathbf{k}}}\left({\xi }_{\mathbf{k}}+\frac{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)\right]\end{array}$$

(164)

Let us consider in what follows the cases at $T=0$ and $T=T_c$.

• At $T=0$.

In order to study the behavior of the gap and of the chemical potential as functions of the coupling, at $T=0$, the equations to solve at mean field level are

$$\begin{array}{ccc}\hfill & & \frac{1}{g}=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{1}{E_{1\mathbf{k}}}\left(1-\frac{h^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)+\frac{1}{E_{2\mathbf{k}}}\left(1+\frac{h^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)\right]\hfill \end{array}$$

(165)

$$\begin{array}{ccc}& & n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[1-\frac{1}{2E_{1\mathbf{k}}}\left({\xi }_{\mathbf{k}}-\frac{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)-\frac{1}{2E_{2\mathbf{k}}}\left({\xi }_{\mathbf{k}}+\frac{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}^2}{\sqrt{{\xi }_{\mathbf{k}}^2{\gamma }_{h\mathbf{k}}^2+{|{\Delta }_0|}^2{h}^2}}\right)\right]\hfill \end{array}$$

(166)

where $n=N/\mathcal{V}$ is the particle density. This set of equations, for $h=0$, has been solved in Ref. [32]. We will focus, instead, on the finite temperature case, in particular at the critical point.

• At T=$T_c$.

If we want to calculate the critical temperature $T_c$, putting by definition ${\Delta }_0\left(T_c\right)=0$, at the mean field level the equation to solve in terms of $T_c$ are the following

$$\begin{array}{ccc}\hfill & & \frac{1}{g}=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}\left(1-\frac{h^2}{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}}\right)+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}\left(1+\frac{h^2}{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}}\right)\right]\hfill \end{array}$$

(167)

$$\begin{array}{ccc}\hfill & & n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left\{1-\frac{1}{2}\left[tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}{2T_c}\right)\right]\right\}\hfill \end{array}$$

(168)

These equations should be solved in order to obtain the critical temperature and the chemical potential, at $T_c$, as functions of the coupling g. The latter can be also written as

$$n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left(1-\frac{sinh({\xi }_{\mathbf{k}}/T_c)}{cosh({\xi }_{\mathbf{k}}/T_c)+cosh({\gamma }_{h\mathbf{k}}/T_c)}\right).$$

(169)

5.2.2. Inclusion of Gaussian Fluctuations at $T=T_c$

The inadequacy of the mean field approach, also with SO coupling, in the strong coupling regime will be clarified later, and its reason is the same as before, based on the fact that the number equation in Equation (168) is that of non-interacting fermions. It is, therefore, necessary to extend the analysis including Gaussian fluctuations around the mean field solution. We will get a second order term in the action, $S_e[\Delta ]=S_e^{\left(0\right)}+S_e^{\left(2\right)}$, in addition to $S_e^{\left(0\right)}$, the action of the non-interacting system but in the presence of spin–orbit and Zeeman terms. As already seen, the following action has already a quadratic term while the trace has to be expanded,

$$S_e[\Delta ]=\int d^{3}r{\int }_0^{\beta }d\tau \frac{|\Delta (\overrightarrow{r},\tau ){|}^2}{g}-\frac{1}{2}\mathrm{Tr}(ln\left({\mathcal{G}}^{-1}\right))+\beta \sum _{\mathbf{k}}{\xi }_{\mathbf{k}}.$$

(170)

The inverse of the Green function can be written in momentum space and in the Matsubara frequency space as in Equation (153). It can be written as

$${\mathcal{G}}^{-1}(\mathsf{k},\mathsf{p})=G^{-1}(\mathsf{k},\mathsf{p})+\widehat{\Delta }(\mathsf{k},\mathsf{p})$$

(171)

where

$$\begin{array}{ccc}\hfill G^{-1}(\mathsf{k},\mathsf{p})& =& \left(\begin{array}{cccc}(i{\omega }_n+{\xi }_{\mathbf{k}}+h)& 0& \gamma \left(\mathbf{k}\right)& 0\\ 0& (i{\omega }_n-{\xi }_{\mathbf{k}}+h)& 0& \gamma \left(\mathbf{k}\right)\\ {\gamma }^{*}\left(\mathbf{k}\right)& 0& (i{\omega }_n+{\xi }_{\mathbf{k}}-h)& 0\\ 0& {\gamma }^{*}\left(\mathbf{k}\right)& 0& (i{\omega }_n-{\xi }_{\mathbf{k}}-h)\end{array}\right){\delta }_{\mathsf{k},\mathsf{p}}\hfill \\ \hfill \widehat{\Delta }(\mathsf{k},\mathsf{p})& =& \left(\begin{array}{cccc}0& -\Delta (\mathsf{p}-\mathsf{k})& 0& 0\\ -{\Delta }^{*}(\mathsf{k}-\mathsf{p})& 0& 0& 0\\ 0& 0& 0& \Delta (\mathsf{p}-\mathsf{k})\\ 0& 0& {\Delta }^{*}(\mathsf{k}-\mathsf{p})& 0\end{array}\right)\hfill \end{array}$$

We can now perform the series expansion of the logarithmic function

$$ln\left({\mathcal{G}}^{-1}\right)=ln(G^{-1}+\widehat{\Delta })=ln\left(G^{-1}(\mathbb{I}+G\widehat{\Delta })\right)=ln\left(G^{-1}\right)+G\widehat{\Delta }-\frac{1}{2}G\widehat{\Delta }G\widehat{\Delta }+O\left({\Delta }^3\right)$$

(172)

The first term gives the action of the non-interacting system in the presence of spin–orbit and Zeeman terms. The second term has a vanishing trace, as well as all the terms with odd powers, $\mathrm{Tr}{\left(G\widehat{\Delta }\right)}^{2n+1}=0$. The action at second order in $\widehat{\Delta }$ is, then, given by

$$S_e[\Delta ]=S_e^{\left(0\right)}+\int d^{3}r{\int }_0^{\beta }d\tau \frac{|\Delta (\overrightarrow{r},\tau ){|}^2}{g}+\frac{1}{4}\mathrm{Tr}\left(G\widehat{\Delta }G\widehat{\Delta }\right)+o\left({\Delta }^3\right)$$

(173)

The inverse of $G^{-1}$, the non-interacting Green function, in momentum space, is given by $G(\mathsf{k},\mathsf{p})=G\left(\mathsf{k}\right){\delta }_{\mathsf{k},\mathsf{p}}$ where

$$G\left(\mathsf{k}\right)=\left(\begin{array}{cccc}\frac{i{\omega }_n+{\xi }_{\mathbf{k}}-h}{{(i{\omega }_n+{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}& 0& \frac{-\gamma \left(\mathbf{k}\right)}{{(i{\omega }_n+{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}& 0\\ 0& \frac{i{\omega }_n-{\xi }_{\mathbf{k}}-h}{{(i{\omega }_n-{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}& 0& \frac{-\gamma \left(\mathbf{k}\right)}{{(i{\omega }_n-{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}\\ \frac{-{\gamma }^{*}\left(\mathbf{k}\right)}{{(i{\omega }_n+{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}& 0& \frac{i{\omega }_n+{\xi }_{\mathbf{k}}+h}{{(i{\omega }_n+{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}& 0\\ 0& \frac{-{\gamma }^{*}\left(\mathbf{k}\right)}{{(i{\omega }_n-{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}& 0& \frac{i{\omega }_n-{\xi }_{\mathbf{k}}+h}{{(i{\omega }_n-{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2}\end{array}\right)$$

(174)

Making the matrix product $G\widehat{\Delta }G\widehat{\Delta }$ and taking the trace we get

$$\frac{1}{4}\mathrm{Tr}\left(G\widehat{\Delta }G\widehat{\Delta }\right)=\frac{1}{2}\mathrm{Tr}\left(G_{22}{\Delta }^{*}{G}_{11}\Delta -G_{24}{\Delta }^{*}{G}_{31}\Delta +G_{44}{\Delta }^{*}{G}_{33}\Delta -G_{42}{\Delta }^{*}{G}_{13}\Delta \right)$$

(175)

where $G_{ij}$ are the matrix elements of G. In momentum space Equation (175) reads

$$\frac{1}{4}\mathrm{Tr}\left(G\widehat{\Delta }G\widehat{\Delta }\right)=\sum _{\mathsf{q}}{\Delta }^{*}\left(\mathsf{q}\right)\chi \left(\mathsf{q}\right)\Delta \left(\mathsf{q}\right)$$

(176)

where we define the function

$$\begin{array}{c}\hfill \chi \left(\mathsf{q}\right)=\frac{1}{2\beta \mathcal{V}}\sum _{\mathsf{k}}[G_{11}\left(\mathsf{k}\right)G_{22}(\mathsf{k}-\mathsf{q})+G_{33}\left(\mathsf{k}\right)G_{44}(\mathsf{k}-\mathsf{q})\\ \hfill -G_{13}\left(\mathsf{k}\right)G_{42}(\mathsf{k}-\mathsf{q})-G_{31}\left(\mathsf{k}\right)G_{24}(\mathsf{k}-\mathsf{q})]\end{array}$$

(177)

The Fourier transform of the first term in Equation (173) is $\int d^{3}r{\int }_0^{\beta }d\tau \frac{|\Delta (\overrightarrow{r},\tau ){|}^2}{g}={\sum }_{\mathsf{q}}\frac{{|\Delta \left(\mathsf{q}\right)|}^2}{g}$, therefore the action, at the second order in $\Delta$, namely at the Gaussian level, can be written as

$$S_e[\Delta ]=S_e^{\left(0\right)}+\sum _{\mathsf{q}}{\Delta }^{*}\left(\mathsf{q}\right){\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)\Delta \left(\mathsf{q}\right)+o\left({\Delta }^3\right)$$

(178)

where the same notation for ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$, used for the case without SO coupling, has been adopted

$${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)=\frac{1}{g}+\chi \left(\mathsf{q}\right)$$

(179)

but where now $\chi \left(\mathsf{q}\right)\equiv \chi (\mathbf{q},i{\nu }_n)$ is more complicated than before. From Equations (174) and (177) we get

$$\chi (\mathbf{q},i{\nu }_n)=\frac{1}{\beta \mathcal{V}}\sum _{\mathbf{k},{\omega }_n}\frac{(i{\omega }_n+{\xi }_{\mathbf{k}})(i{\omega }_n+i{\nu }_m-{\xi }_{\mathbf{k}-\mathbf{q}})+h^2-\mathrm{Re}\left[{\gamma }^{*}\left(\mathbf{k}\right)\gamma \left(\mathbf{k}-\mathbf{q}\right)\right]}{\left({(i{\omega }_n+{\xi }_{\mathbf{k}})}^2-{\gamma }_{h\mathbf{k}}^2\right)\left({(i{\omega }_n-i{\nu }_n-{\xi }_{\mathbf{k}-\mathbf{q}})}^2-{\gamma }_{h\mathbf{k}-\mathbf{q}}^2\right)}$$

(180)

It is convenient to make the shift $\mathbf{k}\to \mathbf{k}+\frac{\mathbf{q}}{2}$ and, after performing the sum over the Matsubara frequencies, we obtain

$$\begin{array}{c}\hfill \chi (\mathbf{q},i{\nu }_n)=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left\{\left(1+C_{\mathbf{k},\mathbf{q}}^h\right)\frac{1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right.\\ \hfill +\left(1+C_{\mathbf{k},\mathbf{q}}^h\right)\frac{1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\\ \hfill +\left(1-C_{\mathbf{k},\mathbf{q}}^h\right)\frac{1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\\ \hfill \left.+\left(1-C_{\mathbf{k},\mathbf{q}}^h\right)\frac{1-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}\right)-n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right\}\end{array}$$

(181)

where

$$\begin{array}{ccc}& & C_{\mathbf{k},\mathbf{q}}^h=\frac{\mathrm{Re}\left[{\gamma }^{*}(\mathbf{k}+\frac{\mathbf{q}}{2})\gamma (\mathbf{k}-\frac{\mathbf{q}}{2})\right]-h^2}{{\gamma }_{h\mathbf{k}+\mathbf{q}/2}{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\hfill \\ \hfill & & =\frac{{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2(k_y^2-\frac{q_y^2}{4})+{(v_{\mathrm{R}}-v_{\mathrm{D}})}^2(k_x^2-\frac{q_x^2}{4})-h^2}{{\gamma }_{h\mathbf{k}+\mathbf{q}/2}{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\hfill \end{array}$$

(182)

with

$${\gamma }_{h\mathbf{k}\pm \mathbf{q}/2}=\sqrt{{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2{\left(k_y\pm q_y/2\right)}^2+{(v_{\mathrm{R}}-v_{\mathrm{D}})}^2{\left(k_x\pm q_x/2\right)}^2+h^2}$$

(183)

For $v_{\mathrm{D}}=0$, namely with only the Rashba coupling, we have

$$C_{\mathbf{k},\mathbf{q}}^h=\frac{|{\mathbf{k}}_{\perp }{|}^2-{|{\mathbf{q}}_{\perp }/2|}^2-{(h/v_{\mathrm{R}})}^2}{\sqrt{[|{\mathbf{k}}_{\perp }-{\mathbf{q}}_{\perp }{/2|}^2+{(h/v_{\mathrm{R}})}^2\left]\right[|{\mathbf{k}}_{\perp }+{\mathbf{q}}_{\perp }/2{|}^2+{(h/v_{\mathrm{R}})}^2]}}$$

(184)

where ${\mathbf{k}}_{\perp }=(k_x,k_y,0)$ are momenta perpendicular to the z-direction, the Zeeman direction.

For $v_{\mathrm{R}}=v_{\mathrm{D}}=v/2$, namely for equal Rashba and Dresselhaus couplings, we have, instead,

$$C_{\mathbf{k},\mathbf{q}}^h=\frac{v^2(k_y^2-q_y^2/4)-h^2}{\sqrt{[v^2{(k_y-q_y/2)}^2+h^2][v^2{(k_y+q_y/2)}^2+h^2]}}$$

(185)

which, for $h=0$, namely without the Zeeman field, it reduces to simply $C_{\mathbf{k},\mathbf{q}}^0=1$.

Employing the symmetry $\mathbf{k}\to -\mathbf{k}$ under the sum we can simplify Equation (182) as it follows

$$\begin{array}{c}\hfill \chi (\mathbf{q},i{\nu }_n)=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left\{\left(1+C_{\mathbf{k},\mathbf{q}}^h\right)\left[\frac{1-2n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right.\right.\\ \hfill \left.+\frac{1-2n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right]\\ \hfill \left.+\left(1-C_{\mathbf{k},\mathbf{q}}^h\right)\frac{2-2n_{\mathrm{F}}\left({\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}^h\right)-2n_{\mathrm{F}}\left({\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}^h\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}^h-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}^h}\right\}\end{array}$$

(186)

or, in terms of the hyperbolic tangent,

$$\begin{array}{c}\hfill \chi (\mathbf{q},i{\nu }_n)=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left\{\left(1+C_{\mathbf{k},\mathbf{q}}^h\right)\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}}{2T}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right.\right.\\ \hfill \left.+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}}{2T}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right]\\ \hfill \left.+\left(1-C_{\mathbf{k},\mathbf{q}}^h\right)\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}}{2T}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}{2T}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right\}\end{array}$$

(187)

We found, therefore, that ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$ is much more complicated that that obtained in the previous standard case without spin–orbit coupling. However, the partition function has the same form as that shown in Equation (79), and, therefore, it is possible to apply the same procedure seen before (in Section 3.2.2) to find the grand canonical potential and the density of particles. This is again due to the fact that the analytical continuation of ${\mathsf{\Gamma }}^{-1}(\mathbf{q},z)$ has only a branch-cut of simple poles in the real axis. The density of particles is, then, given by

$$n=n^{\left(0\right)}+n^{\left(2\right)}$$

(188)

where $n^{\left(0\right)}$ is the mean field contribution as in Equation (168) while the second order term is given by

$$n^{\left(2\right)}=\frac{1}{\beta \mathcal{V}}\frac{\partial }{\partial \mu }\mathrm{Tr}(ln\left({\mathsf{\Gamma }}^{-1}\right))=\frac{1}{\pi \mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }$$

(189)

where the phase shift $\delta (\mathbf{q},\omega )$ is defined as before

$$\delta (\mathbf{q},\omega )=\underset{ϵ\to 0}{lim}\mathrm{Arg}\left(\mathsf{\Gamma }{(\mathbf{q},\omega +iϵ)}^{-1}\right)$$

(190)

after a Wick rotation $i{\nu }_n\to \omega +iϵ$.

5.3. Critical Point

Taking the static and homogeneous limit of Equation (188), after noticing that, at $\mathbf{q}=0$,

$$C_{\mathbf{k},0}^h=1-\frac{2h^2}{{\gamma }_{h\mathbf{k}}^2}$$

(191)

we have that, at $T=T_c$, $\chi (0,0)$ is simply given by

$$\begin{array}{ccc}\hfill \chi (0,0)& =& -\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left\{\left(1-\frac{h^2}{{\gamma }_{h\mathbf{k}}^2}\right)\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}\right]\right.\hfill \\ & & +\left.\frac{h^2}{{\gamma }_{h\mathbf{k}}^2}\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}}\right]\right\}\hfill \end{array}$$

(192)

We find that the equation ${\mathsf{\Gamma }}^{-1}\left(0\right)=0$, namely

$$\frac{1}{g}=-\chi (0,0)$$

(193)

is exactly the same equation reported in Equation (167) after reshuffling the terms. The set of equations that has to be solved, at the Gaussian level, in order to find $T_c$ and $\mu$ at $T_c$ is given by Equation (193), together with Equation (188).

6. BCS-BEC Crossover with Spin-Orbit Coupling

For the same reason discussed before, the introduction of a contact potential gives rise to divergences that can be regularized by the same procedure. For simplicity we will focus on the case without Zeeman term, $h=0$. As shown in Appendix D.1 and in Refs. [44,45], the presence of a SO coupling does not alter the form of the renormalization condition for a contact pseudo-potential. The only caution to be taken is to replace $a_s$ with $a_r$, a scattering length which depends both on the interatomic potential and on the SO term

$$\frac{m}{4\pi a_r}=-\frac{1}{g}+\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{1}{2{ϵ}_{\mathbf{k}}}$$

(194)

therefore, the arguments used to implement the BCS and BEC limits in terms of the scattering length remains unchanged.

6.1. Mean Field

At the mean field level the equations to solve along the crossover are Equations (167) and (168). Implementing Equation (194), for $h=0$, the equations useful to derive the critical temperature $T_c$ and the chemical potential $\mu$ within the mean field theory become

$$\begin{array}{cc}\hfill -\frac{m}{4\pi a_r}=& \frac{1}{4\mathcal{V}}{\sum }_{\mathbf{k}}\left\{\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}-\frac{2}{{ϵ}_{\mathbf{k}}}\right\}\end{array}$$

(195)

$$\begin{array}{cc}\hfill n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}& \left\{1-\frac{1}{2}tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}{2T_c}\right)-\frac{1}{2}tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}{2T_c}\right)\right\}\end{array}$$

(196)

where ${\gamma }_{\mathbf{k}}\equiv {\gamma }_{0\mathbf{k}}$, as reported in Equation (159) with $h=0$. Additionally in this case the number equation is the same as that of free fermions, therefore we expect that the mean field level will fail in correctly reproducing the thermodynamic quantities in the strong coupling limit.

Let us consider Equation (195). Contrary to the case without SO, the mass cannot be eliminated from the equations by simply redefining the momenta, but it can be incorporated in the couplings ${\tilde{v}}_{\mathrm{D}}=\frac{v_{\mathrm{D}}}{v_{\mathrm{F}}}=\frac{m{v}_{\mathrm{D}}}{k_{\mathrm{F}}}$ and similarly for $v_{\mathrm{R}}$. In the continuum limit, with the natural change of variables and substitutions $\overrightarrow{p}=\frac{1}{k_{\mathrm{F}}}\overrightarrow{k}$, $\tilde{\mu }=\frac{\mu }{{ϵ}_{\mathrm{F}}}$, ${\tilde{T}}_c=\frac{T_c}{{ϵ}_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{R}}=\frac{v_{\mathrm{R}}}{v_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{D}}=\frac{v_{\mathrm{D}}}{v_{\mathrm{F}}}$, ${\tilde{\gamma }}_{\mathbf{p}}({\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})=2{\gamma }_{\mathbf{k}}(v_{\mathrm{R}},v_{\mathrm{D}})$, the gap equation can be written as

$$-\frac{m}{4\pi a_r}=\frac{k_{\mathrm{F}}^3}{4{\left(2\pi \right)}^3{ϵ}_{\mathrm{F}}}\int d^{3}p\left\{\frac{tanh\left(\frac{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)}{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}+\frac{tanh\left(\frac{p^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)}{p^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}-\frac{2}{p^2}\right\}$$

(197)

while the number equation, Equation (196), becomes

$$n=k_{\mathrm{F}}^3\int \frac{d^{3}p}{{\left(2\pi \right)}^3}\left[1-\frac{1}{2}tanh\left(\frac{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)-\frac{1}{2}tanh\left(\frac{p^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)\right]$$

(198)

6.1.1. Weak Coupling

At low temperature we expect that the chemical potential of a system of weakly interacting fermions is almost constant

$$\mu \underset{T\to 0}{\approx }{E}_{\mathrm{F}}({ϵ}_{\mathrm{F}},v_{\mathrm{R}},v_{\mathrm{D}})$$

(199)

and therefore, in this regime, $\frac{\mu }{T}\gg 1$. Actually, one can shown perturbativly [32,43] that $\mu \approx {ϵ}_{\mathrm{F}}-m{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2/2$. In this limit, the number equation can be neglected due to the ansatz Equation (199), therefore Equation (197) is enough to self-consistently determine the critical temperature $T_c$ as a function of $a_r$. As we will see, the critical temperature, in the BCS part of the crossover, exhibits an exponential behavior as in the case without SO, however $T_c$ will be considerably improved by turning on the spin–orbit coupling. We will also see that, in this weak coupling limit $T_c$ is almost unaltered by the inclusion of Gaussian fluctuations.

6.1.2. Strong Coupling

In the strong coupling limit, supposing that $\mu /T_c\to -\infty$, which can be verified a posteriori, Equation (197) can be used to determine the chemical potential of the system, in fact the dependence in the equations on the critical temperature disappears making the following approximation

$$tanh\left(\frac{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)\approx 1$$

(200)

so that Equation (197), after introducing the dimensionless coupling $y=\frac{1}{k_{\mathrm{F}}{a}_r}$, becomes simply

$$y=-\frac{2}{{\pi }^2}{\int }_0^{\infty }d{p}_x{\int }_0^{\infty }d{p}_y{\int }_0^{\infty }d{p}_z\left\{\frac{1}{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}+\frac{1}{p^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}-\frac{2}{p^2}\right\}$$

(201)

where

$${\tilde{\gamma }}_{\mathbf{p}}=2\sqrt{{({\tilde{v}}_{\mathrm{R}}+{\tilde{v}}_{\mathrm{D}})}^2{p}_x^2+{({\tilde{v}}_{\mathrm{R}}-{\tilde{v}}_{\mathrm{D}})}^2{p}_y^2}$$

(202)

We can integrate over $p_z$ getting

$$y=\frac{1}{\pi }{\int }_0^{\infty }d{p}_x{\int }_0^{\infty }d{p}_y\left\{\frac{2}{\sqrt{p_x^2+p_y^2}}-\frac{1}{\sqrt{p_x^2+p_y^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}}-\frac{1}{\sqrt{p_x^2+p_y^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}}\right\}$$

(203)

For simplicity, let us consider the case with pure Rashba coupling, so that we have ${\tilde{\gamma }}_{\mathbf{p}}=2{\tilde{v}}_{\mathrm{R}}\sqrt{p_x^2+p_y^2}\equiv 2{\tilde{v}}_{\mathrm{R}}|p_{\perp }|$. The above equation can be written as

$$y=-\frac{1}{2}{\int }_0^{\infty }d{p}_{\perp }{p}_{\perp }\left\{\frac{1}{\sqrt{p_{\perp }^2-2{\tilde{v}}_{\mathrm{R}}{p}_{\perp }-\tilde{\mu }}}+\frac{1}{\sqrt{p_{\perp }^2+2{\tilde{v}}_{\mathrm{R}}{p}_{\perp }-\tilde{\mu }}}-\frac{2}{p_{\perp }}\right\}$$

(204)

With the substitutions $p_{\perp }\to p_{\perp }\pm {\tilde{v}}_{\mathrm{R}}$ in the first two integrals we have

$$y=-\frac{1}{2}\underset{\Lambda \to \infty }{lim}\left\{{\int }_{-{\tilde{v}}_{\mathrm{R}}}^{\Lambda }d{p}_{\perp }\frac{p_{\perp }+{\tilde{v}}_{\mathrm{R}}}{\sqrt{p_{\perp }^2-{\tilde{v}}_{\mathrm{R}}^2-\tilde{\mu }}}+{\int }_{{\tilde{v}}_{\mathrm{R}}}^{\Lambda }d{p}_{\perp }\frac{p_{\perp }-{\tilde{v}}_{\mathrm{R}}}{\sqrt{p_{\perp }^2-{\tilde{v}}_{\mathrm{R}}^2-\tilde{\mu }}}-2{\int }_0^{\Lambda }d{p}_{\perp }\right\}$$

(205)

where $\Lambda$ is the ultraviolet cut-off which can be sent to infinity, $\Lambda \to \infty$. Equation (205) can be rewritten as

$$y=\underset{\Lambda \to \infty }{lim}\left\{\Lambda -{\int }_{{\tilde{v}}_{\mathrm{R}}}^{\Lambda }d{p}_{\perp }\frac{p_{\perp }}{\sqrt{p_{\perp }^2-{\tilde{v}}_{\mathrm{R}}^2-\tilde{\mu }}}\right\}-\frac{{\tilde{v}}_{\mathrm{R}}}{2}{\int }_{-{\tilde{v}}_{\mathrm{R}}}^{{\tilde{v}}_{\mathrm{R}}}d{p}_{\perp }\frac{1}{\sqrt{p_{\perp }^2-{\tilde{v}}_{\mathrm{R}}^2-\tilde{\mu }}}$$

(206)

which gives

$$y=\sqrt{-\tilde{\mu }}-{\tilde{v}}_{\mathrm{R}}arctanh\left(\frac{{\tilde{v}}_{\mathrm{R}}}{\sqrt{-\tilde{\mu }}}\right)\approx \sqrt{-\tilde{\mu }}-\frac{{\tilde{v}}_{\mathrm{R}}^2}{\sqrt{-\tilde{\mu }}}$$

(207)

Taking the square, $y^2\approx -\tilde{\mu }-2{\tilde{v}}_{\mathrm{R}}^2$, and using the definitions $\tilde{\mu }=\mu /{ϵ}_{\mathrm{F}}$ and ${\tilde{v}}_{\mathrm{R}}=v_{\mathrm{R}}/v_{\mathrm{F}}$, reminding that $y=\frac{1}{k_{\mathrm{F}}{a}_r}$, we get

$$\mu =-\frac{1}{2m{a}_r^2}-m{v}_{\mathrm{R}}^2+o\left(a_r\right)$$

(208)

Therefore, in the strong coupling limit (BEC limit), $1/a_r\to \infty$, we found $\mu \simeq -\frac{1}{2m{a}_r^2}$, namely the spin–orbit interaction can be neglected and we get the same expression for the chemical potential as that obtained previously in Equation (102), but with $a_r$ instead of $a_s$. We expect that also in this case the critical temperature in the strong coupling limit at the mean field level is not accurate therefore we will not discuss it here in details.

6.1.3. Full Crossover

Let us consider Equation (195), whose form, in the continuum limit and with rescaled dimensionless parameters, $\overrightarrow{p}=\frac{1}{k_{\mathrm{F}}}\overrightarrow{k}$, $\tilde{\mu }=\frac{\mu }{{ϵ}_{\mathrm{F}}}$, ${\tilde{T}}_c=\frac{T_c}{{ϵ}_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{R}}=\frac{v_{\mathrm{R}}}{v_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{D}}=\frac{v_{\mathrm{D}}}{v_{\mathrm{F}}}$, ${\tilde{\gamma }}_{\mathbf{p}}({\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})=2{\gamma }_{\mathbf{k}}(v_{\mathrm{R}},v_{\mathrm{D}})$, is reported in Equation (197), which can be rewritten in terms of $y=\frac{1}{k_{\mathrm{F}}{a}_r}$ as it follows

$$y=-\int \frac{d^{3}p}{{\left(2\pi \right)}^2}\left[\frac{tanh\left(\frac{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)}{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}+\frac{tanh\left(\frac{p^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)}{p^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}-\frac{2}{p^2}\right]\equiv I_y^{\mathrm{SO}}(\tilde{\mu },{\tilde{T}}_c,{\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})$$

(209)

The number equation, Equation (196), in the continuum limit and in terms of the same dimensionless parameters, is

$$\begin{array}{ccc}\hfill n& =& \frac{k_{\mathrm{F}}^3}{3{\pi }^2}\left\{3{\pi }^2\int \frac{d^{3}p}{{\left(2\pi \right)}^3}\left[1-\frac{1}{2}tanh\left(\frac{p^2-\tilde{\mu }-{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)-\frac{1}{2}tanh\left(\frac{p^2-\tilde{\mu }+{\tilde{\gamma }}_{\mathbf{p}}}{2{\tilde{T}}_c}\right)\right]\right\}\hfill \\ & \equiv & \frac{k_{\mathrm{F}}^3}{3{\pi }^2}{I}_n^{\mathrm{SO}}(\tilde{\mu },{\tilde{T}}_c,{\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})\hfill \end{array}$$

(210)

Reminding that the density of particles is $n=\frac{{\left(2m{ϵ}_{\mathrm{F}}\right)}^{3/2}}{3{\pi }^2}=\frac{k_{\mathrm{F}}^3}{3{\pi }^2}$, the equations to solve are

$$\begin{array}{c}\hfill y=I_y^{\mathrm{SO}}(\tilde{\mu },{\tilde{T}}_c,{\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})\end{array}$$

(211)

$$\begin{array}{c}\hfill 1=I_n^{\mathrm{SO}}(\tilde{\mu },{\tilde{T}}_c,{\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})\end{array}$$

(212)

We can search for the solutions of these two equations, Equations (211) and (212), resorting to a multidimensional root-finding algorithm: for any fixed values of the couplings y, ${\tilde{v}}_{\mathrm{R}}$ and ${\tilde{v}}_{\mathrm{D}}$ we look for the critical temperature $T_c$ and the chemical potential $\mu$ that satisfy those equations simultaneously. In Figure 4 and Figure 5 we report some results for the case with only Rashba coupling.

6.1.4. Special Case: Equal Rashba and Dresselhaus

For $v_{\mathrm{R}}=v_{\mathrm{D}}\equiv v/2$ we have ${\gamma }_{\mathbf{k}}=v|k_y|$, or, in terms of the rescaled quantities $p_y=k_y/k_{\mathrm{F}}$ and $\tilde{v}=v/v_{\mathrm{F}}$,

$${\tilde{\gamma }}_{\mathbf{p}}=2\tilde{v}|p_y|$$

(213)

Both in Equation (209) and in Equation (210) there are functions like $F(p_y^2\pm 2v|p_y|)$ under the integral

$$\frac{1}{2}{\int }_{-\infty }^{+\infty }d{p}_y\left[F\left(p_y^2-2\tilde{v}|p_y|-\tilde{\mu }\right)+F\left(p_y^2+2\tilde{v}|p_y|-\tilde{\mu }\right)\right]$$

(214)

in a symmetric combination so that we can remove the absolute value symbols

$$\frac{1}{2}{\int }_{-\infty }^{+\infty }d{p}_y\left[F\left(p_y^2-2\tilde{v}{p}_y-\tilde{\mu }\right)+F\left(p_y^2+2\tilde{v}{p}_y-\tilde{\mu }\right)\right]$$

(215)

so that after changing the variables $p_y\to p_y\pm \tilde{v}$ we get simply

$${\int }_{-\infty }^{+\infty }d{p}_{y}F\left(p_y^2-{\tilde{v}}^2-\tilde{\mu }\right)$$

(216)

As a general result, for the case with equal Rashba and Dresselhaus couplings at the mean field level the effect of the spin–orbit interaction is just a rigid shift of the chemical potential. The critical temperature $T_c$, is therefore the same as that without SO (see Figure 1) while $\mu$ at finite v is linked to that without SO ($v=0$) in the following way

$$\frac{\mu \left(v\right)}{{ϵ}_{\mathrm{F}}}=\frac{\mu \left(0\right)}{{ϵ}_{\mathrm{F}}}-\frac{v^2}{v_{\mathrm{F}}^2}$$

(217)

which means $\mu \left(v\right)=\mu \left(0\right)-m{v}^2/2$, in agreement with the result obtained by gauge-transforming the fermionic fields [43].

6.2. With Gaussian Fluctuations

Let us now recap the general results with the inclusion of Gaussian fluctuations reporting the full set of equations, Equations (167) and (168) with Equations (188)–(190) together with the regularization Equation (194), useful to derive the chemical potential $\mu$ and the critical temperature $T_c$ as functions of the effective interaction parameter, the inverse scattering length,

$$\frac{m}{4\pi a_r}=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{2}{{ϵ}_{\mathbf{k}}}-\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}\left(1-\frac{h^2}{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}}\right)-\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}\left(1+\frac{h^2}{{\xi }_{\mathbf{k}}{\gamma }_{h\mathbf{k}}}\right)\right]$$

(218)

$$\begin{array}{c}\hfill n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left\{1-\frac{1}{2}\left[tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{h\mathbf{k}}}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{h\mathbf{k}}}{2T_c}\right)\right]\right\}\\ \hfill +\frac{1}{\mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }\frac{d\omega }{\pi }{n}_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }\end{array}$$

(219)

with

$${\gamma }_{h\mathbf{k}}=\sqrt{{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2{k}_y^2+{(v_{\mathrm{R}}-v_{\mathrm{D}})}^2{k}_x^2+h^2}$$

(220)

$$\begin{array}{ccc}\hfill \delta (\mathbf{q},\omega )& =& \mathrm{Arg}\left({\mathsf{\Gamma }}^{-1}(\mathbf{q},\omega +i{0}^{+})\right)\hfill \end{array}$$

(221)

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}^{-1}(\mathbf{q},i{\nu }_n)& =& -\frac{m}{4\pi a_r}+\chi (\mathbf{q},i{\nu }_n)+\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{1}{2{ϵ}_{\mathbf{k}}}\hfill \end{array}$$

(222)

$$\begin{array}{c}\hfill \chi (\mathbf{q},i{\nu }_n)=\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left\{\left(1+C_{\mathbf{k},\mathbf{q}}^h\right)\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}}{2T_c}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right.\right.\\ \hfill \left.+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}}{2T_c}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right]\\ \hfill \left.+\left(1-C_{\mathbf{k},\mathbf{q}}^h\right)\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}+\mathbf{q}/2}}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}{2T_c}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{h\mathbf{k}+\mathbf{q}/2}-{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}\right\}\end{array}$$

(223)

where

$$C_{\mathbf{k},\mathbf{q}}^h=\frac{{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2(k_y^2-\frac{q_y^2}{4})+{(v_{\mathrm{R}}-v_{\mathrm{D}})}^2(k_x^2-\frac{q_x^2}{4})-h^2}{{\gamma }_{h\mathbf{k}+\mathbf{q}/2}{\gamma }_{h\mathbf{k}-\mathbf{q}/2}}$$

(224)

We will consider for simplicity the case with $h=0$ but let us first discuss two special cases.

6.2.1. Special Case: Equal Rashba and Dresselhaus

Let us consider the case where $h=0$ and $v_{\mathrm{R}}=v_{\mathrm{D}}$, namley when Rashba and Dresselhauf couplings are equal. The equations are

$$\frac{m}{4\pi a_r}=-\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}-\frac{2}{{ϵ}_{\mathbf{k}}}\right]$$

(225)

$$\begin{array}{c}\hfill n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left\{1-\frac{1}{2}\left[tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}{2T_c}\right)\right]\right\}\\ \hfill +\frac{1}{\mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }\frac{d\omega }{\pi }{n}_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }\end{array}$$

(226)

where ${\gamma }_{\mathbf{k}}={\gamma }_{0\mathbf{k}}$ and since $v_{\mathrm{R}}=v_{\mathrm{D}}$,

$${\gamma }_{\mathbf{k}}=(v_{\mathrm{R}}+v_{\mathrm{D}})|k_y|\equiv v|k_y|$$

(227)

while $\delta (\mathbf{q},\omega )$ is defined in Equation (221). Since $C_{\mathbf{k},\mathbf{q}}^0=1$, we have that $\chi (\mathbf{q},i{\nu }_n)$ is simply given by

$$\chi (\mathbf{q},i{\nu }_n)=\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}}-{\gamma }_{\mathbf{k}+\mathbf{q}}}{2T_c}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}}-{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}+\mathbf{q}}+{\gamma }_{\mathbf{k}}}+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}}+{\gamma }_{\mathbf{k}+\mathbf{q}}}{2T_c}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}}-{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}+\mathbf{q}}-{\gamma }_{\mathbf{k}}}\right]$$

(228)

As already discussed, at the mean field level, after rescaling $k_y\to k_y\pm mv$ we obtain Equations (93) and (94), namely the mean-field equations without SO, but where the chemical potential is $\mu +m{v}^2/2$, so that Equations (225) and (226) can be simplified as it follows

$$\frac{m}{4\pi a_s}=-\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{tanh\left(\frac{({\xi }_{\mathbf{k}}-m{v}^2/2)}{2T_c}\right)}{{\xi }_{\mathbf{k}}-m{v}^2/2}-\frac{1}{{ϵ}_{\mathbf{k}}}\right]$$

(229)

$$n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left[1-tanh\left(\frac{{\xi }_{\mathbf{k}}-m{v}^2/2}{2T_c}\right)\right]+\frac{1}{\pi \mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }d\omega n_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }$$

(230)

6.2.2. Special Case: Only Zeeman Term

For completness, let us consider the case without spin–orbit couplings but in the presence of only a Zeeman term. In this case ${\gamma }_{h\mathbf{k}}=h$. The equations become simply

$$\frac{m}{4\pi a_r}=-\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left\{\frac{1}{{\xi }_{\mathbf{k}}}\left[tanh\left(\frac{{\xi }_{\mathbf{k}}-h}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}}+h}{2T_c}\right)\right]-\frac{2}{{ϵ}_{\mathbf{k}}}\right\}$$

(231)

$$\begin{array}{c}\hfill n=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\left\{1-\frac{1}{2}\left[tanh\left(\frac{{\xi }_{\mathbf{k}}-h}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}}+h}{2T_c}\right)\right]\right\}\\ \hfill +\frac{1}{\mathcal{V}}\sum _{\mathbf{q}}{\int }_{-\infty }^{+\infty }\frac{d\omega }{\pi }{n}_{\mathrm{B}}\left(\omega \right)\frac{\partial \delta (\mathbf{q},\omega )}{\partial \mu }\end{array}$$

(232)

with $\delta (\mathbf{q},\omega )$ still defined by Equation (221) but where $\chi (\mathbf{q},i{\nu }_n)$, since $C_{\mathbf{k},\mathbf{q}}^h=-1$, is simply given by

$$\chi (\mathbf{q},i{\nu }_n)=\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{tanh\left(\frac{{\xi }_{\mathbf{k}+\mathbf{q}/2}-h}{2T_c}\right)+tanh\left(\frac{{\xi }_{\mathbf{k}-\mathbf{q}/2}+h}{2T_c}\right)}{i{\nu }_n-{\xi }_{\mathbf{k}+\mathbf{q}/2}-{\xi }_{\mathbf{k}-\mathbf{q}/2}}\right]$$

(233)

These are the equations to solve for finding the critical temperature of a polarized Fermi gas [78,79] which should be the same in the presence of a Rabi coupling.

6.2.3. Weak Coupling

In this limit the chemical potential is finite and the critical temperature and the scattering length go to zero, in particular $a_s\to 0^{-}$. As argued in the case without SO coupling, ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$ can be expressed as in Equation (113). As a result, in this limit, the Gaussian correction to the density of particles is proportional to the scattering length and, therefore, is negligible, as shown in Equation (114). Since the chemical potential is almost constant and since also the gap equation is unchanged, in the limit $\frac{1}{k_{\mathrm{F}}{a}_s}\to -\infty$, the critical temperature is expected to be almost the same as in the mean field level, therefore, the Gaussian fluctuations are not so effective at weak coupling, as we will verify by numerical calculations.

6.2.4. Strong Coupling

In the strong coupling limit, $a_r\to 0^{+}$, namely in the BEC limit, with SO coupling we found that $\mu \left(T_c\right)$ diverges with the inverse scattering length $y=\frac{1}{k_{\mathrm{F}}{a}_r}$, Equation (208), getting the same formal result as that without SO coupling, $\mu \to -\frac{1}{2m{a}_r}$. Let us now calculate the critical temperature $T_c$ as a function of the y. In what follows we will consider a generic spin-orbit coupling but without Zeeman term, $h=0$.

In the case without SO term the inverse of the vertex function, ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$, in the strong coupling limit, can be seen as the inverse propagator of free bosons, Equation (126), with an effective mass and a redefined chemical potential, ${\mu }_{\mathrm{B}}=2\mu -E_b$, which undergo the Bose–Einstein condensation below the critical temperature. The calculation in the presence of SO coupling is more difficult to perform due to the absence of spherical symmetry, however we can adopt the method applied in Refs. [25,28,29], doing a second order expansion in the momenta for ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$ in the BEC limit. The two-body physics emerges from ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$, Equation (222), in the so-called vacuum limit, discarding the Fermi distribution, $n_{\mathrm{F}}\to 0$, and putting the chemical potential to zero [25,28,29]. What remains is the inverse of the two-body scattering matrix, $T_{2B}^{-1}(\mathbf{q},i{\nu }_n)\equiv {\mathsf{\Gamma }}^{-1}(\mathbf{q},i{\nu }_n){|}_{n_{\mathrm{F}},\mu =0}$, which, reads

$$\begin{array}{c}\hfill T_{2B}^{-1}(\mathbf{q},i{\nu }_n)=-\frac{m}{4\pi a_r}+\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{2}{{ϵ}_{\mathbf{k}}}+\frac{\left(1+C_{\mathbf{k},\mathbf{q}}^0\right)}{i{\nu }_n-{ϵ}_{\mathbf{k}-\mathbf{q}/2}-{ϵ}_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{\mathbf{k}+\mathbf{q}/2}}\right.\\ \hfill +\frac{\left(1+C_{\mathbf{k},\mathbf{q}}^0\right)}{i{\nu }_n-{ϵ}_{\mathbf{k}-\mathbf{q}/2}-{ϵ}_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{\mathbf{k}+\mathbf{q}/2}}\\ \hfill +\left.\frac{2\left(1-C_{\mathbf{k},\mathbf{q}}^0\right)}{i{\nu }_n-{ϵ}_{\mathbf{k}-\mathbf{q}/2}-{ϵ}_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{\mathbf{k}-\mathbf{q}/2}}\right]\end{array}$$

(234)

• Binding Energy

From the two-body scattering matrix we can evaluate the scattering energies and the existence of a two-particle bound state. We saw in a previous section, Section 2.1, how $T_{2B}\left(z\right)$ has simple poles at the bound states and a branch cut at the scattering states. In order to prove the existence of a bound state in the presence of a spin–orbit term it is enough, therefore, to impose that, at vanishing momentum,

$$T_{2B}^{-1}\left(\mathbf{q}=0,z={\mathcal{E}}_b\right)=0$$

(235)

From Equation (224) we have $C_{\mathbf{k},0}^0=1$, so that Equation (235) implies the following equation

$$-\frac{m}{4\pi a_r}+\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\left\{\frac{1}{{\mathcal{E}}_b-2{ϵ}_{\mathbf{k}}-2{\gamma }_{\mathbf{k}}}+\frac{1}{{\mathcal{E}}_b-2{ϵ}_{\mathbf{k}}+2{\gamma }_{\mathbf{k}}}+\frac{1}{{ϵ}_{\mathbf{k}}}\right\}=0$$

(236)

that should be solved in terms of the bound state energy, measured from the threshold energy for free fermions, as a function of the coupling y. This can be done by simple numerical techniques. In the continuum limit, introducing the following dimensionless quantities $\mathbf{p}=\frac{\mathbf{k}}{k_{\mathrm{F}}}$, ${\tilde{\mathcal{E}}}_b=\frac{{\mathcal{E}}_b}{{ϵ}_{\mathrm{F}}}$, ${\tilde{T}}_c=\frac{T_c}{{ϵ}_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{R}}=\frac{v_{\mathrm{R}}}{v_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{D}}=\frac{v_{\mathrm{D}}}{v_{\mathrm{F}}}$ and ${\tilde{\gamma }}_{\mathbf{p}}({\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})=2{\gamma }_{\mathbf{k}}(v_{\mathrm{R}},v_{\mathrm{D}})$, we get

$$y+\frac{2}{{\pi }^2}{\int }_0^{\infty }d{p}_x{\int }_0^{\infty }d{p}_y{\int }_0^{\infty }d{p}_z\left\{\frac{1}{p^2+{\tilde{\gamma }}_{\mathbf{p}}-{\tilde{\mathcal{E}}}_b/2}+\frac{1}{p^2-{\tilde{\gamma }}_{\mathbf{p}}-{\tilde{\mathcal{E}}}_b/2}-\frac{2}{p^2}\right\}=0$$

(237)

which has the same form of Equation (201). The chemical potential and the bound state energy, $E_b={ϵ}_{\mathrm{F}}{\tilde{E}}_b$, measured at the threshold energy $m{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2$, in terms of dimensionless quantities,

$${\tilde{\mathcal{E}}}_b={\tilde{E}}_b-2{({\tilde{v}}_{\mathrm{R}}+{\tilde{v}}_{\mathrm{D}})}^2$$

(238)

satisfy the same equation, therefore, in the the strong coupling limit, they are connected by the relation $\mu \underset{y\to \infty }{⟶}\frac{{\mathcal{E}}_b}{2}=\frac{E_b}{2}+\frac{m{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2}{2}$. In the deep BEC, neglecting the second term, we have

$$\mu \underset{y\to \infty }{⟶}{E}_b/2$$

(239)

which is the same result obtained in the absence of SO coupling, Equation (123). Taking advantage of Equation (207), in the case of only Rashba coupling ($v_{\mathrm{D}}=0$), we can write an analytical solution

$$y=\sqrt{{\tilde{v}}_{\mathrm{R}}^2-\frac{{\tilde{E}}_b}{2}}-{\tilde{v}}_{\mathrm{R}}arctanh\left({\tilde{v}}_{\mathrm{R}}\sqrt{\frac{2}{2{\tilde{v}}_{\mathrm{R}}^2-{\tilde{E}}_b}}\right)$$

(240)

Inverting Equation (240) we can get $E_b$ as a function of y. For equal Rashba and Dresselhaus couplings, $v_{\mathrm{R}}=v_{\mathrm{D}}$, we have, instead,

$$y=\sqrt{-\frac{{\tilde{E}}_b}{2}}$$

(241)

which is the same result as that without SO coupling. For the more general case one can easily resort to a numerical computation to solve Equation (237), or, after analytically integrating over $p_z$, Equation (203) with ${\tilde{\mathcal{E}}}_b/2$ instead of $\tilde{\mu }$. The results are shown in Figure 6.

For pure Rashba coupling, Figure 6, a bound state exists also in the BCS part of the crossover, where the scattering length is negative. Actually by looking at Equation (240), for any value of y, a negative value of $E_b$ exists. Moreover, in Figure 6 we show that turning on a Dresselhaus term ($v_{\mathrm{D}}\ne 0$), greatly affects the existence of the bound state in the BCS side.

The very important feature is, therefore, that a bound state exists also for negative scattering length, namely in the BCS side. The effect is relevant for a pure Rashba coupling while it is suppressed in the presence of mixed Rashba and Dresselhaus couplings. For equal Rashba and Dresselhaus couplings the bound state exists only in the BEC side, exactly as in the case without SO.

• Effective Masses

At strong coupling two fermions can be treated as a single bosonic particle with a double mass. However in the presence of a spin–orbit interaction we have to reconsider the evaluation of the effective mass. This can be done at low momenta requiring that the dispersion relation of free bosons in the strong coupling limit ($1/k_{\mathrm{F}}{a}_s\to +\infty$) can be written as

$${ϵ}_b\left(\mathbf{q}\right)={\mathcal{E}}_b+\frac{q_x^2}{2m_x}+\frac{q_y^2}{2m_y}+\frac{q_z^2}{2m_z}$$

(242)

where ${\mathcal{E}}_b$, $m_x$ and $m_y$ depend on both the strength of the interaction and the spin–orbit couplings. The effective masses along x and y directions can be generally different in the presence of mixed Rashba and Dresselhaus terms. These masses can be derived imposing that the inverse of the scattering matrix vanishes when measured right at ${ϵ}_b\left(\mathbf{q}\right)$ as in Equation (242)

$$T_{2B}^{-1}\left(\mathbf{q},z={ϵ}_b\left(\mathbf{q}\right)\right)=0$$

(243)

At low energy and momentum this condition is imposed to the second order expansion

$$T_{2B}^{-1}\left(\mathbf{q},{ϵ}_b\left(\mathbf{q}\right)\right)\approx T_{2B}^{-1}\left(0,{\mathcal{E}}_b\right)+{\mathcal{F}}_x({\mathcal{E}}_b,m_x)q_x^2+{\mathcal{F}}_y({\mathcal{E}}_b,m_y)q_y^2+{\mathcal{F}}_z({\mathcal{E}}_b,m_z)q_z^2=0$$

(244)

which implies that

$$T_{2B}^{-1}\left(0,{\mathcal{E}}_b\right)=0,{\mathcal{F}}_x({\mathcal{E}}_b,m_x)=0,{\mathcal{F}}_y({\mathcal{E}}_b,m_y)=0,{\mathcal{F}}_z({\mathcal{E}}_b,m_z)=0$$

(245)

The first equation in Equation (245) is the same encountered before for the bound state energy, Equation (235). By choosing $m_z=2m$ we get immediately that ${\mathcal{F}}_z({\mathcal{E}}_b,2m)=0$ since, in this way, $T_{2B}^{-1}\left(\mathbf{q},{ϵ}_b\left(\mathbf{q}\right)\right)$ loses the dependence on $q_z$. Actually, the two-body scattering matrix reads

$$\begin{array}{ccc}\hfill T_{2B}^{-1}\left(\mathbf{q},{ϵ}_b\left(\mathbf{q}\right)\right)& =& -\frac{m}{4\pi a_r}+\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\frac{2}{{ϵ}_{\mathbf{k}}}\hfill \\ \hfill & & +\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\frac{\left(1+C_{\mathbf{k},\mathbf{q}}^0\right)}{{\mathcal{E}}_b+\frac{q_x^2}{2m_x}+\frac{q_y^2}{2m_y}+\frac{q_z^2}{4m}-{ϵ}_{\mathbf{k}-\mathbf{q}/2}-{ϵ}_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{\mathbf{k}+\mathbf{q}/2}}\hfill \\ & & +\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\frac{\left(1+C_{\mathbf{k},\mathbf{q}}^0\right)}{{\mathcal{E}}_b+\frac{q_x^2}{2m_x}+\frac{q_y^2}{2m_y}+\frac{q_z^2}{4m}-{ϵ}_{\mathbf{k}-\mathbf{q}/2}-{ϵ}_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{\mathbf{k}+\mathbf{q}/2}}\hfill \\ \hfill & & +\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\frac{\left(1-C_{\mathbf{k},\mathbf{q}}^0\right)}{{\mathcal{E}}_b+\frac{q_x^2}{2m_x}+\frac{q_y^2}{2m_y}+\frac{q_z^2}{4m}-{ϵ}_{\mathbf{k}-\mathbf{q}/2}-{ϵ}_{\mathbf{k}+\mathbf{q}/2}+{\gamma }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{\mathbf{k}-\mathbf{q}/2}}\hfill \end{array}$$

(246)

Introducing the following definitions

$$X=\frac{2m}{m_x}-1,Y=\frac{2m}{m_y}-1$$

(247)

the two-body scattering matrix can be written as it follows

$$\begin{array}{ccc}\hfill T_{2B}^{-1}\left(\mathbf{q},{ϵ}_b\left(\mathbf{q}\right)\right)& =& -\frac{m}{4\pi a_r}+\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{1}{2{ϵ}_{\mathbf{k}}}\hfill \\ \hfill & & +\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\frac{\left(1+C_{\mathbf{k},\mathbf{q}}^0\right)}{{\mathcal{E}}_b+\frac{X}{4m}{q}_x^2+\frac{Y}{4m}{q}_y^2-2{ϵ}_{\mathbf{k}}-{\gamma }_{\mathbf{k}-\mathbf{q}/2}-{\gamma }_{\mathbf{k}+\mathbf{q}/2}}\hfill \\ & & +\frac{1}{4\mathcal{V}}\sum _{\mathbf{k}}\frac{\left(1+C_{\mathbf{k},\mathbf{q}}^0\right)}{{\mathcal{E}}_b+\frac{X}{4m}{q}_x^2+\frac{Y}{4m}{q}_y^2-2{ϵ}_{\mathbf{k}}+{\gamma }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{\mathbf{k}+\mathbf{q}/2}}\hfill \\ \hfill & & +\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\frac{\left(1-C_{\mathbf{k},\mathbf{q}}^0\right)}{{\mathcal{E}}_b+\frac{X}{4m}{q}_x^2+\frac{Y}{4m}{q}_y^2-2{ϵ}_{\mathbf{k}}+{\gamma }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{\mathbf{k}-\mathbf{q}/2}}\hfill \end{array}$$

(248)

Expanding for small $\mathbf{q}$ the following quantities

$$\begin{array}{ccc}& & C_{\mathbf{k},\mathbf{q}}^0\approx 1-\frac{{(v_{\mathrm{R}}^2-v_{\mathrm{D}}^2)}^2}{2{\gamma }_{\mathbf{k}}^4}{\left(k_x{q}_y-k_y{q}_x\right)}^2\hfill \end{array}$$

(249)

$$\begin{array}{ccc}& & {\gamma }_{\mathbf{k}-\mathbf{q}/2}+{\gamma }_{\mathbf{k}+\mathbf{q}/2}\approx 2{\gamma }_{\mathbf{k}}+\frac{{(v_{\mathrm{R}}^2-v_{\mathrm{D}}^2)}^2}{4{\gamma }_{\mathbf{k}}^3}{\left(k_x{q}_y-k_y{q}_x\right)}^2\hfill \end{array}$$

(250)

$$\begin{array}{ccc}& & {\gamma }_{\mathbf{k}+\mathbf{q}/2}-{\gamma }_{\mathbf{k}-\mathbf{q}/2}\approx \frac{{(v_{\mathrm{R}}+v_{\mathrm{D}})}^2}{{\gamma }_{\mathbf{k}}}{k}_x{q}_x+\frac{{(v_{\mathrm{R}}-v_{\mathrm{D}})}^2}{{\gamma }_{\mathbf{k}}}{k}_y{q}_y\hfill \end{array}$$

(251)

the two-body scattering matrix, at leading orders, reads

$$\begin{array}{ccc}\hfill & & T_{2B}^{-1}\left(\mathbf{q},{ϵ}_b\left(\mathbf{q}\right)\right)\approx -\frac{m}{4\pi a_r}+\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\left[\frac{1}{{\mathcal{E}}_b-2{ϵ}_{\mathbf{k}}-2{\gamma }_{\mathbf{k}}}+\frac{1}{{\mathcal{E}}_b-2{ϵ}_{\mathbf{k}}+2{\gamma }_{\mathbf{k}}}+\frac{1}{{ϵ}_{\mathbf{k}}}\right]\hfill \\ \hfill & & -\frac{1}{4\mathcal{V}}\left[\frac{X}{m}\sum _{\mathbf{k}}\frac{{({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2+4{\gamma }_{\mathbf{k}}^2}{{\left({({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2-4{\gamma }_{\mathbf{k}}^2\right)}^2}+{\left(v_{\mathrm{R}}^2-v_{\mathrm{D}}^2\right)}^2\sum _{\mathbf{k}}\frac{16k_y^2}{{\left({({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2-4{\gamma }_{\mathbf{k}}^2\right)}^2({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}\right]q_x^2\hfill \\ & & -\frac{1}{4\mathcal{V}}\left[\frac{Y}{m}\sum _{\mathbf{k}}\frac{{({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2+4{\gamma }_{\mathbf{k}}^2}{{\left({({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2-4{\gamma }_{\mathbf{k}}^2\right)}^2}+{\left(v_{\mathrm{R}}^2-v_{\mathrm{D}}^2\right)}^2\sum _{\mathbf{k}}\frac{16k_x^2}{{\left({({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2-4{\gamma }_{\mathbf{k}}^2\right)}^2({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}\right]q_y^2\hfill \\ \hfill & & \equiv T_{2B}^{-1}\left(0,{\mathcal{E}}_b\right)+{\mathcal{F}}_x({\mathcal{E}}_b,m_x)\phantom{|}{q}_x^2+{\mathcal{F}}_y({\mathcal{E}}_b,m_y)\phantom{|}{q}_y^2\hfill \end{array}$$

(252)

The first term is the scattering matrix at $\mathbf{q}=0$, therefore, $T_{2B}^{-1}\left(0,{\mathcal{E}}_b\right)=0$, by definition of the shifted bound state energy ${\mathcal{E}}_b$. Solving ${\mathcal{F}}_x({\mathcal{E}}_b,m_x)=0$ in terms of X and ${\mathcal{F}}_y({\mathcal{E}}_b,m_y)=0$ in terms of Y, and from Equation (247), we can derive the values of the effective masses

$$\begin{array}{c}\hfill \frac{2m}{m_x}=1-m{\left(v_{\mathrm{R}}^2-v_{\mathrm{D}}^2\right)}^2\frac{{\mathcal{B}}_x}{\mathcal{A}}\end{array}$$

(253)

$$\begin{array}{c}\hfill \frac{2m}{m_y}=1-m{\left(v_{\mathrm{R}}^2-v_{\mathrm{D}}^2\right)}^2\frac{{\mathcal{B}}_y}{\mathcal{A}}\end{array}$$

(254)

where

$$\begin{array}{ccc}\hfill & & \mathcal{A}=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{{({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2+4{\gamma }_{\mathbf{k}}^2}{{\left({({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2-4{\gamma }_{\mathbf{k}}^2\right)}^2}\hfill \end{array}$$

(255)

$$\begin{array}{ccc}\hfill & & {\mathcal{B}}_x=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{16k_y^2}{{\left({({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2-4{\gamma }_{\mathbf{k}}^2\right)}^2({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}\hfill \end{array}$$

(256)

$$\begin{array}{ccc}& & {\mathcal{B}}_y=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}}\frac{16k_x^2}{{\left({({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}^2-4{\gamma }_{\mathbf{k}}^2\right)}^2({\mathcal{E}}_b-2{ϵ}_{\mathbf{k}})}\hfill \end{array}$$

(257)

In the specific case of equal Rashba and Dresselhaus couplings ($v_{\mathrm{R}}=\pm v_{\mathrm{D}}$), as already seen, the bound state energy is not distinguishable from the one in absence of SO coupling, which means that a bound state is present only for $y\ge 0$, namely only in the BEC side, where the effective masses are simply $2m/m_x=2m/m_y=1$.

In the continuum limit, rescaling the variables $\mathbf{p}=\frac{\mathbf{k}}{k_{\mathrm{F}}}$, $\tilde{\mu }=\frac{\mu }{{ϵ}_{\mathrm{F}}}$, ${\tilde{\mathcal{E}}}_b=\frac{{\mathcal{E}}_b}{{ϵ}_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{R}}=\frac{v_{\mathrm{R}}}{v_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{D}}=\frac{v_{\mathrm{D}}}{v_{\mathrm{F}}}$, ${\tilde{\gamma }}_{\mathbf{p}}({\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})=2{\gamma }_{\mathbf{k}}(v_{\mathrm{R}},v_{\mathrm{D}})$, Equations (255)–(257) in terms of the dimensionless parameters become

$$\begin{array}{ccc}\hfill & & \mathcal{A}=\frac{k_{\mathrm{F}}^3}{{ϵ}_{\mathrm{F}}^2}\int \frac{d\mathbf{p}}{{\left(2\pi \right)}^3}\frac{{({\tilde{\mathcal{E}}}_b-2p^2)}^2+4{\tilde{\gamma }}_{\mathbf{p}}^2}{{\left({({\tilde{\mathcal{E}}}_b-2p^2)}^2-4{\tilde{\gamma }}_{\mathbf{p}}^2\right)}^2}=\frac{k_{\mathrm{F}}^3}{{ϵ}_{\mathrm{F}}^2}\tilde{\mathcal{A}}\hfill \end{array}$$

(258)

$$\begin{array}{ccc}\hfill & & {\mathcal{B}}_x=\frac{k_{\mathrm{F}}^5}{{ϵ}_{\mathrm{F}}^5}\int \frac{d\mathbf{p}}{{\left(2\pi \right)}^3}\frac{16p_y^2}{{\left({({\tilde{\mathcal{E}}}_b-2p^2)}^2-4{\tilde{\gamma }}_{\mathbf{p}}^2\right)}^2({\tilde{\mathcal{E}}}_b-2p^2)}=\frac{k_{\mathrm{F}}^5}{{ϵ}_{\mathrm{F}}^5}{\tilde{\mathcal{B}}}_x\hfill \end{array}$$

(259)

$$\begin{array}{ccc}& & {\mathcal{B}}_y=\frac{k_{\mathrm{F}}^5}{{ϵ}_{\mathrm{F}}^5}\int \frac{d\mathbf{p}}{{\left(2\pi \right)}^3}\frac{16p_x^2}{{\left({({\tilde{\mathcal{E}}}_b-2p^2)}^2-4{\tilde{\gamma }}_{\mathbf{p}}^2\right)}^2({\tilde{\mathcal{E}}}_b-2p^2)}=\frac{k_{\mathrm{F}}^5}{{ϵ}_{\mathrm{F}}^5}{\tilde{\mathcal{B}}}_y\hfill \end{array}$$

(260)

The effective masses can be written in terms of dimensionless quantities

$$\begin{array}{c}\hfill \frac{2m}{m_x}=1-8{\left({\tilde{v}}_{\mathrm{R}}^2-{\tilde{v}}_{\mathrm{D}}^2\right)}^2\frac{{\tilde{\mathcal{B}}}_x}{\tilde{\mathcal{A}}}\end{array}$$

(261)

$$\begin{array}{c}\hfill \frac{2m}{m_y}=1-8{\left({\tilde{v}}_{\mathrm{R}}^2-{\tilde{v}}_{\mathrm{D}}^2\right)}^2\frac{{\tilde{\mathcal{B}}}_y}{\tilde{\mathcal{A}}}\end{array}$$

(262)

For the case with a pure Rashba coupling, for which $m_x=m_y$, we can perform the integrals $\tilde{\mathcal{A}}$, ${\tilde{\mathcal{B}}}_x$ and ${\tilde{\mathcal{B}}}_y$ analytically, obtaining

$$\begin{array}{c}\tilde{\mathcal{A}}=\frac{1}{2{\pi }^2}{\int }_0^{\infty }d{p}_z{\int }_0^{\infty }d{p}_{\perp }{p}_{\perp }\frac{{({\tilde{\mathcal{E}}}_b-2p_{\perp }^2-2p_z^2)}^2+4{\left(2{\tilde{v}}_{\mathrm{R}}{p}_{\perp }\right)}^2}{{\left[{({\tilde{\mathcal{E}}}_b-2p_{\perp }^2-2p_z^2)}^2-4{\left(2{\tilde{v}}_{\mathrm{R}}{p}_{\perp }\right)}^2\right]}^2}=\frac{-\sqrt{-{\tilde{\mathcal{E}}}_b}}{16\sqrt{2}\pi \left({\tilde{\mathcal{E}}}_b+2{\tilde{v}}_{\mathrm{R}}^2\right)}\end{array}$$

(263)

$$\begin{array}{c}{\tilde{\mathcal{B}}}_x={\tilde{\mathcal{B}}}_y=\frac{1}{2{\pi }^2}{\int }_0^{\infty }d{p}_z{\int }_0^{\infty }d{p}_{\perp }\frac{8p_{\perp }^3}{{\left[{({\tilde{\mathcal{E}}}_b-2p_{\perp }^2-2p_z^2)}^2-4{\left(2{\tilde{v}}_{\mathrm{R}}{p}_{\perp }\right)}^2\right]}^2({\tilde{\mathcal{E}}}_b-2p_{\perp }^2-2p_z^2)}\\ \hfill =\frac{1}{256\sqrt{2}\pi {\tilde{v}}_{\mathrm{R}}^4}\frac{1}{\sqrt{-{\tilde{\mathcal{E}}}_b}}\left[log\left(1+\frac{2{\tilde{v}}_{\mathrm{R}}^2}{{\tilde{\mathcal{E}}}_b}\right)-\frac{2{\tilde{v}}_{\mathrm{R}}^2}{{\tilde{\mathcal{E}}}_b+2{\tilde{v}}_{\mathrm{R}}^2}\right]\end{array}$$

(264)

where, from Equation (238), in the presence of a bound state, we assume that

$${\tilde{\mathcal{E}}}_b\le -2{\tilde{v}}_{\mathrm{R}}^2$$

(265)

The final analytical expressions for the effective masses are, therefore,

$$\frac{2m}{m_x}=\frac{2m}{m_y}=1-\frac{1}{2{\tilde{\mathcal{E}}}_b}\left[\left({\tilde{\mathcal{E}}}_b+2{\tilde{v}}_{\mathrm{R}}^2\right)log\left(1+\frac{2{\tilde{v}}_{\mathrm{R}}^2}{{\tilde{\mathcal{E}}}_b}\right)-2{\tilde{v}}_{\mathrm{R}}^2\right]$$

(266)

which, together with Equation (240) written in terms of ${\tilde{\mathcal{E}}}_b$

$$y=\frac{\sqrt{-{\tilde{\mathcal{E}}}_b}}{\sqrt{2}}-{\tilde{v}}_{\mathrm{R}}arctanh\left({\tilde{v}}_{\mathrm{R}}\frac{\sqrt{2}}{\sqrt{-{\tilde{\mathcal{E}}}_b}}\right)$$

(267)

in order to get their behaviors along the crossover, see Figure 7. This result is in agreement with the T-matrix approach [50]. At the special point $y=\frac{1}{k_{\mathrm{F}}{a}_r}=0$, also called unitary limit, the ratio $2{\tilde{v}}_{\mathrm{R}}^2/{\tilde{\mathcal{E}}}_b\equiv -c_o$ is the solution of the following equation

$$\sqrt{c_o}=tanh(1/\sqrt{c_o})$$

(268)

which is $c_o\approx 0.6948$. As a result $\frac{2m}{m_x}=\frac{2m}{m_y}$ at $y=0$ is the same for any value of ${\tilde{v}}_{\mathrm{R}}$,

$${\left.\frac{2m}{m_x}\right|}_{a_r^{-1}=0}=1-\frac{1}{2}\left[c_o+\left(1-c_o\right)log\left(1-c_o\right)\right]\approx 0.8387$$

(269)

which is the nodal point shown in Figure 7. For the general case of mixed Rashba and Dresselhaus terms the effective masses are given by Equations (261) and (262), together with Equations (258)–(260).

6.2.5. Full Crossover

We can now calculate the critical temperature $T_c$ along the crossover within the bosonic approximation, discussing a possible further improvement useful to refine the calculation in the intermediate regime.

• Bosonic Approximation

We found that the inverse of the scattering matrix, expanding at low momenta around $\mathbf{q}=0$, has the form

$$T_{2B}^{-1}(\mathbf{q},i\nu )\approx i\nu -{ϵ}_b\left(\mathbf{q}\right)=i\nu -{\mathcal{E}}_b-\frac{q_x^2}{2m_x}-\frac{q_y^2}{2m_y}-\frac{q_z^2}{4m}$$

(270)

with ${\mathcal{E}}_b$, $m_x$ and $m_y$ that depend on $y=\frac{1}{k_{\mathrm{F}}{a}_r}$. This result for the two-body problem can be exploited to determine the properties of the vertex function in the strong coupling limit

$${\left.{\mathsf{\Gamma }}^{-1}(\mathbf{q},i\nu )\right|}_{\mu \to -\infty }=T_{2B}^{-1}(\mathbf{q},i\nu +2\mu )$$

(271)

getting, in this regime,

$${\mathsf{\Gamma }}^{-1}(\mathbf{q},i\nu )\approx i\nu +2\mu -{ϵ}_b\left(\mathbf{q}\right)=i\nu -\left(\frac{q_x^2}{2m_x}+\frac{q_y^2}{2m_y}+\frac{q_z^2}{4m}-{\mu }_{\mathrm{B}}\right)$$

(272)

where ${\mu }_{\mathrm{B}}=2\mu -{\mathcal{E}}_b/2\simeq 2\mu -E_b/2$. We find, therefore, that the partition function approaches that of a set of free bosons, with anisotropic masses. It is straightforward to show that the condensation temperature for this kind of systems is actually rescaled by the effective masses

$$T_c=\frac{2\pi }{m_b}{\left(\frac{n_{\mathrm{B}}}{2\zeta (3/2)}\right)}^{\frac{2}{3}}{\left(\frac{m_{\mathrm{B}}}{m_x}\right)}^{\frac{1}{3}}{\left(\frac{m_{\mathrm{B}}}{m_y}\right)}^{\frac{1}{3}}$$

(273)

where $m_{\mathrm{B}}=2m$ and, in the deep BEC, $n_{\mathrm{B}}=\frac{n}{2}=\frac{1}{2}\frac{{\left(2m{ϵ}_{\mathrm{F}}\right)}^{3/2}}{3{\pi }^2}$, therefore

$$\frac{T_c}{{ϵ}_{\mathrm{F}}}\approx 0.218{\left(\frac{2m}{m_x}\right)}^{\frac{1}{3}}{\left(\frac{2m}{m_y}\right)}^{\frac{1}{3}}.$$

(274)

We find that in the BEC limit the effective masses approach those of two fermions,

$$\frac{2m}{m_{\alpha }}\underset{BEC}{\stackrel{}{\to }}1$$

(275)

with $\alpha =x,y$, for any mixture of Rashba and Dresselhaus couplings. For the case of only Rashba (or Dresselhaus) coupling Equation (275) can be read out from Equation (266) in the limit of ${\tilde{\mathcal{E}}}_b\to \infty$.

For two fermions with same momenta in a spin singlet state, composing a boson, the spin–orbit coupling cancels out and a free particle spectrum with mass $m_{\mathrm{B}}=2$ m emerges. From Equation (273), solving in terms of $n_{\mathrm{B}}$, we have

$$n_{\mathrm{B}}=\zeta (3/2){\left(\frac{m{T}_c}{\pi }\right)}^{3/2}{\left(\frac{m_x}{2m}\right)}^{\frac{1}{2}}{\left(\frac{m_y}{2m}\right)}^{\frac{1}{2}}$$

(276)

As done in Equation (127), let us suppose now that the number of bosonic excitations composed by pairs of fermions is $n_{\mathrm{B}}\approx 2n^{\left(2\right)}$, while the rest of fermions remains unpaired so that we have

$$n^{\left(2\right)}\approx 2\zeta (3/2){\left(\frac{m{T}_c}{\pi }\right)}^{3/2}{\left(\frac{m_x}{2m}\right)}^{\frac{1}{2}}{\left(\frac{m_y}{2m}\right)}^{\frac{1}{2}}.$$

(277)

With this approximation the set of equations to solve is the following

$$\frac{m}{4\pi a_r}=-\frac{1}{4\mathsf{\Omega }}\sum _{\mathbf{k}}\left\{\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}+\frac{tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}{2T_c}\right)}{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}-\frac{2}{{ϵ}_{\mathbf{k}}}\right\}$$

(278)

$$\begin{array}{c}\hfill n=\frac{1}{2\mathcal{V}}\sum _{\mathbf{k}}\left\{2-tanh\left(\frac{{\xi }_{\mathbf{k}}-{\gamma }_{\mathbf{k}}}{2T_c}\right)-tanh\left(\frac{{\xi }_{\mathbf{k}}+{\gamma }_{\mathbf{k}}}{2T_c}\right)\right\}\\ \hfill +2\zeta (3/2){\left(\frac{m{T}_c}{\pi }\right)}^{3/2}{\left(\frac{m_x}{2m}\frac{m_y}{2m}\right)}^{1/2}\end{array}$$

(279)

In the continuum limit, from Equations (209) and (210), at fixed density $n=\frac{k_{\mathrm{F}}^3}{3{\pi }^2}$, after rescaling $\mathbf{p}=\frac{1}{k_{\mathrm{F}}}\mathbf{k}$, $\tilde{\mu }=\frac{\mu }{{ϵ}_{\mathrm{F}}}$, ${\tilde{T}}_c=\frac{T_c}{{ϵ}_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{R}}=\frac{v_{\mathrm{R}}}{v_{\mathrm{F}}}$, ${\tilde{v}}_{\mathrm{D}}=\frac{v_{\mathrm{D}}}{v_{\mathrm{F}}}$, ${\tilde{\gamma }}_{\mathbf{p}}({\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})=2{\gamma }_{\mathbf{k}}(v_{\mathrm{R}},v_{\mathrm{D}})$, the above equations can be written in the following form

$$\begin{array}{ccc}& & y=I_y^{\mathrm{SO}}(\tilde{\mu },{\tilde{T}}_c,{\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})\hfill \end{array}$$

(280)

$$\begin{array}{ccc}& & 1=I_n^{\mathrm{SO}}(\tilde{\mu },{\tilde{T}}_c,{\tilde{v}}_{\mathrm{R}},{\tilde{v}}_{\mathrm{D}})+3\zeta (3/2)\sqrt{\frac{\pi }{2}}{\tilde{T}}_c^{3/2}{\left(\frac{m_x}{2m}\frac{m_y}{2m}\right)}^{1/2}\hfill \end{array}$$

(281)

This set of equations can be easily solved numerically. Some results are reported in Figure 8 for the most relevant case with pure Rashba coupling.

These results have to be contrasted with those reported in Figure 2 in the absence of spin–orbit couplings.

• Further Improvements

The existence of a bound state also in the BCS side of the crossover allowed us to extend the results of the effective masses also in this region, although the pole structure of the vertex function $\mathsf{\Gamma }\left(\mathsf{q}\right)$, presented in Equation (272), has been derived in the strong coupling limit. To refine this result it should be sufficient to expand at second order in the momentum $\mathbf{q}$ the full ${\mathsf{\Gamma }}^{-1}\left(\mathsf{q}\right)$ showed in Equations (222) and (224), with $h=0$. The only difference with respect to what done so far is to include the expansion of the hyperbolic tangents of some functions, $f_{\mathbf{q}}=({\xi }_{\mathbf{k}\pm \mathbf{q}/2}\pm {\gamma }_{\mathbf{k}\pm \mathbf{q}/2})/2$, appearing in Equation (224), up to second order in q

$$\begin{array}{ccc}& & tanh\left(f_{\mathbf{q}}/T_c\right)=tanh\left(f_0/T_c\right)+\mathrm{sech}{\left(f_0/T_c\right)}^2\nabla f_{\mathbf{q}}{|}_0·\mathbf{q}/T_c+\hfill \\ \hfill & & +\sum _{ij}\left[\mathrm{sech}{\left(f_0/T_c\right)}^2{\partial }_i{\partial }_j{f}_{\mathbf{q}}{|}_0/2-\mathrm{sech}{\left(f_0/T_c\right)}^{2}tanh\left(f_0/T_c\right){\partial }_i{f}_{\mathbf{q}}{\partial }_j{f}_{\mathbf{q}}{|}_0\right]q_i{q}_j/T_c^2+O\left(q^3\right)\hfill \end{array}$$

(282)

The inclusion of this expansion, with respect to the vacuum limit considered so far, can be seen somehow as the addition of thermal corrections to the zero-temperature results. These corrections are relevant in the intermediate regime, since in the strong coupling limit $f_{\mathbf{q}}/T_c\to \infty$, therefore $tanh(f_{\mathbf{q}}/T_c)\to 1$, recovering the results of the bosonic approximation, while in the weak coupling limit the quantum fluctuations themselves are not relevant at all, as discussed previously.

Alternatively one can calculate the exact Gaussian fluctuations solving numerically the full equations in Equations (218) and (219) together with Equations (221), (222) and (224), as done in the absence of the spin–orbit couplings. However, in the presence of SO the rotational symmetry is broken and the numerics is extremely much more involved and computationally costly.

7. Conclusions

We reviewed the calculations for the critical temperature along the BCS-BEC crossover with and without spin–orbit couplings in the path integral formalism. The final equations, with the inclusion of quantum fluctuations at the Gaussian level have been written in the most general case with Rashba, Dresselhaus and Zeeman terms, even though most of the results presented are done for the most relevant case with only Rashba term. We found that in the weak coupling limit $T_c$ is strongly enhanced, consistently with the increase of the gap [32] which implies an increased binding energy promoting the formation of fermionic pairs, therefore, raising the critical temperature. In the intermediate regime the rapid increase of $T_c$ predicted by the mean-field theory is softened by the evolution of the effective masses, originated by the quantum Gaussian fluctuations, of the composite bosons satisfying $m_x,m_y>2$ m which lower the critical temperature. Intriguingly, either the binding energy and the effective masses have finite values also in the BCS regime. For the Rashba case this finding can be found analytically, in particular, at the unitary limit we found that the masses do not depend on the spin–orbit parameter. Finally, in the deep BEC, we found that $T_c$ goes to the same value as that without SO, consistently with the strong singlet-pair formation.

Before we conclude a last remark is in order. One can expect that the presence of SO coupling could increase the pseudogap regime, a regime of uncondesed fermionic pairs (see for instance [80,81,82] for further details). This consideration is based on the fact that the mean-field pairing temperature is much higher than the critical superfluid temperature $T_c$ because the bosonic excitations are quite relevant already in the BCS regime. As a result the pseudogap regime, supposed to be hosted in between those two temperatures, might be much more extended in the presence of a spin–orbit coupling.