1 Introduction
1.1 Overview of conventional thermal equilibrium BCS-BEC crossover
An [1] ultracold Fermi gas is an extremely dilute gas of charge-neutral Fermi atoms, which is cooled down to \(\mathcal {O}\)(nK). The strength of a pairing interaction between Fermi atoms can experimentally be tuned by adjusting a Feshbach resonance [1,2,3,4,5,6], where two Fermi atoms form a quasi-molecular boson, and it dissociates into two Fermi atoms again. Within the second-order perturbation theory with respect to the Fershbach resonance, one obtains the effective interaction between Fermi atoms, given by [5]
Here, g is a coupling constant of a Feshbach resonance, and \(2\nu\) is the energy difference between the intermediate molecular state (closed channel) and the atomic states before and after the scattering event (open channel). The energy \(2\nu\) is referred to as the threshold energy of a Feshbach resonance. Since the atomic hyperfine states in the closed channel are different from those in the open channel, their Zeeman energies also become different under an external magnetic field B. This naturally leads to the B-dependent threshold energy \(2\nu\). Thus, one can tune the strength of the Feshbach-induced effective interaction \(U_{\text{eff}}^{\text{FR}}\) in Eq. (1), by adjusting the magnitude B of an external magnetic field.
An advantage of this tunable interaction is the realization of the crossover phenomenon from the Bardeen-Cooper-Schrieffer (BCS) state of Cooper pairs to the Bose-Einstein condensation (BEC) of diatomic molecules, as schematically shown in Fig. 1 [7,8,9,10, 3, 11,12,13]. Figure 2 shows the first observation of the superfluid phase transition, as well as the BCS-BEC crossover phenomenon in a \(^{40}{\text{K}}\) Fermi gas. In this experiment, \(^{40}{\text{K}}\) atoms in two different hyperfine states (which are frequently described by pseudospin \(\uparrow\) and \(\downarrow\) in the literature) are trapped in a harmonic potential, and Cooper pairs are formed between them by a Feshbach-induced tunable s-wave pairing interaction. In Fig. 2, \(\Delta B= B-B_{\text{res}}\) is an external magnetic field, measured from the Feshbach resonance field \(B_{\text{res}}\). Physically, \(\Delta B\) is directly related to the strength of a pairing interaction, that is, the decrease of \(\Delta B\) corresponds to the increase of the interaction strength. In particular, the region \(\Delta B > 0\) (\(\Delta B < 0\)) is the weak-coupling BCS (strong-coupling BEC) side. In this experiment, the superfluid state is identified as the parameter region where the number \(N_0\) of condensed Cooper pairs takes a non-zero value. Thus, the superfluid phase transition temperature \(T_{\text{c}}\) exists around the sky-blue area in the T-\(\Delta B\) plane in Fig. 2. We thus find that, starting from the weak-coupling regime (\(\Delta B >0\)), \(T_{\text{c}}\) gradually increases with increasing the strength of pairing interaction, to approach a constant value in the strong-coupling regime (\(\Delta B <0\)), which is known as a typical BCS-BEC crossover behavior of \(T_{\text{c}}\) [7,8,9, 3, 10,11,12,13].
Illustration of BCS-BEC crossover in a thermal equilibrium Fermi gas. Largely overlapping Cooper pairs in the weak-coupling BCS regime gradually shrink with increasing the strength of a pairing interaction, to eventually become the BEC of diatomic molecules
Observed superfluid phase transition in a \(^{40}{\text{K}}\) Fermi gas [3]. The temperature T is normalized by the Fermi temperature \(T_{\text{F}}=0.35{\mu \text{K}}\). \(\Delta B=B-B_{\text{res}}\) is an applied external magnetic field, which is measured from the Feshbach resonance field \(B_{\text{res}}\simeq 202\) G. This axis physically represents the strength of an s-wave pairing interaction associated with a Feshbach resonance. \(N_0/N\) is the number of condensed Fermi atoms, being normalized by the total number N of Fermi atoms. This experiment regards the region with \(N_0>0\) as the superfluid phase, so that the superfluid phase transition temperature \(T_{\text{c}}\) exists around the sky-blue area in the T-\(\Delta B\) plane. Adapted from Ref. [3]
Soon after the realization of the above-mentioned \(^{40}\)K superfluid Fermi gas, the superfluid phase transition, as well as the BCS-BEC crossover phenomenon, has also been realized in a \(^6\)Li Fermi gas [14,15,16]. At present, one can examine superfluid properties in the whole BCS-BEC crossover region by using these two kinds of Fermi gases in cold-atom physics. We also briefly note that the BCS-BEC crossover has recently been discussed in condensed matter physics, such as superconductors like FeSe and ZrNCl [17,18,19,20], as well as astrophysics, such as neutron-star interior [13, 21,22,23].
1.2 Strongly interacting Fermi gas in a driven-dissipative steady state
In cold-atom physics, since the achievement of the superfluid phase transition in \(^{40}\)K and \(^6\)Li Fermi gases, strong-coupling properties in the BCS-BEC crossover region have mainly been studied in the thermal equilibrium case. This is because the usual experimental situation of a trapped Fermi gas is well isolated from the environment and also in the (quasi)equilibrium state. However, a strong interest in non-equilibrium properties of this strongly interacting quantum system has recently emerged, being fueled by the rapid development of non-equilibrium quantum many-body physics [24,25,26,27]. For example, quench dynamics [28,29,30,31] and transport phenomena [32, 33] in a strongly interacting Fermi gas have been experimentally studied. Periodically driven systems [34, 35] and open systems with particle loss [36] have also been realized in an ultracold Fermi gas.
Motivated by such a recent trend in cold-atom physics, in this review, we consider a driven-dissipative Fermi gas, as schematically shown in Fig. 3a [37,38,39,40]. A three-dimensional two-component Fermi gas (main system) with a tunable pairing interaction \(-U\) is coupled with left (L) and right (R) reservoirs, consisting of non-interacting Fermi gases. Both reservoirs are assumed to be huge compared to the main system and are in the thermal equilibrium state at temperature \(T_{\text{env}}\). The main system becomes out of equilibrium when one imposes the chemical potential difference between the reservoirs as \(\mu _{\text{L}}=\mu +\Delta \mu\) and \(\mu _{\text{R}}=\mu -\Delta \mu\), where \(\mu _{\alpha ={\text{L}}, {\text{R}}}\) is the Fermi chemical potential in the \(\alpha\)-resevoir [see Fig. 3b]. Thus, by tuning the pairing interaction strength \(-U\), as well as the chemical potential difference \(\Delta \mu =[\mu _{\text{L}} -\mu _{\text{R}}]/2\) between the reservoirs, we can study non-equilibrium BCS-BEC crossover physics in the main system.
a Model driven-dissipative Fermi gas with a tunable paring interaction \(-U\) [37,38,39,40]. The main system is coupled with the left (\(\alpha ={\text{L}}\)) and the right (\(\alpha ={\text{R}}\)) reservoirs, which are assumed to be in the thermal equilibrium state at the temperature \(T_{\text{env}}\) and the chemical potential \(\mu _{\alpha }\). \(\mathcal {T}_\alpha\) describes the tunneling amplitude between the main system and the \(\alpha\) reservoir. b Schematic energy band structure of the model. We measure the energy from the bottom (\(\varepsilon _{\boldsymbol{p}=0}\)) of the energy band in the main system. In the \(\alpha\) reservoir at \(T_{\text{env}}=0\), the energy band \(\xi ^{\text{res}}_{\boldsymbol{q}}=\boldsymbol{p}^2/(2m) -\mu _{\text{res}}\) is filled up to \(\mu _\alpha\). When the reservoirs have different Fermi levels (\(\Delta \mu \ne 0\)), the main system is driven out of equilibrium due to the pumping and decay of Fermi atoms by the two reservoirs
The model non-equilibrium Fermi gas in Fig. 3a is inspired by the recent transport experiment on a \(^6\text{Li}\) Fermi gas in a two-terminal configuration [41,42,43,44,45]. In this experiment, a two-component Fermi atomic cloud is shaped in a two-terminal setup, where two reservoirs are connected by a mesoscopic channel. By extending this two-terminal setup to the case with multiple junctions, we can prepare a system coupled with multiple reservoirs, just like the model in Fig. 3a. We note that a narrow repulsive potential barrier produced by a Gaussian beam [46,47,48] is also useful to separate a Fermi gas cloud into multiple subsystems. Indeed, a Josephson junction of a \(^6\text{Li}\) Fermi gas is implemented with this technique. By extending the technique, we could divide a Fermi gas into the main system and the two reservoirs shown in Fig. 3a.
To realize the model in Fig. 3a, we need to apply an external magnetic field only to the main system, to adjust the Feshbach-induced pairing interaction \(-U\) there. This could be done by using the technique of the spatial manipulation of interaction strength [49,50,51,52,53,54]. In particular, in a two-component \(^6\text{Li}\) Fermi gas, the combination of the magnetic Feshbach resonance technique and the optical control enables us to adjust a scattering length with high spatial resolution [52,53,54]. Thus, using this technique, one can independently tune the interatomic interaction in the reservoirs and the main system.
1.3 Outline of this article
In this review article, we discuss the non-equilibrium BCS-BEC crossover in the model shown in Fig. 3a, as a paradigmatic example of non-equilibrium quantum many-body phenomenon in strongly interacting Fermi gases. We explain that the driving and the dissipation lead to exotic states that are never obtained in the thermal equilibrium case.
This article is organized as follows. In Section 2, we quickly review the BCS-BEC crossover theory for a thermal equilibrium Fermi gas. We also explain single-particle properties in the crossover region. In Section 3, we deal with the driven-dissipative Fermi gas in Fig. 3a, to consider the non-equilibrium superfluid transition within the mean-field approximation. We show that a non-equilibrium particle energy distribution induces a Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) like unconventional Fermi superfluid in spite of the absence of any spin imbalance. We then proceed to the non-equilibrium BCS-BEC crossover phenomenon in Section 4. Here we extend the thermal equilibrium BCS-BEC crossover theory to the driven-dissipative Fermi gas. We then elucidate non-equilibrium properties of the strongly interacting driven-dissipative Fermi gas. We find that the FFLO-like Fermi superfluid obtained in Section 3 is actually unstable against pairing fluctuations. In Section 5, we present a possible route to stabilize this unconventional Fermi superfluid state.
Throughout this article, we set \(\hbar =k_{\text{B}}=1\), and the system volume is taken to be unity, for simplicity.
2 BCS-BEC crossover in a thermal equilibrium Fermi gas
As a prelude to non-equilibrium BCS-BEC crossover physics in a driven-dissipative Fermi gas, this section reviews single-particle properties of a thermal equilibrium Fermi gas in the BCS-BEC crossover region.
2.1 Strong-coupling theory for a thermal equilibrium Fermi gas
To theoretically describe the observed BCS-BEC crossover behavior of \(T_{\text{c}}\) in Fig. 2, we need to go beyond the mean-field approximation, to include effects of strong pairing fluctuations. This has extensively been attempted by many researchers by various methods, such as quantum Monte Carlo method [55,56,57], renormalization group method [58,59,60,61], functional integral method [62], virial expansion [63, 64], as well as diagrammatic method [8, 9, 11, 65,66,67]. Among them, here we explain the T-matrix approximation (TMA) [68], which is a strong-coupling theory based on the diagrammatic technique.
Although the Feshbach-induced pairing interaction \(U^{\text{FR}}_{\text{eff}}\) in Eq. (1) is quite different from the phonon-mediated one in conventional metallic superconductivity, we can still capture the essence of the BCS-BEC crossover phenomenon, by employing the ordinary BCS Hamiltonian [13], given by
Here, \(a^\dagger _{\boldsymbol{p}, \sigma }\) creates a Fermi atom with momentum \(\boldsymbol{p}\) and (pseudo)spin \(\sigma =\uparrow , \downarrow\). \(\xi _{\boldsymbol{p}} =\boldsymbol{p}^2/(2m) -\mu\) is the kinetic energy measured from the chemical potential \(\mu\), where m is an atomic mass. \(-U\) (\(<0\)) is a contact-type pairing interaction, which is assumed to be tunable by a Feshbach resonance. To remove the ultraviolet divergence coming from the contact-type interaction, we conveniently measure the interaction strength in terms of the s-wave scattering length \(a_s\) [69]. The scattering length \(a_s\) is related to the bare interaction U as
where \(p_{\text{c}}\) is a momentum cutoff. The weak-coupling BCS side and the strong-coupling BEC side are then characterized by \((p_{\text{F}}a_{\text{s}})^{-1} \lesssim 0\) and \((p_{\text{F}}a_{\text{s}})^{-1} \gtrsim 0\), respectively. Here, \(p_{\text{F}}=(3\pi ^2 N)^{1/3}\) is the Fermi momentum of a two-component free gas with N fermions.
In TMA, we perturbatively include effects of pairing interaction \(H_{\text{int}}\) in Eq. (2). In the thermal equilibrium state, this is usually done by using the imaginary-time Matsubara Green’s function technique [70]. However, for later convenience, we employ the real-time Keldysh Green’s function theory [71,72,73] in this article. As we will see later, this formalism naturally allows the application of TMA also to the non-equilibrium case.
We introduce the following \(2\times 2\) matrix single-particle Green’s function:
Here, the superscripts “R”, “A”, and “K” represent the retarded, advanced, and Keldysh components, respectively. These are respectively, defined by [71,72,73]
where \(a_{\boldsymbol{p}, \sigma }(t)\) is the annihilation operator of a Fermi atom in the Heisenberg representation, \([A, B]_\pm = AB \pm BA\), \(\Theta (t)\) is the step function, and the expectation value \(\langle {\cdots }\rangle\) is taken over the BCS Hamiltonian H in Eq. (2). We find from Eqs. (5) and (6) that the retarded Green’s function is related to the advanced one as
In addition, in the thermal equilibrium state, the Keldysh Green’s function is related to the retarded one via the fluctuation-dissipation relation (FDR) [71,72,73] as
Here, \(f(\omega ) = [e^{\omega /T}+1]^{-1}\) is the Fermi-Dirac distribution function. Equations (8) and (9) mean that in the thermal equilibrium case, once we compute the retarded component \(G^{\text{R}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\), the other ones in \(\hat{G}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) are immediately obtained from \(G^{\text{R}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\). However, this is not the case when the system is out of equilibrium. In this case, since the FDR in Eq. (9) no longer holds, the retarded and the Keldysh Green’s functions have to be evaluated independently, as we will explain in Section 4.1.
The \(2\times 2\) matrix Green’s function \(\hat{G}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (4) obey the Dyson equation [71,72,73],
which is diagrammatically drawn as Fig. 4a. Here,
is the bare Green’s function in the case of a non-interacting Fermi gas, where \(\delta\) is an infinitesimally small positive number.
a Dyson equation for \(2\times 2\) matrix TMA Green’s function \(\hat{G}_{\text{TMA}, \sigma }\) (thick solid line). \(\hat{\Sigma }_{\text{TMA}, \sigma }\) is the self-energy correction in TMA. b Particle-particle scattering matrix \(\hat{\Gamma }_0\) in TMA. The wavy line represents the pairing interaction \(-U\). c Illustration of pairing fluctuations described by \(\hat{\Gamma }_0\)
In the Dyson Eq. (10), the self-energy
describes effects of the strong paring interaction \(-U\). The self-energy \(\hat{\Sigma }_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) in TMA is diagrammatically drawn as Fig. 4a, which gives [37, 38]
Here, “\(-\sigma\)” means the opposite component to \(\sigma\). In Eq. (13), \(\Gamma ^{\text{R (K)}}_0(\boldsymbol{q}, \nu )\) is the retarded (Keldysh) component of the particle-particle scattering matrix,
which is given by the series of the ladder-type diagrams shown in Fig. 4b. These ladder diagrams physically describe pairing fluctuations, that is, the sequence of pair formation and dissociation of Fermi atoms, as schematically shown in Fig. 4c. The summation of the ladder diagrams in Fig. 4b gives
where \(n_{\text{B}}(\omega ) = [e^{\omega /T}-1]^{-1}\) is the Bose distribution function. In Eq. (16),
is the lowest order pair correlation function [37, 38]. In obtaining the second line in Eq. (18), we have usesd Eq. (11).
The Dyson Eq. (10) gives the dressed retarded (advanced) Green’s function \(G^{\text{R(A)}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) as
The Keldysh component \(G^{\text{K}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) is then immediately obtained from the FDR in Eq. (9).
In TMA, physical quantities are obtained from the dressed Green’s function \(\hat{G}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (4). The total number N of Fermi atoms is evaluated from the Keldysh component \(G^{\text{K}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\). Noting the definition of the Keldysh component in Eq. (7), we have
The substitution of Eq. (9) into Eq. (20) gives
Here,
are the single-particle spectral function and the density of states, respectively [71,72,73].
One sees from the expression for the spectral function in Eq. (22) how the pairing interaction affects single-particle excitations. The pairing interaction gives rise to the energy shift \(\text{Re}\Sigma ^{\text{R}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\), as well as the broadening of the linewidth \(\text{Im}\Sigma ^{\text{R}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\). We will explain these effects in more detail in Section 2.2.
In the TMA scheme, the strong-coupling effects associated with the pairing interaction are incorporated into the theory by solving the number Eq. (21) to determine the chemical potential \(\mu\) for a given number N of Fermi atoms and the temperature T (\(\ge T_{\text{c}}\)). The superfluid phase transition temperature \(T_{\text{c}}\) is determined by solving the number Eq. (21), together with the Thouless criterion [74]. As shown by Kadanoff and Martin [75, 76], the system experiences superfluid instability, when the retarded particle-particle scattering matrix \(\Gamma ^{\text{R}}_0(\boldsymbol{q}, \nu )\) in Eq. (16) has a pole at \(\boldsymbol{q}=\boldsymbol{q}_{\text{pair}}\) and \(\nu = \mu _{\text{pair}}\), that is,
We note that the momentum \(\boldsymbol{q}_{\text{pair}}\) and the energy \(\mu _{\text{pair}}\) physically describe the center-of-mass momentum and the energy of a Cooper pair, respectively. Since \(\Gamma ^{\text{R}}_0(\boldsymbol{q}, \nu )\) in Eq. (16) is a complex function, Eq. (24) actually consists of two equations,
The latter equation is solved analytically. Substituting \(\Pi ^{\text{R}}_0(\boldsymbol{q}, \nu )\) in Eq. (18) into \(\Gamma ^{\text{R}}_0(\boldsymbol{q}, \nu )\) in Eq. (16), one has
Since \(f(x)+f(-x)=1\), Eq. (27) is satisfied only when \(\mu _{\text{pair}}=0\). Substituting this into Eq. (2526), we have
which is just the well-known Thouless criterion (or the gap equation) [74]. In Eq. (28), \(\boldsymbol{q}_{\text{pair}}\) is chosen so as to obtain the highest \(T_{\text{c}}\). In the thermal equilibrium and spin-balanced (\(N_\uparrow =N_\downarrow\)) case, we obtain \(\boldsymbol{q}_{\text{pair}}=0\) because Cooper pairs must have zero center-of-mass momentum in this case. Setting \(\boldsymbol{q}_{\text{pair}}=0\) in Eq. (28), we have
Figure 5a shows \(T_{\text{c}}\) obtained by solving the TMA coupled Eqs. (21) and (29). For comparison, we also show the results in the case when the Thouless criterion in Eq. (29) is solved for the fixed value of the chemical potential \(\mu =\varepsilon _{\text{F}}\). The latter is just the mean-field approximation, which is also referred to as the Kadanoff-Martin (KM) theory in the literature [75, 76]. Figure 5a indicates that the behavior of the calculated \(T_{\text{c}}\) in TMA agrees well with the experimental result shown in Fig. 2. That is, starting from the weak-coupling BCS regime, \(T_c\) gradually increases with increasing the interaction strength, to approach a constant value in the BEC regime. Since all Fermi atoms form tightly bound molecules in the extreme BEC limit, this value just equals the BEC phase transition temperature
in an ideal gas of \(N_{\text{M}}=N/2\) bosons with the molecular mass \(M_{\text{M}}=2m\). In Eq. (30), \(\zeta (3/2) \simeq 2.612\) is the Riemann zeta function.
a Calculated \(T_{\text{c}}\) in a thermal equilibrium Fermi gas in the BCS-BEC crossover region. “TMA” (solid line) and “KM” (dashed line), respectively, show the results of TMA and the mean-field KM theory. b–d Single-particle properties in the BCS-BEC crossover region. (b1)–(d1) Density of states \(\rho _{\text{TMA},\sigma }(\omega )\) in Eq. (23). (b2)–(d2) Spectral function \(A_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (22). Each panel shows the result at b–d in panel a. In panel (b2), “hole” is the peak line along the hole dispersion \(\omega = -\xi _{\boldsymbol{p}}\)
On the other hand, the mean-field KM theory cannot properly describe the behavior of \(T_{\text{c}}\) in the strong-coupling BEC regime [see the dashed line in Fig. 5a]. This is simply because the KM theory ignores pairing fluctuations, as well as the formation of diatomic molecules in this regime.
2.2 Single-particle properties in the thermal equilibrium BCS-BEC crossover regime
In the BCS-BEC crossover region, a characteristic strong-coupling phenomenon appearing in the single-particle excitations is the pseudogap, where fluctuations in the Cooper channel play an essential role [66, 57, 68, 77,78,79]. To explain this phenomenon, we first recall the non-interacting case, where the single-particle spectral function \(A_{0, \sigma }(\boldsymbol{p}, \omega )\) and the density of states \(\rho _{0, \sigma }(\omega )\) in the non-interacting Fermi gas are, respectively, obtained from \(G^{\text{R}}_{0, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (11) as
As seen from these expressions, the spectral function \(A_{0, \sigma }(\boldsymbol{p}, \omega )\) has the peak line along the particle dispersion \(\omega =\xi _{\boldsymbol{p}}\). The single-particle density of states \(\rho _{0, \sigma }(\omega )\) is a monotonically increasing function of \(\omega\). These simple structures are remarkably modified in the BCS-BEC crossover region due to the presence of strong pairing fluctuations, as shown in Fig. 5b–d. Although the superfluid order parameter vanishes at \(T_{\text{c}}\), Fig. 5(b1) shows that the single-particle density of states \(\rho _{\text{TMA}, \sigma }(\omega )\) still exhibits a dip structure around the Fermi level (\(\omega =0\)). This so-called pseudogap structure [66, 68, 77, 78, 57, 79] becomes more remarkable, as one passes through the BCS-BEC crossover region [see Fig. 5(c1) and 5(d1)].
To quickly understand the role of pairing fluctuations in the pseudogap phenomenon, we approximate the retarded self-energy \(\Sigma ^{\text{R}}_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (1314) in the following manner. Noting that \(\Gamma ^{\text{R, K}}_0(\boldsymbol{q}=0, \nu =0)\) diverges when the Thouless criterion in Eq. (24) is satisfied at \(T_{\text{c}}\), one approximates the self-energy to
Here,
describes effects of pairing fluctuations, which is also referred to as the pseudogap parameter in the literature [68]. Substituting Eq. (33) into the Dyson Eq. (10), we have
Equation (35) indicates that pairing fluctuations described by the pseudogap parameter \(\Delta _{\text{PG}}\) brings about a coupling phenomenon between the particle band \((\omega =\xi _{\boldsymbol{p}})\) and the hole band \((\omega = -\xi _{\boldsymbol{p}})\). Evaluating the single-particle spectral function in Eq. (22) in this so-called pseudogap approximation, one obtains
Equation (36) just has the same form as the single-particle spectral function in the mean-field BCS theory where the superfluid order parameter is replaced by the pseudogap parameter \(\Delta _{\text{PG}}\). Thus, Eq. (36) has the excitation gap \(\Delta _{\text{PG}}\) around \(\omega =0\). From the viewpoint of the above-mentioned particle-hole coupling by pairing fluctuations, the level repulsion between the particle- and hole-band around \(\omega =0\) opens this energy gap [see also Fig. 5(b1)]. While the pseudogap approximation gives a clear single-particle excitation gap as in the BCS superfluid state, quasiparticle lifetime associated with particle-particle scatterings actually rounds this gap structure to some extent. Because of this, this phenomenon appears as a dip in the single-particle density of states, as seen in Fig. 5(b1).
The pseudogap develops with increasing the interaction strength, reflecting the enhancement of pairing fluctuations, as shown in Fig. 5(c1) and (c2). In the strong-coupling BEC regime, the density of states \(\rho _{\text{TMA}, \sigma }(\omega )\), as well as the spectral function \(A_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) has a clear energy gap, rather than the pseudogap. These gapped excitation spectra reflect the formation of tightly-bounded diatomic molecules [68]. The energy gap corresponds to the binding energy \(2|\mu |= 1/(m a_{\text{s}})^2\) of a two-body bound molecule in the strong-coupling BEC limit.
In the current stage of cold Fermi gas physics, it is still difficult to directly observe the single-particle spectral function, as well as the density of states. Regarding this, however, we note that the photoemission spectrum (PES) \(L_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) is observable [80,81,82], which is related to the spectral function \(A_{\text{TMA}, \sigma }(\boldsymbol{p}, \omega )\) as [82]
where \(f(\omega )\) is the Fermi-Dirac distribution function. In a sense, PES may be viewed as the occupied spectral function. Thus, the structural changes of the spectral function in the BCS-BEC crossover region can be observed through the photoemission-type experiment [81, 83].
3 Non-equilibrium superfluid transition in the driven-dissipative Fermi gas: mean-field approach
We now proceed to the non-equilibrium case. In Section 3.1, we present the model Hamiltonian for the driven-dissipative Fermi gas in Fig. 3a. In Section 3.2, we explain how the couplings with the reservoirs affect the single-particle properties of the main system. To apply the mean-field KM theory to the driven-dissipative Fermi gas, we extend the Thouless criterion to the non-equilibrium case in Section 3.3. Using the non-equilibrium Thouless criterion, we study the non-equilibrium superfluid phase transition within the mean-field approximation in Section 3.4.
3.1 Model Hamiltonian for the driven-dissipative Fermi gas
The model driven-dissipative Fermi gas in Fig. 3a is described by the Hamiltonian [37,38,39]
Here,
describes the main system (non-equilibrium interacting Fermi gas) in Fig. 3a, where \(a^\dagger _{\boldsymbol{p}, \sigma }\) is the creation operator of a Fermi atom with pseudospin \(\sigma =\uparrow , \downarrow\) in the main system. \(\varepsilon _{\boldsymbol{p}}=\boldsymbol{p}^2/(2m)\) is the kinetic energy of a Fermi atom with an atomic mass m. \(-U\) is a tunable pairing interaction associated with a Feshbach resonance. As in the thermal equilibrium case discussed in Section 2.1, we measure the interaction strength in terms of the s-wave scattering length \(a_s\), which is related to \(-U\) via Eq. (3).
The left (\(\alpha =\text{L}\)) and right (\(\alpha =\text{R}\)) reservoirs in Fig. 3a are described by
Here, \(c^{\alpha }_{\boldsymbol{q}, \sigma }\) is an anhiration operator of a Fermi atom in the \(\alpha\) reservoir, and \(\xi ^{\text{res}}_{\boldsymbol{p}}= \varepsilon _{\boldsymbol{p}} -\mu _{\text{res}}\) is the single-particle energy in the reservoirs [see Fig. 3b]. We assume that the reservoirs are huge compared to the main system and always in the thermal equilibrium state at the common (environment) temperature \(T_{\text{env}}\). Under this assumption, the Fermi atoms in the reservoirs obey the ordinary Fermi-Dirac distribution function,
The tunnelings of the Fermi atoms between the main system and reservoirs are described by the tunneling Hamiltonian \(H_{\text{mix}}\), given by
where \(\mathcal {T}_{\alpha , \boldsymbol{q}, \boldsymbol{p}}\) is a tunneling matrix element between the main system and the \(\alpha\) (\(=\)L, R) reservoir. For simplicity, we ignore the momentum and \(\alpha\) dependence of the tunneling matrix element, to set \(\mathcal {T}_{\text{L}, \boldsymbol{q}, \boldsymbol{p}}=\mathcal {T}_{\text{R}, \boldsymbol{q}, \boldsymbol{p}} \equiv \mathcal {T}\). In the tunneling Hamiltonian in Eq. (42), the atom tunneling is assumed to occur between randomly distributing spatial positions \(\boldsymbol{R}_i^\alpha\) in the \(\alpha\) reservoir and \(\boldsymbol{r}_i^\alpha\) in the main system (\(i=1,\cdots , N_{\text{t}} \gg 1\)). The introduction of the random tunneling points is just a simple theoretical trick to study the bulk properties of the main system [37,38,39,40]. After taking the spatial average over the randomly distributing tunneling positions, the translational invariance of the main system is recovered.
The factor \(e^{i\mu _\alpha t}\) in Eq. (42) describes the situation where the energy band \(\xi _{\boldsymbol{q}}\) in the \(\alpha\) reservoir is occupied up to their respective Fermi levels \(\mu _\alpha =\mu \pm \Delta \mu\) when \(T_{\text{env}}=0\), as schematically shown in Fig. 3b. Due to this factor, Fermi atoms are injected into and extracted from the main system when \(\mu _{\text{L}} \ne \mu _{\text{R}}\).
We note that the temperature in the main system is not well defined in the non-equilibrium case when \(\Delta \mu \ne 0\). In this case, the superfluid instability in the main system is controlled by the temperature \(T_{\text{env}}\) in the thermal equilibrium reservoirs. To emphasize this, in what follows, we write the superfluid transition temperature in the driven-dissipative Fermi gas as \(T_{\text{env}}^{\text{c}}\).
3.2 Non-equilibrium properties of the driven-dissipative non-interacting Fermi gas
When the main system is in the non-equilibrium steady state due to the couplings with the reservoirs, the Fermi atoms in the main system obey a non-equilibrium energy distribution \(f_{\text{neq}}(\omega )\) that has a different structure from the ordinary Fermi-Dirac distribution function. To obtain the distribution \(f_{\text{neq}}(\omega )\) in the absence of pairing interaction (\(U=0\)), we conveniently introduce the \(2\times 2\) matrix non-equilibrium Green’s function in the main system, given by
which obeys the Dyson equation,
The Dyson Eq. (44) is diagrammatically drawn as Fig. 6. In Eq. (44), \(\hat{G}_{0, \sigma }(\boldsymbol{p}, \omega )\) is the bare Green’s function in the absence of the system-reservoir couplings, given in Eq. (11). We emphasize that the distribution function \(f(\omega )\) in \(\hat{G}_{0, \sigma }(\boldsymbol{p}, \omega )\) has nothing to do with the distribution function \(f_{\text{env}}(\omega )\) in the reservoirs in Eq. (41). We can regard \(f(\omega )\) in \(\hat{G}_{0, \sigma }(\boldsymbol{p}, \omega )\) as the energy distribution in the isolated Fermi gas, before the main system is connected to the reservoirs. As we will see below, the non-equilibrium distribution \(f_{\text{neq}}(\omega )\) in the main system does not depend on \(f(\omega )\), when the system relaxes into the non-equilibrium steady state due to the couplings with the reservoirs.
Dyson equation for \(2\times 2\) matrix non-equilibrium Green’s function \(\hat{G}_{\text{neq}, \sigma }\) (double solid line). The self-energy \(\hat{\Sigma }_{\text{env}, \sigma }\) describes the system-reservoir coupling effects in the second-order Born approximation with respect to the tunneling amplitude \(\mathcal {T}\) (solid square). \(\hat{D}_{\alpha , \sigma }\) (dotted line) denotes the non-interacting Green’s function in the \(\alpha\) (\(=\)L, R) reservoir
In Eq. (44), the self-energy \(\hat{\Sigma }_{\text{env}, \sigma }(\boldsymbol{p}, \omega )\) describes effects of system-reservoir couplings. Within the second-order Born approximation with respect to the tunneling matrix \(\mathcal {T}\), the self-energy is diagrammatically given in Fig. 6. Evaluating this diagrammatic equation, we obtain [38, 39]
For the derivation of Eq. (45), we refer to Ref. [38]. Here,
is the non-interacting \(2\times 2\) matrix Green’s function in the \(\alpha\) (\(=\) L, R) reservoir, and
describes the quasiparticle damping originating from the system-reservoir coupling. In Eq. (47), \(\rho =\rho _\alpha\) is the single-particle density of states in the \(\alpha\) reservoir, where we have ignored the \(\alpha\) dependence of this quantity, for simplicity. We have also ignored the \(\omega\) dependence of \(\rho\), which is sometimes referred to as the wide-band limit approximation in the literature [73]. This approximation is justified in the case when the reservoirs are so huge that the energy dependence of their density of states around \(\omega =\mu\) can be ignored, within the variation of \(\Delta \mu\).
The substitution of Eq. (45) into the Dyson Eq. (44) gives
Here, the Fermi-Dirac function \(f_{\text{env}}(\omega )\) in the reservoirs is given in Eq. (41). As mentioned previously, while the non-equilibrium steady-state Green’s function \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) depends on \(f_{\text{env}}(\omega )\), it is not affected by the distribution function \(f(\omega )\) in the initial state of the main system. This is because the initial memory of the isolated Fermi gas is wiped out due to the coupling with the reservoirs.
The energy distribution function \(f_{\text{neq}}(\omega )\) in the main system is obtained from the Keldysh component of \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (48). Although the FDR, like Eq. (9), does not hold between the retarded and the Keldysh components in the non-equilibrium state (\(\Delta \mu \ne 0\)), they still possess a similar relation. The Keldysh component in Eq. (48) can be written as
where
By analogy with the FDR in the thermal equilibrium case, we can interpret \(f_{\text{neq}}(\omega )\) in Eq. (50) as the non-equilibrium energy distribution function in the driven-dissipative Fermi gas. The distribution function \(f_{\text{neq}}(\omega )\) has a clear two-step structure (at \(\omega =\mu _{\text{L}}\) and \(\omega =\mu _{\text{R}}\)) when \(T_{\text{env}} \ll \Delta \mu\), originating from the different Fermi levels between the left and right reservoirs. We note that such a non-equilibrium energy distribution with the two-step structure has experimentally been observed in mesoscopic wires [84,85,86], as well as carbon nanotubes [87], under a bias voltage V (which corresponds to the chemical potential bias \(\Delta \mu\) in the present model driven-dissipative Fermi gas).
The couplings with the reservoirs affect, not only the energy distribution, but also the spectral function in the main system. The single-particle spectral function \(A_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) is obtained from the retarded component of \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) as
A comparison of \(A_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (51) and \(A_{0, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (31) shows that the system-reservoir couplings give rise to the linewidth \(2\gamma\) in the excitation spectrum. The linewidth physically means the finite lifetime \(\tau\) \((\sim 1/\gamma )\) of the excitation mode \(\varepsilon _{\boldsymbol{p}}\) in the driven-dissipative Fermi gas, due to the atom tunneling between the main system and the reservoirs.
The momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})=\langle {a^\dagger _{\boldsymbol{p}, \sigma } a_{\boldsymbol{p}, \sigma }\rangle }\) is obtained from the Keldysh component \(G^{\text{K}}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) as
Equation (52) indicates that \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) involves information about the single-particle spectral function \(A_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (51), as well as the energy distribution function \(f_{\text{neq}, \sigma }(\omega )\) in Eq. (50).
Figure 7 shows the calculated \(n_{\text{neq}, \sigma }(\boldsymbol{p})\). In the zero bias case \((\Delta \mu =0)\) shown in Fig. 7a, the effects of system-reservoir couplings are dominated by the quasi-particle damping described by \(\gamma\) in Eq. (47). Since the quasi-particle peak in the single-particle spectral function \(A_{\text{neq},\sigma }(\boldsymbol{p}, \omega )\) is broadened by this damping, the smearing of the sharp Fermi edge at \(|\boldsymbol{p}|\simeq p_{\text{F}}\) becomes more remarkably with increasing the value of \(\gamma\), as shown in Fig. 7a.
Calculated non-equilibrium momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) as a function of \(|\boldsymbol{p}|\) in the driven-dissipative Fermi gas. Panels a and c show the effects of the system-reservoir coupling strength \(\gamma\) in the presence and the absence of the chemical potential bias \(\Delta \mu\), respectively. Panel b shows effects of \(\Delta \mu\), when \(\gamma =0.01\mu\). We set \(T_{\text{env}}=0\)
Once non-zero bias \((\Delta \mu \ne 0)\) is imposed, one finds in Fig. 7b that the momentum distribution \(n_{\text{neq},\sigma }(\boldsymbol{p})\) exhibits the two-step structure expected in Eq. (50). However, as shown in Fig. 7c, the two-step structure imprinted on the momentum distribution \(n_{\text{neq},\sigma }(\boldsymbol{p})\) becomes obscure as \(\gamma\) increases due to the broadening of the quasi-particle peak in \(A_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\). Although we do not show it explicitly, the two-step structure is also rounded by thermal excitations as the temperature \(T_{\text{env}}\) increases. Because of these, the conditions for the two-step structure to appear in the momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) are \(\Delta \mu \gg \gamma\) and \(\Delta \mu \gg T_{\text{env}}\).
3.3 Thouless criterion for the driven-dissipative Fermi gas
To extend the Thouless criterion in Eq. (24) to the non-equilibrium steady state, we evaluate the particle-particle scattering matrix \(\hat{\Gamma }_{\text{neq}}(\boldsymbol{q}, \nu )\) in the driven-dissipative Fermi gas. As in the thermal equilibrium case, \(\hat{\Gamma }_{\text{neq}}(\boldsymbol{q}, \nu )\) is also given by the series of the ladder diagrams shown in Fig. 8a. While \(\hat{\Gamma }_0(\boldsymbol{q}, \nu )\) involves the bare Green’s function \(\hat{G}_{0, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (11), \(\hat{\Gamma }_{\text{neq}}(\boldsymbol{q}, \nu )\) involves the non-equilibrium Green’s function \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (48), to take into account the system-reservoir coupling effects. In this sense, \(\hat{\Gamma }_{\text{neq}}(\boldsymbol{q}, \nu )\) physically describes pairing fluctuations in the presence of pumping and decay of Fermi atoms, as schematically shown in Fig. 8b.
a Particle-particle scattering matrix \(\hat{\Gamma }_{\text{neq}}\) in the driven-dissipative Fermi gas. The non-equilibrium Green’s function \(\hat{G}_{\text{neq}, \sigma }\) (double solid line) is diagrammatically given in Fig. 6. \(\hat{\Gamma }_{\text{neq}}\) physically describes paring fluctuations in the presence of pumping and decay of Fermi atoms, as schematically shown in panel (b)
Summing up the ladder diagrams in Fig. 8a, we have [37, 38]
Here,
is the lowest-order pair correlation function in the driven-dissipative Fermi gas. The retarded (advanced) component \(\Pi ^{\text{R(A)}}_{\text{neq}}(\boldsymbol{q}, \nu )\) is obtained by simply replacing \(G^{\text{R, A, K}}_{0, \sigma }(\boldsymbol{p}, \omega )\) with \(G^{\text{R, A, K}}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (18), which reads
In the thermal equilibrium state, the Keldysh component \(\Pi ^{\text{K}}_{\text{neq}}(\boldsymbol{q}, \nu )\) is immediately obtained from the retarded component via the FDR. However, this is not the case for the driven-dissipative Fermi gas. In this non-equilibrium case, the Keldysh component \(\Pi ^{\text{K}}_{\text{neq}}(\boldsymbol{q}, \nu )\) is independently evaluated from the following expression [37, 38]:
As in the thermal equilibrium case, the superfluid phase transition temperature \(T^{\text{c}}_{\text{env}}\) is determined from the retarded particle-particle scattering matrix \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu )\). The non-equilibrium version of the pole Eq. (24) is given by
which actually consists of the following two equations:
In the symmetric coupling case (\(\mathcal {T}_{\text{L}}=\mathcal {T}_{\text{R}}\equiv T\)) we are considering, the latter equation can be solved analytically. When \(\Pi ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu )\) in Eq. (55) is substituted into the retarded component \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu )\) in Eq. (53), Eq. (59) is reduced to
which is satisfied only when \(\mu _{\text{pair}}=2\mu\). We note that in the thermal equilibrium case, the imaginary part of the pole Eq. (24) is satisfied when \(\mu _{\text{pair}}=0\) rather than \(\mu _{\text{pair}}=2\mu\). The difference comes from the different origins of energy in the thermal equilibrium and the present non-equilibrium case. Since we now measure the energy from the bottom of the energy band (\(\varepsilon _{\boldsymbol{p}=0}\)) in the main system [see Fig. 3b], a Cooper pair should have non-zero energy \(\mu _{\text{pair}}=2\mu\). [The factor 2 comes from the fact that a Cooper pair is formed by two Fermi atoms.]
Substituting this solution \(\mu _{\text{pair}}=2\mu\) into Eq. (58), we have
which is the non-equilibrium version of Eq. (28), and is referred to as the non-equilibrium Thouless criterion in what follows. In Eq. (61), \(\boldsymbol{q}_{\text{pair}}\) is determined so as to obtain the highest \(T_{\text{env}}^{\text{c}}\). As noted in Section 2.1, \(\boldsymbol{q}_{\text{pair}}\) physically describes the center-of-mass momentum of a Cooper pair. Thus, when \(\boldsymbol{q}_{\text{pair}}=0\), the main system transitions to the uniform non-equilibrium BCS-type superfluid state, where the Cooper pairs have zero center-of-mass momentum. On the other hand, when \(\boldsymbol{q}_{\text{pair}}\ne 0\), an inhomogeneous superfluid state is realized, being characterized by a spatially oscillating superfluid order parameter (symbolically written as \(\Delta (\boldsymbol{r})= \Delta e^{i\boldsymbol{q}_{\text{pair}}\cdot \boldsymbol{r}}\)). This is analogous to the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state in a superconductor under an external magnetic field [88,89,90].
3.4 Non-equilibrium FFLO superfluid transition in the mean-field approximation
Figure 9a shows the mean-field results of the superfluid phase transition temperature as a function of the system-reservoir coupling strength \(\gamma\) in the zero bias case (\(\Delta \mu =0\)). This figure is obtained by solving Eq. (61) with a fixed value of \(\mu\). (This is an extension of the mean-field KM theory explained in Section 2.1 to the driven-dissipative Fermi gas.) As in the well-known thermal equilibrium case, the superfluid phase is suppressed by thermal fluctuations as the temperature \(T_{\text{env}}\) increases. The superfluid phase is also suppressed as the system-reservoir coupling strength \(\gamma\) increases. This is simply because, as shown in Fig. 7a, the system-reservoir coupling rounds the Fermi edge in the momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) so that \(\gamma\) plays a similar role to the temperature.
a, c Superfluid transition line obtained by solving the non-equilibrium Thouless criterion in Eq. (61) in the weak-coupling BCS regime when \((p_{\text{F}}a_s)^{-1}=-1\). The system experiences the BCS-type (FFLO-type) superfluid instability on the solid (dashed) line. b, d Calculated inverse particle-particle scattering matrix \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu =2\mu )\) as a function of \(|\boldsymbol{q}|\). We note that \(\text{Im}\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu =2\mu )=0\). Each panel shows the results along the paths depicted in panels a and c. The pole of \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu =2\mu )\) (solid circle) is obtained at \(|\boldsymbol{q}|\ne 0\) (\(|\boldsymbol{q}|= 0\)) on the non-equilibrium FFLO (BCS) transition line
In the zero bias case \((\Delta \mu =0)\) shown in Fig. 9a, the pole of \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu )\) always appears at \(\boldsymbol{q}=0\), when \(T_{\text{env}}=T_{\text{env}}^{\text{c}}\) [see Fig. 9b]. As explained previously, this means that the occurrence of the BCS-type uniform superfluid transition at this temperature, being accompanied by BEC of Cooper pairs with zero center-of-mass momentum.
The situation changes when the main system is in the non-equilibrium steady state by the non-zero chemical potential bias \(\Delta \mu \ne 0\). Figure 9c is the phase diagram of the driven-dissipative Fermi gas in terms of the system-reservoir coupling strength \(\gamma\) and the chemical potential bias \(\Delta \mu\), when \(T_{\text{env}}=0\). When \(\Delta \mu\) is very large, the system is in the normal state due to the strong pumping and decay of Fermi atoms. As \(\Delta \mu\) decreases and the non-equilibrium effects are weakened, the system experiences superfluid instability. As shown in Fig. 9(d1), while the pole of \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu )\) appears at \(\boldsymbol{q}=0\) when \(\gamma / \mu \gtrsim 0.04\), it appears at \(|\boldsymbol{q}|=|\boldsymbol{q}_{\text{pair}}|\) (\(> 0\)) when \(\gamma / \mu \lesssim 0.04\). While the former case (\(\boldsymbol{q}=0\)) means the BCS-type superfluid transition at \(T_{\text{env}}^{\text{c}}\), the latter indicates the FFLO-type superfluid transition associated with Cooper pairs with non-zero center-of-mass momentum \(\boldsymbol{q}_{\text{pair}}\ne 0\).
The FFLO state was first proposed in metallic superconductivity under an external magnetic field [88,89,90], and was later discussed in spin-imbalanced (\(N_\uparrow \ne N_\downarrow\)) ultracold Fermi gas [91,92,93,94,95], as well as color superconductivity in quantum chromodynamics [96]. In these systems, the FFLO state is induced by the misalignment of the Fermi surfaces. For example, as schematically shown in Fig. 10a, the Fermi momenta \(p_{\text{F}\sigma }\) of the pseudospin \(\sigma =\uparrow\) and \(\sigma = \downarrow\) Fermi atoms are different from each other in a Fermi gas with the spin imbalance \(N_\uparrow \ne N_\downarrow\). Since Cooper pairs are formed by Fermi atoms near each Fermi surface, Cooper pairs acquire non-zero center-of-mass momentum \(|\boldsymbol{q}_{\text{pair}}| \simeq p_{\text{F}\uparrow } -p_{\text{F}\downarrow }\) due to the Fermi surface mismatch, as schematically shown in Fig. 10b.
Cooper-pair formation in the FFLO state in a thermal equilibrium Fermi gas with spin imbalance (\(N_\uparrow \ne N_\downarrow\)). a Momentum distribution \(n_{\sigma =\uparrow , \downarrow }(\boldsymbol{p})\). Since the system is in the thermal equilibrium state, \(n_{\sigma }(\boldsymbol{p})\) is simply given by the Fermi-Dirac distribution function. However, the spin imbalance leads to \(p_{\text{F}\uparrow }\ne p_{\text{F}\downarrow }\) (where \(p_{\text{F}\sigma }\) is the Fermi momentum in the pseudospin \(\sigma =\uparrow , \downarrow\) component). b The misalignment between \(\sigma =\uparrow\) and \(\downarrow\) Fermi surfaces naturally induces a Cooper pair with non-zero center-of-mass momentum, \(|\boldsymbol{q}_{\text{pair}}| \simeq p_{\text{F} \uparrow } -p_{\text{F} \downarrow }\ne 0\)
However, we emphasize that the present driven-dissipative Fermi gas is a spin-balanced (\(N_\uparrow =N_\downarrow\)) system, so that there is no Fermi surface misalignment in the main system. In this sense, the mechanism of the FFLO-type superfluid in Fig. 9c is different from the conventional one discussed in metallic superconductivity under an external magnetic field.
The origin of this non-equilibrium FFLO superfluid is the non-equilibrium momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) in Eq. (52) [39, 97]. When the momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) has the clear two-step structure [see Fig. 7b], the occupation sharply changes at \(p_{\text{F}1}=\sqrt{2m[\mu +\Delta \mu ]}\) and at \(p_{\text{F}2}=\sqrt{2m[\mu -\Delta \mu ]}\), as schematically shown in Fig. 11a. Then if each edge structure at \(p=p_{\text{F}1}\) and \(p=p_{\text{F}2}\) work like the ordinary Fermi edge at the Fermi surface, each pseudospin component (\(\sigma =\uparrow , \downarrow\)) has two “Fermi surfaces” with different sizes. Then, as schematically shown in Fig. 11b and c, several types of Cooper pairings between Fermi atoms around the “Fermi surfaces” would be possible. Among them, the Cooper pairs depicted in Fig. 11c have non-zero center-of-mass momentum \(|\boldsymbol{q}_{\text{pair}}| \simeq p_{\text{F}1}-p_{\text{F}2}\) \((\ne 0)\), which result in the non-equilibrium FFLO superfluid.
Same figure as Fig. 10 for the driven-dissipative Fermi gas. a The non-equilibrium momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) in Eq. (52). Since the main system is the spin-balanced driven-dissipative Fermi gas, both pseudospin components (\(\sigma =\uparrow , \downarrow\)) obey the same momentum distribution, \(n_{\text{neq}, \uparrow }(\boldsymbol{p})=n_{\text{neq}, \downarrow }(\boldsymbol{p})\). In each component, two Fermi edges imprinted on \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) at \(p=p_{\text{F}1}\) and \(p=p_{\text{F2}}\) work like two “Fermi surfaces” with different sizes. b, c indicate possible Cooper pairings between Fermi atoms around the “Fermi surfaces”. While Cooper pairs have zero center-of-mass momentum in b, they have non-zero center-of-mass momentum in c, which is similar to the FFLO Cooper pair shown in Fig. 10b
4 Non-equilibrium BCS-BEC crossover in the driven-dissipative Fermi gas
In the previous section, we discuss the non-equilibrium superfluid phase transition within the mean-field approximation. In this section, we go beyond the mean-field approximation to study the non-equilibrium BCS-BEC crossover.
4.1 Strong-coupling theory for the driven-dissipative Fermi gas
We extend the thermal-equilibrium TMA explained in Section 2.1 to the case when the system is out of equilibrium. In the non-equilibrium TMA (NETMA) scheme, while \(T_{\text{env}}^{\text{c}}\)-Eq. (61) obtained from the Thouless criterion is still valid, we need to include strong-coupling corrections to the number equation. For this purpose, we introduce the non-equilibrium Green’s function in the main system, given by
This \(2\times 2\) matrix Green’s function obeys the Dyson equation,
which is diagrammatically drawn as Fig. 12. In Eq. (63), \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) is given in Eq. (48). (Note that it already involves the effects of system-reservoir couplings.) The self-energy \(\hat{\Sigma }_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) describes the effects of pairing fluctuations within the TMA level, which is diagrammatically given by Fig. 12. Evaluating the last term in Fig. 12, we obtain [37, 38]
where the particle-particle scattering matrix \(\Gamma ^{\text{R, A, K}}_{\text{neq}}(\boldsymbol{q}, \nu )\) is given in Eq. (53). Substituting Eqs. (64) and (65) into the Dyson Eq. (63), we obtain each component in Eq. (62) as
Dyson equation for NETMA Green’s function \(\hat{G}_{\text{NETMA}, \sigma }\) (thick solid line). The non-equilibrium Green’s function \(\hat{G}_{\text{neq}, \sigma }\) (double solid line) is given in Fig. 6. The self-energy \(\hat{\Sigma }_{\text{NETMA}, \sigma }\) describes the non-equilibrium pairing fluctuation effects. The particle-particle scattering matrix \(\hat{\Gamma }_{\text{neq}}\) is given in Fig. 8a
The number N of Fermi atoms in the main system is computed from the Keldysh component \(G^{\text{K}}_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) as
Here, we have introduced the lesser self-energy [71,72,73],
The number Eq. (68) now includes both the strong-coupling effects associated with the tunable pairing interaction, as well as the system-reservoir coupling effects. As in the thermal equilibrium BCS-BEC crossover theory explained in Section 2.1, we solve the number Eq. (68), together with the non-equilibrium Thouless criterion in Eq. (61), to determine \(\mu\) and \(T_{\text{env}}^{\text{c}}\) for a given parameter set (N, \((p_{\text{F}} a_s)^{-1}\), \(\Delta \mu\)). When \(T_{\text{env}} > T_{\text{env}}^{\text{c}}\), we only solve the number Eq. (68), to determine \(\mu\).
Once \(\mu\) is determined, the single-particle spectral function \(A_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\), the density of states \(\rho _{\text{NETMA}, \sigma }(\omega )\), as well as the photoemission spectrum \(L_{\text{NETMA}, \sigma }(\omega )\) are obtained from the NETMA Green’s function as, respectively, [38]
4.2 Non-equilibrium BCS-BEC crossover
Figure 13 shows the superfluid phase transition temperature \(T_{\text{env}}^{\text{c}}\) as a function of the chemical potential bias \(\Delta \mu\), and the pairing interaction strength \((p_{\text{F}}a_s)^{-1}\), in the non-equilibrium BCS-BEC crossover region obtained by NETMA. In the zero bias case \((\Delta \mu =0)\), the overall behavior of \(T_{\text{env}}^{\text{c}}\) is the same as \(T_{\text{c}}\) in a thermal equilibrium Fermi gas shown in Fig. 5a. Starting from the weak-coupling BCS regime, \(T_{\text{env}}^{\text{c}}\) gradually increases with increasing the interaction strength, to eventually approach a constant value in the strong-coupling BEC regime.
The main system is in the driven-dissipative steady state when the non-zero bias \(\Delta \mu \ne 0\) is imposed. In this non-equilibrium case, \(T_{\text{env}}^{\text{c}}\) is suppressed, as shown in Fig. 13. In the strong-coupling BEC regime, \((p_{\text{F}} a_s)^{-1}\gtrsim 0\), \(T_{\text{env}}^{\text{c}}\) is monotonically suppressed with increasing \(\Delta \mu\). On the other hand, in the weak-coupling BCS side, \((p_{\text{F}} a_s)^{-1}\lesssim 0\), \(T_{\text{env}}^{\text{c}}\) is found to exhibits re-entrant behavior. That is, with decreasing the temperature \(T_{\text{env}}\), although the main system experiences the superfluid phase transition at \(T_{\text{env}}=T_{\text{env}}^{\text{c}}\), the main system again becomes the normal state at low temperatures [see also Fig. 14a].
a Calculated \(T_{\text{env}}^{\text{c}}\) as a function of \(\Delta \mu\), in the driven-dissipative Fermi gas in the BCS region at \((p_{\text{F}}a_s)^{-1}=-0.6\). \(T_{\text{env}}^{\text{c}}\) exhibits re-entrant behavior in this region. b Single-particle density of states \(\rho _{\text{NETMA}, \sigma }(\omega )\) along the paths (b1)–(b3) at panel a. For visibility, we offset each result by 0.5. c Single-particle spectral function \(A_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) at (c1)–(c5) in panel a. In (c4), the broad peak structure “FFLO” originates from FFLO-type pairing fluctuations. The intensity of the spectral function is normalized by \(\varepsilon _{\text{F}}^{-1}\)
Figure 13 indicates that, in NETMA, the non-equilibrium FFLO state, which is obtained in the non-equilibrium mean-field theory as explained in Section 3.4, does not appear in the whole BCS-BEC crossover region. That is, at \(T_{\text{env}}=T_{\text{env}}^{\text{c}}\), the non-equilibrium Thouless criterion in Eq. (61) is always satisfied at \(\boldsymbol{q}_{\text{pair}}=0\) in NETMA. In general, the strong-coupling theory and the mean-field KM theory are expected to give qualitatively the same results in the weak-coupling BCS regime because pairing fluctuations are sufficiently weak there. Indeed, in the thermal equilibrium case, \(T_{\text{c}}\) obtained by TMA approaches the mean-field result in the weak-coupling limit, as shown in Fig. 5a. Nevertheless, in the present case, the non-equilibrium FFLO state predicted in the mean-field theory is not obtained in NETMA even in the weak-coupling BCS regime.
4.2.1 Weak-coupling BCS regime
We discuss the mechanism of the re-entrant behavior of \(T_{\text{env}}^{\text{c}}\) as well as the absence of the non-equilibrium FFLO state, in the weak-coupling BCS regime. Figure 14a shows \(T_{\text{env}}^{\text{c}}\) as a function of \(\Delta \mu\). In Fig. 14b, c, we show the single-particle density of states \(\rho _{\text{NETMA}, \sigma }(\omega )\) along the paths (b1)–(b3) and the spectral function \(A_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) at (c1)–(c5) in Fig. 14a.
When \(\Delta \mu =0\), the pseudogap appears in the density of states, as shown in Fig. 14(b1). We also find from Fig. 14(c1) that the spectral function has the peak along the particle dispersion (\(\omega =\varepsilon _{\boldsymbol{p}}\)), as well as the hole dispersion (\(\omega = -\varepsilon _{-\boldsymbol{p}} +2\mu\)), which is the typical pseudogap structure explained in Section 2.2. Since the pseudogap phenomenon is induced by pairing fluctuations, it should be enhanced (suppressed) as we approach (move away from) \(T_{\text{env}}^{\text{c}}\). Indeed, we see in Fig. 14(b1) that the pseudogap gradually disappears, as one moves away from \(T_{\text{env}}^{\text{c}}\). At (c2) in Fig. 14a, the peak structure along the hole dispersion disappears in the spectral function, as expected.
However, a very different behavior of the pseudogap is obtained when one moves along the path (b2) in Fig. 14a. In this case, although we first approach \(T_{\text{env}}^{\text{c}}\) and then move away from \(T_{\text{env}}^{\text{c}}\), the pseudogap monotonically develops up to the point (c4) in Fig. 14a, as seen in Fig. 14(b2). This tendency is also seen in the single-particle spectral function. Indeed, Figs. 14(c2)–(c4) show that the level repulsion between the particle and the hole dispersions monotonically becomes remarkable, even when one moves from (c3) to (c4) in Fig. 14a.
To understand the reason why the pseudogap phenomenon is remarkably seen at (c4) in Fig. 14a, we plot in Fig. 15b the real part of the particle-particle scattering matrix \(-\text{Re}\big [\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu =2\mu )\big ]\) as a function of \(|\boldsymbol{q}|\), which informs us of the strength of pairing fluctuations. (Note that \(\boldsymbol{q}\) in this quantity physically has the meaning of the center-of-mass momentum of preformed Cooper pairs.) We find from this figure that, when one moves along the path (b) in Fig. 15a, although pairing fluctuations in the BCS-type Cooper channel at \(\boldsymbol{q}=0\) decreases, FFLO-type pairing fluctuations being characterized by \(\boldsymbol{q}\ne 0\) becomes strong.
a Phase diagram of the driven-dissipative Fermi gas at \((p_{\text{F}}a_s)^{-1}=-0.6\). b and c show the real part of the particle-particle scattering matrix \(-\text{Re}\big [\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu =2\mu )\big ]\) as a function of \(|\boldsymbol{q}|\). Each panel shows the result along the path (b) and path (c) in panel a
Since we are considering the spin-balanced case, the ordinary mechanism of the FFLO state does not work here. Instead, the anomalous enhancement of \(-\text{Re}\big [\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu =2\mu )\big ]\) at \(\boldsymbol{q}\ne 0\) seen in Fig. 15b is attributed to the two-step structure imprinted on the non-equilibrium momentum distribution \(n_{\text{neq},\sigma }({\boldsymbol{p}})\), as explained in Section 3.4. Since the thermal excitations are suppressed and the momentum distribution \(n_{\text{neq},\sigma }({\boldsymbol{p}})\) has the clear two-step structure at low temperatures, pairing fluctuations at \(\boldsymbol{q}\ne 0\) are enhanced as \(T_{\text{env}}\) decreases. To conclude, the pronounced pseudogap seen around (c4) in Fig. 14a is induced by strong FFLO-type pairing fluctuations that are enhanced by the two-step structure of the non-equilibrium momentum distribution \(n_{\text{neq},\sigma }({\boldsymbol{p}})\).
Regarding the vanishing FFLO phase transition in the NETMA result, we recall that it has been shown in the thermal equilibrium state that the FFLO state is always unstable against rotational fluctuations in the spatially isotropic system, even in three dimensions [98,99,100,101,102,103]. When the system has a continuous rotational symmetry in space, the FFLO state is infinitely degenerate with respect to the direction of the FFLO-\(\boldsymbol{q}_{\text{pair}}\) vector, as schematically shown in Fig. 16. This infinite degeneracy remarkably enhances FFLO pairing fluctuations, which completely destroy the FFLO long-range order. We note that a similar instability phenomenon is also known in the ordinary BCS superfluid in one and two dimensions, which is sometimes referred to as the Hohenberg-Mermin-Wagner theorem in the literature [104, 105]. However, while the vanishing superfluid long-range order in one and two dimensions is also due to anomalously enhanced pairing fluctuations, in the FFLO case, the instability occurs even in three dimensions. Since the same instability mechanism also works in the present non-equilibrium case [40], the FFLO-type superfluid phase transition seen in the mean-field phase diagram in Fig. 9c completely vanishes in the NETMA phase diagram.
Schematic illustration of a non-equilibrium FFLO Cooper pairs. All superfluid states associated with a Cooper pair depicted in a–d are degenerate due to the continuous rotational symmetry of the gas system
We briefly note that, since the spatial isotropy of the system is crucial for the vanishing FFLO long-range order, the FFLO-type superfluid phase transition would be possible when this symmetry is removed from the present gas system, which we will examine in Section 5.
As explained in Section 3.4, the magnitude of the center-of-mass momentum \(|\boldsymbol{q}_{\text{pair}}|\) of a non-equilibrium FFLO Cooper pair depends on the difference between the “Fermi momenta” of the two effective “Fermi surfaces”. That is, \(|\boldsymbol{q}_{\text{pair}}|\sim p_{\text{F1}}-p_{\text{F2}}=\sqrt{2m[\mu +\Delta \mu ]} -\sqrt{2m[\mu -\Delta \mu ]}\) (see Fig. 11). Indeed, we see in Fig. 15c that the peak position shifts toward \(\boldsymbol{q}=0\) with decreasing \(\Delta \mu\) along the path (c) in Fig. 15a. When we arrive at \(T_{\text{env}}^{\text{c}}\), \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu )\) diverges at \(\boldsymbol{q}=0\) (that is, \(T_{\text{env}}^{\text{c}}\)-Eq. (61) is satisfied at \(\boldsymbol{q}=0\)), and the main system transitions to the uniform non-equilibrium BCS superfluid state.
The strong non-equilibrium FFLO pairing fluctuations also affect the single-particle excitation spectrum. As shown in Fig. 14(c4), the spectral function at (c4) in Fig. 14a exhibits a broad downward spectral peak, which is denoted as “FFLO”. To understand the mechanism of this anomalous spectral structure, we apply the pseudogap approximation to the NETMA self-energy \(\Sigma ^{\text{R}}_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\). As seen from Fig. 15b, \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu )\) has a large intensity at \((\boldsymbol{q}, \nu ) = (\boldsymbol{q}_{\text{pair}}, 2\mu )\) around (c4) in Fig. 14a. Using this, one approximates the self-energy in Eq. (64) to [38]
Here, \(\langle {\cdots }\rangle _{\boldsymbol{q}_{\text{pair}}}\) represents the average over the direction of the \(\boldsymbol{q}_{\text{pair}}\), and
is the pseudogap parameter.
Figure 17 shows the single-particle spectral function \(A_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) computed with the approximated self-energy in Eq. (73). While we set \(\boldsymbol{q}_{\text{pair}}/p_{\text{F}}=0\) in Fig. 17(a1), we take \(\boldsymbol{q}_{\text{pair}}/p_{\text{F}}=0.2\) in Fig. 17(a2). Comparing these two cases, we find that the “FFLO” structure seen in Fig. 17a is induced by non-equilibrium FFLO pairing fluctuations, because the broad peak structure appears only when \(\boldsymbol{q}_{\text{pair}}\ne 0\). The difference between the two cases (\(\boldsymbol{q}_{\text{pair}}=0\) and \(\boldsymbol{q}_{\text{pair}}\ne 0\)) arises from the degrees of freedom with respect to the direction of \(\boldsymbol{q}_{\text{pair}}\) due to the spatial isotropy of the gas system. In the case of BCS-type pairing, where a Cooper pair is formed between Fermi atoms at opposite momenta \(\boldsymbol{p}\) and \(-\boldsymbol{p}\), the excitation modes along the particle dispersion \(\omega = \varepsilon _{\boldsymbol{p}}\) is coupled only to the modes along the hole dispersion \(\omega = -\varepsilon _{-\boldsymbol{p}}+2\mu\). As explained in Section 2.2, the level repulsion between these excitation modes results in the pseudogap. On the other hand, in the case of FFLO-type paring, a Cooper pair is formed between Fermi atoms at momenta \(\boldsymbol{p}\) and \(-\boldsymbol{p}+\boldsymbol{q}_{\text{pair}}\). Then, the excitation modes along the particle dispersion \(\omega = \varepsilon _{\boldsymbol{p}}\) is coupled to the excitation modes along the hole dispersions \(\omega = -\varepsilon _{-\boldsymbol{p}+\boldsymbol{q}_{\text{pair}}} +2\mu\). As schematically shown in Fig. 17b, the latter is a continuous excitation mode between \(\omega = -\big [|\boldsymbol{p}|+|\boldsymbol{q}_{\text{pair}}|\big ]^2/(2m) +2\mu\) and \(\omega = -\big [|\boldsymbol{p}|-|\boldsymbol{q}_{\text{pair}}|\big ]^2/(2m) +2\mu\), because there is no restriction on the direction of \(\boldsymbol{q}_{\text{pair}}\) due to the spatial isotropy of the gas system. The coupling of the continuous spectrum with the excitation modes along the particle dispersion \(\omega =\varepsilon _{\boldsymbol{p}}\) results in the characteristic spectral structure “FFLO” in Fig. 14(c4).
a Single-particle spectral function in the pseudogap approximation. We set (a1) \(\boldsymbol{q}_{\text{pair}}\ne 0\) and (a2) \(\boldsymbol{q}_{\text{pair}}=0\). A broad downward spectral peak “FFLO” in panel (a1) is consistent with Fig. 14(c4). b Mechanism of the broad spectral peak. When a Cooper pair has non-zero center-of-mass momentum \(\varvec{q}_{\text {pair}}\), the excitation modes along the particle dispersion \(\omega = \varepsilon _{\boldsymbol{p}}\) is coupled to the excitation modes along the hole dispersion \(\omega = -\varepsilon _{-\boldsymbol{p}+\boldsymbol{q}_{\text{pair}}} +2\mu\) with various directions of \(\boldsymbol{q}_{\text{pair}}\)
To summarize the discussion, the phase diagram of the driven-dissipative Fermi gas in the weak-coupling BCS regime is obtained as Fig. 18. In this phase diagram, \(T^*_{\text{env}}\) is the so-called pseudogap temperature [68], where the pseudogap appears in the density of states \(\rho _{\text{NETMA}, \sigma }(\omega )\). The pseudogap phase between \(T^{\text{c}}_{\text{env}}\) and \(T^*_{\text{env}}\) can be further divided into two regions:
- BCS-type pseudogap phase (\(\boldsymbol{q}_{\text{pair}}=0\)): the pseudogap is induced by BCS-type pairing fluctuations. In this phase, \(-\text{Re}\big [\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu )\big ]\) has the peak at \(\boldsymbol{q}=0\).
- FFLO-type pseudogap phase (\(\boldsymbol{q}_{\text{pair}}\ne 0\)): the pseudogap is induced by non-equilibrium FFLO-type pairing fluctuations. In this phase, the peak appears at \(\boldsymbol{q}=\boldsymbol{q}_{\text{pair}}\) (\(\ne 0\)) and the spectral function \(A_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) exhibits the broad downward spectral peak.
Phase diagram of the driven-dissipative Fermi gas in the BCS regime, in terms of the temperature \(T_{\text{env}}\) and the chemical potential bias \(\Delta \mu\). \(T_{\text{env}}^{\text{c}}\) (solid line) is the superfluid phase transition temperature and \(T^*_{\text{env}}\) (dashed line) is the pseudogap temperature. The pseudogap phase between \(T_{\text {env}}^{\text {c}}\) and \(T_{\text{env}}^*\) is divided into two regions: one is the pseudogap phase induced by BCS-type pairing fluctuations (green-shaded area). In this phase, \(-\text{Re}\big [\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu )\big ]\) has the peak at \(\boldsymbol{q}=0\). The other pseudogap phase is induced by non-equilibrium FFLO-type pairing fluctuations (blue-shaded area). In this region, the peak appears at \(\boldsymbol{q}=\boldsymbol{q}_{\text{pair}}\) (\(\ne 0\))
Since the non-equilibrium FFLO superfluid state is unstable against its own paring fluctuations, only the BCS-type uniform superfluid state is realized inside the region surrounded by \(T_{\text{env}}^{\text{c}}\) in Fig. 18.
4.2.2 Strong-coupling BEC regime
We next consider the strong-coupling BEC regime. In this regime, as shown in Fig. 19a, \(T_{\text{env}}^{\text{c}}\) does not exhibit re-entrant behavior, but monotonically decreases as the chemical potential bias \(\Delta \mu\) increases. This is simply because the non-equilibrium momentum distribution \(n_{\text{neq}, \sigma }(\boldsymbol{p})\) does not exhibit the two-step structure in this regime, so that, in contrast to the weak-coupling BCS case shown in Fig. 15b, the anomalous enhancement of FFLO-type pairing fluctuations does not occur.
a Calculated \(T_{\text{env}}^{\text{c}}\) in the BEC regime, when \((p_{\text{F}}a_s)^{-1}=0.6\). b Density of states \(\rho _{\text{NETMA}, \sigma }(\omega )\) along the paths (b1) and (b2) in panel a. c Single-particle spectral function \(A_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\). d Photoemission spectrum \(L_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\). Panels (c1)–(c3) and (d1)-(d3) show the results at positions (1)–(3) in panel a
Figure 19b shows how the density of states \(\rho _{\text{NETMA}, \sigma }(\omega )\) varies when one moves along the paths (b1) and (b2) in Fig. 19a. Figure 19c shows the single-particle spectral function \(A_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) at positions (1)–(3) in Fig. 19a. These figures indicate that the single-particle excitation spectra do not change so significantly as in the weak-coupling BCS regime at least within the variations of \(\Delta \mu\) and \(T_{\text{env}}\) along the paths (b1) and (b2) in Fig. 19a.
On the other hand, the non-equilibrium effects appear in the photoemission spectrum \(L_{\text{NETMA}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (72). Comparing Fig. 19(d1) with (d2), the occupancy of the excitation modes along the particle dispersion (upper branch) is higher in (d2). These excitation modes correspond to single-particle excitations of Fermi atoms due to the dissociation of diatomic molecules [68]. Thus, the higher intensity of the upper branch in Fig. 19(d2) means that the diatomic molecules are more destroyed by strong pumping and decay of Fermi atoms in the large \(\Delta \mu\) regime.
We briefly note that in Fig. 19(d3), the occupancy of the upper branch is lower than in (d2). This is because the thermal dissociation of diatomic molecules is suppressed in the low \(T_{\text{env}}\) region.
5 Stable non-equilibrium FFLO superfluid in the driven-dissipative lattice Fermi gas
5.1 Driven-dissipative lattice Fermi gas
As explained in Section 4.2.1, the strong non-equilibrium FFLO pairing fluctuations resulting from the continuous rotational symmetry of the gas system completely destroy the FFLO-type long-range order in the driven-dissipative Fermi gas. However, this also indicates that, if the infinite degeneracy shown in Fig. 16 is lifted and FFLO-type pairing fluctuations can be suppressed, by removing the continuous rotational symmetry of the system, the non-equilibrium FFLO superfluid state may be stabilized.
To explore this possibility, we consider the case when the main system is loaded on a three-dimensional cubic optical lattice [106,107,108]. To model this situation, we describe the main system by the three-dimensional attractive Hubbard model [109], given by
Here, \(-U\) is the on-site s-wave pairing interaction, and the kinetic energy \(\varepsilon _{\boldsymbol{k}}\) of the lattice fermion has the form
with t and \(t'\) being, respectively, the nearest-neighbor and the next-nearest-neighbor (NNN) hopping amplitude. In Eq. (76), the lattice constant a is taken to be unity, for simplicity. In the following discussions, we use the NNN hopping \(t'\) to tune the anisotropy of the Fermi surface. As shown in Fig. 20, the Fermi surface shape gradually deviates from the spherical one as \(t'\) increases. Thus, by adjusting the value of \(t'\), we can control the spatial anisotropy of the main system.
Fermi surface of a free Fermi gas with the band dispersion \(\varepsilon _{\boldsymbol{k}}\) in Eq. (76) and effects of the NNN hopping \(t'\). We take \(n=0.3\)
5.2 Hartree self-energy in the driven-dissipative lattice Fermi gas
In the absence of the optical lattice, the Hartree self-energy (\(\propto U\)) vanishes because \(U\rightarrow 0\) in the limit of \(p_{\text{c}} \rightarrow \infty\), where \(p_{\text{c}}\) is the momentum cutoff in Eq. (3) [13]. On the other hand, in the case of the lattice Fermi gas, \(U \ne 0\) for a given s-wave scattering length \(a_{\text{s}}\), due to the finite bandwidth of the energy band \(\varepsilon _{\boldsymbol{k}}\) in Eq. (76). Thus, the Hartree term gives a non-zero contribution in the lattice system, and we need to incorporate it into the theory.
The Hartree self-energy \(\hat{\Sigma }_{\text{H}, \sigma }(\boldsymbol{k}, \omega )\) is diagrammatically given in Fig. 21a, which gives [40]
where
is the filling fraction of Fermi atoms at each lattice site. Here, \(\tilde{G}^{\text{K}}_{\text{neq}, \sigma }(\boldsymbol{k}, \omega )\) is the Keldysh component of the non-equilibrium Green’s function \(\hat{\tilde{G}}_{\text{neq}, \sigma }(\boldsymbol{k}, \omega )\) in the lattice system, which obeys the Dyson equation [40]
a Dyson equation for the non-equilibrium Green’s function \(\hat{\tilde{G}}_{\text{neq}, \sigma }\) in the presence of the optical lattice. b Particle-particle scattering matrix \(\hat{\tilde{\Gamma }}_{\text{neq}}\) in the driven-dissipative lattice Fermi gas. c Diagrammatic representation of the NENSR Green’s function \(\hat{\tilde{G}}_{\text{NENSR}, \sigma }\) in Eq. (82)
This equation is also diagrammatically drawn as Fig. 21a. In the Dyson Eq. (79), \(\hat{G}_{0, \sigma }(\boldsymbol{k}, \omega )\) is the Green’s function in the isolated non-interacting Fermi gas loaded on the optical lattice, which is obtained by simply replacing \(\varepsilon _{\boldsymbol{p}}=p^2/(2m)\) in Eq. (11) with the dispersion \(\varepsilon _{\boldsymbol{k}}\) in Eq. (76). The self-energy correction \(\hat{\Sigma }_{\text{env}, \sigma }(\boldsymbol{k}, \omega )\) due to the system-reservoir couplings is given in Eq. (45). Solving the Dyson Eq. (79), one has
where the renormalized the kinetic energy \(\tilde{\varepsilon }_{\boldsymbol{k}, \sigma }\) has the form,
Since \(n_{\text{H}, -\sigma }\) in Eq. (81) is determined from Eq. (78), we need to solve self-consistently Eqs. (78) and (80).
5.3 Non-equilibrium Nozières-Schmitt-Rink strong-coupling theory
Although the spatial isotropy can be removed by considering the tight-binding model given in Eq. (76), it also makes numerical calculation harder than the isotropic case, because we cannot simplify computations by using the spatial isotropy of the main system. Thus, we take into account pairing fluctuations in the present lattice system by extending the strong-coupling theory developed by Nozières and Schmitt-Rink (NSR) [65] in the thermal equilibrium state to the non-equilibrium steady state.
The non-equilibrium NSR theory (NENSR) may be viewed as a simplified version of NETMA explained in Section 4.1, because the non-equilibrium Green’s function \(\hat{\tilde{G}}_{\text{NENSR}, \sigma }(\boldsymbol{k}, \omega )\) in NENSR obeys the truncated Dyson equation,
which is diagrammatically drawn as Fig. 21c. [Note that the thick solid line in the last diagram in Fig. 12 is replaced by the double solid line in Fig. 21c.] The NENSR self-energy \(\hat{\tilde{\Sigma }}_{\text{NENSR}, \sigma }(\boldsymbol{k}, \omega )\) is the same as the NETMA self-energy \(\hat{\Sigma }_{\text{NETMA}, \sigma }(\boldsymbol{k}, \omega )\) in Eqs. (64) and (65), except that the ladder diagrams consists of \(\hat{\tilde{G}}_{\text{neq}, \sigma }(\boldsymbol{k}, \omega )\) in Eq. (80) instead of \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) in Eq. (48), as diagrammatically shown in Fig. 21b. Thus, the NENSR self-energy is obtained by replacing \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) and \(\hat{\Gamma }_{\text{neq}}(\boldsymbol{q}, \omega )\) in the NETMA self-energy with \(\hat{\tilde{G}}_{\text{neq}, \sigma }(\boldsymbol{k}, \omega )\) and \(\hat{\tilde{\Gamma }}_{\text{neq}}(\boldsymbol{q}, \omega )\), where
is the particle-particle scattering matrix in the presence of the optical lattice. Here the pair correlation function \(\hat{\tilde{\Pi }}_{\text{neq}}(\boldsymbol{q}, \nu )\) in the lattice system is obtained by replacing \(\hat{G}_{\text{neq}, \sigma }(\boldsymbol{p}, \omega )\) in Eqs. (55) and (56) with \(\hat{\tilde{G}}_{\text{neq}, \sigma }(\boldsymbol{k}, \omega )\).
The equation for the filling fraction n per lattice site, which incorporates the effects of both pairing fluctuations and system-reservoir couplings, is obtained from the Keldysh component \(\tilde{G}^{\text{K}}_{\text{NENSR}, \sigma }(\boldsymbol{k}, \omega )\) of the NENSR Green’s function as
Here, \(n_{\text{H}, \sigma }\) is given in Eq. (78) and
In Eq. (84), \(n_{\text{FL}, \sigma }\) is the fluctuation correction to the filling fraction.
The non-equilibrium Thouless criterion in the NENSR theory reads
As in the NETMA case, the real part and the imaginary part of Eq. (86) gives two equations, and the imaginary part gives the solution \(\mu _{\text{pair}}=2\mu\). Then, substituting this into the real part of Eq. (86), one reaches the following \(T_{\text{env}}^{\text{c}}\) equation for the driven-dissipative lattice Fermi gas:
In the NENSR theory, one solves the \(T_{\text{env}}^{\text{c}}\) Eq. (87), together with the equation for the filling fraction n in Eq. (84), to self-consistently determine \(T_{\text{env}}^{\text{c}}\) and \(\mu (T_{\text{env}}^{\text{c}})\) for a given parameter set \((n, \gamma , \Delta \mu )\). In the \(T_{\text{env}}^{\text{c}}\) Eq. (87), \(\boldsymbol{q}_{\text{pair}}\) is determined so as to obtain the highest \(T_{\text{env}}^{\text{c}}\). The self-consistent solution with \(\boldsymbol{q}_{\text{pair}}=0\) (\(\boldsymbol{q}_{\text{pair}}\ne 0\)) describes the phase transition into the spatially uniform BCS (non-uniform FFLO) type superfluid.
5.4 Stabilization of the non-equilibrium FFLO superfluid
Figure 22a shows the NENSR results on \(T_{\text{env}}^{\text{c}}\) in the driven-dissipative lattice Fermi gas. In contrast to the spatially isotropic case without the optical lattice, where the non-equilibrium FFLO superfluid is completely destroyed by pairing fluctuations as shown in Figs. 14a, 22a and b clearly show that the non-equilibrium FFLO superfluid survives against their pairing fluctuations in the presence of the optical lattice.
a Calculated \(T_{\text{env}}^{\text{c}}\) in the driven-dissipative lattice Fermi gas, as functions of the chemical potential bias \(\Delta \mu\) and the system-reservoir coupling strength \(\gamma\). b Magnitude \(|\boldsymbol{q}_{\text{pair}}|\) of the center-of-mass momentum of a Cooper pair at \(T_{\text{env}}^{\text{c}}\). We take \(t'=0\), \(n=0.3\), and \(U/(6t)=0.8\)
Figure 23(a2) shows the intensity \(-\text{Re}\big [\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu ) \big ]\) of the retarded particle-particle scattering matrix at the solid square in Fig. 23(a1) in the absence of the optical lattice. As explained in Section 4.2.1, this quantity has large intensity around \(|\boldsymbol{q}|=|\boldsymbol{q}_{\text{pair}}|\ne 0\), reflecting the enhancement of non-equilibrium FFLO pairing fluctuations.
(a1) Calculated \(T_{\text{env}}^{\text{c}}\) in the absence of the optical lattice. This is the same plot as Fig. 14a. (a2) Calculated intensity \(-\text{Re}\big [\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu ) \big ]\) of the particle-particle scattering matrix at the solid square in panel (a1). Panel b shows the case in the presence of the optical lattice. In panel b, we set \(U/(6t)=0.8\), \(\gamma /(6t)=0.01\), \(t'=0\), and \(n=0.3\)
Figure 23(b2) shows the result at the solid square in Fig. 23(b1) in the presence of the optical lattice, which is quite different from the gas case shown in Fig. 23(a2). Since the spatial isotropy is explicitly broken by the optical lattice, the ring structure seen in Fig. 23(a2) is not obtained. Instead, Fig. 23(b2) shows that \(-\text{Re}\big [\tilde{\Gamma }^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu ) \big ]\) has peaks at \(\boldsymbol{q}=\boldsymbol{q}^{j, \eta }_{\text{pair}}\) (\(j=x,y,z\) and \(\eta =\pm\)), reflecting the discreate rotational symmetry of the optical lattice potential.
The above difference between the gas case and the lattice case results in a significant difference in the fluctuation correction term \(n_{\text{FL}, \sigma }\) in Eq. (84). In the presence of the optical lattice, noting that the \(\Gamma ^{\text{R}}_{\text{neq}}(\boldsymbol{q}, 2\mu )\) is enhanced around \((\boldsymbol{q}, \nu )= (\boldsymbol{q}_{\text{piar}}^{j, \eta }, 2\mu )\) near the non-equilibrium FFLO transition, we can approximate the NENSR self-energy to
where
is the pseudogap parameter. Substituting Eq. (88) into \(n_{\text{FL}, \sigma }\) in Eq. (84), one has
To evaluate the pseudogap parameter \(\Delta _{\text{PG}}\) in Eq. (89), we also employ the following approximation:
where
with \(N(\mu )\) being the density of states in the main system at \(\omega =\mu\). In obtaining Eq. (91), for simplicity, we have taken the limit \(\gamma \rightarrow 0^+\) in \(\tilde{\Pi }^{\text{R}}_{\text{neq}}(\boldsymbol{q}, \nu )\) as,
and have expanded Eq. (94) around \((\boldsymbol{q}, \nu ) = (\boldsymbol{q}_{\text{pair}}^{j, \eta }, 2\mu )\). Substituting Eq. (91) into Eq. (89), one has
In the same manner, we can evaluate the pseudogap parameter in the spatially isotropic gas system, which reads [40]
with \(q_{\text{c}}\) being a cutoff momentum.
The crucial difference between Eq. (95) and Eq. (96) is that the factor \(1/[|\boldsymbol{q}| - |\boldsymbol{q}_{\text{pair}}|]^2\) in the gas case is replaced by \(1/[\boldsymbol{q} -\boldsymbol{q}^{j, \eta }_{\text{pair}}]^2\) in the lattice Fermi gas. Then, the q integral in the spatially isotropic case in Eq. (96) always diverges as far as \(\boldsymbol{q}_{\text{pair}} \ne 0\), irrespective of the system dimension. As seen from Eq. (90), the fluctuation correction term \(n_{\text{FL}, \sigma }\) in Eq. (90) diverges when the pseudogap parameter \(\Delta ^2_{\text{PG}}\) diverges. Because of this singularity, the non-equilibrium FFLO superfluid is prohibited in the spatially isotropic gas system. On the other hand, replacing \(\boldsymbol{q}-\boldsymbol{q}_{\text{pair}}^{j, \eta }\) by \(\boldsymbol{q}\) in Eq. (95), one finds that the pseudogap parameter \(\Delta ^2_{\text{PG}}\) in the lattice system converges even when \(|\boldsymbol{q}_{\text{pair}}^{j, \eta }|\ne 0\). Thus, in the presence of the optical lattice, the NENSR coupled Eqs. (84) with (86) can be satisfied simultaneously at the non-equilibrium FFLO phase transition, where \(|\boldsymbol{q}_{\text{pair}}|\ne 0\).
The essence of the stabilization of the non-equilibrium FFLO-type long-range order is, strictly speaking, not the optical lattice itself, but rather the resulting anisotropy of the Fermi surface edges in momentum space.
To demonstrate this, we show in Fig. 24 how the stabilization of the non-equilibrium FFLO state is affected by the anisotropy of the Fermi surface. (We recall that, as shown in Fig. 20, the Fermi surface becomes more spherical for larger \(t'\).) As expected, the non-equilibrium FFLO phase transition is gradually suppressed as the Fermi surface shape becomes close to the sphere by increasing the value of \(t'\). In particular, the temperature at the boundary between the non-equilibrium BCS and the FFLO phase (solid circle in Fig. 24), which is referred to as the Lifshitz point in the literature [110, 111], decreases with increasing \(t'\). This is because stronger non-equilibrium FFLO pairing fluctuations are introduced by more spherical Fermi surface edges.
a Calculated \(T_{\text{env}}^{\text{c}}\) and effects of the NNN hopping amplitude \(t'\). b \(|\boldsymbol{q}_{\text{pair}}|\) at \(T_{\text{env}}^{\text{c}}\). We take \(n=0.3\), \(U/(6t)=0.8\), and \(\gamma /(6t)=0.015\). The temperature \(T_{\text{env}}\) is normalized by the value at \(\Delta \mu =0\) (\(\equiv T_{\text{env}}^{\text{c0}}\)). The solid circle shows the Lifshitz point, at which the normal, the non-equilibrium BCS, and the non-equilibrium FFLO phases meet one another
In the current stage of cold Fermi gas physics, a superfluid \(^6\text{Li}\) Fermi gas in an optical lattice has been realized only when the lattice potential is very shallow [112]. Although this experimental situation is somehow different from the simple Hubbard model in Eq. (75), Fig. 24 indicates that the crucial key to stabilize the non-equilibrium FFLO superfluid is not the detailed lattice potential but the resulting anisotropic Fermi surface edge. Thus, if such a shallow optical lattice can still deform the Fermi surface enough to suppress the non-equilibrium FFLO pairing fluctuations, the non-equilibrium FFLO superfluid would be realized there.
6 Summary
Figure 25 summarizes the non-equilibrium BCS-BEC crossover physics in the driven-dissipative Fermi gas. We have calculated the superfluid transition temperature \(T_{\text{env}}^{\text{c}}\) in the non-equilibrium BCS-BEC crossover regime, as shown in Fig. 25a. This is done by extending the T-matrix approximation, developed in thermal equilibrium BCS-BEC crossover physics, to the driven-dissipative non-equilibrium Fermi gas, by employing the Keldysh Green’s function technique.
Summary of this article. a Phase diagram of the driven-dissipative Fermi gas in the non-equilibrium BCS-BEC crossover regime (see Section 4). b In the BEC regime, as the chemical potential bias \(\Delta \mu\) increases, the diatomic molecules are destroyed and the superfluid phase is monotonically suppressed (see Section 4.2.2). c In the BCS regime, the non-equilibrium momentum distribution \(n_{\text{neq},\sigma }(\boldsymbol{p})\) induces FFLO-type Cooper pairs with non-zero center-of-mass momentum \(\boldsymbol{q}_{\text{pair}}\) (see Section 3.4). d In the gas system with continuous rotational symmetry, the non-equilibrium FFLO superfluid is unstable against paring fluctuations due to the infinite degeneracy with respect to the direction of \(\varvec{q}_{\text {pair}}\) (see Section 4.2.1). This instability causes the re-entrant behavior of \(T_{\text{env}}^{\text{c}}\) and the anomalous pseudogap phenomenon. In the presence of the optical lattice, the infinite degeneracy is lifted, so that the non-equilibrium FFLO superfluid can survive against paring fluctuations (see Section 5)
The behavior of \(T_{\text{env}}^{\text{c}}\) is quite different between the weak-coupling BCS and the strong-coupling BEC regime. In the BEC regime, \(T_{\text{env}}^{\text{c}}\) is monotonically suppressed with increasing of the chemical potential bias \(\Delta \mu\). This is because the strong pumping and decay of atoms destroy the diatomic molecules, as schematically shown in Fig. 25b. On the other hand, \(T_{\text{env}}^{\text{c}}\) shows the re-entrant behavior in the BCS regime. This phenomenon is caused by strong FFLO-type pairing fluctuations, being enhanced by the non-equilibrium momentum distribution with the two-step structure, as shown in Fig. 25c. The instability of the non-equilibrium FFLO superfluid against pairing fluctuations results in the re-entrant behavior of \(T_{\text{env}}^{\text{c}}\), as well as the anomalous pseudogap phenomenon.
A possible route to realize the long-range order of this non-equilibrium FFLO state is the removal of the infinite degeneracy with respect to the direction of the center-of-mass momentum \(\boldsymbol{q}_{\text{pair}}\) of the non-equilibrium FFLO Cooper pairs in the isotropic gas system. Indeed, when the main system is loaded on the cubic optical lattice, the lifting of this infinite degeneracy suppresses pairing fluctuations, which stabilizes the long-range order of the non-equilibrium FFLO superfluid in the driven-dissipative lattice Fermi gas.
The recent progress in cold Fermi gas experiments has enabled us to study BCS-BEC crossover physics in a variety of non-equilibrium situations. With this experimental progress, the theoretical understanding of non-equilibrium BCS-BEC crossover physics would become increasingly important in this research field.