Non-Hermitian BCS-BEC crossover of Dirac fermions

We investigate chiral symmetry breaking in a model of Dirac fermions with a complexified coupling constant whose imaginary part represents dissipation. We introduce a chiral chemical potential and observe that for real coupling a relativistic BCS-BEC crossover is realized. We solve the model in the mean-field approximation and construct the phase diagram as a function of the complex coupling. It is found that the dynamical mass increases under dissipation, although the chiral symmetry gets restored if dissipation exceeds a threshold.


Introduction
One of the fundamental tenets of modern quantum mechanics and quantum field theory is that Hamiltonians should be Hermitian. It ensures that energy levels are real and the time evolution is unitary. However, it has been perceived in recent years that there are situations in which physics can be effectively described in terms of non-Hermitian Hamiltonians [1,2]. Such a description proves to be useful for open quantum systems that interact with environments [3][4][5]. In PT -symmetric quantum mechanics, the existence of PT symmetry ensures real energy spectra even when the Hamiltonian is not Hermitian [6][7][8]. This is not just a mathematical possibility but can be realized experimentally [9]. It is not an overstatement to say that the physics of non-Hermitian Hamiltonians is far richer than conventional Hermitian ones; in fact, symmetry classification [10,11] reveals that there are 38 distinct symmetry classes in non-Hermitian systems, while there are only 10 in Hermitian systems. There are a plethora of exotic phenomena that are unique to non-Hermitian systems, such as the non-Hermitian skin effect [12][13][14][15][16] and the non-Hermitian localization transitions [17].
Superfluidity and superconductivity are amongst the most salient and fascinating phenomena in quantum many-body physics [62,63]. Non-Hermitian fermionic superfluidity has been explored in [64][65][66][67][68]. 1 In [67] the authors pointed out that ultracold atomic gases with two-body losses due to inelastic collisions can be naturally described with a complexvalued interaction. They solved the gap equation for fermions on a square/cubic lattice and mapped out the phase diagram in the mean-field approximation. In [68] this analysis was extended to fermions in a continuum model, where the phase diagram across the entire range from weak to strong coupling was obtained. This is a non-Hermitian analog of the well-established BCS-BEC crossover of fermions with s-wave interactions [62,[80][81][82][83][84][85]]when the s-wave scattering length is varied, the system evolves continuously from a weakly interacting BCS regime of loose Cooper pairs to the BEC regime of tightly bound molecules.
The primary goal of the present paper is to generalize the analysis of [67,68] to Dirac fermions and investigate a non-Hermitian relativistic BCS-BEC crossover. 2 What is the motivation of this study? First, Dirac fermions can be realized with ultracold atoms loaded on an optical lattice [86][87][88] and hence the experimental protocol proposed in [67] to produce a complex coupling can in principle be applied to this case as well. Second, a relativistic BCS-BEC crossover is believed to take place in QCD and QCD-like theories at finite density [89][90][91][92][93][94][95][96][97][98] (see [99] for a review). While the interaction strength and the density can be separately varied in nonrelativistic systems, this is not the case in QCD -actually we face a density-induced crossover: at low density, quark matter is strongly coupled and tightly bound diquarks condense, whereas at high density the coupling is weak due to asymptotic freedom and a BCS-type description is justified [100]. The study of such a crossover (including the possibility of phase transitions at intermediate densities) is potentially relevant to compact star phenomenology and heavy-ion collision experiments.
In this paper we investigate χSB in the NJL model with a complex coupling. Generally, χSB is considered to be a strong coupling phenomenon and a weakly coupled BCS picture does not apply. However, at finite chiral chemical potential [21,[101][102][103][104][105][106][107][108], χSB occurs at an arbitrarily weak coupling due to the fact that a chiral chemical potential induces a nonzero density of states for fermions at low energy and serves as a catalyst of χSB [107,108]. By tuning both the coupling strength and the chiral chemical potential, we probe the entire range of the BCS-BEC crossover for χSB and construct a complete phase diagram as a function of the complex four-fermion coupling. A novel mechanism for emergence and disappearance of complex saddles of the action is also illustrated.
2 Although [49] solved a zero-dimensional model of Dirac fermions with a complex four-fermion coupling, its higher-dimensional analog has not been thoroughly studied yet.
a chirally imbalanced matter in gauge theories is intrinsically unstable due to the axial anomaly [109][110][111], our model has no coupling to gauge fields and there is no instability due to anomalies. This paper is organized as follows. In section 2 the model is defined and the thermodynamic potential is derived. In section 3 the phase diagram for real coupling is presented. The detrimental effect of the baryon chemical potential on χSB is illustrated. In section 4 we turn on an imaginary part of the coupling, solve the gap equation numerically, map out the phase diagram, and determine the boundary between the normal phase, a metastable χSB phase, and a stable χSB phase. The results are then compared with those for nonrelativistic fermions [67,68]. We conclude in section 5.

The NJL model
We consider a model with the partition function where the Euclidean Lagrangian is given by µ is the quark chemical potential and µ 5 is the chiral (or axial) chemical potential. We set the current mass to zero. Eq. (2.1) is the same model as in [107,108] except that here we have nonzero µ.
It is straightforward to introduce N f flavors of quarks and let N f → ∞ with G ∝ 1/N f , which allows us to rigorously justify a saddle point analysis, but we shall stick to N f = 1 for simplicity of exposition.
With the Hubbard-Stratonovich transformation, quarks can be readily integrated out and yields The bosonic fields are related to fermionic observables as σ = G ψψ and π = G ψiγ 5 ψ .
In the mean-field approximation, we have where a momentum cutoff Λ was introduced to remove UV divergences. Let us define dimensionless variables which leads to a dimensionless action S ≡ S/Λ 4 given by In the following, we will assume 0 < µ 5 < 1 so that the Fermi surface stays inside the domain of integration. In generic open quantum systems, temperature is not well defined, and we treat t as a formal parameter used to define the path integral for the partition function, as in [67]. We will set t to zero in the ensuing analysis.

Phase diagram for real coupling
Let us begin with a discussion for real coupling g > 0. In the zero-temperature limit t → +0, (2.8) reduces to Numerical minimization of S allows us to determine the dynamical mass M as a function of g and µ 5 . Our result for µ = 0 is presented in Figure 1. As one can see from the left plot, while there is a nonzero critical coupling ≈ 20 at µ 5 = 0, it goes away for µ 5 > 0: χSB occurs for any nonzero coupling. This is the catalysis effect emphasized in e.g., [107,108]. The dynamical mass increases monotonically with µ 5 . What is the distinction between the BCS regime and the BEC regime? In nonrelativistic systems the unitarity limit marks a crisp boundary, but what about relativistic systems? The locus of a "transition" between the two regimes is not uniquely defined. As suggested in [99], one popular criterion is to see the dispersion relation E(p) of quasiparticles. If E(p) takes a minimum at |p| ≈ p F , it is in the BCS regime, and if E(p) is a monotonically increasing function of |p|, it is in the BEC regime. This works pretty well in dense QCD where E(p) = ( p 2 + M 2 − µ) 2 + ∆ 2 (with ∆ the superconducting gap) experiences such a transition when the dynamical mass M is equal to µ. In the current model, however, the dispersion (2.6) takes a minimum at |p| = µ 5 regardless of the dynamical mass. Yet another criterion is to compare the interparticle distance and the size of Cooper pairs. If the wave function of Cooper pairs extends beyond the average interparticle distance, it is in the BCS regime, otherwise in the BEC regime. In [113,114] the size of Cooper pairs in relativistic color superconductors was computed as a function of the quark density and such a BCS-BEC-type crossover was indicated. Unfortunately, a similar analysis is difficult for our model because the interaction is pointlike and the gap has no momentum dependence. As a rule of thumb, let us take the inverse of M as the size of a Cooper pair, and take the inverse of µ 5 as the average inter-quark distance. Then the region with 1/M > 1/ µ 5 (1/M < 1/ µ 5 ) corresponds to the BCS (BEC) regime, respectively. According to this crude estimate we labeled each regime in the right plot of Figure 1.
Next we proceed to the analysis for nonzero µ. The phase diagrams are shown in Figure 2. The main observation here is that for any µ 5 , there is a critical µ beyond which the chiral symmetry is restored. This is because µ induces a mismatch of Fermi surfaces and disrupts Cooper pairing. The critical value of µ is known as the Chandrasekhar-Clogston limit [115,116]. Analogous situations arise in both condensed matter [117,118] and QCD [119][120][121][122][123][124]. When the ordinary isotropic Cooper pairing is hampered, a nonstandard pairing that breaks translation symmetry is likely to set in, though it is beyond the scope of this paper.

Phase diagram for complex coupling
Finally we complexify the coupling constant. Throughout this section we set µ = 0 to simplify the ensuing numerical analysis. Eq. (2.8) reduces to where g ∈ C and M ∈ C. When M is purely imaginary, (x ± µ 5 ) 2 + M 2 may be exactly on the branch cut of the complex square root, which makes the zero-temperature limit t → +0 ill-defined. When M is not purely imaginary, we have for t → +0 This integral can be performed with the formula where tanh −1 (z) = 1 2 log 1+z 1−z is the inverse function of tanh(z). This way we obtain Since S is a function of M 2 , the trivial vacuum M = 0 is always a solution to ∂S/∂M = 0. We are interested in the gap equation for M = 0, which reads (4.6) We have varied g on the complex plane and numerically searched for a solution to (4.6) for each g. It turned out that there was no solution for Re g < 0, indicating that chiral symmetry is unbroken for Re g < 0. This is natural because Re g < 0 is a repulsive interaction. Furthermore, the phase structures for Im g > 0 and Im g < 0 are symmetric about the real axis of g. Hence we will assume Re g > 0 and Im g > 0 in the following. By monitoring the magnitude of the gradient ∂S/∂M we found an interesting mechanism that changes the number of saddle points of S. In Figure 3 Figure 3(b), a new saddle point is suddenly born out of the imaginary axis of M . So there are now two saddle points. When the imaginary part of g is further increased, as shown in Figure 3(c), the old saddle is absorbed into the imaginary axis and we are left with a single saddle. In this fashion the number of saddles (i.e., the solutions to (4.6)) can jump abruptly.
When there are multiple saddles, the dominant one is definitely the one that has the lowest value of Re S. Following [67,68], we define three phases as below.
• Normal phase: (4.6) has no solution, i.e., M = 0 is the only saddle of S. • Metastable χSB phase: there are solutions to (4.6), but their Re S are higher than that for M = 0. • Stable χSB phase: there are solutions to (4.6) whose Re S are lower than that for M = 0.
In Figure 4 (left) we display the phase diagram for µ 5 = 0.5 on the complex g plane. In the vicinity of the real axis we have a stable χSB phase. As Im g increases, we are driven   into a metastable χSB phase via a quantum phase transition. At small Re g the metastable saddle goes away and the chiral symmetry is completely restored. Interestingly enough, the normal phase at small Re g sharply ends at Im g 22.5; this means that a very strong dissipation can trigger χSB (albeit a metastable one). The global phase structure we found here is quite similar to the one for nonrelativistic fermions on a lattice [67], though we note that the phase diagram at small Re g was not presented in [67] because of the limitation of numerical calculations.
To gain more insights, in Figure 4 (right) we plot the magnitude of the gap |M |. When there are multiple saddles, we took the one that has the lowest Re S. Notice that nothing dramatic happens at the boundary between the metastable χSB phase and the stable χSB phase. It is worth noting that the gap magnitude tends to be enhanced by dissipation. As a crude guide we drew the boundary |M | ∼ µ 5 between the BCS regime and the BEC regime. We emphasize that this is a rule of thumb and a more rigorous characterization of the crossover region in non-Hermitian superfluids is left as an open problem. We point out that for Im g 20 the gap reaches the UV cutoff scale (|M | ∼ 1); this implies that all the calculations for Im g 20 are extremely sensitive to the regularization scheme used, and hence one has to be careful about physical interpretations. Figure 5 shows the real and imaginary parts of M . The left panel shows that Re M grows monotonically with Re g and is largely independent of Im g. Note that Re M approaches zero along the boundary with the normal phase. This means that, when one moves out of the normal phase, a nontrivial solution to the gap equation emerges out of the imaginary axis of M . The right panel shows that Im M grows monotonically with Im g.
It is instructive to compare our findings with preceding works.
• Ref. [67] found an enhancement of superfluidity by dissipation, and attributed it to the continuous quantum Zeno effect which suppresses tunneling and reinforces onsite molecule formation. This effect is unique to a lattice system. Nevertheless we observed a similar enhancement in a continuum model; we speculate that the presence of a hard momentum cutoff plays a role similar to that of a lattice. • Ref. [68] solved the BCS-BEC crossover for a complex scattering length and found that the superfluid phase is stable even in the limit of strong dissipation, which is surprising. The difference of [68] from [67] and this paper may be due to the fact that [68] considered a complex-valued chemical potential. 3 • Ref. [66] reports that a non-Hermitian perturbation (an imaginary magnetic field) enhances superfluidity. • Ref. [34] solved the NJL model supplemented with a non-Hermitian bilinear term that preserves chiral symmetry. It was found that, as the non-Hermitian coupling grows, the dynamical mass first rises and then drops to zero.
Finally we turn to the quasiparticle excitation spectra E(p) = (|p|/Λ − µ 5 ) 2 + M 2 which is generally complex for complex g. The real and (scaled) imaginary parts of E(p) are displayed in Figure 6. At weak coupling, the imaginary part of E(p) is sharply concentrated around the Fermi level (left panel). By contrast, at stronger coupling, the imaginary part has a much broader support (right panel). These are results for weak dissipation (Im g = 0.5). If dissipation is stronger, a more peculiar thing can occur: at the boundary between the normal phase and the metastable χSB phase, M is purely imaginary (cf. Figure 5), hence E(p) is purely imaginary for µ 5 − |M | ≤ |p|/Λ ≤ µ 5 + |M | and is real otherwise. The bounds |p|/Λ = µ 5 ± |M | are an example of the so-called exceptional points [127] where the dimensionality of the eigenspace of the Hamiltonian decreases.

Conclusions and outlook
In summary, we have investigated χSB of Dirac fermions in four dimensions at finite chiral chemical potential for a complex-valued four-fermion coupling. We varied the coupling over the entire range from the weakly bound BCS regime to the tightly bound BEC regime, and numerically constructed the phase diagram. Our primary result is that the imaginary coupling tends to enhance χSB up to a certain threshold; when the imaginary coupling exceeds the threshold, the chiral symmetry is restored, although a nontrivial solution to the gap equation continues to exist. We illustrated how complex saddles of the action come into existence and go away through the imaginary axis of the complex gap plane along which the action has a branch cut. We also worked out the complex energy spectra of quasiparticles. Our results can in principle be tested in experiments using ultracold atomic gases and other materials that host Dirac fermions [128][129][130]. 4 However, it may be difficult to draw implications for quark matter in compact stars from this work, because (i) observable signals from compact stars are extremely scarce, and (ii) the pointlike four-fermion interaction is a very crude approximation to the non-Abelian gauge interaction between quarks.
There are miscellaneous future directions in which this work can be extended. A partial list is given below.
• To incorporate the competition between the chiral condensate and the diquark condensate (i.e., a competition between a Dirac mass and a Majorana mass) along the lines of [131,132] • To investigate collective fluctuations around the saddle points and test the validity of the mean-field approximation for non-Hermitian χSB • To analytically prove numerical findings in this work, e.g., the existence of a nontrivial solution to the gap equation for large Im g • To examine dissipative χSB under an external magnetic field that catalyzes χSB [133,134] • To study U(1) A vortices • To study the effect of nonzero µ on non-Hermitian χSB • To analyze the structure of Cooper pairs as a function of the complex coupling, and provide a more precise description of the non-Hermitian BCS-BEC crossover • To use the renormalization group to improve the mean-field analysis • To clarify how to apply the Lefschetz thimble approach [58] to non-Hermitian χSB where the complex action and its gradient have branch cuts • To use the results of this paper to benchmark algorithms (such as the complex Langevin method [135]) for simulating complex-action theories • To extend this work to lower spatial dimensions • To test the conclusions of this paper with other low-energy effective models such as holography [136] (see [137] for a non-Hermitian extension)