Theory of Non-Hermitian Fermionic Superfluidity with a Complex-Valued Interaction

TL;DR

This paper develops a mean-field theory for non-Hermitian fermionic superfluids with complex interactions, revealing unique phase transitions and the effects of dissipation on superfluid properties.

Contribution

It introduces a novel NH gap equation and uncovers unconventional phase transitions driven by non-Hermitian effects in fermionic superfluids.

Findings

Superfluidity exhibits reentrant behavior with increasing dissipation.

Superfluid gap is enhanced by dissipation in the strong interaction regime.

Unique phase transitions are linked to exceptional points in the quasiparticle spectrum.

Abstract

Motivated by recent experimental advances in ultracold atoms, we analyze a non-Hermitian (NH) BCS Hamiltonian with a complex-valued interaction arising from inelastic scattering between fermions. We develop a mean-field theory to obtain a NH gap equation for order parameters, which are different from the standard BCS ones due to the inequivalence of left and right eigenstates in the NH physics. We find unconventional phase transitions unique to NH systems: superfluidity shows reentrant behavior with increasing dissipation, as a consequence of non-diagonalizable exceptional points, lines, and surfaces in the quasiparticle Hamiltonian for weak attractive interactions. For strong attractive interactions, the superfluid gap never collapses but is enhanced by dissipation due to an interplay between the BCS-BEC crossover and the quantum Zeno effect. Our results lay the groundwork for studies of fermionic superfluidity subject to inelastic collisions.

Figures

Equations

$$\dot{\rho }$$

$$=-i(H_{\mathrm{eff}}\rho -\rho H_{\mathrm{eff}}^{†})+\gamma \sum _i{L}_i\rho L_i^{†},$$

$$H_{\mathrm{eff}}=\sum _{\mathbf{k}\sigma }{\xi }_{\mathbf{k}}{c}_{\mathbf{k}\sigma }^{†}{c}_{\mathbf{k}\sigma }-U\sum _i{c}_{i\uparrow }^{†}{c}_{i\downarrow }^{†}{c}_{i\downarrow }{c}_{i\uparrow },$$

$$Z=\sum _n{e}^{-\beta E_n}=\sum _n{}_L⟨E_n|e^{-\beta H_{\mathrm{eff}}}|E_n{⟩}_R.$$

$$S_{\mathrm{eff}}(\bar{\Delta },\Delta )=-\sum _{{\omega }_n,\mathbf{k}}\log ({\omega }_n^2+{\xi }_{\mathbf{k}}^2+\bar{\Delta }\Delta )+\frac{\beta N}{U}\bar{\Delta }\Delta .$$

$$\frac{N}{U}=\sum _{\mathbf{k}}\frac{1}{2\sqrt{{\xi }_{\mathbf{k}}^2+\bar{\Delta }\Delta }}\tanh \frac{\beta \sqrt{{\xi }_{\mathbf{k}}^2+\bar{\Delta }\Delta }}{2},$$

$$\Delta =-\frac{U}{N}\sum _{\mathbf{k}}{}_L⟨c_{-\mathbf{k}\downarrow }{c}_{\mathbf{k}\uparrow }{⟩}_R,$$

$$\bar{\Delta }=-\frac{U}{N}\sum _{\mathbf{k}}{}_L⟨c_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}{⟩}_R,$$

$$H_{\mathrm{eff}}=\sum _{\mathbf{k}\sigma }{\xi }_{\mathbf{k}}{c}_{\mathbf{k}\sigma }^{†}{c}_{\mathbf{k}\sigma }-\frac{U}{N}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }}{c}_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}{c}_{-{\mathbf{k}}^{\prime }\downarrow }{c}_{{\mathbf{k}}^{\prime }\uparrow }.$$

$$H_{\mathrm{MF}}=\sum _{\mathbf{k}}\left(\begin{array}{cc}{c}_{\mathbf{k}\uparrow }^{†}& c_{-\mathbf{k}\downarrow }\end{array}\right)\left(\begin{array}{cc}{\xi }_{\mathbf{k}}& \Delta \\ \bar{\Delta }& -{\xi }_{\mathbf{k}}\end{array}\right)\left(\begin{array}{c}{c}_{\mathbf{k}\uparrow }\\ c_{-\mathbf{k}\downarrow }^{†}\end{array}\right),$$

$$|\mathrm{BCS}{⟩}_R=\prod _{\mathbf{k}}\left(u_{\mathbf{k}}+v_{\mathbf{k}}{c}_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}\right)|0⟩,$$

$$|\mathrm{BCS}{⟩}_L=\prod _{\mathbf{k}}\left(u_{\mathbf{k}}^{\ast }+{\bar{v}}_{\mathbf{k}}^{\ast }{c}_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}\right)|0⟩,$$

$$\Delta (\theta )={\Delta }_0{e}^{i\theta },$$

$$\bar{\Delta }(\theta )={\Delta }_0{e}^{-i\theta },$$

$$Z=\int \mathcal{D}\bar{c}\mathcal{D}c\exp \{-S(\bar{c},c)\},$$

$$S(\bar{c},c)={\int }_0^{\beta }d\tau \left(\sum _{\mathbf{k}\sigma }{\bar{c}}_{\mathbf{k}\sigma }(\tau )({\partial }_{\tau }+{\xi }_{\mathbf{k}})c_{\mathbf{k}\sigma }(\tau )-U\sum _i{\bar{c}}_{i\uparrow }(\tau ){\bar{c}}_{i\downarrow }(\tau )c_{i\downarrow }(\tau )c_{i\uparrow }(\tau )\right),$$

$$Z=\int \mathcal{D}\bar{c}\mathcal{D}c\mathcal{D}\bar{\Delta }\mathcal{D}\Delta e^{-S(\bar{\Delta },\Delta ,\bar{c},c)},$$

$$S(\bar{\Delta },\Delta ,\bar{c},c)={\int }_0^{\beta }d\tau \left[\sum _{\mathbf{k}\sigma }{\bar{c}}_{\mathbf{k}\sigma }(\tau )({\partial }_{\tau }+{\xi }_{\mathbf{k}})c_{\mathbf{k}\sigma }(\tau )+\sum _i\left(\frac{{\bar{\Delta }}_i(\tau ){\Delta }_i(\tau )}{U}+{\bar{\Delta }}_i(\tau )c_{i\downarrow }(\tau )c_{i\uparrow }(\tau )+{\Delta }_i(\tau ){\bar{c}}_{i\uparrow }(\tau ){\bar{c}}_{i\downarrow }(\tau )\right)\right],$$

$$S(\bar{\Delta },\Delta ,\bar{c},c)=$$

$$+\sum _{\mathbf{k},{\mathbf{k}}^{\prime \prime },{\omega }_n,{\Omega }_l}\frac{1}{\sqrt{\beta N}}{\bar{\Delta }}_{{\mathbf{k}}^{\prime \prime }}({\Omega }_l)c_{{\mathbf{k}}^{\prime \prime }-\mathbf{k}\downarrow }({\Omega }_l-{\omega }_n)c_{\mathbf{k}\uparrow }({\omega }_n)$$

$$+\sum _{\mathbf{k},{\mathbf{k}}^{\prime \prime },{\omega }_n,{\Omega }_l}\frac{1}{\sqrt{\beta N}}{\Delta }_{{\mathbf{k}}^{\prime \prime }}({\Omega }_l){\bar{c}}_{\mathbf{k}\uparrow }({\omega }_n){\bar{c}}_{{\mathbf{k}}^{\prime \prime }-\mathbf{k}\downarrow }({\Omega }_l-{\omega }_n),$$

$$S(\bar{c},c,\bar{\Delta },\Delta )=\frac{\beta N}{U}\bar{\Delta }\Delta +\sum _{\mathbf{k},{\omega }_n}\left(\begin{array}{cc}{\bar{c}}_{\mathbf{k}\uparrow }({\omega }_n)& c_{-\mathbf{k}\downarrow }(-{\omega }_n)\end{array}\right)\left(\begin{array}{cc}-i{\omega }_n+{\xi }_{\mathbf{k}}& \Delta \\ \bar{\Delta }& -i{\omega }_n-{\xi }_{\mathbf{k}}\end{array}\right)\left(\begin{array}{c}{c}_{\mathbf{k}\uparrow }({\omega }_n)\\ {\bar{c}}_{-\mathbf{k}\downarrow }(-{\omega }_n)\end{array}\right),$$

$$Z=\int \mathcal{D}\bar{\Delta }\mathcal{D}\Delta e^{-S_{\mathrm{eff}}(\Delta ,\bar{\Delta })},$$

$$S_{\mathrm{eff}}(\bar{\Delta },\Delta )$$

$$=-\sum _{{\omega }_n,\mathbf{k}}\ln ({\omega }_n^2+{\xi }_{\mathbf{k}}^2+\bar{\Delta }\Delta )+\frac{\beta N}{U}\bar{\Delta }\Delta ,$$

$$e^{-S_{\mathrm{eff}}(\Delta ,\bar{\Delta })}=\int \mathcal{D}\bar{c}\mathcal{D}c{e}^{-S(\bar{\Delta },\Delta ,\bar{c},c)}$$

$$\Delta =-\frac{U}{N}\sum _{\mathbf{k}}{}_L⟨c_{-\mathbf{k}\downarrow }{c}_{\mathbf{k}\uparrow }{⟩}_R,$$

$$\bar{\Delta }=-\frac{U}{N}\sum _{\mathbf{k}}{}_L⟨c_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}{⟩}_R.$$

$$\sum _{\mathbf{k}}{}_L⟨c_{-\mathbf{k}\downarrow }{c}_{\mathbf{k}\uparrow }{⟩}_R=\frac{\int \mathcal{D}\bar{c}\mathcal{D}c{\sum }_{\mathbf{k},{\omega }_n}{c}_{-\mathbf{k}\downarrow }(-{\omega }_n)c_{\mathbf{k}\uparrow }({\omega }_n)e^{-S}}{\int \mathcal{D}\bar{c}\mathcal{D}c{e}^{-S}},$$

$$\sum _{\mathbf{k}}{}_L⟨c_{\mathbf{k}\uparrow }^{†}{c}_{-\mathbf{k}\downarrow }^{†}{⟩}_R=\frac{\int \mathcal{D}\bar{c}\mathcal{D}c{\sum }_{\mathbf{k},{\omega }_n}{\bar{c}}_{\mathbf{k}\uparrow }({\omega }_n){\bar{c}}_{-\mathbf{k}\downarrow }(-{\omega }_n)e^{-S}}{\int \mathcal{D}\bar{c}\mathcal{D}c{e}^{-S}}.$$

$$H_{\mathrm{MF}}=\sum _{\mathbf{k}}\left(\begin{array}{cc}{c}_{\mathbf{k}\uparrow }^{†}& c_{-\mathbf{k}\downarrow }\end{array}\right)\left(\begin{array}{cc}{\xi }_{\mathbf{k}}& \Delta \\ \bar{\Delta }& -{\xi }_{\mathbf{k}}\end{array}\right)\left(\begin{array}{c}{c}_{\mathbf{k}\uparrow }\\ c_{-\mathbf{k}\downarrow }^{†}\end{array}\right)=\sum _{\mathbf{k}}\left(\begin{array}{cc}{\bar{\gamma }}_{\mathbf{k}\uparrow }& {\gamma }_{-\mathbf{k}\downarrow }\end{array}\right)\left(\begin{array}{cc}{E}_{\mathbf{k}}& 0\\ 0& -E_{\mathbf{k}}\end{array}\right)\left(\begin{array}{c}{\gamma }_{\mathbf{k}\uparrow }\\ {\bar{\gamma }}_{-\mathbf{k}\downarrow }\end{array}\right),$$

$${\bar{\gamma }}_{\mathbf{k}\uparrow }=u_{\mathbf{k}}{c}_{\mathbf{k}\uparrow }^{†}-{\bar{v}}_{\mathbf{k}}{c}_{-\mathbf{k}\downarrow },$$

$${\bar{\gamma }}_{-\mathbf{k}\downarrow }={\bar{v}}_{\mathbf{k}}{c}_{\mathbf{k}\uparrow }+u_{\mathbf{k}}{c}_{-\mathbf{k}\downarrow }^{†},$$

$${\gamma }_{\mathbf{k}\uparrow }=u_{\mathbf{k}}{c}_{\mathbf{k}\uparrow }-v_{\mathbf{k}}{c}_{-\mathbf{k}\downarrow }^{†},$$

$${\gamma }_{-\mathbf{k}\downarrow }=v_{\mathbf{k}}{c}_{\mathbf{k}\uparrow }^{†}+u_{\mathbf{k}}{c}_{-\mathbf{k}\downarrow },$$

$$u_{\mathbf{k}}=\sqrt{\frac{E_{\mathbf{k}}+{\xi }_{\mathbf{k}}}{2E_{\mathbf{k}}}},v_{\mathbf{k}}=-\sqrt{\frac{(E_{\mathbf{k}}-{\xi }_{\mathbf{k}})}{2E_{\mathbf{k}}}}\frac{\sqrt{\Delta }}{\sqrt{\bar{\Delta }}},{\bar{v}}_{\mathbf{k}}=-\sqrt{\frac{(E_{\mathbf{k}}-{\xi }_{\mathbf{k}})}{2E_{\mathbf{k}}}}\frac{\sqrt{\bar{\Delta }}}{\sqrt{\Delta }}.$$

$$\rho =\sum _M{\rho }^{(M)}.$$

Full text

Introduction.– In recent years, non-Hermitian (NH) quantum systems have been actively studied both experimentally and theoretically Bender and Boettcher (1998); Bender (2007); Nakagawa et al. (2018); Lourenço et al. (2018); Ashida et al. (2016, 2017); Kawabata et al. (2017); Rudner and Levitov (2009); Lee (2016); Shen et al. (2018); Kunst et al. (2018); Yao et al. (2018); Kawabata et al. (2018a); Gong et al. (2018); Kawabata et al. (2019, 2018b); Zhou and Lee (2019); Lee and Chan (2014); Ghatak and Das (2018); Zhou and Cui (2018); Chtchelkatchev et al. (2012); Guo et al. (2009); Rüter et al. (2010); Feng et al. (2011); Regensburger et al. (2012); Zeuner et al. (2015); Li et al. (2019); Xiao et al. (2017); Zhou et al. (2018); Xiao et al. (2018); Dürr et al. (2009); García-Ripoll et al. (2009); Daley et al. (2009); Kantian et al. (2009); Yoshida et al. (2018); Yao and Wang (2018). NH quantum systems arise when the system undergoes dissipation to an environment Daley (2014); Dalibard et al. (1992). It has been revealed that non-Hermiticity drastically alters the properties of a number of quantum phenomena that have been established in the Hermitian physics, such as quantum phase transitions Bender and Boettcher (1998); Bender (2007); Nakagawa et al. (2018); Lourenço et al. (2018), quantum critical behavior Ashida et al. (2016, 2017); Kawabata et al. (2017), topological phases Rudner and Levitov (2009); Lee (2016); Shen et al. (2018); Kunst et al. (2018); Yao et al. (2018); Kawabata et al. (2018a); Gong et al. (2018); Kawabata et al. (2019, 2018b); Zhou and Lee (2019), and magnetism Lee and Chan (2014). Such theoretical predictions have been confirmed experimentally by using optical systems and ultracold atoms Guo et al. (2009); Rüter et al. (2010); Feng et al. (2011); Regensburger et al. (2012); Zeuner et al. (2015); Li et al. (2019); Xiao et al. (2017); Zhou et al. (2018); Xiao et al. (2018). However, since most of the previous studies dealt with single-particle physics, exploration of many-body physics in NH systems is still in its infancy Nakagawa et al. (2018); Lourenço et al. (2018); Ashida et al. (2016, 2017); Dürr et al. (2009); García-Ripoll et al. (2009); Daley et al. (2009); Kantian et al. (2009).

Fermionic superfuidity is one of the most striking quantum many-body phenomena, which has been a subject of intensive investigation in condensed matter physics Bardeen et al. (1957). More recently, ultracold atomic systems have opened a new arena to study fermionic superfluidity Regal et al. (2004); Zwierlein et al. (2004); Bartenstein et al. (2004), where they are subject to losses due to inelastic collisions. For example, if we consider a superfluid mediated by the orbital Feshbach resonance Zhang et al. (2015); Iskin (2016); Xu et al. (2016); He et al. (2016), which controls an interaction between the ground state and an excited state of an atom Höfer et al. (2015); Pagano et al. (2015); Cappellini et al. (2019); Darkwah Oppong et al. (2019), loss inevitably occurs due to inelastic processes between different orbitals. Such inelastic two-body losses cause the decay of eigenstates of the Hamiltonian and may be described by complex-valued interactions, thus providing an ideal platform to study NH fermionic superfluids. Despite its growing importance, however, theory for NH fermionic superfluidity has not been established yet Ghatak and Das (2018); Zhou and Cui (2018); Chtchelkatchev et al. (2012).

In this Letter, we demonstrate how fermionic superfluidity in ultracold atoms is modified under inelastic collisions, by generalizing the standard BCS theory to a situation in which fermions interact with each other via a complex-valued attraction. We elucidate that the non-Hermiticity alters several fundamental properties of superfluidity; for example, the order parameters of particles and holes are not necessarily complex conjugate to each other in the NH physics, and the Bogoliubov quasiparticles obey neither Fermi nor Bose statistics since eigenstates are, in general, not orthogonal to each other.

Furthermore, we find that the non-Hermiticity leads to unique quantum phase transitions in superfluids. For a weak interaction, the real part of the superfluid gap is first suppressed and then quenched with increasing dissipation. Remarkably, superfluidity is restored beyond a certain strength of dissipation and the superfluid gap is even enhanced afterwards with increasing dissipation. We show that these phase transitions emerge from exceptional points, lines and surfaces that are unique to the NH physics, where the Hamiltonian cannot be diagonalized Heiss (2012); Berry (2004). For a strong interaction, superfluidity is not suppressed and never breaks down because fermions are paired to form molecules on each site, thereby avoiding intersite decoherence. Our finding can experimentally be tested in various ultracold atomic species under inelastic collisions such as 173Yb, 40K, and 6Li Cappellini et al. (2019); Darkwah Oppong et al. (2019); Stewart et al. (2008); Chin et al. (2004, 2006); Stöferle et al. (2006); Mitra et al. (2018).

Model.– We consider ultracold fermionic atoms with an attractive interaction in a three-dimensional optical lattice. When atoms undergo inelastic collisions, the scattered atoms are lost from the system since a large internal energy is converted to the kinetic energy. An atomic gas undergoing two-body losses due to inelastic collisions is described by a quantum master equation Daley (2014)

[TABLE]

where $L_{i}$ is a Lindblad operator that describes a loss at site $i$ with rate $\gamma$, and $\rho$ is the density matrix of the atomic gas. When the quantum-jump term, which is the last term in Eq. (1), is negligible, the system is described by an effective NH Hamiltonian $H_{\mathrm{eff}}=H-\frac{i}{2}\gamma\sum_{i}L_{i}^{\dagger}{L}_{i}$. Such a situation is realized when we consider the dynamics over a sufficiently short time compared with the inverse loss rate $1/\gamma$ Dürr et al. (2009), which characterizes the timescale where the effect of quantum jumps becomes significant. In this case, the lowest real part of the eigenspectrum gives the effective ground state, and the imaginary part of energy corresponds to a decay rate of each eigenstate. The two-body loss is described by $L_{i}=c_{i\downarrow}c_{i\uparrow}$, giving a NH BCS Hamiltonian

[TABLE]

with a complex-valued interaction $U=U_{1}+i\gamma/2$, where $U_{1},\gamma>0$. Here, $\xi_{\bm{k}}\equiv\epsilon_{\bm{k}}-\mu$, $\epsilon_{\bm{k}}$ is the energy dispersion, $\mu$ is the chemical potential, and $c_{\bm{k}\sigma}$ and $c_{i\sigma}$ denote annihilation operators of a spin-$\sigma$ fermion with momentum $\bm{k}$ and at site $i$, respectively. In this Letter, we formulate a mean-field theory from $H_{\mathrm{eff}}$ and elucidate how unconventional properties of superfluidity emerge in NH BCS systems.

Formulation of the NH mean-field theory.– We first clarify how the standard BCS mean-field theory is changed due to non-Hermiticity by formulating it with a path-integral approach. We start with a partition function defined as

[TABLE]

Here, ${}_{L}\langle E_{n}|$ and $|E_{n}\rangle_{R}$ are left and right eigenstates of $H_{\mathrm{eff}}$ with eigenenergy $E_{n}$, and they satisfy an orthonormal relation ${}_{L}\langle E_{n}|E_{m}\rangle_{R}=\delta_{nm}$ Brody (2013). We note that, as temperature is not well defined in generic open quantum systems, we only consider the infinite limit of $\beta$ to elucidate the physics of the ground state and calculate the excitation spectrum. Thus, $\beta$ is a parameter used to formulate a path integral and should not be regarded as the temperature of the system. We use a path-integral representation of the partition function (3) to perform the Hubbard-Stratonovich transformation with auxiliary fields $\Delta$, $\bar{\Delta}$ and then integrate out the fermionic degrees of freedom to obtain $Z=\int\mathcal{D}\bar{\Delta}\mathcal{D}{\Delta}e^{-S_{\mathrm{eff}}(\Delta,\bar{\Delta})}$, where $S_{\mathrm{eff}}$ is the effective action given by Sup

[TABLE]

Here, $N$ denotes the number of lattice sites, and $\omega_{n}$ is the Matsubara frequency of fermions. The saddle point condition for the partition function, $\partial{S}_{\mathrm{eff}}/\partial\Delta=\partial{S}_{\mathrm{eff}}/\partial\bar{\Delta}=0$, yields the NH gap equation

[TABLE]

if there exists a nontrivial solution other than $\Delta=\bar{\Delta}=0$, where we set $\beta$ to infinity to obtain an effective ground state. The chemical potential $\mu$ is determined so that the mean particle number in the non-Hermitian ensemble (3) is equal to the particle number of the density-matrix sector of interest Sup .

In NH physics, four distinct types of order parameters can be defined according to whether left and right eigenstates are assigned to the bra or ket vectors in the expectation value. Importantly, the expectation value of an operator $A$ calculated from the non-Hermitian ensemble (3) should correspond to $~{}_{L}\langle{A}\rangle_{R}\equiv\sum_{n}~{}_{L}\langle E_{n}|A|E_{n}\rangle_{R}e^{-\beta{E}_{n}}/Z$. Thus, the order parameters corresponding to the superfluid gap are given by Sup

[TABLE]

indicating that $\bar{\Delta}\neq\Delta^{*}$ since $|E_{n}\rangle_{L}\neq|E_{n}\rangle_{R}$. As discussed below, this leads to various intriguing consequences on the properties of a NH superfluid.

To elucidate the effect of non-Hermiticity, let us apply the mean-field decoupling to the NH BCS Hamiltonian with the simplest $s$-wave pairing interaction

[TABLE]

Substituting $c_{\bm{k}\uparrow}^{\dagger}{c}_{-\bm{k}\downarrow}^{\dagger}={}_{L}\langle c_{\bm{k}\uparrow}^{\dagger}{c}_{-\bm{k}\downarrow}^{\dagger}\rangle_{R}+\delta({c}_{\bm{k}\uparrow}^{\dagger}{c}_{-\bm{k}\downarrow}^{\dagger})$ and $c_{-\bm{k}\downarrow}c_{\bm{k}\uparrow}={}_{L}\langle c_{-\bm{k}\downarrow}c_{\bm{k}\uparrow}\rangle_{R}+\delta({c}_{-\bm{k}\downarrow}c_{\bm{k}\uparrow})$ into Eq. (8) and neglecting the second-order terms in $\delta$, we obtain the mean-field Hamiltonian

[TABLE]

which is diagonalized as $H_{\mathrm{MF}}=\sum_{\bm{k}}E_{\bm{k}}(\bar{\gamma}_{\bm{k}\uparrow}\gamma_{\bm{k}\uparrow}+\bar{\gamma}_{-\bm{k}\downarrow}\gamma_{-\bm{k}\downarrow})-\sum_{\bm{k}}E_{\bm{k}}$. The quasiparticle operators $\bar{\gamma}_{\bm{k}\sigma}$, $\gamma_{\bm{k}\sigma}$ and the corresponding energy $E_{\bm{k}}$ are given by $\bar{\gamma}_{\bm{k}\uparrow}=u_{\bm{k}}c_{\bm{k}\uparrow}^{\dagger}-\bar{v}_{\bm{k}}{c}_{-\bm{k}\downarrow}$, $\gamma_{\bm{k}\uparrow}=u_{\bm{k}}c_{\bm{k}\uparrow}-v_{\bm{k}}c_{-\bm{k}\downarrow}^{\dagger}$ (and similar equations hold for $\bar{\gamma}_{\bm{k}\downarrow}$, $\gamma_{\bm{k}\downarrow}$) and $E_{\bm{k}}=\sqrt{\xi_{\bm{k}}^{2}+\bar{\Delta}\Delta}$, respectively, where the coefficients satisfy $u_{\bm{k}}^{2}+v_{\bm{k}}\bar{v}_{\bm{k}}=1$ (for their explicit forms, see Supplemental Material Sup ). In the Hermitian limit, $\bar{\gamma}_{\bm{k}\sigma}$ and $\bar{v}_{\bm{k}}$ respectively reduce to $\gamma^{\dagger}_{\bm{k}\sigma}$ and $v^{*}_{\bm{k}}$ which describe the Bogoliubov quaiparticles. Here, we note that the mean-field Hamiltonian is non-Hermitian since $\Delta^{*}\neq\bar{\Delta}$, and the quasiparticle operators $\gamma_{\bm{k}\sigma}$ and $\bar{\gamma}_{\bm{k}\sigma}$ are not Hermitian conjugate to each other. Therefore, the Hamiltonian cannot be diagonalized via a unitary transformation.

As a result, the right and left ground states are defined by $\gamma_{\bm{k}\sigma}|\mathrm{BCS}\rangle_{R}=0$ and $\bar{\gamma}_{\bm{k}\sigma}^{\dagger}|\mathrm{BCS}\rangle_{L}=0$, respectively, where

[TABLE]

and $\ket{0}$ is the vacuum for fermions. They satisfy ${}_{L}\langle\mathrm{BCS}|\mathrm{BCS}\rangle_{R}=1$ and reproduce the ordinary BCS ground state in the Hermitian limit. We thus obtain $H\bar{\gamma}_{\bm{k}\sigma}|\mathrm{BCS}\rangle_{R}=E_{\bm{k}}\bar{\gamma}_{\bm{k}\sigma}\mathrm{BCS}\rangle_{R}$ and $H^{\dagger}\gamma_{\bm{k}\sigma}^{\dagger}|\mathrm{BCS}\rangle_{L}=E_{\bm{k}}^{*}\gamma_{\bm{k}\sigma}^{\dagger}|\mathrm{BCS}\rangle_{L}$, which imply that $\bar{\gamma}_{\bm{k}\sigma}$ and $\gamma_{\bm{k}\sigma}^{\dagger}$ create the right and left eigenstates, respectively, when acted on the ground state. Here, we have shifted the ground state energy to zero. Using Eqs. (6), (7), (10), and (11), we obtain the $\beta\to\infty$ limit of the NH gap equation as $N/U=\sum_{\bm{k}}1/2E_{\bm{k}}$, which is solved self-consistently. We note that the quasiparticle operators satisfy an anticommutation relation $\{\gamma_{\bm{k}\sigma},\bar{\gamma}_{\bm{k}^{\prime}\sigma^{\prime}}\}=\delta_{\bm{k}\bm{k}^{\prime}}\delta_{\sigma\sigma^{\prime}}$, although these quasiparticles obey neither Fermi nor Bose statistics due to $\gamma_{\bm{k}\sigma}^{\dagger}\neq\bar{\gamma}_{\bm{k}\sigma}$, reflecting non-Hermiticity of the mean-field Hamiltonian.

Here, we point out an important relation between the order parameters $\Delta$ and $\bar{\Delta}$. In the Hermitian case, they are complex conjugate to each other and we can choose a gauge where $\Delta$ is real without loss of generality. This is equivalent to requiring $H_{\mathrm{MF}}^{\dagger}=H_{\mathrm{MF}}^{*}$ in the matrix representation in the Fock-state basis in the Hilbert space. Now we consider the NH case. As in the Hermitian case, the NH BCS Hamiltonian (8) satisfies a symmetry relation $H^{\dagger}=H^{\ast}$ under the matrix representation in terms of Fock states, indicating that the left eigenstates are obtained through complex conjugation of the right ones. The NH BCS Hamiltonian has the U(1) symmetry as in the Hermitian case and this is not affected by the complex nature of the interaction. Then, when a superfluid is formed, its ground states become degenerate due to spontaneous U(1) symmetry breaking. The BCS ground states (10) and (11) are consistent with these properties if

[TABLE]

where $\Delta_{0}\in\mathbb{C}$ and $\theta$ is the U(1) phase. By choosing a special gauge for which $H_{\mathrm{MF}}^{\dagger}=H_{\mathrm{MF}}^{*}$ is satisfied, we have $\Delta=\bar{\Delta}$. Here, we note that the relation (12) is specific to the NH BCS Hamiltonian (8) and may be changed depending on symmetry of a NH Hamiltonian.

Quantum phase transitions of a NH superfluid.– We solve the gap equation at $\beta\to\infty$ numerically (for analytic solutions in the case of constant density of states, see Supplemental Material Sup ). Figure 1 shows the superfluid order parameter $\Delta_{0}$. Here, for simplicity, we consider a system with particle-hole symmetry, and set the chemical potential measured from the Fermi energy to zero Sup . For small $U_{1}$, $\mathrm{Re}\Delta_{0}$ is suppressed by dissipation $\gamma$ and then vanishes, indicating a breakdown of superfluidity. Remarkably, as $\gamma$ increases, the superfluid solution reappears, and the gap size $\mathrm{Re}\Delta_{0}$ is enhanced due to dissipation, eventually exceeding the value in the Hermitian limit. On the other hand, for strong attractive interaction, $\mathrm{Re}\Delta_{0}$ is not suppressed, but rather enhanced due to dissipation.

The qualitative difference between the cases of weak and strong attractions is explained by an interplay between the BCS-BEC crossover Leggett (1980a, b); Nozières and Schmitt-Rink (1985); Diener et al. (2008) and dissipation. In the strong-dissipation limit, the behavior of the system is governed by the continuous quantum Zeno effect (QZE) Syassen et al. (2008); Mark et al. (2012); Barontini et al. (2013); Zhu et al. (2014); Tomita et al. (2017); Yan et al. (2013), which suppresses tunneling to neighboring sites, leading to localization of particles. In the case of weak attraction, the inter-site coherence of Cooper pairs is suppressed by dissipation and the superfluidity is destroyed. However, localization due to the QZE facilitates formation of on-site molecules of fermions for strong dissipation and consequently superfluidity reappears. In fact, the solution of the gap equation approaches $\Delta_{0}=U/2$ in the strong-dissipation limit, as shown by the dashed lines in Fig. 1. This is consistent with the fact that the physics is dominated by the on-site interaction under the QZE, supporting our mean-field analysis. On the other hand, under strong attraction, fermions form bosonic molecules almost at single sites and thus the molecules can survive under dissipation. In this case, the effect of dissipation is to give rise to an effective one-body loss of molecules and the remaining molecules can undergo Bose-Einstein condensation.

We here point out that the breakdown and restoration of superfluidity present clear signatures of the emergence of exceptional points, where the Hamiltonian cannot be diagonalized Heiss (2012); Berry (2004). In fact, when $\mathrm{Re}\Delta_{0}=0$, the mean-field Hamiltonian $H_{\mathrm{MF}}$ cannot be diagonalized for $\xi_{\bm{k}}=\pm\mathrm{Im}\Delta_{0}$. Figure 2(a) shows the real part of the energy spectrum of quasiparticles in two dimensions. The regions where the orange and blue surfaces merge form exceptional points, lines (Fig. 2(c)), and surfaces (Fig. 2(d)) in one-, two-, and three-dimensional systems, respectively. Such characteristic behavior has its origin in a parity-particle-hole ($CP$) symmetry $CPH_{\mathrm{MF}}(\bm{k})(CP)^{-1}=-H_{\mathrm{MF}}(\bm{k})$ of the mean-field Hamiltonian $H_{\mathrm{MF}}(\bm{k})=\epsilon_{\bm{k}}\sigma_{z}+i\mathrm{Im}\Delta_{0}\sigma_{x}$, where $CP=\sigma_{x}K$, $\sigma_{x,z}$ are the Pauli matrices, and $K$ is complex conjugation Yoshida et al. (2019); Budich et al. (2019); Okugawa and Yokoyama (2019); Zhou et al. (2019). In fact, as a function of the momentum, $H_{\mathrm{MF}}(\bm{k})$ exhibits spontaneous breaking of the $CP$ symmetry at the exceptional points, whose dimensionality is indeed protected by the symmetry constraint Sup . Thus, the quantum phase transitions of the NH superfluid cannot be classified into the conventional first- or second-order phase transitions in Hermitian systems, but are attributed to the emergence of exceptional points unique to non-Hermiticity. We note here that the two exceptional points corresponding to the breakdown and restoration of superfluidity merge and disappear as the strength of attraction increases. Intriguingly, as a remnant of the merged exceptional points, the real part of the superfluid gap for intermediate strengths of $U_{1}$ shows a characteristic minimum at a certain strength of dissipation Sup ; the gap is first suppressed by dissipation, but enhanced again by the QZE as the dissipation is further increased.

Furthermore, the emergence of exceptional manifolds leads to some intriguing dynamics in the NH superfluid. In Fig. 2(b), the imaginary part of the quasiparticle energy takes a positive finite value only in between the exceptional lines or surfaces, amplifying quasiparticle distribution in the particular region of the Brillouin zone through the time evolution. The characteristic structure in Fig. 2(b), which can be used as a smoking gun of the non-Hermiticity, can be observed as long as $\mathrm{Im}\Delta_{0}>0$ even in a region away from the breakdown and restoration points.

We note that the nontrivial solution of the NH gap equation may give a metastable superfluid, which corresponds to a local minimum of the real part of the energy. Whether the superfluid is metastable or not is decided from comparison of ground-state energies between the superfluid state and the normal state, as detailed in Supplemental Material Sup .

From these results, we obtain a phase diagram of the NH BCS model as shown in Fig. 3. In the blue region, the superfluid state is an effective ground state of the NH BCS Hamiltonian. When the dissipation is increased, the superfluid state remains stable if the attraction is sufficiently strong. When the system enters the red region, the superfluid state becomes metastable with respect to the real part of the energy. The metastable superfluid undergoes an unconventional quantum phase transition due to exceptional points in the case of weak attractions, leading to the disappearance of the superfluid state in the yellow region. Here we remark that a similar phase diagram is obtained from an exact solution of a one-dimensional NH Hubbard model Nakagawa et al. . Although a similar formulation to obtain the gap equation can be made for a continuum system, the reentrant superfluidity is unique to the lattice system since the localization due to the QZE cannot occur without a lattice. On the other hand, the breakdown of the superfluidity and the related phase transitions can occur in the continuum system.

Towards experimental realization.– The NH quantum phase transitions can be observed by controlling a two-body loss rate. Since the superfluid is metastable in the red region, it can be realized by slowly increasing dissipation from the blue region. For weak attraction $U_{1}$, the superfluid undergoes an unconventional phase transition to the normal state due to the exceptional points. To observe the reappearance of the superfluidity under large $\gamma$, we may first prepare a metastable superfluid at large $U_{1}$ and $\gamma$, and then decrease $U_{1}$. Finally, the metastablity of the superfluid can be confirmed through comparison of the results between fast and slow increases of the dissipation from $\gamma=0$. For a detailed discussion including the relevant timescales, see Supplemental Material Sup .

These NH superfluids are expected to be realized with ultracold atoms under inelastic collisions. For example, a superfluid of 173Yb atoms with an orbital Feshbach resonance Zhang et al. (2015) offers one such candidate since it is inevitably accompanied by two-body losses as observed experimentally Höfer et al. (2015); Pagano et al. (2015); Cappellini et al. (2019). The effect of non-Hermiticity on the superfluid gap can be observed by spectroscopy with Raman transitions between hyperfine levels or a clock transition Cappellini et al. (2019); Darkwah Oppong et al. (2019). Furthermore, concerning the control of two-body loss rates, introducing dissipation with photoassociation techniques Tomita et al. (2017) may also enable the realization of NH superfluids with 40K and 6Li Stewart et al. (2008); Chin et al. (2004, 2006); Stöferle et al. (2006); Mitra et al. (2018). Although the strength of dissipation is usually fixed by scattering properties of atoms, dissipation engineering using photoassociation techniques will be a feasible method for realizing the NH fermionic superfluidity.

Conclusions.– We have investigated how the BCS superfluidity is extended to NH quantum systems under inelastic interactions. We have elucidated some remarkable features unique to the NH fermionic superfluidity, such as exotic Bogoliubov quasiparticles which belong to neither fermions nor bosons and found unconventional quantum phase transitions unique to non-Hermiticity. In particular, for weak attraction, it has been revealed that the superfluidity breaks down with increasing dissipation but shows reentrant behavior as dissipation is further increased. Remarkably, these phase transitions are accompanied by distinctive features of the non-Hermiticity, i.e. the emergence of exceptional points, lines and surfaces in the quasiparticle Hamiltonian for one-, two- and three-dimensions. These characteristic features will play a decisive role in detecting NH phase transitions in experiments. On the other hand, for strong attraction, the superfluid state is not suppressed but enhanced due to the confinement of molecules to single sites via the QZE. While we have focused on a conventional $s$-wave superfluid, $p$-wave, $d$-wave and other exotic superfluids in NH systems will also be relevant for experiments, and merit future investigation.

Appendix A Non-Hermitian mean-field theory of fermionic superfluidity

The partition function is written in the path-integral representation as

[TABLE]

where $c$ and $\bar{c}$ are Grassmann variables, and $\tau$ is the imaginary time. After performing the Hubbard-Stratonovich transformation, the partition function is rewritten as

[TABLE]

[TABLE]

where $\bar{\Delta}_{i}(\tau)$ and $\Delta_{i}(\tau)$ are auxiliary bosonic fields, which are not necessarily complex conjugate to each other in the saddle-point approximation. Substituting $c_{i\sigma}(\tau)=\frac{1}{\sqrt{\beta{N}}}\sum_{\bm{k},\omega_{n}}e^{i(\bm{r}_{i}\cdot\bm{k}-\omega_{n}\tau)}c_{\bm{k}\sigma}(\omega_{n})$, $\bar{c}_{i\sigma}(\tau)=\frac{1}{\sqrt{\beta{N}}}\sum_{\bm{k},\omega_{n}}e^{-i(\bm{r}_{i}\cdot\bm{k}-\omega_{n}\tau)}\bar{c}_{\bm{k}\sigma}(\omega_{n})$, $\Delta_{i}(\tau)=\frac{1}{\sqrt{\beta{N}}}\sum_{\bm{k},\Omega_{l}}e^{i(\bm{r}_{i}\cdot\bm{k}-\Omega_{l}\tau)}\Delta_{\bm{k}}(\Omega_{l})$ and $\bar{\Delta}_{i}(\tau)=\frac{1}{\sqrt{\beta{N}}}\sum_{\bm{k},\Omega_{l}}e^{-i(\bm{r}_{i}\cdot\bm{k}-\Omega_{l}\tau)}\bar{\Delta}_{\bm{k}}(\Omega_{l})$ into Eq. (S3), we have

[TABLE]

where $\omega_{n}$ and $\Omega_{l}$ are the Matsubara frequencies for fermions and bosons, respectively. In the mean-field theory, we ignore the spatial and temporal fluctuations of $\Delta_{\bm{k}^{\prime\prime}}(\Omega_{l})$ and $\bar{\Delta}_{\bm{k}^{\prime\prime}}(\Omega_{l})$ and thus set $\bm{k}^{\prime\prime}=0$ and $\Omega_{l}=0$. Then the action is described as

[TABLE]

where $\Delta=\frac{1}{\sqrt{\beta{N}}}\Delta_{\bm{k}=0}(0)$ and $\bar{\Delta}=\frac{1}{\sqrt{\beta{N}}}\bar{\Delta}_{\bm{k}=0}(0)$. Integrating out the fermionic degrees of freedom by using the formula $\int\exp\left[-\sum_{i,j}\bar{c}_{i}A_{ij}c_{j}\right]\prod_{i=1}^{n}\mathcal{D}\bar{c}_{i}\mathcal{D}c_{i}=\mathrm{det}A$, we obtain

[TABLE]

where

[TABLE]

is the effective action. By differentiating

[TABLE]

with respect to $\Delta$ and $\bar{\Delta}$ and using the saddle point condition $\partial{S}_{\mathrm{eff}}/\partial\Delta=\partial{S}_{\mathrm{eff}}/\partial\bar{\Delta}=0$, the order parameters are given by

[TABLE]

Here, we have used the expectation values defined as

[TABLE]

By using Eqs. (S9) and (S10), we decouple the NH Hamiltonian (8) and obtain

[TABLE]

where the quasiparticle operators are given by

[TABLE]

and the coefficients are given by

[TABLE]

Appendix B Chemical potential in non-Hermitian systems

Here we clarify how the chemical potential is determined in non-Hermitian quantum systems. A density matrix $\rho$ of an atomic gas is decomposed into sectors, each of which has a definite particle number $M$ as

[TABLE]

In the short-time dynamics of the atomic gas, the time evolution of the density matrix $\rho^{(M)}$ is described by the NH Hamiltonian $H_{\mathrm{eff}}$ restricted to the $M$-particle Hilbert space Dürr et al. (2009). An effective energy spectrum in the $M$-particle Hilbert space can be extracted by tuning the chemical potential so that the mean particle number in the non-Hermitian ensemble (3) in the main text is

[TABLE]

where $\hat{N}=\sum_{\bm{k},\sigma}c_{\bm{k}\sigma}^{\dagger}c_{\bm{k}\sigma}$ is the particle-number operator. In the $\beta\to\infty$ limit, Eq. (S16) is calculated as

[TABLE]

where the BCS ground states (10) and (11) in the main text are used. For a cubic or square lattice considered in the main text, a half-filling condition $M=N$ is achieved by setting $\mu=0$ in Eq. (S17).

We note that the expectation value of particle numbers taken by the right BCS eigenstates does not coincide with Eq. (S17):

[TABLE]

since the BCS states do not conserve the particle number. However, this is not a contradiction since in the $M$-particle Hilbert space, the number-conserving states

[TABLE]

are realized instead of the BCS states Leggett and Sols (1991); Leggett (2006). Here, $\alpha_{\bm{k}}\equiv u_{\bm{k}}/v_{\bm{k}}(\theta=0)$. For these states, the particle number does not depend on the choice of an expectation value:

[TABLE]

Appendix C Calculation of the condensation energy of the superfluid state

We here discuss how to obtain the condensation energy of the superfluid state for consideration of its stability. Using Eq. (S7), the condensation energy is given by the difference in energy between the superfluid and normal states as

[TABLE]

Using the integration contour in Fig. S1, we can calculate the sum over the Matsubara frequency as

[TABLE]

where $f(z)=(e^{\beta{z}}+1)^{-1}$ is the Fermi distribution function. We here note that the integrand in Eq. (S23) has a branch cut on the lines connecting the branch points of Eq. (S22) as shown in Fig. S1. We thus obtain

[TABLE]

In the $\beta\to\infty$ limit, the condensation energy is given by

[TABLE]

Figure S2 shows the condensation energy for the parameters corresponding to those in Fig. 1 in the main text. For $U_{1}/t=1.8$, we find that the superfluid state is energetically stable only in the weak-dissipation region. For a region where the real part of the condensation energy is positive, the nontrivial solution of the NH gap equation gives a local minimum of the real part of energy, leading to a metastable superfluid solution. We also find that the nontrivial solution of the gap equation in the strong-dissipation regime gives a metastable state because the energy is higher than that of the trivial solution. On the other hand, for $U_{1}/t=10$, the superfluid is always energetically stable, and we can observe the steady enhancement of the superfluid gap.

Appendix D Numerical results in an intermediate regime $U_{1}/t=2.5$

We have conducted numerical calculations for an intermediate strength of the interaction $U_{1}/t=2.5$ as shown in Fig. S3. We see from Fig. S3(a) that the real part of the superfluid gap is suppressed by dissipation but does not vanish. Moreover, it is even enhanced and eventually exceeds the value of the Hermitian limit as the dissipation increases, which is attributed to the confinement of Cooper pairs to individual sites due to the QZE. On the other hand, from Fig. S3(b), the superfluid state under sufficiently strong dissipation becomes metastable due to the positive condensation energy. While the superfluid does not exhibit an unconventional phase transition with exceptional points in this case, a remnant of the phase transition in a weak attraction can be seen as the minimum of $\mathrm{Re}\Delta_{0}$ in Fig. S3(a).

Appendix E Spontaneous $CP$ symmetry breaking as an origin of exceptional points

At the breakdown and restoration points of the superfluid gap, the gap is pure imaginary and the mean-field Hamiltonian $H_{\mathrm{MF}}(\bm{k})=\epsilon_{\bm{k}}\sigma_{z}+i\mathrm{Im}\Delta_{0}\sigma_{x}$ satisfies the following $CP$ symmetry Okugawa and Yokoyama (2019)

[TABLE]

where the operator $CP$ is given by $CP=\sigma_{x}K$. This is similar to parity-time ($\mathcal{PT}$) symmetry Bender and Boettcher (1998); Bender (2007), which is expressed for a general Hamiltonian $H(\bm{k})$ as

[TABLE]

and they are indeed equivalent to each other as can be seen if we multiply the Hamiltonian by $i$ Kawabata et al. (2019). Consequently, as the $\mathcal{PT}$ symmetry dictates that the eigenspectrum of the Hamiltonian appears as real or complex-conjugate pairs Bender (2007), the $CP$ symmetry dictates that the eigenspectrum appears as pure imaginary or anti-complex-conjugate pairs (Fig. S4). One can easily confirm that the $CP$ symmetry is unbroken if the eigenspectrum is pure imaginary; otherwise it is spontaneously broken.

A generic two-band Hamiltonian with the $CP$ symmetry can be written as

[TABLE]

where $a(\bm{k}),b_{z}(\bm{k}),d_{x}(\bm{k}),d_{y}(\bm{k})\in\mathbb{R}$ ($\sigma_{0}$ is the $2\times 2$ identity matrix). Since its eigenvalues are given by $ia(\bm{k})\pm\sqrt{b_{z}(\bm{k})^{2}-d_{x}(\bm{k})^{2}-d_{y}(\bm{k})^{2}}$, the spontaneous $CP$-symmetry-breaking transition occurs when

[TABLE]

which is accompanied by the emergence of exceptional points. We note that only a single condition (S29) is needed for the exceptional point, while two conditions are needed for an exceptional point in a system without any symmetry Yoshida et al. (2019); Budich et al. (2019); Okugawa and Yokoyama (2019); Zhou et al. (2019). Because of this symmetry constraint, the exceptional points form a $(d-1)$-dimensional surface in $d$-dimensional systems. In our case, we have $a(\bm{k})=d_{y}(\bm{k})=0,b_{z}(\bm{k})=\xi_{\bm{k}}$, and $d_{x}(\bm{k})=\mathrm{Im}\Delta_{0}$ from the mean-field Hamiltonian $H_{\mathrm{MF}}(\bm{k})$. The quasiparticle energy spectra and exceptional points depicted in Fig. 2 in the main text are consistent with the above general argument.

Appendix F Analytical calculation of the gap equation and the condensation energy at a constant density of states

The gap equation (Eq. (5) in the main text) can be solved analytically if the density of states $\rho_{0}$ is constant. Under this assumption, the gap equation in the $\beta\to\infty$ limit reads

[TABLE]

where $\omega_{D}=1/2\rho_{0}$ is the energy cutoff. Here, we set a branch cut of $\sqrt{z}$ and $\log{z}$ to $z\in(-\infty,0)$ and assume that $\mathrm{Re}\Delta_{0}$ is positive without loss of generality. Results obtained below do not depend on the choice of the branch cut. Then, the argument is restricted to $\mathrm{Arg}\frac{\sqrt{\omega_{D}^{2}+\Delta_{0}^{2}}+\omega_{D}}{\Delta_{0}}\in(-\frac{\pi}{2},\frac{\pi}{2})$ and the gap equation has nontrivial solutions if and only if $U_{1}$ satisfies $\rho_{0}U_{1}<\frac{1}{\pi}$ and $\gamma$ satisfies $0\leq\gamma\leq\gamma^{c_{1}}$ or $\gamma\geq\gamma^{c_{2}}$, where $\gamma^{c_{1}}$ and $\gamma^{c_{2}}$ are determined from $(\rho_{0}\pi{U}_{1})^{2}+(\rho_{0}\pi\gamma/2-1)^{2}=1$ as

[TABLE]

Under the above conditions, the nontrivial solution of the gap equation at $\beta\to\infty$ is given by

[TABLE]

This gives the behavior consistent with Fig. 1 in the main text.

Next, we calculate the condensation energy (S25) for a constant density of states. In this case, Eq. (S25) can be rewritten as

[TABLE]

where $\alpha=\frac{\Delta_{0}}{\omega_{D}}$. This gives the behavior consistent with Fig. S2.

Finally, we show a phase diagram in Fig. S5 obtained from the analytical solutions of the gap equation with a constant density of states. In Fig. S5, the yellow region corresponds to the normal state, where the gap equation has only a trivial solution $\Delta_{0}=0$. This phase is surrounded by the red region, where the gap equation has the nontrivial solution (S33) which corresponds to the metastable superfluid state, since the nontrivial solution gives a local minimum of the real part of energy due to a positive condensation energy. When the attractive interaction $U_{1}$ is sufficiently strong, the system is in the superfluid phase (blue region in Fig. S5), where the nontrivial solution of the gap equation gives an effective ground state. The phase diagram is qualitatively consistent with Fig. 3 in the main text.

Appendix G Details of experimental setups for probing the non-Hermitian superfluid

Here we discuss an experimental setup for probing the NH superfluidity. As mentioned in the main text, the NH dynamics is realized when we neglect the quantum-jump term in the master equation (1). To fulfill this condition, the loss rate $\gamma$ should be much smaller than the energy scale that governs the timescale for relaxation towards a quasi-equilibrium state. In superfluid states, the thermalization proceeds in a timescale of hopping of atoms to neighboring sites Chin et al. (2006). Thus, for a BCS superfluid in which $\gamma,U_{1}\ll t$ is satisfied, the NH breakdown of superfluid may be observed. We note that, as inferred from Figs. 3 and S5, the breakdown of superfluid due to exceptional points can be induced by small dissipation in the BCS regime, which justifies the NH dynamics. On the other hand, for the reentrant superfluidity, which is another unique phenomenon in the NH superfluid, dissipation larger than the hopping is required (see Figs. 3 and S5). This seems to invalidate the assumption for the ignorance of the quantum-jump term. However, as mentioned in the main text, the physics behind the reentrant superfluidity is the QZE, which suppresses the hopping and facilitates the formation of on-site molecules. The manifestation of the QZE can be observed by using the following protocol. Hence, a large dissipation is not incompatible with the ignorance of the quantum-jump term.

The experimental protocol is illustrated in Fig. S6. We first prepare a superfluid state for large attraction $U_{1}$, where atom pairs are confined to each lattice site. Then, we introduce large dissipation $\gamma/t\sim 5$ using photoassociation techniques Tomita et al. (2017). Finally, we decrease the attraction $U_{1}$ by tuning the magnetic field for a Feshbach resonance, which can be operated fast on the order of few $10\mu{s}$ Chin et al. (2006). If we do not have dissipation, atoms in the on-site molecular pairs tunnel to neighboring sites at a hopping rate $t$, which gives a timescale of the order of 100 $\mu$s Chin et al. (2006). However, under a large dissipation in the quantum Zeno regime, the tunneling of atoms is suppressed by the QZE and thus the dissociation of molecules after ramping down the attraction $U_{1}$ is delayed. Consequently, even after a timescale of $1/t$, atoms surviving in the system will still form on-site molecules. Such QZE-assisted molecules can be regarded as a signature of the reentrant superfluidity. In ultracold atoms, molecules (double occupancy) can be detected by measuring binding energies, for example, using radio-frequency spectroscopy Stöferle et al. (2006), Raman spectroscopy, or clock-laser spectroscopy Cappellini et al. (2019). Ramping down the lattice depth will give a shorter timescale of hopping, and for the region of the reentrant superfluidity, the above methods can be used to determine the phase boundaries due to exceptional points.

Finally, we note that the NH dynamics is faithfully realized when we perform a quantum-gas-microscopy measurement of an atom number and then postselect the measurement outcome which does not contain any atom loss Ashida et al. (2016, 2017). Since the quantum-gas microscopy for the attractive Hubbard model was already realized Mitra et al. (2018), this method can also be used for an unambiguous observation of NH superfluidity.