$\eta$ Pairing of Light-Emitting Fermions:
Nonequilibrium Pairing Mechanism at High Temperatures

In a seminal paper Yang (1989), Yang showed that exact eigenstates of the Fermi-Hubbard model can be constructed through the $\eta$-pairing mechanism. Although the Hubbard model in two and three dimensions is not exactly solvable due to the competing nature of hopping and interaction, a hidden $\eta$ symmetry of the model permits $\eta$-pairing eigenstates Yang and Zhang (1990); Pernici (1990). The $\eta$-pairing states exhibit an off-diagonal long-range order Yang (1962), giving a rare example of fermionic superfluidity that is an exact eigenstate. However, since they cannot be the ground state of the Hubbard model except for the limiting case of an infinite attractive interaction Yang (1989); Singh and Scalettar (1991), $\eta$ pairing requires nonequilibrium situations such as adiabatic state preparation Rosch et al. (2008); Kantian et al. (2010), laser irradiation Kaneko et al. (2019, 2020); Werner et al. (2019); Li et al. (2020), periodic driving Kitamura and Aoki (2016); Peronaci et al. (2020); Cook and Clark (2020); Tindall et al. (2021a, b), and tailored dissipative processes Diehl et al. (2008); Kraus et al. (2008); Bernier et al. (2013); Tindall et al. (2019); Zhang and Song (2020a, b).

Control of superfluidity with light lies at the forefront of nonequilibrium quantum many-body physics Fausti et al. (2011); Kaiser et al. (2014); Hu et al. (2014); Mitrano et al. (2016); Cantaluppi et al. (2018); Suzuki et al. (2019); Buzzi et al. (2020); Budden et al. (2020). Since irradiation of light generally heats up a sample, how to create quantum coherence of a superfluid by the light-matter coupling is a highly nontrivial question with several interesting possibilities Kaneko et al. (2019, 2020); Werner et al. (2019); Li et al. (2020); Fausti et al. (2011); Kaiser et al. (2014); Hu et al. (2014); Mitrano et al. (2016); Cantaluppi et al. (2018); Suzuki et al. (2019); Buzzi et al. (2020); Budden et al. (2020). Importantly, a pairing mechanism in nonequilibrium quantum matter may be fundamentally different from that of the BCS theory Bardeen et al. (1957) because of its nonequilibrium nature.

In this Letter, we propose a nonequilibrium $\eta$-pairing mechanism induced by spontaneous emission of light from strongly correlated fermionic atoms. We consider a cold-atom Fermi-Hubbard system in an optical lattice in which a spin-up state of an atom can decay to a spin-down state via a coupling to an excited state that decays by spontaneously emitting a photon, as illustrated in Fig. 1. We show that the nonequilibrium steady states of this light-matter-coupled system are given by the $\eta$-pairing states. Here, the $\eta$-pairing states are stabilized due to double protection of fermion pairs from possible decay processes. First, as shown in Fig. 1(a), the Pauli exclusion principle prohibits on-site pairs (i.e., doublons) from radiative decay Sandner et al. (2011). Second, complete destructive interference between doublon-decay processes prevents doublons from breaking up into single particles, thereby maintaining the $\eta$-pairing states as a coherent steady state. This pairing mechanism is inherently of nonequilibrium nature and makes a marked contrast to the BCS mechanism Bardeen et al. (1957) in which an attractive interaction between fermions and the Fermi-surface instability play essential roles. Remarkably, unlike previous proposals Bernier et al. (2013); Tindall et al. (2019), the $\eta$-pairing steady state is observed even when the dynamics starts from an infinite-temperature initial state, because the entropy of the system can be decreased by dissipative processes of spontaneous emission which is described by non-Hermitian Lindblad operators. This indicates that one can create coherent atom pairs at an initial temperature far above the BCS critical temperature. We propose an experimental implementation of the model with alkaline-earth-like fermionic atoms in an optical lattice and discuss how to detect the $\eta$-pairing state through measurement of a momentum distribution of doublons.

Model.– We generalize a quantum optical master equation Breuer and Petruccione (2007), which has been used to describe the dynamics of noninteracting atoms subject to spontaneous emission of photons, to interacting fermionic atoms Lehmberg (1970); Pichler et al. (2010); Sarkar et al. (2014). We consider spin-$1/2$ fermionic atoms in a $d$-dimensional cubic optical lattice. The Hamiltonian of our system is described by the Hubbard model $H$ when the lattice potential is sufficiently deep and the single-band approximation is valid Esslinger (2010):

where $c_{j\sigma}$ is the annihilation operator of a spin-$\sigma$ fermion at a lattice position $\bm{R}_{j}\in\{0,\cdots,L-1\}^{d}$, $n_{j\sigma}\equiv c_{j\sigma}^{\dagger}c_{j\sigma}$ is the density operator, $t$ is the tunneling amplitude, and $U$ is the on-site interaction strength. Here $\delta$ is the energy difference between the internal spin states, which can be set to zero without loss of generality sup . We assume that atoms tunnel only between nearest-neighbor sites. We impose the periodic boundary condition to ensure the translational invariance of the model, but the following results can be generalized to the case of the open boundary condition.

We consider a situation in which an atom in a spin-up state can decay to a spin-down state via spontaneous emission of a photon [see Fig. 1(a)]. Such a decay process can be induced by weakly hybridizing the spin-up state with an excited level of the atom [see Fig. 1(b)]. A detailed discussion will be given later and in Supplemental Material sup with a concrete experimental implementation. Within the Born-Markov approximation, we can trace out the photon degrees of freedom and obtain the Lindblad-type quantum master equation that describes the time evolution of the density matrix $\rho$ of the atomic gas as Breuer and Petruccione (2007)

where $\tau$ is the time, $\mathcal{L}$ is the Liouvillian superoperator, and the Lindblad operator $L_{j}=\sqrt{\gamma}c_{j\downarrow}^{\dagger}c_{j\uparrow}$ describes a spontaneous emission process with rate $\gamma>0$ at site $j$.

$\eta$-pairing steady state.– The Hubbard Hamiltonian (1) possesses the $\eta$ SU(2) symmetry generated by $\eta^{x}=(\eta^{+}+\eta^{-})/2,\eta^{y}=(\eta^{+}-\eta^{-})/2i$, and $\eta^{z}=\frac{1}{2}\sum_{j}(n_{j\uparrow}+n_{j\downarrow}-1)$, where $\eta^{+}=\sum_{j}e^{i\bm{Q}\cdot\bm{R}_{j}}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}$ creates a doublon with crystal momentum $\bm{Q}=(\pi,\cdots,\pi)$, and $\eta^{-}=(\eta^{+})^{\dagger}$ Yang (1989); Yang and Zhang (1990). Thus, if a state $\ket{\Psi}$ is an eigenstate of $H$, so is $\eta^{+}\ket{\Psi}$. For example, Yang’s $\eta$-pairing state $\ket{\psi_{N}}\equiv(\eta^{+})^{N/2}\ket{0}$, where $\ket{0}$ is the vacuum of fermions, is an $N$-particle eigenstate of $H$ Yang (1989). The $\eta$-pairing state can be regarded as a condensate of doublons at momentum $\bm{Q}$ and exhibits an off-diagonal long-range order Yang (1989, 1962). We also introduce a family of eigenstates $\ket{\psi_{2N_{p}};\bm{k}_{1},\cdots,\bm{k}_{N_{\mathrm{ex}}}}\equiv(\eta^{+})^{N_{p}}c_{\bm{k}_{1}\downarrow}^{\dagger}\cdots c_{\bm{k}_{N_{\mathrm{ex}}}\downarrow}^{\dagger}\ket{0}$, where $c_{\bm{k}\sigma}=\frac{1}{\sqrt{N_{s}}}\sum_{j}c_{j\sigma}e^{-i\bm{k}\cdot\bm{R}_{j}}$ is the annihilation operator of a single-particle eigenstate with crystal momentum $\bm{k}$, and $N_{s}$ denotes the number of lattice sites. This state also shows the $\eta$-pairing correlations but contains excess spin-down particles that do not form pairs.

Since $L_{j}\ket{\psi}=0$ for $\ket{\psi}=\ket{\psi_{N}}$ and $\ket{\psi}=\ket{\psi_{2N_{p}};\bm{k}_{1},\cdots,\bm{k}_{N_{\mathrm{ex}}}}$, we find that the $\eta$-pairing states $\rho=\ket{\psi}\bra{\psi}$ are steady states of the quantum master equation (2): $\mathcal{L}\ket{\psi}\bra{\psi}=0$. Equivalently, no photon is emitted from the $\eta$-pairing states, implying that they are many-body dark states Diehl et al. (2008); Kraus et al. (2008). The steadiness of the $\eta$-pairing states is physically understood as follows. First, on-site fermion pairs (doublons) are stabilized since a radiative decay of a spin-up particle is Pauli-blocked by a spin-down particle sitting at the same site [see Fig. 1(a)] Sandner et al. (2011). Second, although doublons may have varying center-of-mass momenta and decay via hopping processes, $\eta$ pairs never decay since they constitute eigenstates of the Hubbard Hamiltonian. Physically, this is understood as complete destructive interference between different decay processes of doublons, as can be seen from $-t\sum_{\langle i,j\rangle}\sum_{\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+\mathrm{H.c.})\ket{\psi_{N}}=0$ sup . The above pairing mechanism is clearly distinct from the BCS mechanism which arises from the Fermi-surface instability due to an attractive interaction between fermions.

Symmetries, conserved quantity, and exact Liouvillian eigenmodes.– Since the system has many dark states, a general steady state is a statistical mixture of the $\eta$-pairing states $\ket{\psi_{N}}$ and $\ket{\psi_{2N_{p}};\bm{k}_{1},\cdots,\bm{k}_{N_{\mathrm{ex}}}}$. The steady-state distribution of the $\eta$-pairing states depends on the initial state and the model parameters. To seek an optimal condition, we exploit the fact that the Lindblad operators commute with the generators of the $\eta$ SU(2) symmetry: $[L_{j},\eta^{+}]=[L_{j},\eta^{-}]=0$. This property and the $\eta$ SU(2) symmetry of the Hamiltonian lead to the conservation of

where $\langle\mathcal{O}\rangle=\mathrm{Tr}[\mathcal{O}\rho]$ Buča and Prosen (2012); Albert and Jiang (2014); Tindall et al. (2019). Thus, to achieve large $\eta$-pairing correlations in a steady state, we should prepare an initial state having sufficiently large $C$. This implies that high-energy or high-temperature states rather than low-energy ones of the repulsive Hubbard model are appropriate for an initial state, since the former involve a large number of doublons which give a large $\sum_{j}\langle n_{j\uparrow}n_{j\downarrow}\rangle$ in Eq. (3), whereas the latter do not. The conservation of $C$ is similar to the situation in Ref. Tindall et al. (2019); however, we emphasize that the physical mechanism of the $\eta$ pairing is directly linked to the Fermi statistics and the destructive interference, but not to heating discussed in Ref. Tindall et al. (2019).

The time scale of the formation dynamics of the $\eta$-pairing state is governed by the eigenvalues of the Liouvillian $\mathcal{L}$. Remarkably, some eigenvalues and the corresponding eigenmodes of the Liouvillian can exactly be obtained from the procedure in Refs. Torres (2014); Nakagawa et al. (2021). To see this, we decompose the Liouvillian as $\mathcal{L}=\mathcal{K}+\mathcal{J}$, where $\mathcal{K}\rho\equiv-i(H_{\mathrm{eff}}\rho-\rho H_{\mathrm{eff}}^{\dagger})$ describes an effective non-Hermitian Hamiltonian dynamics with $H_{\mathrm{eff}}\equiv H-(i/2)\sum_{j}L_{j}^{\dagger}L_{j}$, and $\mathcal{J}\rho\equiv\sum_{j}L_{j}\rho L_{j}^{\dagger}$ is the quantum-jump term. The non-Hermitian Hamiltonian $H_{\mathrm{eff}}$ commutes with the total magnetization operator $S_{z}=\sum_{j}(n_{j\uparrow}-n_{j\downarrow})/2$, while the quantum-jump term decreases the magnetization by one. It follows that in the basis that diagonalizes $\mathcal{K}$ the Liouvillian $\mathcal{L}$ is a triangular-matrix form in which the off-diagonal elements are determined from $\mathcal{J}$. Thus, any Liouvillian eigenvalue $\lambda$ is the same as that of $\mathcal{K}$ and can be calculated from a pair of eigenvalues $E_{a}$ and $E_{b}$ of the non-Hermitian Hubbard model $H_{\mathrm{eff}}$ as $\lambda=-i(E_{a}-E_{b}^{*})$ sup . Such a feature does not exist in models with other Lindblad operators used in Refs. Diehl et al. (2008); Kraus et al. (2008); Bernier et al. (2013); Tindall et al. (2019). In particular, for the one-dimensional ($d=1$) case, the present model is exactly solvable by the Bethe ansatz method Nakagawa et al. (2021); Essler et al. (2010).

The Liouvillian excitation spectrum of the $\eta$-pairing steady states involves single-particle and collective excitations as in conventional BCS superfluids. For example, single-particle excitations used to show the metastability of the $\eta$-pairing state Yang (1989) provide exact Liouvillian eigenmodes with degenerate eigenvalues $\lambda=-\gamma$ in any spatial dimensions sup . Such single-particle modes decay on the time scale of $1/\gamma$. The late-stage dynamics near the steady state is governed by collective modes of the $\eta$ SU(2) pseudospins. In fact, the Bethe ansatz solution shows the dispersion relation of the collective modes as $\lambda(P)=2\mathrm{Im}\sqrt{16t^{2}\cos^{2}(P/2)+(U+i\gamma/2)^{2}}-\gamma\simeq 2\mathrm{Im}[J](1+\cos P)$ for $P\simeq\pi$, where $P$ is the momentum of the excitation and $J=4t^{2}/(U+i\gamma/2)$ is an effective exchange coupling between pseudospins, giving a gapless Liouvillian spectrum that implies a power-law tail in the relaxation dynamics in the thermodynamic limit Cai and Barthel (2013); sup . Note that the dispersion relation of the collective modes near $P=\pi$ is quadratic rather than linear in contrast to the Nambu-Goldstone mode of a conventional charge-neutral superfluid. This property highlights the difference between $\eta$-pairing and BCS superfluids; the relevant symmetry of the former is SU(2), while that of the latter is U(1).

Moreover, by the Shiba transformation Shiba (1972) $Vc_{j\uparrow}V^{\dagger}=c_{j\uparrow},Vc_{j\downarrow}V^{\dagger}=e^{i\bm{Q}\cdot\bm{R}_{j}}c_{j\downarrow}^{\dagger}$, the model discussed in this Letter can be mapped to a dissipative Hubbard model subject to two-body particle loss Nakagawa et al. (2020, 2021) described by a Lindblad operator $e^{i\bm{Q}\cdot\bm{R}_{j}}VL_{j}V^{\dagger}=\sqrt{\gamma}c_{j\downarrow}c_{j\uparrow}$. Thus, the Liouvillian spectrum of the model described in Eqs. (1) and (2) is equivalent to that of a dissipative Hubbard model subject to two-body particle loss after changing the sign of the interaction strength. The $\eta$-pairing states of the model correspond to ferromagnetic steady states discussed in Refs. Nakagawa et al. (2020, 2021) through this mapping, while the physical mechanisms of their stabilization are different. This correspondence extends the Shiba symmetry between the repulsive and attractive Hubbard models to the symmetry between spontaneous emission and particle loss, making a nontrivial connection between different dissipative processes of many-body systems.

Dynamics of pair correlations.– We numerically solve the quantum master equation (2) in the one-dimensional case by the quantum-trajectory method Dalibard et al. (1992); Carmichael (1993); Dum et al. (1992); Daley (2014). We take 1000 quantum trajectories, where the time evolution is calculated by exact diagonalization of the non-Hermitian Hamiltonian $H_{\mathrm{eff}}$. We consider an initial state $\rho(0)=\ket{\psi_{0}}\bra{\psi_{0}}$ of an 8-site system, where $\ket{\psi_{0}}=\prod_{j=0,2,4,6}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\ket{0}$ is a density-wave state of doublons. Figure 2 (a) shows the dynamics of the pair correlation function $C(j,0;\tau)\equiv\mathrm{Tr}[c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}c_{0\downarrow}c_{0\uparrow}\rho(\tau)]$. In the steady state, the amplitude of the pair correlation function does not depend on sites and its sign alternates from site to site, indicating the formation of $\eta$ pairs with the center-of-mass momentum $\bm{Q}$. The imaginary part of the pair correlation function is negligible in the steady state, while the real part does not significantly depend on the model parameters sup . Figure 2(b) shows the dynamics of the density of doublons $n_{d}\equiv(1/N_{s})\sum_{j}\langle n_{j\uparrow}n_{j\downarrow}\rangle$. After a long time, the density of doublons reaches a nonzero steady-state value. Since the number of doublons in the initial state is $\mathcal{O}(N_{s})$ and $C\geq 0$, the amplitude of the pair correlations in the steady states is estimated to be $\mathcal{O}(1/N_{s})$, which can be observed as a long-range order in finite-size systems such as trapped atoms Tindall et al. (2021a); Mazurenko et al. (2017). The dynamics in Fig. 2 exhibits a two-step relaxation. Since the decay rate $\gamma$ is sufficiently larger than the hopping amplitude, the pair correlations initially increase over the time scale of hopping which breaks the doublons in the initial state, and the slow relaxation at a late stage near the steady state is governed by the collective modes with an exchange coupling $\mathrm{Im}[J]=-0.16$.

We also study the dynamics from an infinite-temperature initial state $\rho(0)=I/d_{\mathrm{H}}$, where $I$ and $d_{\mathrm{H}}$ are the identity operator and the dimension of the Hilbert space of the 8-particle sector, respectively. To numerically simulate the dynamics with the quantum-trajectory method, we sample a random pure initial state in each quantum trajectory. Figures 2 (c) and (d) show the dynamics of the pair correlations and that of the density of doublons, respectively. Notably, the $\eta$-pairing correlations develop even if the dynamics starts from the infinite-temperature initial state. This behavior is in marked contrast to the heating-induced $\eta$-pairing mechanism in Ref. Tindall et al. (2019), where the $\eta$ paring cannot develop from the infinite-temperature state because of the Hermiticity of Lindblad operators.

Doublon momentum distribution.– The formation of an $\eta$-pairing state can be detected from the momentum distribution of doublons, which can be measured with time-of-flight techniques after doublons are converted into molecules by a Feshbach resonance Regal et al. (2004); Zwierlein et al. (2004); Stöferle et al. (2006); Altman and Vishwanath (2005); Perali et al. (2005); Matyjaśkiewicz et al. (2008). The momentum distribution of doublons is defined as $R_{\bm{k}}\equiv\langle d_{\bm{k}}^{\dagger}d_{\bm{k}}\rangle/\sum_{\bm{k}}\langle d_{\bm{k}}^{\dagger}d_{\bm{k}}\rangle$, where $d_{\bm{k}}\equiv(1/\sqrt{N_{s}})\sum_{j}c_{j\downarrow}c_{j\uparrow}e^{-i\bm{k}\cdot\bm{R}_{j}}$ is the annihilation operator of a doublon with momentum $\bm{k}$. Here, we have normalized the momentum distribution by the total number of doublons such that $0\leq R_{\bm{k}}\leq 1$ is satisfied, since the number of doublons is not conserved during the dynamics. Figure 3 shows the dynamics of the doublon momentum distribution for the two initial states as in Fig. 2. The build-up dynamics of the distribution at $\bm{k}=\pm\bm{Q}$ clearly signals the formation of $\eta$ pairing.

In Fig. 3, one can see that the doublons for $\bm{k}\neq\pm\bm{Q}$ are still populated after a long time. The residual doublon distribution is not due to imperfection of the $\eta$-pairing formation but arises from a fundamental constraint due to the Fermi statistics. In fact, the doublon momentum distribution can be bounded from above as $\langle d_{\bm{k}}^{\dagger}d_{\bm{k}}\rangle\leq\Lambda_{2}/2$ sup , where $\Lambda_{2}$ is the maximum eigenvalue of the two-particle reduced density matrix $(\rho_{2})_{i\sigma_{1}j\sigma_{2};k\sigma_{3}l\sigma_{4}}\equiv\langle c_{j\sigma_{2}}^{\dagger}c_{i\sigma_{1}}^{\dagger}c_{k\sigma_{3}}c_{l\sigma_{4}}\rangle$ Yang (1962). Using the bound of $\Lambda_{2}$ from the Pauli exclusion principle Yang (1962), we obtain

where the bound is saturated by Yang’s $\eta$-pairing state $\ket{\psi_{N}}$ for $\bm{k}=\bm{Q}$ sup . Moreover, provided that the model does not have a dark state other than $\ket{\psi_{N}}$ and $\ket{\psi_{N_{p}};\bm{k}_{1},\cdots,\bm{k}_{N_{\mathrm{ex}}}}$, we can derive a bound on $R_{\bm{Q}}$ in the steady state as sup

The bound (5) is consistent with the numerical results in Fig. 3. If we take the thermodynamic limit $N_{s}\to\infty$ while keeping the density $n\equiv N/(2N_{s})$ constant, the inequality (5) reduces to $R_{\bm{Q}}\leq 1-n$. Thus, the normalized doublon momentum distribution at momentum $\pm\bm{Q}$ cannot reach unity for a system with nonzero density because of the Pauli exclusion principle, indicating a fundamental distinction between doublons and bosons sup .

Experimental implementation.– The validity of the Hubbard model for the description of ultracold fermionic atoms in an optical lattice has been tested by a number of experiments Esslinger (2010). In particular, it has been experimentally verified in Ref. Gall et al. (2020) that a cold-atom quantum simulator of the Hubbard model possesses the Shiba symmetry Shiba (1972); Ho et al. (2009) between the repulsive and attractive Hubbard models, which leads to the $\eta$ SU(2) symmetry when combined with the spin SU(2) symmetry Yang and Zhang (1990). The Lindblad operator $L_{j}$ for spontaneous decay from a spin-up state to a spin-down state can be implemented by weakly coupling the spin-up state to an excited state that can decay to the spin-down state [see Fig. 1(b)] Sandner et al. (2011); sup . With large detuning of the laser that couples the spin-up state and the excited state, the latter can be adiabatically eliminated Pichler et al. (2010); Sarkar et al. (2014); Reiter and Sørensen (2012), and the model (2) is obtained sup . To avoid recoil of atoms that induces transitions to neighboring sites or higher bands, a sufficiently deep optical lattice should be used sup ; Pichler et al. (2010); Sarkar et al. (2014).

In this implementation, a spontaneous decay from the excited state to the spin-up state leads to an additional Lindblad operator $L_{j}^{(d)}=\sqrt{\gamma_{d}}n_{j\uparrow}$, which causes dephasing of atoms. Since the inverse of the dephasing rate $\gamma_{d}$ determines the lifetime of the $\eta$-pairing states as confirmed numerically sup , the condition $\gamma_{d}\ll\gamma,t$ needs to be met. A possible candidate that achieves this condition is alkaline-earth-like fermionic atoms such as ¹⁷¹Yb, ¹⁷³Yb, and ⁸⁷Sr sup . We can implement the Fermi-Hubbard model by trapping spin-polarized atoms in the ground ${}^{1}S_{0}$ state and in the metastable ${}^{3}P_{0}$ state in a magic-wavelength optical lattice Gorshkov et al. (2010). We can also realize an effective spontaneous decay by coupling the ${}^{3}P_{0}$ state to a particular hyperfine state in the ${}^{1}P_{1}$ state that can quickly decay to the ${}^{1}S_{0}$ state Sandner et al. (2011). Since the spontaneous decay from the ${}^{1}P_{1}$ state to the ${}^{3}P_{0}$ state is negligible in the experimental time scale, the condition $\gamma_{d}\ll\gamma,t$ can well be satisfied.

Summary.– In this Letter, we have proposed a nonequilibrium pairing mechanism for fermionic atoms using spontaneous emission of light. The nonequilibrium $\eta$ pairing has been shown to be formed through a non-BCS mechanism, where the Pauli blocking of radiative decay of on-site fermion pairs and the destructive interference between doublon-decay processes play vital roles. We have shown that $\eta$-pairing correlations in the steady state are supported by high density of doublons in the initial state, which indicates that coherent $\eta$ pairs can be created at temperatures comparable to or even higher than the Fermi temperature, as long as the single-band approximation is valid. Since charge-neutral atoms do not have the dynamical instability of $\eta$-pairing states due to coupling with dynamical electromagnetic fields as shown recently Tsuji et al. (2021), the scheme proposed in this Letter for ultracold atoms serves as a promising candidate for observation of the $\eta$-pairing state well above the BCS critical temperature.

Supplemental Material for
“$\eta$ Pairing of Light-Emitting Fermions: Nonequilibrium Pairing Mechanism at High Temperatures”