Stable non-equilibrium Fulde-Ferrell-Larkin-Ovchinnikov state in a spin-imbalanced driven-dissipative Fermi gas loaded on a three-dimensional cubic optical lattice

The realization of unconventional superfluid states is one of the most exciting challenges in the current stage of cold Fermi gas physics. At present, although various non-$s$-wave pairing states have been discovered in metallic superconductivity, as well as in liquid ³He, the simplest isotropic $s$-wave superfluid state has only been realized in ⁴⁰K Jin2004 and ⁶Li Zwierlein2004 ; Kinast2004 ; Bartenstein2004 Fermi gases. Since the high tunability is an advantage of ultracold Fermi gases Chin2010 ; Barontini2013 ; Labouvie2015 ; Labouvie2016 ; Tomita2017 ; Chong2018 , once this challenge is achieved, one would be able to examine its various superfluid properties in a wide parameter region. Indeed, in ⁴⁰K and ⁶Li superfluid Fermi gases, a tunable pairing interaction associated with a Feshbach resonance Chin2010 has enabled systematic studies about how the weak-coupling Bardeen-Cooper-Schrieffer (BCS) type superfluid discussed in metallic superconductivity continuously changes to the Bose-Einstein condensation (BEC) of tightly bound molecules, with increasing the $s$-wave interaction strength Levin2005 ; Zwerger2012 ; Strinati2018 ; Ohashi2020 .

Among various candidates discussed in cold Fermi gas physics, the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state Fulde1964 ; Larkin1964 ; Takada1969 ; Matsuda2007 is a promising one. This unconventional Fermi superfluid is characterized by a spatially oscillating superfluid order parameter $\Delta({\bm{r}})$, which is symbolically written as

$$\Delta({\bm{r}})=\Delta e^{i{\bm{Q}}_{\rm FF}\cdot{\bm{r}}},$$

(1)

where ${\bm{Q}}_{\rm FF}$ physically describes the center-of-mass-momentum of a FFLO Cooper pair. Although the FFLO state was originally proposed in the context of metallic superconductivity under an external magnetic field Fulde1964 ; Larkin1964 ; Takada1969 ; Matsuda2007 , it has also recently been discussed in ultracold Fermi gases Strinati2018 ; Hu2006 ; Parish2007 ; Liao2010 ; Chevy2010 ; Kinnunen2018 . Cooper pairs in the FFLO state are formed between Fermi surfaces with different sizes as shown in Fig. 1(a), giving ${\bm{Q}}_{\rm FF}\neq 0$. Recently, the observation of the FFLO state has been reported in several superconducting materials, such as CeCoIn₅ Bianchi2002 ; Bianchi2003 ; Kumagai2006 , CeCu₂Si₂ Kitagawa2018 , KFe₂As₂ Cho2017 , FeSe Ok2020 ; Kasahara2020 , and $\kappa$-(BEDT-TTF)₂Cu(NCS)₂ Singleton2000 ; Wright2011 ; Mayaffre2014 . Thus, the realization of a FFLO superfluid Fermi gas would be important, in order for cold atom physics to catch up with this recent exciting progress in condensed matter physics.

At a glance, ultracold Fermi gases seem suitable for the FFLO state: (1) Although the FFLO state is known to be weak against impurity scatterings, one can prepare a very clean Fermi gas. (2) The splitting of Fermi surfaces between $\uparrow$-(pseudo)spin and $\downarrow$-(pseudo)spin components can immediately be prepared in a spin-imbalanced Fermi gas.

However, in spite of these advantages, the realization of the FFLO superfluid Fermi gas has not been reported yet. This seems because the current cold Fermi gas experiments are facing the following two serious difficulties: (i) In the presence of spin imbalance, the system undergoes the phase separation into the BCS superfluid region and the normal-state region of unpaired excess atoms, before reaching the desired FFLO phase transition Sheehy2006 ; He2008 ; Shin2008 . (ii) The spatial isotropy of the gas cloud anomalously enhances FFLO pairing fluctuations, which completely destroy the FFLO long-range order, even in three dimensions Shimahara1998 ; Ohashi2002 ; Radzihovsky2009 ; Radzihovsky2011 ; Yin2014 ; Jakubczyk2017 ; Wang2018 ; Zdybel2021 ; Kawamura2020_JLTP ; Kawamura2020 .

Regarding these problems, we have recently proposed the following two ideas Kawamura2020_JLTP ; Kawamura2020 ; Kawamura2022_2 : For (i), instead of using a spin-imbalanced Fermi gas, we proposed to use a spin-balanced driven dissipative Fermi gas, being coupled with two reservoirs as schematically shown in Fig. 2(a) Kawamura2020 . This is a non-equilibrium open system, where losses of particles and energy are continuously compensated by the two reservoirs Yamaguchi2012 ; Hanai2016 ; Hanai2017 ; Hanai2018 , and is known to exhibit various interesting phenomena that cannot be examined in the thermal equilibrium state Cross1993 ; Hoyle2006 ; Cross2009 . In this non-equilibrium system, we showed that, under a certain condition, the momentum distribution of Fermi atoms has a two-edge structure, originating from the chemical potential difference between the two reservoirs, as illustrated in Fig. 2(b). These edges work like two ‘Fermi surfaces’ with different sizes, which enhances the FFLO pair correlation without spin imbalance. Indeed, we showed within the non-equilibrium mean-field theory that the FFLO superfluid state is obtained Kawamura2022 , where Cooper pairs are formed between Fermi atoms near the two edges, as schematically shown in Fig. 1(b). However, we also found that the difficulty (ii) also exists in this non-equilibrium case, so that this desired mean-field solution is immediately destroyed, once FFLO pairing fluctuations are taken into account Kawamura2020_JLTP ; Kawamura2020 .

For the difficulty (ii), we clarified in a thermal-equilibrium spin-imbalanced Fermi gas that the destruction of the FFLO long-range order by FFLO pairing fluctuations can be avoided, when the spatial isotropy of the gas cloud is removed by loading the system on a three-dimensional cubic optical lattice Kawamura2022_2 . However, we also found that, as far as we deal with a spin-imbalanced Fermi gas, the stabilized FFLO state competes with the above-mentioned phase separation, so that careful parameter tuning is still needed.

In this paper, by combining these two ideas, we explore a possible route to reach the FFLO superfluid phase transition in ultracold Fermi gases, without facing the phase separation, as well as the destruction by FFLO pairing fluctuations. For this purpose, we again consider the model driven-dissipative two-component Fermi gas shown in Fig. 2(a), but this time we impose a three-dimensional optical lattice potential to the system. To include pairing fluctuations, we extend the strong-coupling theory developed in the thermal-equilibrium state by Nozières and Schmitt-Rink (NSR) Nozieres1985 , to the case when the system is out of equilibrium. Using this, we examine whether or not the combined two-step structure of the Fermi momentum distribution with the background optical lattice can stabilize the non-equilibrium FFLO (NFFLO) state, overwhelming the above-mentioned difficulties.

We note that the NFFLO state discussed in this paper is, strictly speaking, somehow different from the ordinary thermal-equilibrium FFLO state in the spin-imbalanced system: As illustrated in Fig. 1(a), the ordinary FFLO Cooper pairs are formed between $\uparrow$-spin particles near a larger Fermi surface and $\downarrow$-spin particles near a smaller one. In the NFFLO state, on the other hand, since the momentum distribution of each spin component has two edges [see Fig. 1(b)], Cooper pairs are formed between, not only $\uparrow$-spin particles near a larger Fermi surface edge and $\downarrow$-spin particles near a smaller Fermi surface edge, but also $\uparrow$-spin particles near the smaller Fermi surface edge and $\downarrow$-spin particles near the larger Fermi surface edge.

We also examine how these two kinds of FFLO states are related to each other, by considering the case with spin imbalance. In the thermal-equilibrium state, the spin imbalance causes the so-called Zeeman splitting between the $\uparrow$-spin and $\downarrow$-spin Fermi surfaces. When the system becomes out of equilibrium by adjusting the chemical potential difference between the two reservoirs, the spin imbalance further splits each edge structure in the momentum distribution into two, so that the system looks as if there are four Fermi surfaces [see Fig. 1(c)]. Including these ‘multiple Fermi surfaces’ within the framework of the non-equilibrium mean-field BCS theory, we draw a superfluid phase diagram with respect to the environmental temperature, the chemical potential difference, and a fictitious magnetic field to tune the spin imbalance.

Before ending this section, we quickly summarize recent progress in the driven-dissipative system. Recent experimental progress has enabled us to examine, not only classical, but also quantum many-body driven-dissipative systems, such as exciton-polaritons Kasprzak2006 ; Carusotto2013 ; Sanvitto2016 , superconducting circuits Houck2012 ; Ma2019 , optical cavities Brennecke2013 ; Vaidya2018 , as well as ultracold atomic gases Barontini2013 ; Labouvie2015 ; Labouvie2016 ; Tomita2017 ; Chong2018 . At present, the same driven-dissipative ultracold Fermi gas as the model shown in Fig. 2(a) has not been realized, the realization would be possible within the current experimental technology, by extending the recent transport experiment on a ${}^{6}{\rm Li}$ Fermi gas in a two-terminal configuration Brantut2012 ; Krinner2015 ; Husmann2015 ; Nagy2016 ; Krinner2017 , or employing a tilted triple-well optical trap Graefe2006 ; Mossmann2006 ; Liu2007 . For a more detailed experimental proposal, see Ref. Kawamura2022 .

This paper is organized as follows. In Sec. II, we explain how to extend the mean-field BCS theory, as well as the strong-coupling NSR theory, developed in the thermal equilibrium state to the non-equilibrium steady state of a driven-dissipative Fermi gas. In this extension, we take into account the effects of a background optical lattice, as well as spin imbalance. Using these theories, we examine the NFFLO phase transition and effects of pairing fluctuations in Sec. III. We also show how the optical lattice stabilizes the NFFLO long-range order there. In Sec. IV. we consider a driven-dissipative lattice Fermi gas with spin imbalance. In the phase diagram with respect to the environmental temperature, the chemical potential difference between the reservoirs, and a fictitious magnetic field to adjust the spin imbalance, we identify the region where the non-equilibrium BCS, NFFLO, and ordinary FFLO states appear. Throughout this paper, we set $\hbar=k_{\rm B}=1$, and the system volume $V$ is taken to be unity, for simplicity.

In this section, we deal with the spin-balanced case, by setting $h=0$ in Eqs. (8) and (9). The spin-imbalanced case is separately discussed in Sec. IV.

The upper panels in Fig. 6 show $T_{\rm env}^{\rm c}$ in a driven-dissipative spin-balanced Fermi gas, loaded on the cubic optical lattice. In this figure, we distinguish between the NBCS and NFFLO phase transitions from whether $|{\bm{Q}}_{\rm FF}|$ equals zero or not in the lower panels in Fig. 6. In contrast to the case with no optical lattice (where the mean-field NFFLO solution is completely destroyed by NFFLO pairing fluctuations Kawamura2020_JLTP ; Kawamura2020 , as shown in Fig. 7), Figs. 6(a1) and (b1) show that the NFFLO long-range order survives against pairing fluctuations in the lattice system. (We summarize the NMF and the NNSR theories in the absence of the optical lattice in Appendix B.)

We find from the comparison of Fig. 6(a1) with Fig. 6(b1) that pairing fluctuations tend to decrease $T_{\rm env}^{\rm c}$. However, apart from this difference, these figures also indicate that, once the NFFLO state is stabilized by the optical lattice, the NMF theory (which completely ignores pairing fluctuations) already captures the essential behavior of $T_{\rm env}^{\rm c}$ as a function of $\delta\mu$ and $\gamma$, at least when the filling fraction equals $n=n_{\uparrow}+n_{\downarrow}=0.3$. Indeed, as an example, when we extract $T_{\rm env}^{\rm c}$ at $\gamma/(6t)=0.005$ from Figs. 6(a1) and (b1), the NNSR result is found to be very similar to the NMF result, as shown in Fig. 8. Furthermore, when Fig. 8 is replotted with respect to the scaled variables $T_{\rm env}^{\rm c}/T_{\rm env}^{\rm c0}$ and $\delta\mu/T_{\rm env}^{\rm c0}$ (where $T_{\rm env}^{\rm c0}$ is the superfluid phase transition temperature at $\delta\mu=0$), the scaled NMF and NFFLO results almost coincide with each other, as shown in the inset in Fig. 8. That is, although pairing fluctuations remarkably damage the mean-field NFFLO solution in the absence of optical lattice, their effects are not so crucial in a lattice Fermi gas, at least when $n=0.3$.

Figure 9(a1) shows the intensity $-{\rm Re}\big{[}\Gamma^{\rm R}(\bm{q},\nu=2\mu)\big{]}$ of the retarded particle-particle scattering matrix near the region where the re-entrant phenomenon of the NBCS phase transition occurs [solid square in Fig. 9(a2)], in the absence of the optical lattice. Since this quantity physically describes pairing fluctuations with the center-of-mass momentum $\bm{q}$, the fact of the large intensity around $|{\bm{q}}|=|{\bm{Q}}|\neq 0$ means the enhancement of NFFLO pairing fluctuations there. We also point out that $|{\bm{Q}}|$ is directly related to the size difference between two ‘Fermi surface’ edges that are imprinted on the momentum distribution of Fermi atoms by the two reservoirs [see Fig. 9(a3)]. This means that these edges really work like two Fermi surfaces with different sizes, as in the spin-imbalanced case.

Figure 9(b1) shows the results in the presence of the three-dimensional cubic optical lattice, which is quite different from Fig. 9(a1). [Figure 9(b1) is obtained at the solid square in Fig. 9(b2).]: Since the spatial isotropy is broken by the optical lattice, the ring structure seen in Fig. 9(a1) is not obtained, but the $-{\rm Re}\big{[}\Gamma^{\rm R}(\bm{q},\nu=2\mu)\big{]}$ exhibits four peaks, reflecting the four-fold rotational symmetry of the cubic lattice. However, for example, the peak at ${\bm{q}}={\bm{Q}}^{+}_{x}$ in Fig. 9(b1) is still related to the size difference between the two Fermi surface edges, as shown in Fig. 9(b3). That is, these edges also work like two Fermi surfaces, to enhance NFFLO pairing fluctuations around ${\bm{q}}={\bm{Q}}^{+}_{x}$, as well as the other equivalent peaks in Fig. 9(b1).

The above-mentioned difference seen in Figs. 9(a1) and (b1) makes a significant difference in the NNSR fluctuation correction terms in Eqs. (38) and (68): In the presence of the optical lattice, noting that the particle-particle scattering matrix $\hat{\Gamma}(\bm{q},\nu)$ in Eq. (28) is enhanced around $(\bm{q},\nu)=(\bm{Q}^{\eta}_{j},2\mu)$ near the NFFLO phase transition (where ${\bm{Q}}^{\eta=\pm}_{j=x,y,z}$ represent the four peak positions in Fig. 9(b1), as well as other two peak positions existing along the $q_{z}$ axis), we approximate the self-energy in Eqs. (36) and (37) to, near $T_{\rm env}^{\rm c}$,

	$\displaystyle\hat{\Sigma}_{{\rm NNSR},\sigma}(\bm{k},\omega)$
	$\displaystyle\simeq\scalebox{0.92}{$\displaystyle-\Delta_{\rm pg}^{2}\sum_{\eta=\pm}\sum_{j=x,y,z}\begin{pmatrix}G^{\rm A}_{{\rm NMF},-\sigma}&G^{\rm K}_{{\rm NMF},-\sigma}\\[4.0pt] 0&G^{\rm R}_{{\rm NMF},-\sigma}\end{pmatrix}(\bm{Q}^{\eta}_{j}-\bm{k},2\mu-\omega)$}$
	$\displaystyle=\Delta_{\rm pg}^{2}\sum_{\eta=\pm}\sum_{j=x,y,z}\hat{G}^{*}_{{\rm NMF},-\sigma}(\bm{Q}^{\eta}_{j}-\bm{k},2\mu-\omega).$		(40)

Here, the so-called pseudogap parameter,

$$\Delta_{\rm pg}^{2}=\frac{i}{2}\sum_{\bm{q}}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi}\Gamma^{\rm K}(\bm{q},\nu),$$

(41)

physically describes the strength of pairing fluctuations Chen2005 ; Tsuchiya2009 . Evaluating the fluctuation correction $n_{{\rm FL},\sigma}$ involved in Eq. (38) by using Eq. (40), one has

	$\displaystyle n_{{\rm FL},\sigma}$	$\displaystyle=\frac{i\Delta_{\rm pg}^{2}}{2}\sum_{\eta=\pm}\sum_{j=x,y,z}\sum_{\bm{k}}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\Big{[}\hat{G}_{{\rm NMF},\sigma}(\bm{k},\omega)$
		$\displaystyle\hskip 14.22636pt\times\hat{G}^{*}_{{\rm NMF},-\sigma}(\bm{Q}_{j}^{\eta}-\bm{k},2\mu-\omega)\hat{G}_{{\rm NMF},\sigma}(\bm{k},\omega)\Big{]}^{\rm K}.$		(42)

To evaluate the pseudogap parameter $\Delta_{\rm pg}^{2}$, we approximate $\Gamma^{\rm R}(\bm{q},\nu)$ to

$$\Gamma^{\rm R}(\bm{q},\nu)\simeq\sum_{\eta=\pm}\sum_{j=x,y,z}\frac{-U}{C\big{[}\bm{q}-\bm{Q}_{j}^{\eta}\big{]}^{2}-i\lambda\big{[}\nu-2\mu\big{]}},$$

(43)

where we have assumed that $T_{\rm env}^{\rm c}$ satisfies Eq. (33), and

	$\displaystyle C=\frac{U}{2}\left.\nabla_{\bm{q}}^{2}\Pi^{\rm R}(\bm{q},2\mu)\right|_{\bm{q}=\bm{Q}_{j}^{\pm}},$		(44)
	$\displaystyle\lambda=\frac{\pi U}{8T_{\rm env}}N(\mu){\rm sech}^{2}\left(\frac{\delta\mu}{2T_{\rm env}}\right),$		(45)

with $N(\mu)$ being the density of states in the main system at $\omega=\mu$. In obtaining Eq. (43), for simplicity, we have taken the limit $\gamma\to+0$ in $\Pi^{\rm R}(\bm{q},\nu)$ in Eq. (29), giving

	$\displaystyle\lim_{\gamma\to+0}\Pi^{\rm R}(\bm{q},\nu)$
	$\displaystyle=\sum_{\bm{p}}\frac{1-f(\tilde{\varepsilon}_{\bm{p}+\bm{q}/2,\uparrow}-\mu-\delta\mu)-f(\tilde{\varepsilon}_{\bm{p}-\bm{q}/2,\downarrow}-\mu+\delta\mu)}{\nu_{+}-\tilde{\varepsilon}_{\bm{p}+\bm{q}/2,\uparrow}-\tilde{\varepsilon}_{\bm{p}-\bm{q}/2,\downarrow}},$		(46)

and have expanded it around $(\bm{q},\nu)=(\bm{Q}_{j=x,y,z}^{\eta=\pm},2\mu)$. Using Eqs. (28) and (43), one reaches

	$\displaystyle\Delta_{\rm pg}^{2}$	$\displaystyle=\frac{i}{2}\sum_{\bm{q}}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi}\left|\Gamma^{\rm R}(\bm{q},\nu)\right|^{2}\Pi^{\rm K}(\bm{q},\nu)$
		$\displaystyle\simeq\sum_{\eta=\pm}\sum_{j=x,y,z}\frac{iU^{2}\Pi^{\rm K}(\bm{Q}_{j}^{\eta},2\mu)}{2}$
		$\displaystyle\hskip 22.76228pt\times\sum_{\bm{q}}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi}\frac{1}{C^{2}\big{[}\bm{q}-\bm{Q}_{j}^{\eta}\big{]}^{4}+\lambda^{2}\big{[}\nu-2\mu\big{]}^{2}}$
		$\displaystyle=\sum_{\eta=\pm}\sum_{j=x,y,z}\frac{iU^{2}\Pi^{\rm K}(\bm{Q}_{j}^{\eta},2\mu)}{4\lambda C}\sum_{\bm{q}}\frac{1}{\big{[}\bm{q}-\bm{Q}_{j}^{\eta}\big{]}^{2}}.$		(47)

In deriving the second line, we have employed the same approximation as that used in deriving Eq. (40). Replacing $\bm{q}-\bm{Q}_{j}^{\eta}$ by $\bm{q}$ in Eq. (47), one finds that $\Delta_{\rm pg}^{2}$ converges in three dimension, irrespective of the value of $\bm{Q}_{j}^{\eta}$. This immediately concludes the convergence of $n_{{\rm FL},\sigma}$ [which is proportional to $\Delta_{\rm pg}^{2}$, see Eq. (42)]. Thus, the NNSR coupled equations (33) and (38) can be satisfied simultaneously at the NFFLO phase transition (where ${\bm{Q}}_{j}^{\eta}\neq 0$). We briefly note that the six NFFLO vectors ${\bm{Q}}_{j=x,y,z}^{\eta=\pm}$ have the same magnitude, being equal to $|{\bm{Q}}_{\rm FF}|$ in Fig. 6(b2).

A quite different phenomenon occurs in the absence of the optical lattice: In this spatially isotropic case, when the intensity $-{\rm Re}\big{[}\Gamma^{\rm R}(\bm{q},\nu=2\mu)\big{]}$ of the retarded particle-particle scattering matrix exhibits a ring structure as seen in Fig. 9(a1), the pseudogap parameter $\Delta_{\rm pg}^{2}$ can be approximated to

$$\Delta_{\rm pg}^{2}\simeq\frac{iU^{2}\Pi^{\rm K}(\bm{Q},2\mu)}{4\lambda C}\int_{0}^{q_{\rm c}}\frac{q^{2}dq}{2\pi^{2}}\frac{1}{\big{[}|\bm{q}|-|\bm{Q}|\big{]}^{2}},$$

(48)

where $q_{\rm c}$ is a cutoff momentum. (For the derivation, see Appendix C.) Comparing Eq. (48) with Eq. (47), the factor $1/\big{[}\bm{q}-\bm{Q}_{j}^{\eta}\big{]}^{2}$ in the lattice case is now replaced by $1/\big{[}|\bm{q}|-|\bm{Q}|\big{]}^{2}$, reflecting the isotropic edge positions shown in Fig. 9(a3). Then, the $q$-integration in the pseudogap parameter $\Delta_{\rm pg}$ in Eq. (48) always diverges as far as $\bm{Q}\neq 0$, even in three dimensions. This means that the divergence of the correction term $N_{{\rm FL},\sigma}$ in the NNSR number equation (68). Because of this singularity, the NNSR coupled equations (60) and (68) are never satisfied simultaneously, which prohibits the NFFLO phase transition.

We note that the essence of the stabilization mechanism of the NFFLO state is, strictly speaking, not the detailed lattice potential itself, but the resulting anisotropy of the ‘Fermi surface’ edges shown in Fig. 9(b3). Indeed, as seen in the left panels in Fig. 10, when one deforms the shape of these edges to be more spherical by increasing the value of the next nearest neighbor hopping $t^{\prime}$ (see Fig. 11), the NFFLO region obtained in the NNSR theory shrinks. Since the NMF result is not sensitive to $t^{\prime}$ (see the right panels in Fig. 10), the suppression of the NFFLO region seen in Fig. 10(a1) is found to be due to stronger NFFLO pairing fluctuations by more spherical ‘Fermi surface’ edges.

We have shown in Sec. III and in our recent paper Kawamura2022 that the removal of the spatial isotropy of a Fermi gas by a three-dimensional cubic optical lattice is a promising route to stabilize both the thermal equilibrium and non-equilibrium FFLO states against pairing fluctuations. In this section, we examine how these FFLO states are related to each other, by considering a spin-imbalanced driven-dissipative lattice Fermi gas. As shown in Fig. 8, once the NFFLO state is stabilized by the optical lattice, the essential behavior of the phase transition temperature $T_{\rm env}^{\rm c}$ can be captured by the simpler NMF theory at $n=0.3$. Keeping this in mind, in this section, we treat the spin-imbalanced system at this filling fraction within the NMF scheme.

Figure 12(a) shows the phase diagram of a driven-dissipative lattice Fermi gas in the non-equilibrium steady state, with respect to the environmental temperature $T_{\rm env}$, the chemical potential bias $\delta\mu$, and the fictitious magnetic field $h$ to adjust the spin imbalance of the main system. In this figure, the $T_{\rm env}$-$\delta\mu$ plane at $h=0$ describes the spin-balanced non-equilibrium steady state discussed in Sec. III, where the two Fermi surface edges produced by the two reservoirs lead to the NFFLO phase transition in the region of large chemical potential bias $\delta\mu$. On the other hand, the $T_{\rm env}$-$h$ plane at $\delta\mu=0$ corresponds to the spin-imbalanced thermal equilibrium state, where the Zeeman splitting of $\uparrow$-spin and $\downarrow$-spin Fermi surfaces brings about the ordinary FFLO phase transition in the region of high magnetic field $h$ Fulde1964 ; Larkin1964 ; Takada1969 ; Matsuda2007 . Except for these limits, the system with $\delta\mu\neq 0$ and $h\neq 0$ has four Fermi surface like edges in the momentum distribution of Fermi atoms, as illustrated in Fig. 1(c). To be more correct, for example, the two edges at ${\rm FS}^{\rm L}$ and ${\rm FS}^{\rm R}$ in Fig. 11 ($\delta\mu\neq 0$ and $h=0$), respectively, split into ${\rm FS}_{\uparrow}^{\rm L}$ and ${\rm FS}_{\downarrow}^{\rm L}$, and ${\rm FS}_{\uparrow}^{\rm R}$ and ${\rm FS}_{\downarrow}^{\rm R}$, as shown in Fig. 13(a1). In what follows, we simply call these four edges Fermi surfaces, unless any confusion may occur.

We first focus on the phase diagram at $T_{\rm env}=0$, which is explicitly shown in Fig. 12(b). To see the role of the four Fermi surfaces ${\rm FS}_{\sigma=\uparrow,\downarrow}^{\rm L,R}$, we plot in Fig. 13(b) the inverse $\Gamma^{\rm R}(\bm{q}=(q_{x},0,0),\nu=2\mu)^{-1}$ of the retarded particle-particle scattering matrix at the phase boundaries (A1)-(A4) in Fig. 12(b) note.QFF . While $-{\rm Re}\big{[}\Gamma^{\rm R}(\bm{q},\nu=2\mu)^{-1}\big{]}$ has a single minimum in the spin-balanced case ($h=0$), it has two dips in the presence of spin imbalance ($h\neq 0$), which physically means the enhancement of pairing fluctuations around these dip momenta. (The same enhancement can also be seen in the $-q_{x}$ direction, as well as $\pm q_{y}$ and $\pm q_{z}$ directions because of the four-fold rotational symmetry of the cubic lattice.)