Critical temperature in the BCS-BEC crossover with spin-orbit coupling

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.


I. 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 arXiv:2104.00570v1 [cond-mat.quant-gas] 1 Apr 2021 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. The discovery of the high-temperature superconductors 5 , 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 6 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. λ 2 th = 2π 2 mk B T ∼ n −2/3 . The degenerate temperatures, for the densities mentioned before, are smaller than a µK. These are the conditions achieved in the first experimental realization of a Bose-Einstein condensation 7 . 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 8,9 in bosonic systems. Soon after these techniques have been extended to create and manipulate ultracold gases for fermionic systems [10][11][12][13] 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 14,15 , studying also the thermodynamical properties 16,17 . 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 [18][19][20][21] .
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 spinorbit effects with the so-called Rashba 22 and Dresselhaus 23 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.  The aim of this paper is to review some basic concepts on the BCS-BEC crossover, 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 in the intermediate region of the crossover and in the strong coupling regime.

A. 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 an the effective charge Ze of the nucleus. This central force is known to lead to a Bohr energy spectrum with n = 1, 2, ..., where α = e 2 4π 0 c ≈ 1 137 is the so-called fine structure constant and m r ≈ 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 operatorsL 2 andL 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,Ĥ so = αso 2L ·Ŝ, therefore, only the total angular momentumĴ =L+Ŝ is conserved. This additional term causes the known correction to the spectrum, called fine structure, shifting the degeneracy in the energy spectrum by ∆E J = αso 2 [J(J + 1) − L(L + 1) − S(S + 1)]. 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, ∆E J = 0. However, also the coupling between the spins of the electronsŜ and of the nucleusÎ can remove the degeneracies and it is generally approximated by a dipole-dipole interaction This term, again, does not lead to the conservation of bothÎ andŜ separately but only of the total spin of the atomF =Î +Ŝ, 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 operatorsF 2 andF z . The effect of this contribution is a small correction to the Bohr spectrum, usually called hyperfine structure Since S = 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 ± 1 2 , so the hyperfine correction is given by 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 where the Zeeman couplings are inversely proportional to the mass of the particles, therefore γ N γ 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Î z andŜ z . Diagonalizing the Hamiltonian in Eq (5) we can find the corrections to the Bohr spectrum. For example, in the case of an hydrogen atom I = 1 2 and we get 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 traps 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, with a confining force F(r) = −∇V (r). 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 oneor 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 fardetuning, the confinement is independent from the internal state of the atom. When an atom is put in a laser beam, the electric field E(r) = E 0 (r)e iωt induces an electric dipole moment p(r) = α ω E 0 (r)e iωt , where α ω 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 which is simply directly proportional to the intensity I(r) 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 F(r) = −∇V (r) ∝ ∇I(r) . 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. 68,69 . The trapped 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. 70 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 µ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 71 . 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.

B. 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 18 two states (|F, F z ) of 87 Rb were chosen with pseudo-spins | ↑ = |1, 0 and | ↓ = |1, −1 . In order to understand schematically how it works, let us suppose that the atom is in an initial state | ↓ , 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 | ↑ . The inverse process is equally allowed. Due to Doppler effect, this coupling will be necessary momentum-sensitive and, as shown in Refs. 18,19 , the effective generated term is equivalent to a spin-orbit coupling of a spinfull particle moving in a static electric field, let us say, along the z axis. An effective additional term in the Hamiltonian can read where B SO is the magnetic field seen by the moving particle. This is called a Rashba coupling, in analogy to the same coupling occurring in crystals. With the same technique other interactions can be made in laboratory such as a Dresselhaus term or a Weyl term H W = λ σ · k. The two lasers are detuned by a frequency 2δ from the Raman resonance and couple the states | ↑ and | ↓ with the Raman strength 2Ω. This generates additional terms, Ωσ z and δσ y . We will neglect the detunig term and will consider the Raman coupling as an effective artificial Zeeman term 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.

II. 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 pseudopotential, 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. [72][73][74] . 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 = m 1 m 2 m 1 +m 2 for the reduced mass, in the relative frame the Schödinger equation reads 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, In this case Eq. (12) becomes a free Schrödinger equation whose solution is given by the composition of an incoming plane wave with momentum p and a scattering state with momentum p . It is relevant to look for elastic collisions, that can be studied fixing the scattering energy to E = 2 p = p 2 m , equal to the incoming wave energy. The solution is given by The homogeneous solution (V = 0) is simply |ψ p = |p , therefore, it is possible to write an implicit solution, known as the Lippmann-Schwinger equation with an infinitesimally small positive real number. In the coordinate representation, it reads ψ p (r) = r|p + r| 1 At long distances, (see Appendix A for details) Eq. (16) can be written as after defining q ≡ pr = p r r, and the scattering amplitude 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, Due to the spherically symmetric potential, the scattering amplitude depends only on the modulus p and the angle θ = arccos( q·p p 2 ) = arccos(r ·p). It is convenient, therefore, to perform a spherical wave expansion with P the Legendre polynomials. In order to preserve the normalization of the wavefunction, one has to impose the following condition for the coefficients f (p) (see Appendix B) which defines the phase shift δ (p) and where we rescaled, here and in what follows, p → p for simplicity, to get rid of . This relation imposes that when the relative momentum vanishes, also the phase shift should be zero, δ (0) = 0. The scattering amplitude, at fixed angular momentum, can be rewritten as In the radial Schrödinger equation (see Appendix C), for = 0, besides the potential V (r) 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 ( = 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 ≈ r * Λ th . For ultracold gases the temperatures are typically T ∼ nK and, since the thermal wavelength Λ th ∼ T −1/2 , it is allowed to consider ≈ 0, namely, the scattering amplitude is dominated by the s-wave contribution It is expected that scattering processes modified by the presence of V (r) for = 0 are greatly suppressed at sufficiently low energy. Moreover, also the = 0 contribution is affected by the low energy hypothesis. As shown in Ref. 73 , it is possible to define the scattering length at fixed considering the low energy limit which comes from the low energy limit of the phase shift 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 ∼ 1/Λ) it is possible to expand the phase δ 0 (p) = −a s p + O(p 2 ), so defining the scattering length as Actually only a s = a =0 has the dimension of a length, contrary to the other terms a with > 0. Considering the relation in Eq. (22) and the following expansion the scattering amplitude is then given, at low energy, by the expression 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 swave scattering length, and r eff , an effective range proportional to r * . We will not treat r 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 σ(p) = σ (p), we have σ =0 (p) = 8π p 2 (2 + 1) sin 2 (δ (p)) ∼

A. 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 operatorT 2B , called T scattering matrix, By this definition it is possible to rewrite the scattering amplitude as a function ofT 2B In the same way also the Lippmann-Schwinger equation, introduced earlier, can be written aŝ whose solution, in terms ofT 2B , is the followinĝ This solution can be expressed as an expansion in the potentialV Given a complete set of eigenstates ofĤ, inserting a completeness relation, one can writê |ψp ψp| z−2 pV , where α < 0 are the energies of the bound states, |ψ p the scattering states and Ω is a volume. It is then clear that the two-body scattering matrix T 2B (z) 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 Eq. (28) can be also rewritten in terms of the scattering energy E = 2 p 2 m so that At low energies (E → 0 − ) it is found that, for positive scattering length, a s > 0, this quantity has a simple pole at 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 Fano-Feshbach resonance method one can control by a magnetic field the value of a s . We refer to Refs. 10-13 for a detailed description of this technique.

B. Renormalization of the contact potential
Let us consider the contact potential in the position representation The action of this operator on the eigenstates of the momentum operator is We will want this potential to reproduce the scattering amplitude of a realistic potential. Exploiting the Lippmann-Schwinger equation for T 2B , written in Eq. (31), we have We now impose that, for the contact potential, the zero energy scattering amplitude should give the same result of a realistic potential. From Eq. (34) we have while from Eq. (38), using Eq. (37), we get This relation may seem inconsistent at first sight due to divergence of the last term. However, as we will see in what follows, this result will turn 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.

With spin-orbit coupling
As shown in Refs. 43,44 , in the presence of a spin-orbit coupling, the equation useful to cure the ultraviolet divergences is formally the same as Eq. (40), (see Eq. (D26)), 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).

III. 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 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 Eq. (41) is associated to the grand canonical partition function in the path integral formalism which is where the action is give by It is possible to take advantage of the properties of Gaussian integrals (see for instance Ref. 74 ) to handle the interacting term introducing the complex auxiliary field ∆(r, τ ) The grand canonical partition function can be rewritten, therefore, as it follows with the use of a new action, dependent now also on an auxiliary field The physical meaning of ∆ will be clear from the saddle point equation. It is more convenient to introduce the Nambu spinors so that the action of the system can be written as with G −1 , the inverse of the interacting Green function, defined as Actually in order to put S[ψ † , ψ, ∆] in a quadratic form in the spinor representation we exploited the anticommutation relations of the Grassmann fields, getting also the last term in Eq. (48). Defining z = e −β k ξ k , this term can be put outside the functional integral. The grand canonical partition function is the following Performing the standard Gaussian integral over the fermionic degrees of freedom one gets which defines the effective action depending only on the auxiliary field

A. Gap equation
Starting from this action it is possible to write down the saddle point equation and look for an homogeneous solution. With ∆( r, τ ) = ∆ 0 the effective action is given by where V is the volume. In the momentum space and in the Matsubara frequencies the matrix G −1 defined as in Eq. (49) is which has poles in real frequencies at ±E k where The minimum of the spectrum is at k = µ with the minimum energy |∆ 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|∆ 0 |. The value of this minimum energy can be obtained imposing δSe δ∆ 0 = 0, getting Performing the standard sum over the Matsubara frequencies and going to the continuum limit 1 V k → 1 (2π) 3 d 3 k, one gets the so called gap equation The contact potential is a great simplification in the model, however it gives rise to ultraviolet divergences. Actually Eq. (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 ω D of the underlying lattice. Typically, in classical superconductors ω D f = µ and a finite gap can exists only for particles near the Fermi energy. Exploiting the approximation of constant density of states (DOS), where ν 0 is the DOS at the Fermi energy, the gap equation becomes Solving the gap equation Eq. (58) one gets, at T = 0, Supposing that at a critical temperature, T = T c , the thermal fluctuations spoil the superconductivity, ∆ 0 (T c ) = 0, one gets with γ ≈ 0.5772, the Euler-Mascheroni number. Finally, for T < T c the gap is

B. 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 ∆(x) = ∆ 0 + δ∆(x), 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

Mean field
At the mean field level the action is proportional to the grand canonical potential from which is possible to derive the equation for the mean number of particles Summing over the Matsubara frequencies one gets Let us consider the cases at T = 0 and T = T c . At T=0. In order to study the behaviors of the gap ∆ 0 and of the chemical potential µ as functions of the coupling, at T = 0, at the mean field level, the equations to solve are where n = N/V, the particle density.
On the other hand if we want to calculate the critical temperature T c , putting by definition ∆ 0 (T c ) = 0, at the mean field level the equation to solve in terms of T c are

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 Eq. (52), in momentum space, can be written as Expanding the logarithm for small ∆ one gets The linear term in∆ has a vanishing trace, therefore, up to second order the action is given by Making explicit the quadratic term in ∆ we havê All the calculation can be done in momentum space where k = (k, iω n ) and q = (q, iν m ) are the four-momenta andĜ(k) is given bŷ whose components fulfill the property G 11 (k) = −G 22 (−k) that can be exploited to rewrite the quadratic term. Calling for simplicity G 11 (k) = G(k), we have where we defined the function χ(q) which, after performing the summation over the Matsubara frequencies iω n , reads where n f is the Fermi distribution. The action at second order in ∆, corresponding to the inclusion of the Gaussian fluctuation around mean field solution, is where we introduced the function The partition function at the Gaussian level is, therefore, with the definition Z 0 = e −S (0) e , which can be written as the grand canonical potential is shown to be The number equation, at the Gaussian level, is given by where The first term in Eq. (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 Γ −1 (q, z) , has poles in the complex plane at z = z j (j = 0, 1, 2 . . . ) and a branch cut on the real axis, the sum over Matsubara frequencies can be rewritten as In our case the integrand has poles only on the cut, so that this formula reduces to where n b is the Bose distribution. Since Γ(q, ω) has to be real on the real axis, it can be written as a modulus times a phase. Defining a phase shift δ(q) as in Ref. 4 , we can write At the Gaussian level, therefore, the density of particles can written in the following form 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 Eq. (57) has an ultraviolet divergence that has to be cured.

IV. 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 Sec. II B we showed that an effective contact potential, V (r) = −gδ(r), 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 Eq. (40), reported here for convenience Let us substitute the coupling g from Eq. (89) in the gap equation Eq. (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 The divergence in the gap equation is exactly canceled out by the divergence in Eq. (89). The limits of strong (g → ∞) and weak coupling (g → 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 Λ so that Eq. (89) becomes In the limit g → 0, namely in the BCS limit, the cut-off Λ does not play any role and eventually can be sent to infinity, therefore the scattering length is negative and 1 as → −∞. In the limit g → ∞, namely in the BEC limit, we get, instead, Sending the cut-off to infinity the limit of strong interaction can be captured for 1 as → +∞. 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.

A. Mean field theory
At the mean field level one has to simultaneously solve the gap and the number equations found previously, We will discuss in what follows the weak (BCS) and the strong (BEC) interacting limits.

Weak coupling
Let us consider the number equation Eq. (93), which, in the continuum limit, reminding that k = k 2 2m and introducing the density of state, becomes For T c very small one can resort to the well-known Sommerfield expansion, getting for the density of particle at low temperature , which is the typical experimental situation, this equation of state implies , and solving for µ we find an explicit relationship between the chemical potential and the temperature The chemical potential is, therefore, almost constant, µ f . Thermal corrections occur only at second order. Exploiting this result, we can solve the gap equation Eq. (92), as done in Ref. 74 . With constant chemical potential, Eq. (92) 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 = /µ, Eq. (92) becomes At µ = f and supposing T c f one gets π kfas = −2 ln( 8e γ f e 2 πTc ), with γ 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 a critical temperature exponentially small in the weak coupling limit, in agreement with the BCS theory and consistently with the hypothesis T c f .

Strong coupling
In the strong coupling limit, we will solve first the gap equation Eq. (92) to obtain the chemical potential µ which is expected to be very large and negative, as in the case of free fermions at high temperatures, and that |µ(T c )| T c . Taking advantage of this ansatz, whose validity will be proved a posteriori, we have The relation between the chemical potential and the scattering length, in the strong coupling limit, is, therefore, given by where we used Eq. (35) (here = 1), reminding, from Sec. II A, that for positive scattering length the two-body problem supports the existence of a bound state with energy E b . We will use Eq. (101) in the number equation Eq. (93) to derive the critical temperature. After introducing the density of state and after a rescaling, in the limit 1 as → +∞, we have At fixed density of particles along the crossover, , therefore in the strong coupling limit the leading term for the critical temperature at the mean field level is given by 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 ∆(T c ) = 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 / f 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 µ p = 2µ 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.

Along the crossover
Let us consider Eq. (92) and Eq. (93) in the continuum limit. Defining y = 1 kfas andμ = µ 2Tc , working at fixed density, we can write These integrals can be calculated numerically for different values ofμ, building the lists of values for µ f , y and Tc f , y . The graphs reported in Fig. 1 are obtained in this way, confirming the analytical results at mean field level.

B. Beyond mean field: Gaussian fluctuations
As seen previously, the second order expansion of the action is given by Eq. (77), which is with Γ −1 (q) given by Eq. (78). After regularizing the contact interaction through Eq. (89), and after a shift of the momenta k → k + q/2 in order, for simplicity, to get rid of the angular dependence in the denominator on the angle defined by k · q = kq cos θ, we can write The sum is invariant under k → −k, therefore the numerator [1 − n f (ξ k+q/2 ) − n f (ξ k−q/2 )] can be written as [1 − 2n f (ξ k+q/2 )] = tanh ξ k+q/2 2Tc , getting simply In the continuum limit and writing explicitly the spectra ξ k±q/2 = |k±q/2| 2 2m we have The dependence on cos θ is only in the hyperbolic tangent and we can integrate over it, using where A(k, q) = ln cosh We have, at the end, the following expression Let us consider the weak and the strong coupling limits.

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 → 0 − . Without entering into the details we have where f (q, µ) some complex function of the four-momentum q (momentum and frequency) and chemical potential µ. In the weak coupling limit, approximately Γ −1 (q) does not depend on critical temperature since, for small T c we have A(k, q) ∼ 1 2Tc (k+q/2) 2 2m − µ . Calculating explicitly the second order correction to the density of particles, we have ∂µ ∝ a s (113) 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, µ ≈ f , and since also the gap equation is unchanged, in the limit 1 k F as → −∞, 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.

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 Γ −1 . As we have seen, in the BEC limit, βµ → −∞ and 1/k f a s → ∞. For a large and negative chemical potential Inserting this result in Eq. (111), after integration, we get Taking advantage of the definition of the bound state where R(q) = 4π which is simply R(q) ≈ 8π m 2 as for small fluctuations in the deep BEC limit. The quantity µ 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 Rescaling the fields φ(q) = ∆(q) √

R(q)
and defining m b = 2m the partition function at the Gaussian level is equivalent to that for free bosons with mass m b (equal to the mass of two fermions) Now we can take advantage of the well-known expression for the critical temperature for the Bose-Einsten condensation which is given by If in the deep BEC regime all fermions form pairs, therefore n b = n/2. At fixed fermionic density According to the BEC theory, at the critical point, the effective chemical potential approaches zero from negative values, therefore which is exactly the solution of the gap equation in the strong coupling limit, Eq. (101).

Full crossover
As shown in Sec. III B 2, the set of equations to solve at second order is with the phase shift δ(q, ω) = −Arg Γ −1 (q, ω + i ) . We need, therefore, to find the second term of the number equation, n (2) , possibly resorting also to some approximations. a. 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 which describes a system of free bosons, with µ b = 2µ − E b and m 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 (2) in the number equation Eq. (124) with the density of free bosons whose inverse propagator is given by Eq. (125). It is well known that the critical temperature T c of free bosons is given by Eq. (120). Let us suppose that the number of bosonic excitations composed by pairs of fermions is n b = 2n (2) , while the rest of fermions remains unpaired. In this approximation, using Eq. (120), we have, therefore, As a result, we have to solve Eq. (123) together with the number equation which, in the bosonic approximation, reads After defining the following quantities in the continuum limit, Eq. (123) and Eq. (127), at fixed density n = k 3 f 3π 2 can be written as Making the same substitutions in the vertex function, after a Wick rotation in the frequencies (ν n → ω +i , with > 0 an infinitesimal quantity), we can define Γ −1 (q, ω +i ) = Γ −1 (Q, ν +i ), with the rescaled real frequency ν = ω 4Tc , We remind that, looking at Eq. (88), the interesting part of the function Γ −1 (Q, ν +i ) is actually its argument, therefore, its prefactor is irrelevant. Within the range of integration the integrand has a pole in p 0 = μ + ν − Q 2 /4, if p 0 real and positive. Defining the function we can write with the prefactor C = − m For (μ + ν − Q 2 /4) > 0, namely for p 0 real and positive, Γ −1 has an imaginary part. We can write the vertex function Γ −1 and its argumentδ as it follows where the real and imaginary parts of Γ −1 are given by 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 Eq. (88), which by rescaling all the parameters, reads Writing the number of particles in terms of f , the final number equation can be rewritten as We have, therefore, to solve this equation together with the gap equation Eq. (103). 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, 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.

V. 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 Sec. I B. Let us consider the following Hamiltonian H = H 0 + H I , in the momentum space, made by a single-particle term and an interacting term where σ x = 0 1 1 0 and σ y = 0 −i i 0 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 r (σ x k y − σ y k x ), Eq. (9), and the Dresselhaus term, v d (σ x k y + σ y k x ), Eq. (10). We included also an effective artificial Zeeman term hσ z , Eq. (11), which can be tuned by a Raman laser 18 . This Hamiltonian is associated, in the path-integral formalism, by the Grassmann variables ψ and ψ † , to the partition function with the action where ψ s (k) and ψ † s (k) are Grassmann variables depending on the four-momentum k = (ω n , k). As already seen, the term with four fermions can be decoupled by means of an Hubbard-Stratonovich transformation, introducing the auxiliary field ∆(q) in this way the partition function takes the form Z = Dψ † Dψ D∆ * D∆ e −S[ψ † ,ψ,∆] where, after defining the new action reads It is now convenient to introduce the following Nambu-Jona-Lasinio multispinors so that Eq. (148) can be written as where the sum β k ξ k comes from interchanging the fermionic fields and the corresponding equal-time limiting procedure in order to write the action in a matrix form 74 , namely In Eq. (150) we also introduced the following inverse Green function 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 obtaining an effective action depending only on the pairing function A. Gap equation As done in Sec. III, we first search for the homogeneous pairing which minimizes the action, selecting only the contribution at q = 0 in momentum space. Calling ∆(q = 0) = ∆ 0 the determinant appearing in Eq. (153) is simply where the energies are and where (147), and ξ k = k 2 /2m − µ. v r and v d are the Rashba and Dresselhaus velocities and h is the Zeeman field. The effective action, therefore, takes the form The saddle point equation sets the mean field solution and is given by the condition ∂Se ∂∆ 0 = 0 After summing over the Matsubara frequencies the gap equation, in the presence of spin-orbit interaction and Zeeman term, reads (160)

B. 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 Eq. (62). The zeroth order is equivalent to the mean field approximation. We will then include the fluctuations at the Gaussian level.

Mean field
The effective action, at mean field level, is the action in Eq. (158), evaluated at the homogeneous solution of the gap equation, corresponding to the saddle point approximation. The associated grand canonical potential is simply Ω[∆ 0 , ] = βS e [∆ 0 ]. The mean number of particles, defined by N = − ∂Ω ∂µ is, therefore, given by Summing over the Matsubara frequencies, the number equation, at the mean-field level, reads (162) Let us consider 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 where n = N/V is the particle density. This set of equations, for h = 0, has been solved in Ref. 31 . 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 ∆ 0 (T c ) = 0, at the mean field level the equation to solve in terms of T c are the following 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 . (167)

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 Eq. (166) 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 [∆] = S (0) e + S (2) e , in addition to S (0) e , 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, The inverse of the Green function can be written in momentum space and in the Matsubara frequency space as in Eq. (151). It can be written as We can now perform the series expansion of the logarithmic function 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, Tr(G∆) 2n+1 = 0. The action at second order in∆ is, then, given by The inverse of G −1 , the non-interacting Green function, in momentum space, is given by Making the matrix product G∆G∆ and taking the trace we get where G ij are the matrix elements of G. In momentum space Eq. (175) reads where we define the function , therefore the action, at the second order in ∆, namely at the Gaussian level, can be written as where the same notation for Γ −1 (q), used for the case without SO coupling, has been adopted but where now χ(q) is more complicated than before. From Eqs. (174) and (177) we get It is convenient to make the shift k → k + q 2 and, after performing the sum over the Matsubara frequencies, we obtain For v d = 0, namely with only the Rashba coupling, we have where k ⊥ = (k x , k y , 0) are momenta perpendicular to the z-direction, the Zeeman direction.
For v r = v d = v/2, namely for equal Rashba and Dresselhaus couplings, we have, instead, which, for h = 0, namely without the Zeeman field, it reduces to simply C 0 k,q = 1. Employing the symmetry k → −k under the sum we can simplify Eq. (181) as it follows or, in terms of the hyperbolic tangent, We found, therefore, that Γ −1 (q) 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 Eq (79), and, therefore, it is possible to apply the same procedure seen before (in Sec. III B 2) to find the grand canonical potential and the density of particles. This is again due to the fact that the analytical continuation of Γ −1 (q, z) has only a branch-cut of simple poles in the real axis. The density of particles is, then, given by where n (0) is the mean field contribution as in Eq. (166) while the second order term is given by where the phase shift δ(q, ω) is defined as before after a Wick rotation iν n → ω + i .

C. Critical point
Taking the static and homogeneous limit of Eq. (187), after noticing that, at q = 0, we have that, at T = T c , χ(0, 0) is simply given by We find that the equation Γ −1 (0) = 0, namely is exactly the same equation reported in Eq. (165) after reshuffling the terms. The set of equations that has to be solved, at the Gaussian level, in order to find T c and µ at T c is given by Eq. (193), together with Eq. (188).

VI. 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 and in Refs. 43,44 , the presence of a SO coupling does not alter the form of the renormalization condition for a contact pseudopotential. 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 therefore, the arguments used to implement the BCS and BEC limits in terms of the scattering length remains unchanged.

A. Mean field
At the mean field level the equations to solve along the crossover are Eqs. (165), (166). Implementing Eq. (194), for h = 0, the equations useful to derive the critical temperature T c and the chemical potential µ within the mean field theory become where γ k ≡ γ 0k , as reported in Eq. (157) with h = 0. Also 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 Eq. (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 v d = vd vf = mvd kf and similarly for v r . In the continuum limit, with the natural change of variables and substitutions p = 1

the gap equation can be written as
while the number equation, Eq. (196), becomes

Weak coupling
At low temperature we expect that the chemical potential of a system of weakly interacting fermions is almost constant and therefore, in this regime, µ T 1. Actually, one can shown perturbativly 31,42 that µ ≈ f − m(v r + v d ) 2 /2. In this limit, the number equation can be neglected due to the ansatz Eq. (199), therefore Eq. (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.

Strong coupling
In the strong coupling limit, supposing that µ/T c → −∞, which can be verified a posteriori, Eq. (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 so that Eq. (197), after introducing the dimensionless coupling y = 1 kfar , becomes simply We can integrate over p z getting For simplicity, let us consider the case with pure Rashba coupling, so thatγ p = 2ṽ r p 2 x + p 2 y ≡ 2ṽ r |p ⊥ |. The above equation can be written as With the substitutions p ⊥ → p ⊥ ±ṽ r in the first two integrals we have where Λ is the ultraviolet cut-off which can be sent to infinity, Λ → ∞. Eq. (205) can be rewritten as which gives Taking the square, y 2 ≈ −μ − 2ṽ 2 r , and using the definitionsμ = µ/ f andṽ r = v r /v f , reminding that y = 1 kfar , we get Therefore, in the strong coupling limit (BEC limit), 1/a r → ∞, we found µ − 1 2ma 2 r , namely the spin-orbit interaction can be neglected and we get the same expression for the chemical potential as that obtained previously in Eq. (101), 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.

Full crossover
Let us consider Eq. (195), whose form, in the continuum limit and with rescaled dimensionless Eq. (197), which can be rewritten in terms of y = 1 kfar as it follows The number equation, Eq. (196), in the continuum limit and in terms of the same dimensionless parameters, is Reminding that the density of particles is n = (2m f) 3/2 3π 2 = k 3 f 3π 2 , the equations to solve are We can search for the solutions of these two equations, Eqs. multidimensional root-finding algorithm: for any fixed values of the couplings y,ṽ r andṽ d we look for the critical temperature T c and the chemical potential µ that satisfy those equations simultaneously. In Figs. 4 and 5 we report some results for the case with only Rashba coupling. Special case: equal Rashba and Dresselhaus couplings. For v r = v d ≡ v/2 we have γ k = v|k y |, or, in terms of the rescaled quantities p y = k y /k f andṽ = v/v f , in a symmetric combination so that we can remove the absolute value symbols so that after changing the variables p y → p y ±ṽ we get simply +∞ −∞ dp y F p 2 y −ṽ 2 −μ 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 Fig. 1) while µ at finite v is linked to that without SO (v = 0) in the following way We will consider for simplicity the case with h = 0 but let us first discuss two special cases.

Special case: equal Rashba and Dresselhaus
Let us consider the case where h = 0 and v r = v d , namley when Rashba and Dresselhauf couplings are equal. The equations are (226) where γ k = γ 0k and since v r = v d , while δ(q, ω) is defined in Eq. (221). Since C 0 k,q = 1, we have that χ(q, iν n ) is simply given by As already discussed, at the mean field level, after rescaling k y → k y ± mv we obtain Eqs. (92) and (93), namely the mean-field equations without SO, but where the chemical potential is µ + mv 2 /2, so that Eqs. (225), (226) can be simplified as it follows

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 γ hk = h. The equations become simply with δ(q, ω) still defined by Eq. (221) but where χ(q, iν n ), since C h k,q = −1, is simply given by These are the equations to solve for finding the critical temperature of a polarized Fermi gas 75,76 which should be the same in the presence of a Rabi coupling.

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 → 0 − . As argued in the case without SO coupling, Γ −1 (q) can be expressed as in Eq. (112). 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 Eq. (113). Since the chemical potential is almost constant and since also the gap equation is unchanged, in the limit 1 k F as → −∞, 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.

Strong coupling
In the strong coupling limit, a r → 0 + , namely in the BEC limit, with SO coupling we found that µ(T c ) diverges with the inverse scattering length y = 1 kfar , Eq. (208), getting the same formal result as that without SO coupling, µ → − 1 2mar . 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, Γ −1 (q), in the strong coupling limit, can be seen as the inverse propagator of free bosons, Eq. (125), with an effective mass and a redefined chemical potential, µ b = 2µ − 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. 24,27,28 , doing a second order expansion in the momenta for Γ −1 (q) in the BEC limit. The two-body physics emerges from Γ −1 (q), Eq. (222), in the so-called vacuum limit, discarding the Fermi distribution, n f → 0, and putting the chemical potential to zero 24,27,28 . What remains is the inverse of the two-body scattering matrix, T −1 2B (q, iν n ) ≡ Γ −1 (q, iν n )| nf,µ=0 , which, reads Bound state. 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, Sec. II A, how T 2B (z) 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, From Eq. (224) we have C 0 k,0 = 1, so that Eq. (235) implies the following equation 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 p = k which has the same form of Eq. (201). The chemical potential and the bound state energy, E b = fẼb , measured at the threshold energy m(v r + v d ) 2 , in terms of dimensionless quantities, satisfy the same equation, therefore, in the the strong coupling limit, they are connected by the In the deep BEC, neglecting the second term, we have which is the same result obtained in the absence of SO coupling, Eq. (122). Taking advantage of Eq. (207), in the case of only Rashba coupling (v d = 0), we can write an analytical solution Inverting Eq. (240) we can get E b as a function of y. For equal Rashba and Dresselhaus couplings, v r = v d , we have, instead, 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 Eq. (237), or, after analytically integrating over p z , Eq. (203) withẼ b /2 instead ofμ. The results are shown in Fig. 6. For pure Rashba coupling, Fig. 6, a bound state exists also in the BCS part of the crossover, where the scattering length is negative. Actually by looking at Eq. (240), for any value of y, a negative value of E b exists. Moreover, in Fig. 6 we show that turning on a Dresselhaus term (v d = 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 f a s → +∞) can be written as where 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 (q) as in Eq. (242) At low energy and momentum this condition is imposed to the second order expansion in q which implies that The first equation in Eq. (245) is the same encountered before for the bound state energy, Eq. (235). By choosing m z = 2m we get immediately that F z (E b , 2m) = 0 since, in this way, T −1 2B (q, b (q)) loses the dependence on q z . Actually the two-body scattering matrix reads Introducing the following definitions the two-body scattering matrix can be written as it follows Expanding for small q the following quantities the two-body scattering matrix, at leading orders, reads The first term is the scattering matrix at q = 0, therefore, T −1 2B (0, E b ) = 0, by definition of the shifted bound state energy E b . Solving F x (E b , m x ) = 0 in terms of X and F y (E b , m y ) = 0 in terms of Y , and from Eq. (247), we can derive the values of the effective masses In the specific case of equal Rashba and Dresselhaus couplings (v r = ±v 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 ≥ 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 (256) and (257) in terms of the dimensionless parameters become The effective masses can be written in terms of dimensionless quantities For the case with a pure Rashba coupling, for which m x = m y , we can perform the integralsÃ, B x andB y analytically, obtaining where, from Eq. (238), in the presence of a bound state, we assume that The final analytical expressions for the effective masses are, therefore, which, together with Eq. (240) written in terms ofẼ b in order to get their behaviors along the crossover, see Fig. 7. This result is in agreement with the T -matrix approach 50 . At the special point y = 1 kfar = 0, also called unitary limit, the ratio 2ṽ 2 r /Ẽ b ≡ −c o is the solution of the following equation which is c o ≈ 0.6948. As a result 2m mx = 2m my at y = 0 is the same for any value ofṽ r , 2m m x a −1 which is the nodal point shown in Fig. 7. For the general case of mixed Rashba and Dresselhaus terms the effective masses are given by Eqs. (261), (262), together with Eqs. (258)-(260).

Full crossover
a. Bosonic approximation. We found that the inverse of the scattering matrix, expanding at low momenta around q = 0, has the form with E b , m x and m y that depend on y = 1 kfar . This result for the two-body problem can be exploited to determine the properties of the vertex function in the strong coupling limit getting, in this regime, 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 where m b = 2m and, in the deep BEC, We find that in the BEC limit the effective masses approach those of two fermions, 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 b = 2m emerges. From Eq. (273), solving in terms of n b , we have As done in Eq. (126), let us suppose now that the number of bosonic excitations composed by pairs of fermions is n b ≈ 2n (2) , while the rest of fermions remains unpaired so that we have With this approximation the set of equations to solve is the following In the continuum limit, from Eqs. (209), (210), at fixed density n = k 3 f 3π 2 , after rescaling p = 1 kf k, , the above equations can be written in the following form This set of equations can be easily solved numerically. Some results are reported in Fig. 8 for the most relevant case with pure Rashba coupling. These results have to be contrasted with those reported in Fig. 2 in the absence of spin-orbit couplings. b. 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 Γ(q), presented in Eq. (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 q the full Γ −1 (q) showed in Eqs. (222), (223), 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 q = (ξ k±q/2 ± γ k±q/2 )/2, appearing in Eq. (223), up to second order in q 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 q /T c → ∞, therefore tanh(f q /T c ) → 1, recovering the results in the bosonic approximation, while in the weak coupling limit the quantum fluctuation themselves are not relevant at all, as discussed previously. Alternatively one can calculate the exact Gaussian fluctuations solving numerically the full equations in Eqs. (218), (219) together with Eqs. (221), (222), (223), 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 computational costly.

VII. 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 31 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 > 2m 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.
which, inserted in Eq. (A1), gives Eq. (16). In the large r limit, using the expansion 1 |r−r | = where q = pr/r and p → p, to get rid of in the expression, and where we drop, for simplicity, the prefactor 1/(2π) 3 . For elastic scattering on a short-range spherically symmetric potential, the scattering amplitude can depend now only on the modulus p and the angle θ = arccos( q·p p 2 ). It is so correct to perform a spherical wave expansion of the function f (p, θ) as in Eq. (20). Let us expand the plane wave also on this basis, e ip·r = ∞ =0 (2 + 1)i  (pr)P (cos θ) The  (pr) are the usual spherical Bessel functions, whose large distance expansion is i  (pr) = i e ipr− π 2 − e −ipr+ π 2 2ipr + O(r −2 ) = e ipr − e −i(pr− π) 2ipr + O(r −2 ) (B3) The scattering wavefunction is then a superposition of incoming and outgoing spherical waves ψ p (r) = 1 2ipr ∞ =0 (2 + 1) −e −i(pr− π) + (1 + 2ipf (p))e ipr P (cos θ) + O(r −2 ) (B4) We found, then, that the complex amplitude of the outgoing spherical wave is modified by the presence of the interaction, however, to preserve normalization, its modulus has to be constant.
In the presence of particles with spin, instead, it is necessary to take into account the spin part of the wavefunction. For example, for fermions, if the scattering happens in a singlet (antisymmetric) state |S , we should impose that the wavefuction is symmetric in the position space |ψ p (r) = e ip·r + e −ip·r + f (p, θ) e ipr r + f (p, π − θ) e −ipr r |S (B12) By these considerations the cross sections for singlet (σ S ) and triplet (σ T ) spin-1/2 fermions result σ S (p) = 8π p 2 even (2 + 1) sin 2 (δ (p)) , σ T (p) = 8π p 2 odd.
where a s is defined in Eq. (26), therefore the corresponding wavefunction is The same form for the wavefunction can be obtained imposing the so-called Bethe-Peierls boundary conditions 1 rψ p (r) One should evaluate the scattering length for a realistic interatomic potential, however in many cases this requires difficult calculations. For this reason it is often convenient to use a simpler pseudo-potential that can reproduce the same scattering amplitude at low energies. Fortunately, the relevant experimental parameter in ultracold atomic physics, tunable by the so-called Feshbach resonance technique, is the scattering length a s . spin spaces. For clarity different brackets are associated to different spaces, |· ∈ H, |·) ∈ H r , |· ∈ H s 1 ⊗ H s 2 , so that the position dependent states of the spin is such that (r|ψ = |ψ(r) . The same study on two-body problem, seen before, should be done also in this case, always for a spherically symmetric and spin-independent potentialÛ (X (1) −X (2) ), where X (i) are the coordinates of the particles. The total two-body Hamiltonian is then 1b +Ĥ (2) 1b +Û (X (1) −X (2) ) (D2) This can be easily separated in the Hamiltonian of the center of mass and Hamiltonian of the relative motion. The latter Hamiltonian iŝ with the definitions of the relative quantities, p = P (1) −P (2) 2 , K = P (1) + P (2) , r = X (1) − X (2) , c = M (1) − M (2) , B(K) = Z (1) + Z (2) + λK · M (1) + M (2) , m = m 1 m 2 m 1 +m 2 . The total momentum K is conserved during the scattering process and it can be treated as a constant. A basis of eigenstates for the free Hamiltonian is given by the tensor product of eigenstates of the momentum operator and of spin eigenstates |ψ (0) k,α = |k) ⊗ |k, α with α = 0, 1, 2, 3, labeling the four configurations of two spins, also momentum dependent, which in relative position representation are (r|ψ (0) k,α = |ψ (0) k,α (r) = 1 (2π) 3/2 e ikr |k, α . Since the momentum part of the eigenfunction is known (plane waves as in the previous case) the eigenproblem can be restricted exclusively to the spin part. In the relative coordinates the HamiltonianĤ 0 can be decomposed in the following formĤ 0 =Ĥ 1 2×2 + β(p, ξ)σ x + γ(p, ξ)σ y + η(p, ξ)σ z = The coefficients of the Pauli matrices depend on the specific parameters of the model and linearly by the momentum; for example, η should have the form η(p, ξ) = v(ξ)·p + q(ξ), for some v and q model dependent set of parameters. The eigenvalues ofĤ (1) 1b are given by
For convention | ↑ 1 is the eigenvector associated to E + and | ↓ 1 to E − . The same problem is identically solved for H (2) 1b with the exception of a sign which has to be considered and so the eigenvalues will have the same form but with a different sign in the spin-orbit part of the spectrum. The eigenvectors | ↑ 2 and | ↓ 2 are identically assigned and the energies are With a proper unitary transformation it is always possible to put the Hamiltonian H 0 in a diagonal form The eigenvalues of the total Hamiltonian can be obtained building a basis of eigenstates; with a simple change in the notation |m i = |m i , with m = ↑, ↓, this can be done building the basis |m 1 , m 2 = |m 1 ⊗ |m 2 (D10) The associated eigenvalues are thereforê The four states of spin are labeled by α = 0, 1, 2, 3 and, for simplicity, we will indicate with t = (α, k), which implicitly depends on ξ, the specific eigenstate of the free Hamiltonian with eigenvalue E t . To obtain a renormalization relation similar to Eq. (40) it is necessary to write the Lippmann-Schwinger equation for the stationary scattering state |ψ t (r) = |ψ t |r ) = α d 3 k e ik·(r−r ) E − E t + i0 + |k, α k, α|.
(D14) The breaking of the rotational symmetry induced by the SO term makes hard the evaluation of the Green's function and its asymptotic behavior at large distances, therefore, the treatment done before is not easily adaptable to this case. However, as proved in Ref. 44 , the scattering