Dirac spectrum and the BEC-BCS crossover in QCD at nonzero isospin asymmetry

For large isospin asymmetries, perturbation theory predicts the QCD ground state to be a superfluid phase of $u$ and $\bar{d}$ Cooper pairs. This phase, which is denoted as the BCS phase, is expected to be smoothly connected to the standard phase with Bose-Einstein condensation (BEC) of charged pions at $\mu_I\ge m_\pi/2$ by an analytic crossover. A first hint for the existence of the BCS phase, which is likely characterised by the presence of both, deconfinement and charged pion condensation, is coming from the lattice observation that the deconfinement crossover smoothly penetrates into the BEC phase. To further scrutinize the existence of the BCS phase, in this proceedings article we investigate the complex spectrum of the massive Dirac operator in 2+1-flavor QCD at nonzero temperature and isospin chemical potential. The spectral density near the origin is related to the BCS gap via a generalization of the Banks-Casher relation to the case of complex Dirac eigenvalues (derived for the zero-temperature, high-density limits of QCD at nonzero isospin chemical potential).


Introduction
Quantum chromodynamics (QCD) is established as the fundamental theory governing nuclear matter and hadrons. It holds that hadrons are made from quarks, their antimatter siblings, and gluons that carry the strong/color force binding quarks to each other while being themselves color charged objects. The fundamental laws of QCD are elegantly concise, however, understanding the structural complexity of hadrons in terms of quarks and gluons governed by those laws remains an open challenge.
In the highest energy RHIC and LHC collisions, strongly interacting matter at the relevant energies contains almost as many antiquarks as quarks, which, borrowing condensed matter physics nomenclature, one could call "undoped strongly interacting matter". However, there exist many physical settings, like non-central heavy-ion collisions, the structure of compact stars and the evolution of the early Universe, where instead strongly interacting matter is doped with e.g. an excess of down quarks over up quarks. In the physical systems mentioned above this translates into an excess of neutrons over protons or positively charged pions over negatively charged pions. To understand these systems, we must map the phase diagram of QCD as a function of both temperature and "doping", which we can express in terms of (negative) isospin density n I = n u − n d or, equivalently, in the grand canonical approach to QCD, in terms of isospin chemical potential µ I = (µ u − µ d )/2.
While systems of isospin-asymmetric matter are typically characterized by nonzero baryon density too, that is they also carry an excess of matter over antimatter encoded in a nonzero baryon chemical potential µ B , switching on µ B would hinder direct lattice simulations due to the complex action problem. As a first step it is then certainly useful to study the QCD phase diagram in the (T, µ I ) plane at µ B = 0, which has the advantage of being fully accessible to standard lattice Monte Carlo techniques in principle. As anticipated by perturbation theory and model calculations [1,2] ( Fig. 1(a)), lattice simulations found [3,4] ( Fig. 1(b)), an interesting and complex structure with at least three phases.
The main questions we address in this proceedings revolves around the existence of the BCS phase and the location of its boundaries. The existence of a BCS phase is expected because perturbation theory, which is applicable in the limit |µ I | Λ QCD , predicts that the attractive gluon interaction forms pseudoscalar Cooper pairs of u andd quarks at zero temperature [1]. Model calculations also confirmed the existence of a BCS phase at nonzero temperature (see e.g. [2]). The transition between the BEC phase and the BCS phase is expected to be an analytic crossover, given that the symmetry breaking pattern is the same. Lattice simulations also show large values for the Polyakov loop within the BEC phase [3]. Those can be considered to be a hint for a superconducting ground state with deconfined quarks, that is for the BCS phase. Here, we propose to look for further signatures for the existence and location of the BCS phase in the complex spectrum of the Dirac operator, which can be related to the BCS gap in a Banks-Casher type relation [5] (cf. Eq. (6)).

Simulation setup and observables
The fermion matrix M ud within the action S ud =ψM ud ψ for the light quarks ψ = (u, d) in Euclidean spacetime and in the continuum, reads (1) Here A µ is the gluon field and τ a are the Pauli matrices. Note that besides the terms including the isospin chemical potential µ I and the light quark mass m ud , in M ud there also is an explicit symmetry breaking term including the parameter λ, referred to as pionic source, that couples to the charged pion field π ± =ūγ 5 d −dγ 5 u. This unphysical term is needed to enable the observation of the spontaneous breaking of the continuous U τ 3 (1) symmetry and will act as a regulator for the simulations in the BEC phase [6][7][8]. Physical results are obtained by taking the λ → 0 limit.
For our measurements we consider 2+1-flavor QCD with µ I > 0 and λ > 0 as already simulated to map out the phase diagram shown in Fig 1(b) [3]. The lattices considered so far are N 3 s × N t lattices with N t = 6 at various temperatures T = 1/(N t a). The Dirac operator is discretized employing the staggered formulation and the rooting procedure. The partition function of this system is given in terms of the path integral over the gluon link variables U µ = exp(iaA µ ), where β = 6/g 2 is the inverse gauge coupling, S G the tree-level Symanzik improved gluon action, M ud the light quark matrix in the basis of the up and down quarks and M s the strange quark matrix, The argument of / D indicates the chemical potential µ I and η 5 the staggered equivalent of γ 5 . The positivity of the integrand of Z can be shown. In particular, both determinants in the measure of the path integral (2) are positive. The quark masses are tuned to their physical values along the line of constant physics (LCP) from Ref. [9], with the pion mass m π ≈ 135 MeV.
In the above setup our main observable is the spectrum of complex eigenvalues of the massless Dirac operator / D(µ I ). For the up quark, the eigenproblem reads where the eigenvalues ν n are complex numbers. The eigenproblem for the down quark can be obtained from Eq. (4) using chiral symmetry, i.e. / D(µ I )η 5 + η 5 / D(µ I ) = 0, and hermiticity, i.e. η 5 / D(µ I )η 5 = / D(−µ I ) † and reads The fact that [ / D(µ I ), / D † (µ I )] = 0, i.e. that / D(µ I ) is not a normal operator, entails that its left and right eigenvectors do not coincide. However, following Eqs. (4) and (5), for each eigenvalue in the up quark sector there is a complex conjugate pair in the down quark sector.  Our choice of observable is motivated by the extension of the Banks-Casher relation to the case of complex Dirac eigenvalues derived in Ref. [5] for the zero-temperature, high-density limits of QCD at nonzero isospin chemical potential. The derived Banks-Casher type relation for massless quarks reads This relation gives us a prescription on how to obtain information on the BCS gap ∆ from the density of the complex Dirac eigenvalues extrapolated at the origin ρ(0). In Ref. [5] the main idea for the extension of the Banks-Casher relation for T = 0 and |µ I | Λ QCD is to write down the partition function Z (M) as a function of the quark mass matrix M, both in the fundamental QCD-like theory and in the corresponding effective theory. Taking suitable derivatives then yields an expression proportional to ρ(0) in the fundamental theory and ∆ 2 in the effective theory. The Banks-Casher-type relation is obtained by identifying these results which leads to Eq. (6). A similar relation is also expected to hold at nonzero quark masses and temperatures.

Results
To solve the eigenproblem of Eq. (4) we employed the Scalable Library for Eigenvalue Problem Computations (SLEPc) [10], which is a software package for the solution of large sparse eigenproblems on parallel computers. The solver used for the eigenvalue problem is a Krylov-Schur solver whose implementation within SLEPc is suited for non-Hermitian problems. We compute about 150 eigenvalues of the non-hermitian Dirac operator, which are the closest (in modulo) to the origin.
Since the simulations are carried out for physical pion masses, away from the chiral limit (i.e. m ud = 0), we try to extrapolate the density ρ(ν) to m ud + i · 0 rather than to zero neglecting, at first, possible corrections due to non-zero masses and temperatures. We evaluate ρ(m ud ) by using kernel density estimation (KDE), a non-parametric way to estimate the multivariate probability density function from  . For the latter case the shaded areas correspond to the range of values P ren. takes within the green BEC boundary in Fig. 1(b) at the two considered temperatures.
the measured spectrum. Such technique is implemented in the python library scikit-learn [11], which we employ for the analysis. It can be observed, by inspecting the contour plots in Fig. 2, how only for µ I large enough, i.e. within the BEC phase, the spectrum is wide enough in the real direction to encompass the red dot in Fig. 2 at m ud resulting in ρ(m ud ) = 0. At µ I < m π /2 the eigenvalues are, instead, clustered along the imaginary axis and ρ(m ud ) = 0. At the largest simulated µ I values there is a tendency ρ(m ud ) → 0 due to the drift of the eigenvalues away from the real axis. However it should be noted that the impact of cutoff effects for larger and larger µ I values remains to be assessed by a systematic comparison with results on finer lattices.
Quantitative results for the spectral density are shown in Fig. 3. It is interesting to match the µ I -and Tdependence of ρ(m ud ) with the location of the boundary of the BEC phase as determined by the onset of the pion condensate Σ π and with the location of the deconfinement crossover within the BEC phase as hinted for by a specific value of the renormalized Polyakov loop, that is P ren. = 1. This is done both in Fig. 3(b) and in Fig. 4 by using results for Σ π and P ren. obtained in the same setup in Ref [3].
What can be observed is that the signal for the extrapolated spectral density seems to become nonzero around µ BEC I (T), that is at the location of the BEC phase boundary for the considered temperature. However, results also show that the extrapolated spectral density drops to zero again at larger values of µ I . Notice that lattice artefacts are expected to suppress ρ(m ud ), just as they do with Σ π [6][7][8]. Disentangling the signal for the BCS-BEC crossover from discretization errors at large µ I is therefore difficult and a more systematic study is certainly needed to draw realistic conclusions. As one can see from Fig. 5, the extrapolated spectral density still shows significant volume effects as well as a dependence of the results on the pionic source λ.

Discussion and conclusions
The presented results clearly show that the extrapolated spectral density is sensitive to the BEC boundary. However, to be able to draw conclusions on whether, and in which µ I range, there is sensitivity to the BEC-BCS crossover as well, a more systematic analysis is needed. Such analysis will allow us to establish the expected quantitative connection between the measured density and the BCS gap. Larger volumes, finer lattice spacings and a λ → 0 extrapolation must be considered, and this is ongoing work. Finer lattices, in particular, will help us identifying lattice artefacts due to cutoff effects at large µ I . Moreover, given that the Banks-Casher relation that we intend to use as a prescription to connect the spectral density with the BCS gap is strictly valid only for T = 0 and in the |µ I | Λ QCD limit, a generalization of this relation away from this limit is desired. In addition, we might have to consider larger isospin chemical potentials and smaller temperatures.