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

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.

Fermionic superfuidity is one of the most striking quantum many-body phenomena, which has been a subject of intensive investigation in condensed matter physics [38].More recently, ultracold atomic systems have opened a new arena to study fermionic superfluidity [39][40][41], where they are subject to losses due to inelastic collisions.For example, if we consider a superfluid mediated by the orbital Feshbach resonance [42][43][44][45], which controls an interaction between the ground state and an excited state of an atom [46][47][48][49], loss inevitably occurs due to inelastic processes between different orbitals.Such inelastic twobody 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 [19][20][21].
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 no longer 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 [50,51].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 173 Yb, 40 K, and 6 Li [48,49,[52][53][54][55][56].
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 [36] where L i is a Lindblad operator that describes a loss at site i with rate γ, and ρ 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 Such a situation is realized when we consider the dynamics over a sufficiently short time compared with the inverse loss rate 1/γ [31], 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 with a complex-valued interaction U = U 1 + iγ/2, where U 1 , γ > 0. Here, ξ k ≡ k − µ, k is the energy dispersion, µ is the chemical potential, and c kσ and c iσ denote annihilation operators of a spin-σ fermion with momentum k and at site i, respectively.In this Letter, we formulate a mean-field theory from H eff and elucidate how unconventional properties of superfluidity emerge in NH BCS systems.
Formulation of the NH mean-field theory.-Wefirst 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 Here, L E n | and |E n R are left and right eigenstates of H eff with eigenenergy E n , and they satisfy an orthonormal relation L E n |E m R = δ nm [57].We note that, as temperature is not well defined in generic open quantum systems, we only consider the infinite limit of β to elucidate the physics of the ground state and calculate the excitation spectrum.Thus, β 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 ∆, ∆ and then integrate out the fermionic degrees of freedom to obtain Z = D ∆D∆e −S eff (∆, ∆) , where S eff is the effective action given by [58] Here, N denotes the number of lattice sites, and ω n is the Matsubara frequency of fermions.The saddle point condition for the partition function, ∂S eff /∂∆ = ∂S eff /∂ ∆ = 0, yields the NH gap equation if there exists a nontrivial solution other than ∆ = ∆ = 0, where we set β to infinity to obtain an effective ground state.The chemical potential µ 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 [58].
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 A R ≡ n L E n |A|E n R e −βEn /Z.Thus, the order parameters corresponding to the superfluid gap are given by [58] indicating that ∆ = ∆ * since |E n L = |E n 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 8) and neglecting the second-order terms in δ, we obtain the mean-field Hamiltonian which is diagonalized as The quasiparticle operators γkσ , γ kσ and the corresponding energy E k are given by γk↑ = u k c † k↑ − vk c −k↓ , γ k↑ = u k c k↑ − v k c † −k↓ (and similar equations hold for γk↓ , γ k↓ ) and E k = ξ 2 k + ∆∆, respectively, where the coefficients satisfy u 2 k + v k vk = 1 (for their explicit forms, see Supplemental Material [58]).In the Hermitian limit, γkσ and vk respectively reduce to γ † kσ and v * k which describe the Bogoliubov quaiparticles.Here, we note that the mean-field Hamiltonian is non-Hermitian since ∆ * = ∆, and the quasiparticle operators γ kσ and γkσ 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 γ kσ |BCS R = 0 and γ † kσ |BCS L = 0, respectively, where and |0 is the vacuum for fermions.They satisfy L BCS|BCS R = 1 and reproduce the ordinary BCS ground state in the Hermitian limit.We thus obtain which imply that γkσ and γ † kσ 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 β → ∞ limit of the NH gap equation as N/U = k 1/2E k , which is solved self-consistently.We note that the quasiparticle operators satisfy an anticommutation relation {γ kσ , γk σ } = δ kk δ σσ , although these quasiparticles obey neither Fermi nor Bose statistics due to γ † kσ = γkσ , reflecting non-Hermiticity of the mean-field Hamiltonian.
Here, we point out an important relation between the order parameters ∆ and ∆.In the Hermitian case, they are complex conjugate to each other and we can choose a gauge where ∆ is real without loss of generality.This is equivalent to requiring H † MF = H * 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 † = H * 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 where ∆ 0 ∈ C and θ is the U(1) phase.By choosing a special gauge for which Quantum phase transitions of a NH superfluid.-Wesolve the gap equation at β → ∞ numerically (for analytic solutions in the case of constant density of states, see Supplemental Material [58]).Figure 1 shows the superfluid order parameter ∆ 0 .Here, for simplicity, we consider a system with particle-hole symmetry, and set the chemical potential measured from the Fermi energy to zero [58].For small U 1 , Re∆ 0 is suppressed by dissipation γ and then vanishes, indicating a breakdown of superfluidity.Remarkably, as γ increases, the superfluid solution reappears, and the gap size Re∆ 0 is enhanced due to dissipation, eventually exceeding the value in the Hermitian limit.On the other hand, for strong attractive interaction, Re∆ 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 [59][60][61][62] and dissipation.
In the strong-dissipation limit, the behavior of the system is governed by the continuous quantum Zeno effect (QZE) [63][64][65][66][67][68], 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 ∆ 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 [50,51].In fact, when Re∆ 0 = 0, the mean-field Hamiltonian H MF cannot be diagonalized for ξ k = ±Im∆ 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 threedimensional systems, respectively.Such characteristic behavior has its origin in a parity-particle-hole (CP ) symmetry CP H MF (k)(CP ) −1 = −H MF (k) of the meanfield Hamiltonian H MF (k) = k σ z + iIm∆ 0 σ x , where CP = σ x K, σ x,z are the Pauli matrices, and K is complex conjugation [69][70][71][72].In fact, as a function of the momentum, H MF (k) exhibits spontaneous breaking of the CP symmetry at the exceptional points, whose dimensionality is indeed protected by the symmetry constraint [58].Thus, the quantum phase transitions of the NH superfluid cannot be classified into the conventional firstor 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 [58]; 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 ergy 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 Im∆ 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 [58].
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 [73].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.-TheNH quantum phase transitions can be observed by controlling a twobody 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 γ, we may first prepare a metastable superfluid at large U 1 and γ, 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 γ = 0.For a detailed discussion including the relevant timescales, see Supplemental Material [58].
These NH superfluids are expected to be realized with ultracold atoms under inelastic collisions.For example, a superfluid of 173 Yb atoms with an orbital Feshbach resonance [42] offers one such candidate since it is inevitably accompanied by two-body losses as observed experimentally [46][47][48].The effect of non-Hermiticity on the superfluid gap can be observed by spectroscopy with Raman transitions between hyperfine levels or a clock transition [48,49].Furthermore, concerning the control of twobody loss rates, introducing dissipation with photoassociation techniques [67] may also enable the realization of NH superfluids with 40 K and 6 Li [52][53][54][55][56].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.

Supplemental Material for "Theory of Non-Hermitian Fermionic Superfluidity with a Complex-Valued
Interaction" Non-Hermitian mean-field theory of fermionic superfluidity The partition function is written in the path-integral representation as where c and c are Grassmann variables, and τ is the imaginary time.After performing the Hubbard-Stratonovich transformation, the partition function is rewritten as where ∆i (τ ) and ∆ i (τ ) are auxiliary bosonic fields, which are not necessarily complex conjugate to each other in the saddle-point approximation. Substituting where ω n and Ω l are the Matsubara frequencies for fermions and bosons, respectively.In the mean-field theory, we ignore the spatial and temporal fluctuations of ∆ k (Ω l ) and ∆k (Ω l ) and thus set k = 0 and Ω l = 0. Then the action is described as where ∆ = 1 √ βN ∆ k=0 (0) and ∆ = 1 √ βN ∆k=0 (0).Integrating out the fermionic degrees of freedom by using the formula exp − i,j ci A ij c j n i=1 Dc i Dc i = detA, we obtain where is the effective action.By differentiating with respect to ∆ and ∆ and using the saddle point condition ∂S eff /∂∆ = ∂S eff /∂ ∆ = 0, the order parameters are given by Here, we have used the expectation values defined as By using Eqs.(S9) and (S10), we decouple the NH Hamiltonian ( 8) and obtain where the quasiparticle operators are given by and the coefficients are given by Chemical potential in non-Hermitian systems Here we clarify how the chemical potential is determined in non-Hermitian quantum systems.A density matrix ρ of an atomic gas is decomposed into sectors, each of which has a definite particle number M as (S15) In the short-time dynamics of the atomic gas, the time evolution of the density matrix ρ (M ) is described by the NH Hamiltonian H eff restricted to the M -particle Hilbert space [31].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 where N = k,σ c † kσ c kσ is the particle-number operator.In the β → ∞ limit, Eq. ( S16) is calculated as 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 µ = 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): 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 are realized instead of the BCS states [74,75].Here, α k ≡ u k /v k (θ = 0).For these states, the particle number does not depend on the choice of an expectation value: 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 Using the integration contour in Fig. S1, we can calculate the sum over the Matsubara frequency as where f (z) = (e β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 In the β → ∞ limit, the condensation energy is given by 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.
Condensation energy (S25) as a function of γ/t.We here assume the cubic lattice, and the interaction and the chemical potential are set to the same values as in Fig. 1 in the main text, respectively.The dashed lines show the asymptotic behavior in the strong-dissipation limit.The inset in (a) shows an enlarged view of the weak-dissipation regime.
Numerical results in an intermediate regime U1/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 Re∆ 0 in Fig. S3(a).The chemical potential measured from the Fermi energy is set to zero and a cubic lattice is assumed as in Fig. 1.The dashed lines show the asymptotic behavior in the strong-dissipation limit.The inset in (b) shows an enlarged view in the weak-dissipation regime.

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 MF (k) = k σ z + iIm∆ 0 σ x satisfies the following CP symmetry [71] CP H MF (k)(CP ) where the operator CP is given by CP = σ x K.This is similar to parity-time (PT ) symmetry [1,2], which is expressed for a general Hamiltonian H(k) as and they are indeed equivalent to each other as can be seen if we multiply the Hamiltonian by i [15].Consequently, as the PT symmetry dictates that the eigenspectrum of the Hamiltonian appears as real or complex-conjugate pairs [2], 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 where a(k), b z (k), d x (k), d y (k) ∈ R (σ 0 is the 2 × 2 identity matrix).Since its eigenvalues are given by ia(k , the spontaneous CP -symmetry-breaking transition occurs when 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 [69][70][71][72].Because of this symmetry constraint, the exceptional points form a (d − 1)-dimensional surface in d-dimensional systems.In our case, we have a(k) = d y (k) = 0, b z (k) = ξ k , and d x (k) = Im∆ 0 from the mean-field Hamiltonian H MF (k).The quasiparticle energy spectra and exceptional points depicted in Fig. 2 in the main text are consistent with the above general argument.
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 ρ 0 is constant.Under this assumption, the gap equation in the β → ∞ limit reads where ω D = 1/2ρ 0 is the energy cutoff.Here, we set a branch cut of √ z and log z to z ∈ (−∞, 0) and assume that Re∆ 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 Arg (ρ 0 πU 1 ) 2 + (ρ 0 πγ/2 − 1) 2 = 1 as Under the above conditions, the nontrivial solution of the gap equation at β → ∞ is given by 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 where α = ∆0 ω 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 ∆ 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.

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 γ 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 [54].Thus, for a BCS superfluid in which γ, U 1 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 γ/t ∼ 5 using photoassociation techniques [67].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µs [54].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 µs [54].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 [55], Raman spectroscopy, or clock-laser spectroscopy [48].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 [5,6].Since the quantum-gas microscopy for the attractive Hubbard model was already realized [56], this method can also be used for an unambiguous observation of NH superfluidity.

A B C
 Normal  Metastable superfluid  Superfluid FIG.S6.Experimental protocol for observing the reentrant superfluidity.We first prepare on-site pairing superfluidity for large attraction U1 (A) and then introduce large dissipation γ using a sudden application of photoassocitiation (B).The system will reach the reentrant superfluid region by decreasing U1 (C).

FIG. 1 .
FIG. 1. Numerical solution of ∆0 as a function of γ/t obtained from the NH gap equation (5) at β → ∞ and µ = 0 for (a) U1/t = 1.8 and (b) U1/t = 10.We assume a cubic lattice with energy dispersion k = −2t(cos kx +cos ky +cos kz), where t is the hopping amplitude.The dashed lines denote the asymptotic behavior in the strong-dissipation limit.The insets show (a) an enlarged view near the origin (weak dissipation) and (b) that of the real part.

FIG. 2 .
FIG. 2. (a) Real and (b) imaginary parts of the quasiparticle energy spectrum E k = ± 2 k + ∆ 2 0 (orange for positive sign, blue for negative sign) at critical points.(c) Exceptional lines in two dimensions.(d) Exceptional surfaces in three dimensions.The energy dispersion is for a square lattice k = −2(cos kx + cos ky) in (a), (b) and (c), and for a cubic lattice k = −2(cos kx + cos ky + cos kz) in (d).The gap is set to ∆0 = 0.19i for (a), (b) and (c), and ∆0 = 0.4i for (d).
FIG. 3. Phase diagram of the NH BCS model at β → ∞ and µ = 0.The yellow region corresponds to the normal state.The red region shows the metastable superfluid state for which a nontrivial solution of the gap equation gives a local energy minimum.The blue region shows the superfluid state corresponding to a nontrivial solution of the gap equation that gives an effective ground state.A region with small U1 is not shown because of the limitation of numerical calculations.

FIG. S1 .
FIG. S1.Integration contour in the calculation of the condensation energy.

×
FIG. S3.(a) Superfluid gap obtained from the NH gap equation at β → ∞ and (b) the condensation energy (S25) for U1 = 2.5.The chemical potential measured from the Fermi energy is set to zero and a cubic lattice is assumed as in Fig.1.The dashed lines show the asymptotic behavior in the strong-dissipation limit.The inset in (b) shows an enlarged view in the weak-dissipation regime.
FIG. S4.Schematic diagrams for the eigenspectrum of a Hamiltonian that has PT symmetry (left) or CP symmetry (right).The dashed lines indicate the points of spontaneous symmetry breaking.

2 )
FIG. S5.Phase diagram obtained from Eqs. (S33) and (S34).The yellow region corresponds to the normal state where the gap equation has only a trivial solution ∆0 = 0.The red indicates a metastable superfluid state which corresponds to nontrivial solutions of the gap equation but is not an effective ground state.The blue region corresponds to a superfluid phase where a nontrivial solution of the gap equation gives an effective ground state.