Proposed Fermi-surface reservoir-engineering and application to realizing unconventional Fermi superfluids in a driven-dissipative non-equilibrium Fermi gas

The model driven-dissipative two-component Fermi gas shown in Fig. 1(a) is described by the Hamiltonian,

$$H=H_{\rm sys}+H_{\rm env}+H_{\rm mix},$$

(1)

where each term has the form, in the Nambu representation Schrieffer ; Nambu1960 ,

	$\displaystyle H_{\rm sys}$	$\displaystyle=-\int d\bm{r}\Psi^{\dagger}(\bm{r}){\nabla_{\bm{r}}^{2}\over 2m}\tau_{3}\Psi(\bm{r})$
		$\displaystyle\hskip 11.38092pt-U\int d\bm{r}\Psi^{\dagger}({\bm{r}})\tau_{+}\Psi({\bm{r}})\Psi^{\dagger}({\bm{r}})\tau_{-}\Psi({\bm{r}}),$		(2)
	$\displaystyle H_{\rm env}$	$\displaystyle=\sum_{\alpha={\rm L,R}}\int d{\bm{R}}\Phi_{\alpha}^{\dagger}(\bm{R})\left[-{\nabla_{\bm{R}}^{2}\over 2m}-\omega_{\alpha}\right]\tau_{3}\Phi_{\alpha}(\bm{R}),$		(3)
	$\displaystyle H_{\rm mix}$	$\displaystyle=\sum_{\alpha={\rm L},{\rm R}}\sum_{i=1}^{N_{\rm t}}\left[\Lambda_{\alpha}\Phi_{\alpha}^{\dagger}(\bm{R}^{\alpha}_{i})e^{i\tau_{3}\mu_{\alpha}t}\tau_{3}\Psi(\bm{r}^{\alpha}_{i})+{\rm H.c.}\right].$		(4)

In these Hamiltonians,

	$\displaystyle\Psi({\bm{r}})$	$\displaystyle=\left(\begin{array}[]{c}\psi_{\uparrow}({\bm{r}})\\ \psi_{\downarrow}^{\dagger}({\bm{r}})\\ \end{array}\right),$		(7)
	$\displaystyle\Phi_{\alpha}({\bm{R}})$	$\displaystyle=\left(\begin{array}[]{c}\phi_{\alpha,\uparrow}({\bm{R}})\\ \phi_{\alpha,\downarrow}^{\dagger}({\bm{R}})\end{array}\right),$		(10)

represent the two-component Nambu fields describing fermions in the main system and the $\alpha={\rm L,R}$ reservoir, respectively (where $\psi_{\sigma}({\bm{r}})$ ($\phi_{\alpha,\sigma}({\bm{R}})$) is the annihilation operator of a fermion with pseudo-spin $\sigma=\uparrow,\downarrow$ and particle mass $m$ in the main system ($\alpha$-reservoir)). The corresponding Pauli matrices $\tau_{i}$ ($i=1,2,3$), as well as $\tau_{\pm}=[\tau_{1}\pm i\tau_{2}]/2$, act on the particle-hole space.

$H_{\rm sys}$ in Eq. (2) describes the main system in Fig. 1(a), consisting of a two-component Fermi gas with an $s$-wave pairing interaction $-U~{}(<0)$. Because $H_{\rm sys}$ involves the ultraviolet divergence, as usual in cold Fermi gas physics, we remove this singularity by measuring the interaction strength in terms of the $s$-wave scattering length $a_{s}$ Randeria1995 . It is related to the pairing interaction $-U$ as

$${4\pi a_{s}\over m}=-{U\over\displaystyle 1-U\sum_{\bm{p}}{1\over 2\varepsilon_{\bm{p}}}},$$

(11)

where $\varepsilon_{\bm{p}}={\bm{p}}^{2}/(2m)$ is the kinetic energy of a Fermi atom with an atomic $m$. In this paper, we only deal with the weak-coupling regime, by setting $(p_{\rm F}a_{s})^{-1}=-1$. Here, $p_{\rm F}=\sqrt{2m\mu}$, where $\mu\equiv[\mu_{\rm R}+\mu_{\rm L}]/2~{}(>0)$ is the averaged chemical potential between the two reservoirs (where $\mu_{\alpha}$ is the chemical potential of the $\alpha$-reservoir in Fig. 1(a)).

The two reservoirs ($\alpha={\rm L,R}$) are described by the Hamiltonian $H_{\rm env}$ in Eq. (3). As schematically shown in Fig. 3, $-\omega_{\alpha={\rm L},{\rm R}}$ in Eq. (3) gives the bottom of the energy band in the $\alpha$-reservoir, when the energy is measured from the bottom ($\varepsilon_{\bm{p}=0}=0$) of the energy band in the main system. We assume that both reservoirs are huge compared to the main system (which is satisfied by taking $\omega_{\alpha}$ to be sufficiently large compared to $\mu$), and are always in the thermal equilibrium state at the common environment temperature $T_{\rm env}$ noteT . The particle occupation in each reservoir obeys the ordinary Fermi distribution function at $T_{\rm env}$,

$$f(\omega)=\frac{1}{e^{\omega/T_{\rm env}}+1}.$$

(12)

The coupling between the main system and the reservoirs is described by $H_{\rm mix}$ in Eq. (4) with the coupling strength $\Lambda_{\alpha={\rm L},{\rm R}}$. For simplicity, we set $\Lambda_{\rm L}=\Lambda_{\rm R}\equiv\Lambda$ in what follows. The particle tunneling is assumed to occur between randomly distributing spatial positions $\bm{R}_{i}^{\alpha}$ in the $\alpha$-reservoir and ${\bm{r}}^{\alpha}_{i}$ in the main system ($i=1,\cdots N_{\rm t}\gg 1$). As we will see later (see Eq. (40) and the paragraphs nearby), the random tunneling points introduced here assures the supply/decay of the particle from/to the reservoir to be uniform noteRandom . This mimics the experimental situations that we propose below, where the injection of particles are also approximately uniform (to be discussed later). This phenomenological modeling parameters $\Lambda$ and $N_{t}$ would only appear in the expression of the linewidth $\gamma$ [see Eq. (46)], which we interpret as a fitting parameter that would be determined experimentally Hanai2016 ; Hanai2017 ; Hanai2018 .

In Eq. (4), the factor $\exp(i\tau_{3}\mu_{\alpha}t)$ describes the situation that the energy band in the $\alpha$-reservoir is filled up to $\mu_{\alpha}$ (at $T_{\rm env}=0$), as schematically shown in Fig. 3 Kawamura2020 ; Kawamura2020_JLTP . By imposing the chemical-potential bias $\mu_{\rm L}=\mu+\delta\mu$ and $\mu_{\rm R}=\mu-\delta\mu$ ($\delta\mu>0$) between the two reservoirs, we realize the non-equilibrium steady state in the main system. In this paper, we fix the value of the average chemical potential $\mu$, and tune the non-equilibrium situation of the main system by adjusting the tunneling matrix element $\Lambda$, as well as the chemical-potential bias $\delta\mu$ noteP .

At the end of this subsection, we comment on this model and its feasibility in experiments. We expect the proposed setup can experimentally be realized, both in ultracold Fermi systems and electron systems. In the former systems, several groups have succeeded in splitting a cloud of a trapped Fermi gas into two Brantut2012 ; Krinner2015 ; Husmann2015 ; Krinner2017 ; Husmann_thesis ; Kanasz-Nagy2016 . This system is known to be well described by two systems coupled via the tunneling between them Uchino2017 ; Yao2018 ; Sekino2020 ; Furutani2020 , similarly to our Hamiltonian (Eq. (1)). Thus, we strongly believe that our system composed of three systems (the main system and the left and right reservoir) can be realized by extending this two-terminal setup into three, and applying a magnetic field to tune the strength of a pairing interaction associated with a Feshbach resonance only to the main system by using the techniques developed in Bauer2009 ; Yamazaki2010 ; Fu2013 ; Jagannathan2016 ; Arunkumar2018 ; Arunkumar2019 . In this setup, the tunneling process between the system and reservoirs would be overwhelmingly complicated: When a particle is injected from the reservoir to the system, they relax, decohere, and diffuse via inelastic collision processes. However, since this process seems to occur very quickly in the two-terminal configuration in cold atomic systems, the pumping can be regarded as being effectively uniform (which justifies our modeling with random tunneling points). This is supported by the verification of the Landauer formula Krinner2015 and the a.c. and d.c. Josephson effects Levy2007 , which uses the above assumption of fast dissipation over space. In Appendix C, we further argue that it is possible to achieve the parameter regime necessary to realize the unconventional states (i.e. $\gamma/\mu<0.01$) within the current experimental techniques.

In electron systems, a metal (that turns into a superconductor at low temperature) under a strong voltage bias ($eV\gg k_{\rm B}T$) is also a promising experimental setup to realize our proposed non-equilibrium FF-type superconducting state. In fact, a non-equilibrium quasiparticle distribution with the two-step structure analogous to the one depicted in Fig. 1(a) has been observed in voltage-biased mesoscopic wires Pothier1997 ; Poither1997_2 ; Anthore2003 , as well as in carbon nanotubes Chen2009 . The non-equilibrium quasiparticle distribution in a voltage-biased mesoscopic superconducting wire has also been investigated theoretically based on the quasiclassical theory, which has shown that a non-equilibrium distribution with the two-step structure is realized near the center of the wire, when the wire is sufficiently long compared to the superconducting coherence length Keizer2006 ; Vercruyssen2012 . Thus, we expect that the exotic superconducting state induced by the non-equilibrium quasiparticle distribution can also be observed in voltage-biased mesoscopic superconducting wires connected to normal-metal electrodes Boogaard2004 ; Li2011 ; Vercruyssen2012 or thin superconducting films sandwiched between electrodes Ohtomo2004 ; Reyren2007 ; Ueno2008 . In Appendix C, we provide further arguments on the feasibility of the realization of our proposal.

To deal with the non-equilibrium superfluid state in the main system, we conveniently employ the $4\times 4$ matrix Nambu-Keldysh Green’s function Rammer2007 ; Zagoskin2014 , given by

$$\hat{\mathcal{G}}(1,2)=\left(\begin{array}[]{cc}\mathcal{G}^{\rm R}(1,2)&\mathcal{G}^{\rm K}(1,2)\\ 0&\mathcal{G}^{\rm A}(1,2)\end{array}\right)=\left(\begin{array}[]{cc}-i\Theta(t_{1}-t_{2})\langle[\Psi(1)\ \raise 1.29167pt\hbox{$\diamond$}\kern-3.99994pt\lower 3.01385pt\hbox{$,$}\ \Psi^{\dagger}(2)]_{+}\rangle&-i\langle[\Psi(1)\ \raise 1.29167pt\hbox{$\diamond$}\kern-3.99994pt\lower 3.01385pt\hbox{$,$}\ \Psi^{\dagger}(2)]_{-}\rangle\\ 0&i\Theta(t_{2}-t_{1})\langle[\Psi(1)\ \raise 1.29167pt\hbox{$\diamond$}\kern-3.99994pt\lower 3.01385pt\hbox{$,$}\ \Psi^{\dagger}(2)]_{+}\rangle\end{array}\right),$$

(13)

where $\Theta(t)$ is the step function, and the abbreviated notations $1\equiv({\bm{r}}_{1},t_{1})$ and $2\equiv({\bm{r}}_{2},t_{2})$ are used. In Eq. (13), $\mathcal{G}^{\rm R}$, $\mathcal{G}^{\rm A}$, and $\mathcal{G}^{\rm K}$ are the $2\times 2$ matrix retarded, advanced, and Keldysh Green’s functions in the two-component Nambu space, respectively. In Eq. (13),

$$[\Psi(1)\ \raise 1.29167pt\hbox{$\diamond$}\kern-3.99994pt\lower 3.01385pt\hbox{$,$}\ \Psi^{\dagger}(2)]_{\pm}=\Psi(1)\diamond\Psi^{\dagger}(2)\pm\Psi^{\dagger}(2)\diamond\Psi(1),$$

(14)

where “$\diamond$” denotes the operation,

	$\displaystyle\Psi(1)\diamond\Psi^{\dagger}(2)=\left(\begin{array}[]{cc}\psi_{\uparrow}(1)\psi^{\dagger}_{\uparrow}(2)&\psi_{\uparrow}(1)\psi_{\downarrow}(2)\\ \psi_{\downarrow}^{\dagger}(1)\psi_{\uparrow}^{\dagger}(2)&\psi_{\downarrow}^{\dagger}(1)\psi_{\downarrow}(2)\end{array}\right),$		(17)
	$\displaystyle\Psi^{\dagger}(2)\diamond\Psi(1)=\left(\begin{array}[]{cc}\psi^{\dagger}_{\uparrow}(2)\psi_{\uparrow}(1)&\psi_{\downarrow}(2)\psi_{\uparrow}(1)\\ \psi^{\dagger}_{\uparrow}(2)\psi^{\dagger}_{\downarrow}(1)&\psi_{\downarrow}(2)\psi^{\dagger}_{\downarrow}(1)\end{array}\right).$		(20)

For later convenience, we also introduce the $4\times 4$ matrix lesser Green’s function $\mathcal{G}^{<}(1,2)$, which is related to $\mathcal{G}^{\rm R,K,A}(1,2)$ as

	$\displaystyle\mathcal{G}^{<}(1,2)$	$\displaystyle=i\braket{\Psi^{\dagger}(2)\diamond\Psi(1)}$
		$\displaystyle=\frac{1}{2}\big{[}\mathcal{G}^{\rm K}(1,2)-\mathcal{G}^{\rm R}(1,2)+\mathcal{G}^{\rm A}(1,2)\big{]}.$		(21)

We briefly note that the diagonal and off-diagonal components of $\mathcal{G}^{<}$ are related to the particle density and the pair amplitude, respectively.

In the Nambu-Keldysh scheme, effects of the pairing interaction $-U$ and the system-reservoir couplings $\Lambda_{\alpha={\rm L,R}}$ can be summarized by the $4\times 4$ matrix self-energy correction,

$$\hat{\Sigma}(1,2)=\left(\begin{array}[]{cc}\Sigma^{\rm R}(1,2)&\Sigma^{\rm K}(1,2)\\[4.0pt] 0&\Sigma^{\rm A}(1,2)\end{array}\right)=\hat{\Sigma}_{\rm int}(1,2)+\hat{\Sigma}_{\rm env}(1,2),$$

(22)

which appears in the non-equilibrium Nambu-Keldysh Dyson equation Rammer2007 ; Zagoskin2014 ,

$$\hat{\mathcal{G}}(1,2)=\hat{\mathcal{G}}_{0}(1,2)+\big{[}\hat{\mathcal{G}}_{0}\circ\hat{\Sigma}\circ\hat{\mathcal{G}}\big{]}(1,2).$$

(23)

Here,

	$\displaystyle\big{[}A\circ B\big{]}(1,2)$	$\displaystyle=\scalebox{0.92}{$\displaystyle\int d\bm{r}_{3}\int_{-\infty}^{\infty}dt_{3}A(\bm{r}_{1},t_{1},\bm{r}_{3},t_{3})B(\bm{r}_{3},t_{3},\bm{r}_{2},t_{2})$}$
		$\displaystyle=\int d3A(1,3)B(3,2),$		(24)

and

	$\displaystyle\hat{\mathcal{G}}_{0}(1,2)$	$\displaystyle=\sum_{\bm{p}}\int{d\omega\over 2\pi}e^{i{\bm{p}}\cdot({\bm{r}}_{1}-{\bm{r}}_{2})-i\omega(t_{1}-t_{2})}\hat{\mathcal{G}}_{0}({\bm{p}},\omega)$
		$\displaystyle=\sum_{\bm{p}}\int{d\omega\over 2\pi}e^{i{\bm{p}}\cdot({\bm{r}}_{1}-{\bm{r}}_{2})-i\omega(t_{1}-t_{2})}$
		$\displaystyle\hskip 11.38092pt\times\left(\begin{array}[]{cc}\frac{1}{\omega_{+}-\varepsilon_{\bm{p}}\tau_{3}}&-2\pi i\delta(\omega-\varepsilon_{\bm{p}})\big{[}1-2f_{\rm ini}(\omega)\big{]}\\ 0&\frac{1}{\omega_{-}-\varepsilon_{\bm{p}}\tau_{3}}\end{array}\right)$		(27)

is the bare Green’s function in the initial thermal equilibrium state at $t=-\infty$, where the system-reservoir couplings $\Lambda_{\alpha={\rm L,R}}$, as well as the pairing interaction $-U$, were absent. (Note that $\hat{\mathcal{G}}_{0}(1,2)$ only depends on the relative coordinate as $\hat{\mathcal{G}}_{0}(1,2)=\hat{\mathcal{G}}_{0}(1-2)$.) In Eq. (27), $\omega_{\pm}=\omega\pm i\delta$ (where $\delta$ is an infinitesimally small positive number), and $f_{\rm ini}(\omega)=1/[e^{\omega/T_{\rm ini}}+1]$ is the Fermi distribution function with $T_{\rm ini}$ being the initial temperature of the main system at $t=-\infty$. Although the Dyson equation (23) looks like depending on the initial state through $\hat{\mathcal{G}}_{0}(1,2)$, we will soon find that the dressed Green’s function $\hat{\mathcal{G}}(1,2)$ in the final non-equilibrium steady state, which we are interested in, actually loses the initial memory Kawamura2020 ; Hanai2016 ; Hanai2017 ; Hanai2018 ; Kawamura2020_JLTP ; Szymanska2006 ; Szymanska2007 ; Yamaguchi2012 ; Yamaguchi2015 .

In Eq. (22), $\hat{\Sigma}_{\rm int}$ and $\hat{\Sigma}_{\rm env}$ describe effects of the pairing interaction $-U$ and the system-reservoir couplings $\Lambda_{\alpha={\rm L,R}}~{}(=\Lambda)$, respectively. In the mean-field BCS approximation Hanai2016 ; Hanai2017 ; Hanai2018 , the former is diagrammatically drawn as Fig. 4(a), which gives

	$\displaystyle\hat{\Sigma}_{\rm int}(1,2)=iU\sum_{s=\pm}\sum_{\alpha=1,2}\big{(}\tau_{s}\otimes\eta^{+}_{\alpha}\big{)}{\rm Tr}_{\rm N}{\rm Tr}_{\rm K}\big{[}\big{(}\tau_{-s}\otimes\eta^{-}_{\alpha}\big{)}\hat{\mathcal{G}}(1,2)\big{]}\delta(1-2)$
	$\displaystyle=\frac{iU}{2}\sum_{s=\pm}\left(\begin{array}[]{cc}\tau_{s}{\rm Tr}_{\rm N}\big{[}\tau_{-s}\mathcal{G}^{\rm K}(1,2)\big{]}&\tau_{s}{\rm Tr}_{\rm N}\big{[}\tau_{-s}\mathcal{G}^{\rm R}(1,2)+\tau_{-s}\mathcal{G}^{\rm A}(1,2)\big{]}\\[4.0pt] \tau_{s}{\rm Tr}_{\rm N}\big{[}\tau_{-s}\mathcal{G}^{\rm R}(1,2)+\tau_{-s}\mathcal{G}^{\rm A}(1,2)\big{]}&\tau_{s}{\rm Tr}_{\rm N}\big{[}\tau_{-s}\mathcal{G}^{\rm K}(1,2)\big{]}\end{array}\right)\delta(1-2).$		(30)

Here,

$$\eta_{\alpha}^{+}=\frac{1}{\sqrt{2}}\sigma_{2-\alpha},{\hskip 14.22636pt}\eta_{\alpha}^{-}=\frac{1}{\sqrt{2}}\sigma_{\alpha-1}$$

(31)

are vertex matrices Kawamura2020 , where $\sigma_{i=1,2,3}$ are the Pauli matrices acting on the Keldysh space. ${\rm Tr}_{\rm N}$ and ${\rm Tr}_{\rm K}$ stand for taking the trace over the Nambu and the Keldysh space, respectively.

We introduce the superfluid order parameter $\Delta(\bm{r}_{1},t_{1})$, which is related to the off-diagonal component of $\mathcal{G}^{\rm K}$ and $\mathcal{G}^{<}$ as

$$\Delta(\bm{r}_{1},t_{1})\equiv U\braket{\psi_{\downarrow}(1)\psi_{\uparrow}(1)}=-\frac{iU}{2}\mathcal{G}^{\rm K}(1,1)_{12}=-iU\mathcal{G}^{<}(1,1)_{12}.$$

(32)

Then, Eq. (30) can be simply written as

$\displaystyle\hat{\Sigma}_{\rm int}(1,2)=$

$\displaystyle\left(\begin{array}[]{cc}-\Delta(1)\tau_{+}-\Delta^{*}(1)\tau_{-}&0\\[4.0pt] 0&-\Delta(1)\tau_{+}-\Delta^{*}(1)\tau_{-}\end{array}\right)\delta(1-2).$

(35)

We note that the off-diagonal components of $\hat{\Sigma}_{\rm int}(1,2)$ identically vanish, because $\mathcal{G}^{\rm R(A)}(1,1)_{12}=\mathcal{G}^{\rm R(A)}(1,1)_{21}=0$.

For the self-energy correction $\hat{\Sigma}_{\rm env}$ in Eq. (22), we take into account the system-reservoir couplings within the second-order Born approximation, as diagrammatically shown in Fig. 4(b). Evaluation of this diagram gives

	$\displaystyle\hat{\Sigma}_{\rm env}({\bm{p}},{\bm{p}}^{\prime},t_{1},t_{2})$	$\displaystyle=$	$\displaystyle\int d{\bm{r}}_{1}\int d{\bm{r}}_{2}\hat{\Sigma}_{\rm env}(1,2)e^{-i({\bm{p}}\cdot{\bm{r}}_{1}+{\bm{p}}^{\prime}\cdot{\bm{r}}_{2})}$		(36)
		$\displaystyle=$	$\displaystyle|\Lambda|^{2}\sum_{\alpha={\rm L,R}}\sum_{i,j}^{N_{\rm t}}\hat{\mathcal{D}}_{\alpha}(\bm{R}_{i}^{\alpha}-\bm{R}_{j}^{\alpha},t_{1}-t_{2})e^{-i\tau_{3}\mu_{\alpha}(t_{1}-t_{2})}e^{-i({\bm{p}}\cdot{\bm{r}}_{i}+{\bm{p}}^{\prime}\cdot{\bm{r}}_{j})},$

where

	$\displaystyle\hat{\mathcal{D}}_{\alpha}(\bm{R}_{i}^{\alpha}-\bm{R}_{j}^{\alpha},t_{1}-t_{2})$	$\displaystyle=\sum_{\bm{q}}\int{d\omega\over 2\pi}e^{i{\bm{q}}\cdot({\bm{R}}_{i}^{\alpha}-{\bm{R}}_{j}^{\alpha})-i\omega(t_{1}-t_{2})}\hat{\mathcal{D}}_{0}({\bm{q}},\omega)$
		$\displaystyle=\sum_{\bm{q}}\int{d\omega\over 2\pi}e^{i{\bm{q}}\cdot({\bm{R}}_{i}^{\alpha}-{\bm{R}}_{j}^{\alpha})-i\omega(t_{1}-t_{2})}\left(\begin{array}[]{cc}\frac{1}{\omega_{+}-\xi^{\alpha}_{\bm{q}}\tau_{3}}&-2\pi i\delta(\omega-\xi^{\alpha}_{\bm{q}})\tanh\left(\frac{\omega}{2T_{\rm env}}\right)\\ 0&\frac{1}{\omega_{-}-\xi^{\alpha}_{\bm{q}}\tau_{3}}\end{array}\right)$		(39)

is the non-interacting Nambu-Keldysh Green’s function in the $\alpha$-reservoir, with $\xi_{\bm{q}}^{\alpha}=\varepsilon_{\bm{q}}-\omega_{\alpha}$ being the kinetic energy measured from the bottom energy $-\omega_{\alpha}$ of the band in the $\alpha$-reservoir. (Note that the reservoirs are assumed to be in thermal equilibrium.) When one takes the spatial averages over the randomly distributing tunneling positions $\bm{R}_{i}^{\alpha}$ and $\bm{r}_{i}^{\alpha}$ in Eq. (36), the resulting self-energy $\braket{{\hat{\Sigma}}_{\rm env}({\bm{p}},{\bm{p}}^{\prime},t_{1},t_{2})}_{\rm av}$ recovers its the translational invariance as Hanai2016 ; Hanai2017 ; Hanai2018 ; Kawamura2020

$$\braket{{\hat{\Sigma}}_{\rm env}({\bm{p}},{\bm{p}}^{\prime},t_{1},t_{2})}_{\rm av}=\hat{\Sigma}_{\rm env}({\bm{p}},t_{1},t_{2})\delta_{{\bm{p}},{\bm{p}}^{\prime}},$$

(40)

where

$\displaystyle\hat{\Sigma}_{\rm env}({\bm{p}},t_{1},t_{2})=N_{t}|\Lambda|^{2}\sum_{{\bm{q}},\alpha={\rm L,R}}\hat{\mathcal{D}}_{\alpha}({\bm{q}},t_{1}-t_{2})e^{-i\mu_{\alpha}(t_{1}-t_{2})\tau_{3}}.$

(41)

Carrying out the Fourier transformation with respect to the relative time $t_{1}-t_{2}$, we have

$$\hat{\Sigma}_{\rm env}(\bm{p},\omega)=N_{\rm t}|\Lambda|^{2}\sum_{{\bm{q}},\alpha={\rm L,R}}\hat{\mathcal{D}}_{\alpha}\big{(}\bm{q},\omega-\mu_{\alpha}\tau_{3}\big{)}.$$

(42)

For simplicity, we employ the so-called wide-band limit approximation Stefanucci2013 , that is, we assume white reservoirs with the constant density of states $\rho_{\alpha}(\omega)\equiv\rho$. This approximation is justified when the reservoirs are so huge that the energy dependence of the reservoir density of states around the Fermi level can be ignored Kawamura2020 , which is just the situation we are considering ($\mu\ll\omega_{\alpha}$). Then, replacing $\bm{q}$ summation in Eq. (42) by the $\xi^{\alpha}$ integration, one has

	$\displaystyle\hat{\Sigma}_{\rm env}(\bm{p},\omega)$
	$\displaystyle=\scalebox{0.93}{$\displaystyle\left(\begin{array}[]{cc}-2i\gamma\tau_{0}&-2i\gamma\left[\tanh\left(\frac{\omega-\tau_{3}\mu_{\rm L}}{2T_{\rm env}}\right)+\tanh\left(\frac{\omega-\tau_{3}\mu_{\rm R}}{2T_{\rm env}}\right)\right]\tau_{0}\\ 0&2i\gamma\tau_{0}\end{array}\right).$}$		(45)

Here,

$$\gamma=\pi N_{\rm t}\rho|\Lambda|^{2}$$

(46)

is the quasi-particle damping rate, and $\tau_{0}$ is the $2\times 2$ unit matrix acting on the particle-hole Nambu space.

In this paper, we explore stable non-equilibrium superfluid steady states, having the following type of the order parameter:

$$\Delta(\bm{r},t)=\Delta_{0}e^{i\bm{Q}\cdot\bm{r}}e^{-2i\mu t}.$$

(47)

Without loss of generality, we take $\Delta_{0}>0$. When ${\bm{Q}}=0$, Eq. (47) describes the BCS-type uniform superfluid, which has been discussed in exciton(-polariton) systems Szymanska2006 ; Hanai2016 ; Hanai2017 ; Hanai2018 ; Hanai2019 ; Yamaguchi2012 ; Yamaguchi2015 . When ${\bm{Q}}\neq 0$, Eq. (47) has the same form as the order parameter in the Fulde-Ferrell (FF) superfluid state, discussed in superconductivity under an external magnetic field Fulde1964 ; Larkin1964 ; Takada1969 ; Shimahara1994 ; Matsuda2007 , as well as in a spin-polarized Fermi gas Hu2006 ; Parish2007 ; Liao2010 ; Chevy2010 ; Kinnunen2018 ; Strinati2018 . Although the Larkin-Ovchinnikov type solution Larkin1964 , $\Delta(\bm{r},t)=\Delta_{0}\cos(\bm{Q}\cdot\bm{r})e^{-2i\mu t}$, is also conceivable in our model, leaving this possibility as our future study, we only deal with the FF-type solution in this paper. Regarding this, we emphasize that the main system has no spin imbalance.

To treat the superfluid order parameter $\Delta(\bm{r},t)$ in Eq. (47), it is convenient to formally remove time and spatial dependence from it, which is achieved by employing the following gauge transformation Rammer2007 :

	$\displaystyle\hat{\tilde{\mathcal{G}}}(1,2)$	$\displaystyle=\left(\begin{array}[]{cc}{\tilde{\mathcal{G}}}^{\rm R}(1,2)&{\tilde{\mathcal{G}}}^{\rm K}(1,2)\\ 0&{\tilde{\mathcal{G}}}^{\rm A}(1,2)\end{array}\right)$		(50)
		$\displaystyle\equiv e^{-i\chi(1)\tau_{3}\otimes\sigma_{0}}\hat{\mathcal{G}}(1,2)e^{i\chi(2)\tau_{3}\otimes\sigma_{0}},$		(51)
	$\displaystyle\hat{\tilde{\Sigma}}(1,2)$	$\displaystyle=\left(\begin{array}[]{cc}{\tilde{\Sigma}}^{\rm R}(1,2)&{\tilde{\Sigma}}^{\rm K}(1,2)\\ 0&{\tilde{\Sigma}}^{\rm A}(1,2)\end{array}\right)$		(54)
		$\displaystyle\equiv e^{-i\chi(1)\tau_{3}\otimes\sigma_{0}}\hat{\Sigma}(1,2)e^{i\chi(2)\tau_{3}\otimes\sigma_{0}},$		(55)

where

$$\chi(\bm{r},t)={1\over 2}\bm{Q}\cdot\bm{r}-\mu t,$$

(56)

and $\sigma_{0}$ is the $2\times 2$ unit matrix acting on the Keldysh space. The Nambu-Keldysh Dyson equation (23) after this manipulation is given by

$$\hat{\tilde{\mathcal{G}}}(1,2)=\hat{\tilde{\mathcal{G}}}_{0}(1,2)+\big{[}\hat{\tilde{\mathcal{G}}}_{0}\circ\hat{\tilde{\Sigma}}\circ\hat{\tilde{\mathcal{G}}}\big{]}(1,2).$$

(57)

Here,

$$\hat{\tilde{\mathcal{G}}}_{0}(1,2)=\hat{\tilde{\mathcal{G}}}_{0}(1-2)=e^{i[\mu(t_{1}-t_{2})-{\bm{Q}\over 2}\cdot(\bm{r}_{1}-\bm{r}_{2})]\tau_{3}}\hat{\mathcal{G}}_{0}(1-2).$$

(58)

is the gauge-transformed bare Green’s function. In the energy and momentum space, Eq. (57) has the form,

	$\displaystyle\hat{\tilde{\mathcal{G}}}(\bm{p},\omega)$	$\displaystyle=\left(\begin{array}[]{cc}\tilde{\mathcal{G}}^{\rm R}({\bm{p}},\omega)&\tilde{\mathcal{G}}^{\rm K}({\bm{p}},\omega)\\ 0&\tilde{\mathcal{G}}^{\rm A}({\bm{p}},\omega)\end{array}\right)$		(61)
		$\displaystyle=\hat{\tilde{\mathcal{G}}}_{0}(\bm{p},\omega)+\hat{\tilde{\mathcal{G}}}_{0}(\bm{p},\omega)\hat{\tilde{\Sigma}}(\bm{p},\omega)\hat{\tilde{\mathcal{G}}}(\bm{p},\omega),$		(62)

where

	$\displaystyle\hat{\tilde{\mathcal{G}}}_{0}(\bm{p},\omega)$	$\displaystyle=\int d{\bm{r}}\int dte^{i[\omega t-{\bm{p}}\cdot{\bm{r}}]}\hat{\tilde{\mathcal{G}}}_{0}({\bm{r}},t)=\hat{\mathcal{G}}_{0}(\bm{p}+{\bm{Q}\over 2}\tau_{3},\omega+\mu\tau_{3})$
		$\displaystyle=\left(\begin{array}[]{cc}\frac{1}{\omega_{+}-\xi_{\bm{p}+(\bm{Q}/2)\tau_{3}}\tau_{3}}&-2\pi i\delta(\omega-\xi_{\bm{p}+(\bm{Q}/2)\tau_{3}}\tau_{3})\big{[}1-2f_{\rm ini}(\omega+\mu\tau_{3})\big{]}\\ 0&\frac{1}{\omega_{-}-\xi_{\bm{p}+(\bm{Q}/2)\tau_{3}}\tau_{3}}\end{array}\right),$		(65)

with $\xi_{\bm{p}}=\varepsilon_{\bm{p}}-\mu$. The interaction component $\hat{\tilde{\Sigma}}_{\rm int}$ of the self-energy $\hat{\tilde{\Sigma}}=\hat{\tilde{\Sigma}}_{\rm int}+\hat{\tilde{\Sigma}}_{\rm env}$ in Eq. (57) is obtained from Eq. (35) as

$$\hat{\tilde{\Sigma}}_{\rm int}(1,2)=\hat{\tilde{\Sigma}}_{\rm int}(1-2)=\left(\begin{array}[]{cc}-\Delta_{0}\tau_{1}&0\\ 0&-\Delta_{0}\tau_{1}\end{array}\right)\delta(1-2).$$

(66)

In the energy and momentum space, this self-energy has the form,

$$\hat{\tilde{\Sigma}}_{\rm int}(\bm{p},\omega)=\left(\begin{array}[]{cc}-\Delta_{0}\tau_{1}&0\\ 0&-\Delta_{0}\tau_{1}\end{array}\right).$$

(67)

In the same manner, the self-energy $\hat{\tilde{\Sigma}}_{\rm env}(1,2)$ coming from the system-reservoir couplings can also be obtained from Eq. (45). Thus, it has the form, in the energy and momentum space,

	$\displaystyle\hat{\tilde{\Sigma}}_{\rm env}(\bm{p},\omega)=\hat{\Sigma}_{\rm env}(\bm{p}+{\bm{Q}\over 2}\tau_{3},\omega+\mu\tau_{3})$
	$\displaystyle=\left(\begin{array}[]{cc}-2i\gamma\tau_{0}&-2i\gamma\left[\tanh\left(\frac{\omega-\delta\mu}{2T_{\rm env}}\right)+\tanh\left(\frac{\omega+\delta\mu}{2T_{\rm env}}\right)\right]\tau_{0}\\ 0&2i\gamma\tau_{0}\end{array}\right).$		(70)

The self-energy $\hat{\tilde{\Sigma}}(\bm{p},\omega)$ involved in the Dyson equation (62) is then given by the sum of Eqs. (67) and (70).

Solving the Dyson equation (62), we obtain

	$\displaystyle\tilde{\mathcal{G}}^{\rm R}(\bm{p},\omega)$	$\displaystyle=\sum_{\eta=\pm}\frac{1}{\omega+2i\gamma-E^{0,\eta}_{\bm{p},\bm{Q}}}\Xi^{0,\eta}_{\bm{p},\bm{Q}},$		(71)
	$\displaystyle\tilde{\mathcal{G}}^{\rm A}(\bm{p},\omega)$	$\displaystyle=\sum_{\eta=\pm}\frac{1}{\omega-2i\gamma-E^{0,\eta}_{\bm{p},\bm{Q}}}\Xi^{0,\eta}_{\bm{p},\bm{Q}},$		(72)
	$\displaystyle\tilde{\mathcal{G}}^{\rm K}(\bm{p},\omega)$	$\displaystyle=\sum_{\eta=\pm}\frac{4i\gamma[1-2F(\omega)]}{[\omega-\eta E^{0,\eta}_{\bm{p},\bm{Q}}]^{2}+4\gamma^{2}}\Xi^{0,\eta}_{\bm{p},\bm{Q}},$		(73)

where

	$\displaystyle E^{0,\pm}_{\bm{p},\bm{Q}}=\sqrt{\bigl{(}\xi^{\rm(s)}_{\bm{p},\bm{Q}}\bigr{)}^{2}+\Delta_{0}^{2}}\pm\xi^{\rm(a)}_{\bm{p},\bm{Q}}\equiv E^{0}_{\bm{p},\bm{Q}}\pm\xi^{\rm(a)}_{\bm{p},\bm{Q}},$		(74)
	$\displaystyle\Xi^{0,\pm}_{\bm{p},\bm{Q}}={1\over 2}\left[\tau_{0}\pm{\xi_{{\bm{p}},{\bm{Q}}}^{\rm(s)}\over E^{0}_{\bm{p},\bm{Q}}}\tau_{3}\mp{\Delta_{0}\over E^{0}_{\bm{p},\bm{Q}}}\tau_{1}\right],$		(75)

with $\xi^{\rm(s)}_{\bm{p},\bm{Q}}=[\xi_{\bm{p}+\bm{Q}/2}+\xi_{-\bm{p}+\bm{Q}/2}]/2$, and $\xi^{\rm(a)}_{\bm{p},\bm{Q}}=[\xi_{\bm{p}+\bm{Q}/2}-\xi_{-\bm{p}+\bm{Q}/2}]/2$. In Eq. (73),

$$F(\omega)={1\over 2}[f(\omega+\delta\mu)+f(\omega-\delta\mu)]$$

(76)

works as the non-equilibrium distribution function in the main system (although $T_{\rm env}$ is used in $f(\omega)$, see Eq. (12)). When we set $\omega=\xi_{\bm{p}}$, the resulting momentum distribution $F(\xi_{\bm{p}})$ has two Fermi-surface-like edges at $p_{\rm F1}=\sqrt{2m\mu_{\rm R}}$ and $p_{\rm F2}=\sqrt{2m\mu_{\rm L}}$, as schematically drawn in Fig. 1(a).

The superfluid order parameter $\Delta({\bm{r}},t)$ in Eq. (47) is self-consistently determined from Eq. (32). To evaluate this equation, we note that the gauge-transformed lesser Green’s function $\tilde{\mathcal{G}}^{<}(\bm{p},\omega)$ in the energy and momentum space is obtained from Eqs. (21) and (71)-(73) as

	$\displaystyle\tilde{\mathcal{G}}^{<}(\bm{p},\omega)$	$\displaystyle=\frac{1}{2}\big{[}\tilde{\mathcal{G}}^{\rm K}(\bm{p},\omega)-\tilde{\mathcal{G}}^{\rm R}(\bm{p},\omega)+\tilde{\mathcal{G}}^{\rm A}(\bm{p},\omega)\big{]}$
		$\displaystyle=\sum_{\eta=\pm}\frac{4i\gamma F(\omega)}{[\omega-\eta E^{0,\eta}_{\bm{p},\bm{Q}}]^{2}+4\gamma^{2}}\Xi^{0,\eta}_{\bm{p},\bm{Q}}.$		(77)

Using Eq. (77) in evaluating Eq. (32), we obtain the gap equation,

	$\displaystyle 1$	$\displaystyle=U\sum_{\bm{p}}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}$
		$\displaystyle\hskip 17.07182pt\times\frac{4\gamma\big{[}\omega-\xi^{\rm(a)}_{\bm{p},\bm{Q}}\big{]}\big{[}1-2F(\omega)\big{]}}{\bigl{[}(\omega-E^{0+}_{\bm{p},\bm{Q}})^{2}+4\gamma^{2}\bigr{]}\bigl{[}(\omega+E^{0-}_{\bm{p},\bm{Q}})^{2}+4\gamma^{2}\bigr{]}}.$		(78)

To quickly see the relation to the ordinary BCS gap equation, we take the thermal equilibrium and uniform limit, by setting $(\gamma,\delta\mu,{\bm{Q}})\to(+0,0,0)$. Then, Eq. (78) is reduced to

$$1=U\sum_{\bm{p}}\frac{1}{2E_{\bm{p}}}\tanh\left(\frac{E_{\bm{p}}}{2T_{\rm env}}\right),$$

(79)

where $E_{\bm{p}}=E_{{\bm{p}},{\bm{Q}}=0}^{0}=\sqrt{\xi_{\bm{p}}^{2}+\Delta_{0}^{2}}$. Equation (79) is just the same form as the ordinary BCS gap equation Schrieffer , when one interprets $\mu$ and $T_{\rm env}$ as the Fermi chemical potential and the temperature in the main system, respectively. Thus, Eq. (78) may be viewed as a non-equilibrium extension of the BCS gap equation.

Because $\Delta({\bm{r}},t)$ in Eq. (47) involves two parameters $\Delta_{0}$ and ${\bm{Q}}$, we actually need one more equation to completely fix these parameters note . For this purpose, we impose the vanishing condition for the net current ${\bm{J}}_{\rm net}$ in the main system. In the Nambu-Keldysh formalism, it is given by

$$\bm{J}_{\rm net}=\sum_{\sigma=\uparrow,\downarrow}\sum_{\bm{p}}\bigl{[}\bm{p}+\bm{Q}/2\bigr{]}n_{\bm{p},\sigma}=0,$$

(80)

where the Fermi momentum distribution $n_{\bm{p},\sigma}$ in the pseudo-spin $\sigma$ component is related to the diagonal component of the lesser Green’s function $\tilde{\mathcal{G}}^{<}(\bm{p},\omega)$ in Eq. (77) as

	$\displaystyle n_{\bm{p},\uparrow}$	$\displaystyle=-i\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\tilde{\mathcal{G}}^{<}(\bm{p},\omega)_{11},$		(81)
	$\displaystyle n_{\bm{p},\downarrow}$	$\displaystyle=1-i\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\tilde{\mathcal{G}}^{<}(-\bm{p},\omega)_{22}.$		(82)

We solve the gap equation (78) under the condition in Eq. (80), to self-consistently determine $(\Delta_{0},{\bm{Q}})$ for a given set $(\gamma,\delta\mu,T_{\rm env})$ of environment parameters.

We note that the net current $\bm{J}_{\rm net}$ is the current flowing inside the main system, but not to be confused with the current flowing from the reservoir to the main system (the latter is always nonzero unless we consider the chemical equilibrium case $\mu_{\rm L}=\mu_{\rm R}$). Although we assume here that the former is zero, strictly speaking, there is no a priori way of determining it since it ultimately depends on the choice of a boundary condition noteBloch . Hence, one may also consider more general cases with ${\bm{J}}_{\rm net}\neq 0$. Even in such a case, however, the essential physics remains unchanged: The non-equilibrium FF-type superfluid state with a non-zero net current can be realized as a stable non-equilibrium steady state (unless the current exceeds the Landau velocity, which would destroy all superfluid states). Thus, we choose the simplest case of the vanishing-current boundary condition with Eq. (80) for the sake of brevity.

We also note that the normal state ($\Delta_{0}=\bm{Q}=0$) always satisfies the gap equation (78) and the vanishing current condition in Eq. (80). However, as pointed out in our previous work Kawamura2020 ; Kawamura2020_JLTP , the normal state becomes unstable, when the particle-particle scattering vertex $\chi(\bm{Q},\nu)$ develops a pole at $\nu=2\mu$. Here, ${\bm{Q}}$ and $\nu$ are the center-of-mass momentum and the total energy of two particles participating in the Cooper channel, respectively. Within the random phase approximation in terms of $-U$, one has Kawamura2020 ; Kawamura2020_JLTP , in NESS,

$$\chi(\bm{Q},\nu=2\mu)={-U\over\displaystyle 1-{U\over 4\pi}\sum_{\eta,\zeta=\pm}{1\over\xi^{\rm s}_{\bm{p},\bm{Q}}}{\rm Tan}^{-1}\left(\frac{\xi_{\bm{p},\bm{Q}}^{\eta,\zeta}}{2\gamma}\right)},$$

(83)

where $\xi_{{\bm{p}},{\bm{Q}}}^{\eta,\zeta}=\xi_{\bm{p}+\eta\bm{Q}/2}+\zeta\delta\mu$. We determine the region where the normal state is stable on the phase diagram from the condition that $\chi(\bm{Q},\nu=2\mu)<0$ for any $\bm{Q}$, that is, $\chi(\bm{Q},\nu)$ has no pole at $\nu=2\mu$.