**1. 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 *nI* = *nu* − *nd* 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, that is, they also carry an excess of matter over antimatter encoded in a nonzero baryon chemical potential *μ<sup>B</sup>*, 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

(*<sup>T</sup>*, *μ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] (Figure 1a), lattice simulations found [3,4] (Figure 1b) an interesting and complex structure with at least three phases.

The main questions we address revolve around the existence of the BCS phase and the location of its boundaries. The existence of a BCS phase is expected due to perturbation theory, which is applicable in the limit |*μI*| <sup>Λ</sup>*QCD* and predicts that the attractive gluon interaction forms pseudoscalar Cooper pairs of *u* and ¯*d* 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. Equation (6)).

**Figure 1.** Conjectured (**a**) [1,2] and measured (**b**) [3,4] phase diagram of quantum chromodynamics (QCD) at pure isospin chemical potential.

#### **2. Simulation Setup and Observables**

The fermion matrix M*ud* within the action *Sud* = *ψ*¯M*ud ψ* for the light quarks *ψ* = (*<sup>u</sup>*, *d*) in Euclidean spacetime and in the continuum, reads

$$\mathcal{M}\_{\rm ud} = \gamma\_{\mu} (\partial\_{\mu} + iA\_{\mu}) \, \mathbb{1} + m\_{\rm ud} \mathbb{1} + \mu\_{I} \gamma\_{4} \, \tau\_{3} + i\lambda \gamma\_{5} \, \tau\_{2} \,. \tag{1}$$

Here *<sup>A</sup>μ* is the gluon field and *τa* are the Pauli matrices. Note that in addition to the terms, including the isospin chemical potential *μI* and the light quark mass *mud*, in M*ud* there is also an explicit symmetry breaking term including the parameter *λ*, referred to as pionic source, that couples to the charged pion field *π*<sup>±</sup> = *<sup>u</sup>*¯*γ*5*d* − ¯*dγ*5*u*. This unphysical term is needed to enable the observation of the spontaneous breaking of the continuous U*<sup>τ</sup>*3 (1) symmetry and will act as a regulator for the simulations in the BEC phase [6–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 Figure 1b [3]. The lattices considered so far are *N*<sup>3</sup> *s* × *Nt* lattices with *Nt* = 6 at various temperatures *T* = 1/(*Nta*). 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 *<sup>U</sup>μ* = exp(*iaAμ*),

$$\mathcal{Z} = \int \mathcal{D}l I\_{\mu} \, e^{-\oint \mathcal{S}\_{\mathcal{G}}} \left( \det \mathcal{M}\_{\text{ud}} \right)^{1/4} \left( \det \mathcal{M}\_{\text{s}} \right)^{1/4} \,, \tag{2}$$

where *β* = 6/*g*<sup>2</sup> is the inverse gauge coupling, *SG* 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,

$$
\mathcal{M}\_{\rm nd} = \begin{pmatrix}
\mathcal{D}(\mu\_I) + m\_{\rm ud} & \lambda \eta\_5 \\
\end{pmatrix}, \qquad \mathcal{M}\_{\rm s} = \mathcal{D}(0) + m\_{\rm s}.\tag{3}
$$

The argumen<sup>t</sup> 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 in Equation (2) are positive. The quark masses are tuned to their physical values along the line of constant physics (LCP) from [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 *<sup>D</sup>*/(*μI*). For the up quark, the eigenproblem reads

$$D(\mu\_I)\,\psi\_n = \nu\_n\,\psi\_n.\tag{4}$$

where the eigenvalues *νn* are complex numbers. The eigenproblem for the down quark can be obtained from Equation (4) using chiral symmetry, i.e., *<sup>D</sup>*/(*μI*)*η*5 + *η*5*D*/(*μI*) = 0, and hermiticity, i.e., *η*5*D*/(*μI*)*η*5 = *<sup>D</sup>*/(−*μI*)† and reads

$$
\hat{\psi}\_n^\dagger \mathcal{D}(-\mu\_I) = \hat{\psi}\_n^\dagger \,\nu\_n^\*, \qquad \hat{\psi}\_n = \eta\_5 \psi\_n. \tag{5}
$$

The fact that [*D*/(*μI*), *<sup>D</sup>*/†(*μI*)] = 0, i.e., that *<sup>D</sup>*/(*μI*) is not a normal operator, entails that its left and right eigenvectors do not coincide. However, following Equations (4) and (5), for each eigenvalue in the up quark sector there is a complex conjugate pair in the down quark sector.

**Figure 2.** Contour plots of the complex spectrum of the Dirac operator as obtained for a *Ns* = 24, *Nt* = 6 lattice at *λ*/*mud* ∼ 0.29, and *T* = 155 MeV for isospin chemical potentials *μI*/*<sup>m</sup>π* = 0.61 (**a**), *μI*/*<sup>m</sup>π* = 0.91 (**b**), and *μI*/*<sup>m</sup>π* = 2.30 (**c**). The red dot indicates *mud* + *i* · 0.

Our choice of observable is motivated by the extension of the Banks-Casher relation to the case of complex Dirac eigenvalues derived in [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

$$
\Delta^2 = \frac{2\pi^3}{9} \rho(0). \tag{6}
$$

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 [5], the main idea for the extension of the Banks-Casher relation for *T* = 0 and |*μI*| <sup>Λ</sup>*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 Δ<sup>2</sup> in the effective theory. The Banks-Casher-type relation is obtained by identifying these results which leads to Equation (6). A similar relation is also expected to hold at nonzero quark masses and temperatures.
