Continuum Goldstone spectrum of two-color QCD at finite density with staggered quarks

We carry out lattice simulations of two-color QCD and spectroscopy at finite density with two flavors of rooted-staggered quarks and a diquark source term. As in a previous four-flavor study, for small values of the inverse gauge coupling we observe a Goldstone spectrum which reflects the symmetry-breaking pattern of a Gaussian symplectic chiral random-matrix ensemble (GSE) with Dyson index $\beta_D=4$, which corresponds to any-color QCD with adjoint quarks in the continuum instead of QC$_2$D wih fundamental quarks. We show that this unphysical behavior occurs only inside of the bulk phase of $SU(2)$ gauge theory, where the density of $Z_2$ monopoles is high. Using an improved gauge action and a somewhat larger inverse coupling to suppress these monopoles, we demonstrate that the continuum Goldstone spectrum of two-color QCD, corresponding to a Gaussian orthogonal ensemble (GOE) with Dyson index $\beta_D=1$, is recovered also with rooted-staggered quarks once simulations are performed away from the bulk phase. We further demonstrate how this change of random-matrix ensemble is reflected in the distribution of eigenvalues of the Dirac operator. By computing the unfolded level spacings inside and outside of the bulk phase, we demonstrate that, starting with the low-lying eigenmodes which determine the infrared physics, the distribution of eigenmodes continuously changes from the GSE to the GOE one as monopoles are suppressed.


I. INTRODUCTION
The QCD phase diagram continues to be subject of intense theoretical and experimental studies. The region of high baryon density at relatively low temperatures is of particular relevance for the inner cores of neutron stars, and at somewhat higher temperatures for neutron-star mergers. It is probed experimentally in the beam-energy scan at RHIC and the future heavy-ion programs at J-PARC, NICA and FAIR. In this regime one usually expects a chiral first-order transition ending a critical point, but inhomogeneous phases or more exotic states of matter like a quarkyonic phase have been proposed to occur as well. At even higher densities, beyond reach of current experiments and astrophysical observations, asymptotic freedom and the attractive perturbative interactions between quarks close to the Fermi surface entail the formation of Cooper pairs and color superconductivity.
Unfortunately, QCD at high densities remains inaccessible to stochastic integration methods, since the fermion determinant becomes complex at finite chemical potential µ. This leads to an insurmountably hard fermion-sign problem precisely where a finite baryon density starts to build up in the ground state. The problem does not arise on the other hand in certain QCD-like theories, e.g. two-color QCD (QC 2 D) or G 2 -QCD [2,3], which show chiral symmetry breaking, confinement and asymptotic freedom as well. These theories can thus be approached with standard Monte-Carlo techniques on the lattice and provide therefore interesting testbeds to develop and test algorithms for QCD at finite density.
Beside this technical aspect, QCD-like theories are interesting in their own right. On the lattice, QC 2 D has been studied with staggered [1,[4][5][6][7][8][9][10][11][12] and Wilson fermions [13][14][15][16][17][18]. In contrast to QCD, the color-singlet baryons are diquarks and hence bosonic in two-color QCD, for example, while fermionic baryons do not exist in its spectrum. The physics of the bosonic diquark baryons qualitatively resembles QCD at finite isospin density with pion condensation [19] and is by now fairly well understood [20][21][22]. There are firm predictions for diquark condensation when µ reaches half the pion mass m π from chiral effective field theory and random matrix theory [6,[23][24][25][26][27][28][29], and model studies of the BEC-BCS crossover inside the condensed phase [30][31][32][33]. In QC 2 D the lightest diquarks play a dual role as two-color baryons and pseudo-Goldstone bosons of the dynamical breaking of an extended chiral symmetry. When they condense, they are expected to form a superfluid which changes in nature from a Bose-Einstein condensate of tightly bound diquarks to a BCS-like pairing of quarks as chiral symmetry gets gradually restored with increasing density.
The phase diagram of the quark-meson-diquark model for QC 2 D from the functional renormalization group [32] is shown in Figure 1. The qualitative features resemble lattice results [16,17], especially when the Polyakov-loop variable is included in the effective model description [33]. Evidence of the BEC-BCS crossover inside the diquarkcondensation phase was also provided from lattice simulations [9] albeit still close to the bulk phase of SU (2) (see below).
Anti-unitary symmetries of the Dirac operators in QCD-like theories without a fermion-sign problem are both, a blessing and a curse for the rooted-staggered fermion formulation. On one hand, phase ambiguities when rooting a complex determinant [34] do not occur.  [32,33]: half-value of the chiral condensate (red), and second-order phase boundary for diquark condensation (solid blue) with rough indication of the BEC-BCS crossover (dotted blue).
In fact, even for a single staggered fermion the determinant remains positive at finite µ, whereas this requires two flavors of continuum Dirac fermions. This is because of the missing C 2 = −1 from the charge conjugation matrix C for Dirac spinors in the anti-unitary symmetries of staggered Dirac operators. On the other hand, it implies that the corresponding Gaussian chiral random matrix ensembles get swapped, staggered fermions reflect the behavior of the Gaussian symplectic ensemble (GSE) when the continuum Dirac fermions show that of the Gaussian orthogonal ensemble (GOE) and vice versa. In particular, the staggered Dirac operator of fundamental quarks in QC 2 D has the GSE Dyson index β D = 4, while for continuum or Wilson fermions it is the GOE one, β D = 1. For adjoint quarks in any-color QCD (or fundamental quarks in G 2 -QCD) it is just the other way round. This is why the sector of positive fermion determinant of QC 2 D with adjoint quarks was studied as a replacement for the continuum theory within the correct random matrix ensemble in the early days [5].
Here we have addressed the following question: with the full SU (4) taste symmetry in the continuum limit, it should be possible to define a standard (tasteless) charge conjugation from that for staggered quarks [35] as well. Does this imply that the correct symmetry breaking pattern, corresponding to the random matrix ensemble of the continuum two-color Dirac operator is recovered also with staggered quarks in the continuum limit?
The answer seems to be positive. The particular evidence for this that we provide is the behavior of the Goldstone pion inside the diquark condensation phase which shows the characteristic change indicative of the change of the Dyson index as the continuum limit is approached: In the continuum, the extended SU (2N f ) chiral symmetry is dynamically broken down to the compact sym-plectic Sp(N f ) with fundamental quarks and Dyson index β D = 1 [24]. 1 For N f = 2 on the bosonic level this amounts to the simple vector-like breaking of SO(6) → SO(5) with coset S 5 and five Goldstone bosons, three pions and a scalar (anti-)diquark pair. The exact chiral symmetry extending the usual U (1) e × U (1) o of the staggered action in the two-color case is U (2) on the other hand [4]. With fundamental quarks it breaks down to U (1) V . Up to the Goldstone pion of the broken U (1) it therefore resembles the N f = 1 case of continuum quarks with Dyson index β D = 4 as in adjoint QCD [24] or G 2 -QCD [36] which is SU (2) → U (1). 2 The effective field theory prediction for the Goldstone spectrum with the U (2) → U (1) chiral symmetry breaking of the fundamental staggered two-color action was explicitly worked out in [1]. Most important for our purposes is the behavior of the Goldstone pion inside the diquark condensation phase at µ ≥ m π /2 where it resembles that of the symmetric pion branch P S of the Dyson index β D = 4 case in the continuum [24], with a mass that decreases ∼ m 2 π /2µ, although this branch strictly speaking only exist for N f ≥ 2 there.
We will demonstrate that the behavior of this Goldstone pion branch in the continuum limit indeed changes to that of the symmetric P S mode in the β D = 1 case of the continuum theory which increases ∼ 2µ inside the diquark condensation phase [24]. Again this mode exists in the continuum only for N f ≥ 2, starting with a multiplicity of 3 for the three degenerate pions in the two-flavor case. This again indicates that the taste symmetry needs to be at least partially restored to achieve this.
Here it is also important to note that the previous results with staggered quarks in QC 2 D that led to the β D = 4 Goldstone spectrum were obtained inside a socalled bulk phase on the lattice, where the lattice spacing is almost independent of the inverse gauge coupling. The large effects of this on spectroscopy and thermodynamics at vanishing chemical potential have been investigated in Ref. [37]. The influence of bulk effects on simulations at finite density has not yet been discussed and is thus the main focus of our study. Indeed, we will find that inside the bulk phase, the Goldstone spectrum for Dyson index β D = 4 as for adjoint continuum quarks is reproduced, while outside the bulk phase we observe the correct Goldstone spectrum of two-color QCD, corresponding to the Dyson index β D = 1 of the continuum Dirac operator, also with rooted-staggered quarks. We demonstrate that this change of Dyson index is reflected in the eigenvalue statistics of the Dirac operator. By computing the unfolded level spacings, we demonstrate that the distribu-tion of eigenmodes is completely dominated by the GSE in the bulk phase, but obtains larger contributions from the GOE as one leaves the bulk phase. This change turns out to be continuous, and builds up starting with the low-lying eigenmodes. As infrared physics is controlled by these lower levels, we observe that one correctly reproduces the continuum theory Goldstone spectrum even when the higher eigenmodes are still dominated by the GSE.
This paper is organized as follows: In Section II we introduce the lattice action, including a diquark source term, and the simulation parameters used in this work. In Section III we study the quark-number density and chiral and diquark condensates for small values of the gauge coupling (which turns out to be inside the bulk phase) and compare our results to the predictions of leading order chiral perturbation theory. Section IV explains the bulk phase and introduces its order parameter, the Z 2 -monopole density, ending with a discussion on new lattice parameters to suppress bulk effects. In Section V we then repeat the study of Section III with our new set of parameters and discuss lattice discretization and finite volume effects, such as the effective quenching of the theory in the saturated regime and problems with additive renormalization of the chiral condensate at finite µ. In Section VI after providing a detailed description of chiral symmetry-breaking pattern of staggered fermions, we present our main result, which is the numerical measurement and comparison of the Goldstone spectrum at non-vanishing chemical potential, inside and outside the bulk phase. And finally, in Section VII we present the unfolded level spacings of the Dirac operator and demonstrate that change of Goldstone spectrum is accompanied by a change of eigenvalue distribution. We end with our conclusion and an outlook in Section VIII.

II. LATTICE SETUP
For the low-temperature scan of the Goldstone spectrum of two-color QCD, we use standard rooted staggered fermions, with the Dirac operator at non-vanishing baryon chemical potential µ. To study spontaneous symmetry breaking and competing order on a finite lattice, we add a diquark source λ corresponding to a Majorana mass term in the Lagrangian, which explicitly breaks the chiral U (2) symmetry of the massless staggered action down to a U (1) as well but in a direction different from that of the Dirac mass term [4][5][6][7], Physical results are then retrieved in the λ → 0 limit. The diquark condensate is obtained from The staggered fermion action is conveniently expressed in a Nambu-Gorkov basis with an inverse Nambu-Gorkov propagator Grassmann integration produces the square root of the determinant of A in the path integral measure. This is seen most easily when considering (χ,χ T ) as a set of independent real Grassmann variables, whose Gaussian integral results in a (positive) Pfaffian which agrees with In fact, a direct use of hybrid Monte-Carlo (HMC) based on a Gaussian pseudo-fermion integral over (AA † ) −1 would even produce (det A) 2 in the measure. However, because is block diagonal with both blocks having the same determinant, we can use size-half pseudo-fermion fields (in the Nambu-Gorkov space) to remove this further doubling. On the other hand, we can not use size-half fields in the even-odd staggered lattice at finite µ. This means that without any rooting, we compute det A and describe eight fermion species instead of the usual four staggered tastes. Therefore we use standard rooting techniques to compute (det A) N f /8 for each continuum flavor, i.e. we take the fourth root to simulate with N f = 2 as in Ref. [9]. The rooting is achieved by a rational approximation of the fermion matrix in the pseudofermion action of the HMC algorithm. The diquark source explicitly breaks baryon number conservation, or more precisely the U (1) V of the staggered action, and hence in addition to the usual ψψ contractions in the calculation of the correlation functions, we now also have ψψ and ψψ contractions. These correspond to the diagonal terms of the propagator G = A −1 obtained from Eq. (5), In this paper we compare results with N f = 2 staggered flavors at β = 1.5 on a 12 3 × 24 lattice with am = 0.025, which turns out to be deep inside the bulk phase, to N f = 2 staggered flavors at β = 1.7 on a 16 3 × 32 lattice with am = 0.01 and an improved gauge action so that this is just outside the bulk phase. The parameters of the simulations at β = 1.5 with the unimproved gauge action correspond to those used in Ref. [1] for N f = 4.

III. EFFECTIVE FIELD THEORY PREDICTIONS
Kogut et al. have studied the symmetries of QCDlike theories at finite baryon density with pseudoreal quarks in the fundamental representation using chiral effective Lagrangians [1,24]. A linear sigma model for the symmetry-breaking pattern and Goldstone spectrum of the staggered two-color action was used to describe the data in [1]. For the purpose of illustrating the basic features of the diquark-condensation transition at µ = µ c = m π /2 here we fit our data to the somewhat simpler form of the leading-order chiral perturbation theory (χPT) predictions from the non-linear sigma model [24]. This describes the rotation of the vacuum alignment from the chiral qq into the diquark qq condensate at a fixed With explicit diquark source λ, the rotation angle α(µ) is obtained from such that a non-zero value α 0 is obtained already at µ = 0 which depends on the relative size of the Majorana and Dirac quark masses, i.e. tan α 0 = λ/m. Chiral and diquark condensate, and quark-number density n as functions of µ are then given by As a first test we have performed simulations with the parameters of Ref. [1], i.e. a lattice gauge coupling of β = 1.5, quark mass am = 0.025, and diquark source aλ = 0.0025 on a 12 3 ×24 lattice and the standard Wilson plaquette action, so as to reproduce their results with the square root of the determinant in (6) With the fourth root for N f = 2 continuum flavors and the same lattice parameters we have then obtained the results shown in Figure 2. They are fitted to the leadingorder forms from chiral perturbation theory in Eq. (10) with (9) for the vacuum alignment angle α(µ). From these fits we obtain the fit parameters F, G, and the critical chemical potential, which results as aµ c = 0.1889 (5). Within the errors this agrees with our spectroscopic result from the pion correlator at λ = 0.0025, µ = 0 which Fit of lattice data to leading-order χPT form of chiral (purple) and diquark (orange) condensates, and quarknumber density (brown) from Eqs. (9) and (10), with N f = 2, β = 1.5, am = 0.025, aλ = 0.0025 on a 12 3 × 24 lattice with unimproved Wilson gauge action.
yields am π /2 = 0.1887 (6). The agreement between fits and data in Figure 2 is nearly perfect up to aµ ∼ 0.3. The deviations at larger chemical potentials were attributed to a µ-dependence of the total condensate Σ c in Ref. [1] where linear sigma model fits were therefore used instead. These imply qq ∼ µ and hence n ∼ µ 3 at large µ, as predicted for a BCS-like pairing in QCD at large isospin density already in [19]. The corresponding rise in the isospin density beyond the χPT prediction was observed in lattice QCD simulations [20], and it was traced to the BEC-BCS crossover in a functional renormalization group study of the quark-meson model as an effective theory with linearly realized chiral symmetry and order-parameter fluctuations beyond mean field [21].
The same interpretation of this rise in the diquark density, beyond the χPT prediction, as an indication of the BEC-BCS crossover in QC 2 D, was also adopted in Ref. [9]. Because these results were obtained on rather coarse lattices, it is therefore important to verify that they are not qualitatively affected by strong discretization artifacts such as the Z 2 monopoles in the bulk phase of SU (2).

IV. THE BULK TRANSITION
Most lattice gauge theories, as for instance SU (2), SU (3) or G 2 gauge theory, exhibit a bulk phase in the strong-coupling regime, characterized by the presence of unphysical lattice artifacts such as electric vortices and magnetic monopoles [38][39][40][41][42][43]. These dominate the ultraviolet behaviour, such that the lattice spacing is nearly independent of the coupling constant, and taking a continuum limit is not possible. In the physical weak-coupling regime the short distance physics is gov-erned by asymptotic freedom and the continuum limit is approached by β → ∞. Depending on which gauge group and which representation, both regions are either separated by a true phase boundary or a cross-over, the later being the case for the fundamental representation of SU (2) considered here (cf. Figure 1 in Ref. [42]). The bulk transition is almost independent of the lattice size and persists also in the presence of fermions [44,45].
An order parameter for the strong-coupling to weakcoupling transition is the Z 2 monopole density [40] where C runs over all elementary cubes of the lattice. z is sensitive to preferred signs of the plaquettes on the faces of these cubes, which are aligned below the bulk transition: In the bulk phase z is non-vanishing, while it vanishes in the physical weak-coupling regime.
With β = 1.5 and the unimproved Wilson gauge action, we have found the Z 2 monopole density at aµ = 0 to be z = 0.884040(95), and thus we expect bulk effects to be dominant in this regime. We have also confirmed that z is only very weakly affected by the inclusion of dynamical quarks. For a fixed inverse gauge coupling β, the monopole density can be significantly reduced with Symanzik's gauge action [46] (we employ the tree-level improved variant here). This is illustrated in Figure 3, which shows the β-dependence of z for the improved and unimproved actions, with and without dynamical fermions in each case. In principle, one would like to suppress Z 2 monopoles as much as possible. In practice, finite volume effects become increasingly severe at larger β due to a smaller physical lattice spacing a, and one is forced to make a compromise. This is illustrated in Figure 4, which shows the β-dependence of various meson masses obtained in a previous study [37]. For small physical volumes (aN ) 4 (larger β) the meson masses degenerate, signaling an explicit breaking of chiral-symmetry by the finite system size, while for small inverse couplings β bulk effects are dominant. For our simulations of the continuum physics we therefore chose β = 1.7 with N s = 16 and N t = 32, as a compromise between small a and small β, where m π /m ρ = 0.58 (5). Using the improved action the Z 2 monopole density at β = 1.7 and aµ = 0 is z = 0.27340(66). Although there is still a substantial amount of monopoles on the lattice, this choice of parameters pushes the simulations on the weak-coupling side of the bulk crossover, as our spectroscopic results discussed below clearly demonstrate.

V. LEAVING THE BULK PHASE
Simulating at finite µ with the improved action and our new choice of lattice parameters (β = 1.7, am = 0.01, 16 3 × 32), we carry out a study of the µ dependence of different observables outside of the bulk phase. We first observe that we can again fit the quark-number density and diquark condensate to leading-order χPT predictions from Eqs. (9) and (10) (note here that the same expressions are predicted by χPT for both the Gaussian orthogonal and Gaussian symplectic ensembles [1,24]). We simulate with three different values of the explicit diquark source aλ = 0.0050; 0.0025; 0.0010 and apply the fits directly at finite λ. Results are shown in Fig. 5, together with extrapolations to λ → 0.
We find that our results for aλ = 0.0010 already agree within one standard deviation with the limit of vanishing diquark source. Attempting to extract the critical chemical potential for diquark condensation from fits to aλ = 0.0010 yields aµ c = 0.172(21) however, which is slightly overestimated compared to our spectroscopic result (aµ c = 0.1456 (28)) from the pion correlator at aλ = 0.0010, µ = 0. We find that a consistent value (aµ c = 0.1356(86)) is obtained from a χPT fit to the λ → 0 extrapolation of qq . We also observe significant deviations from the χPT predictions at around aµ ∼ 0.3, which we interpret as signaling the onset of the BEC-BCS crossover. We thus conclude that the behavior of qq and n is not qualitatively different outside of the bulk phase.
On the other hand, from an observed decrease of the chiral condensate above µ c (shown in Figure 6), which we did not observe at β = 1.5, we infer the presence of UVdivergence, such that renormalization is required. As discussed in Ref. [47], it is possible to renormalize the chiral condensate at finite temperature using the chiral susceptibility χ mq , as both contain the same UV-divergent term c U V , viz. Since the UV divergence originates mainly from the connected chiral susceptibility χ con (also shown in Figure 6, we neglect λ-dependent contributions here), a renormalized condensate can be defined as Σ = qq mq − m q χ con . We observe that both the condensate and χ con exhibit a similar decrease at large µ and thus conclude that the UV-divergence c U V is µ-dependent. At the chiral transition the disconnected susceptibility χ dis contains a singular contribution. At the diquark condensation transition, we find that χ con has a singular part as well, as its peak height is bounded from above by a finite volume (see Figure 7). This singularity will dominate over c U V /a at finite a in the infinite volume limit. The chiral condensate on the other hand, at zero temperature, must remain independent of µ for µ < µ c . It does not have such a singular contribution and it would be unphysical to introduce one with the connected susceptibility subtraction. At any rate, this would introduce a µ-dependence below µ c and hence a Silver-Blaze problem. Therefore, a different (µ-dependent) subtraction of the chiral condensate is required. Likewise, it is impossi-  ble to remove the UV-divergence by subtracting a heavy quark condensate like qq mq − mq m q qq m q , since the pion mass and thus the position of the diquark onset strongly depend on the quark mass.
Measuring the µ dependence of the Z 2 monopole density and the quark number density, we observe that both quantities saturate at large µ (see Figure 8 and 9). In these figures, µ has been normalized with the critical chemical potential µ c = m π /2. With increasing µ the Z 2 monopole density approaches its quenched value, while at the same point the quark number density saturates. We conclude that in this high chemical potential regime the lattice is fully occupied with fermions, such that the system effectively becomes quenched.
Finally, we observe that the Polyakov loop is rather insensitive to the chemical potential with staggered quarks, and in fact coincides with its value in the quenched limit (see Figure 9). This is in contrast to lattice simulations of two-color QCD with Wilson fermion [10,16] and G 2 -QCD with Wilson fermions [2], where the Polyakov loop shows a peak around half filling. Also, in a previous effective Polyakov loop model study for QCD-like theories [48,49] with heavy Wilson quarks it has been seen that the Polyakov loop expectation value has a peak at the inflection point of the quark number density.
To explain this discrepancy, we consider that in twodimensional two-color QCD, where large temporal extends of the lattice are feasible, the peak vanishes in the limit of N t → ∞ [50], suggesting that the non-vanishing Polyakov loop might be an effect of the residual temperature due to the finite lattice volume. This is also in agreement with recent lattice simulations with staggered fermions at zero temperature but larger inverse gauge coupling, leading to a larger residual temperature [51].
Here the Polyakov loop also increases with increasing chemical potential. For Wilson fermions, a larger lattice spacing might lead to the excitement of heavy doublers beyond some aµ, such that the free energy becomes finite. As these heavy doublers are not present in the staggered formalism, this behaviour is not observed here at comparable lattice spacings.

VI. CHIRAL SYMMETRY BREAKING PATTERN AND THE GOLDSTONE SPECTRUM
Having established a set of parameters (β = 1.7, 16 3 × 32, am = 0.01) with which we expect to reproduce the continuum physics, we now turn to the primary focus of this work which is to study the Goldstone spectrum. We begin by reviewing the symmetry-breaking channels of the staggered action of two-color lattice QCD and discussing the associated Goldstone modes and correlation functions, which will then be compared inside and outside of the bulk phase. These issues were previously discussed in Ref. [4]. We present a compact summary here to keep this paper self-contained.
For this purpose, it is convenient to introduce a new basis for the fermion fields, given bȳ which separates even and odd sites. The kinetic part of the staggered action (4) then reads [4] S kin = n∈Λ ,ν η ν (n) 2 X e (n) e µδν,4 0 0 e −µδν,4 U ν (n)X o (n +ν) −X e (n) where the sum runs over even sites only. It follows that in the limit m = λ = µ = 0 the fermion action is invariant under The original U (1) e × U (1) o symmetry of the N f = 1 staggered action for two-color QCD is therefore enlarged to U (2) in this limit [1,4,7]. Applying the same basis transformation to the mass and the diquark source terms one obtains where we understand Pauli matrices σ i to act in the basis (13) and τ i to act on color indices. Hence, the condensates are indistinguishable at µ = 0 as they are connected [4]. The Goldstone modes are derived by applying infinitesimal U (2) rotations toχχ and χχ, where the coefficient of O(δ) is then identified as the Goldstone mode [4]. The results are shown in Table I. Both condensates leave one generator of U (2) unbroken and hence induce the same symmetry-breaking pattern U (2) → U (1) at µ = 0. Since at µ = 0 the symmetry is reduced from U (2) to U (1) e × U (1) o , one is left with two generators {1, σ 3 } of U (2), which correspond to the staggered U (1) and baryon number conservation, respectively.
The above can be generalized to N f > 1 staggered fermions. There one has a U (1) × U (1) symmetry for each flavor, leading to U (N f )×U (N f ), which is extended to U (2N f ) at µ = 0. The same Goldstone modes listed in Table I appear also for N f > 1 (and in particular for the N f = 2 case considered in this paper), but with different multiplicities. The full symmetry-breaking pattern is summarized in Fig. 10. It can be seen that anycolor QCD with quarks in the adjoint representation in the continuum exhibits the same pattern of symmetry breaking, with an additional breaking of U (1) A due to the axial anomaly [24]. The Goldstone spectrum consists of two meson modes, the (pseudoscalar) pion π and the scalar meson f 0 , and two diquark modes, a scalar diquark qq and a pseudoscalar diquark qq. In Table II, we show the employed interpolating operators for the f 0 and π modes taken from Ref. [52]. It is important to realize that channel 1 not only contains the desired scalar meson, but also an excited pion. However, these two states can be separated during the fitting procedure as they have opposite parity. The groundstate pion is exclusively contained in channel 2. The interpolating operators from [1] are employed for the (pseudo-)scalar diquark modes, shown in Table III These modes furthermore contain contributions from their corresponding anti-diquarks. Nevertheless, the anti-diquark modes become less important with increasing chemical potential as the propagation of particles is favored (e aµ ) over the propagation of anti-particles (e −aµ ) at non-vanishing chemical potential.
Channel Operator J P C States 1χχ 0 ++ f0  To extract the ground state masses of the particle states, we employ the zero-momentum projected corre-lations functions of the form Note that only the connected contributions are considered here. For the different channels shown in Tables II  and III we obtain: • Channel 1 -Scalar Meson • Channel 2 -Pion / Pseudoscalar Meson • Channel 3 -Scalar Diquark • Channel 4 -Pseudoscalar Diquark We used the notation G = D † D + λ 2 −1 D † here, which corresponds to only the off-diagonal terms in Eq. (7). The diagonal terms produce corrections of order O(λ 2 ), which we refrain from listing here explicitly as their contributions are negligibly small for the values of λ considered in this work. We did in fact include these corrections for the results shown in Figs. 11 and 13. The results in Fig. 12 were obtained without them.

A. Goldstone spectrum in the bulk phase
We now study the (pseudo) Goldstone spectrum on a 12 3 × 24 lattice at β = 1.5 with quark mass am = 0.025 and diquark source λ = 0.0025. As the combined condensate qq 2 + qq 2 rotates from a chiral to a diquark condensate with increasing chemical potential, we assume that the Goldstone modes corresponding to a given U (2) generator mix in a similar way, i.e. rotate into each other with the same rotation angle α(µ), as described by Eq. (9). Hence, we introduce the two combined modes qq/f 0 : 1 2 χ T τ 2 χ +χτ 2χ T cos α +χχ sin α π/ qq:χ χ cos α + 1 2 χ T τ 2 χ +χτ 2 χ T sin α . Note that this simple mode mixing holds only to leading order χPT. Note also that in addition to the terms of order O(λ 2 ) discussed above, the correlators of the combined modesqq/f 0 and π/ qq contain additional terms of order O(λ) from the mixed parts (i.e. the terms ∼ sin α cos α), as diagonal elements of the full propagator in Nambu-Gorkov space (7) contribute to these. We again refrain from listing these terms here explicitly but included them in our simulations.
To obtain the masses of the (pseudo) Goldstone modes, we measure the zero-momentum projected connected correlation functions and extract their masses from fits to cosh(·) (which combines exponential decays forwards and backwards in Euclidean time). For π/ qq a factor (−1) t is inserted to account for negative parity. In the case of theqq/f 0 mode we have to use the fitting function as the correlator contains a contribution from the scalar diquark mode qq in addition to the mixing of the scalar anti-diquarkqq with the scalar meson f 0 . At µ = 0 the combined modeqq/f 0 is a massive pseudo Goldstone mode for all values of µ. Including a nonvanishing quark mass m = 0, also the combined mode π/ qq becomes a pseudo Goldstone mode as the U (1) A symmetry gets broken. The only true Goldstone mode is given by the scalar diquark mode qq in the limit λ → 0 and for µ > µ c . Generally for µ < µ c and λ → 0, the pion mass m π stays constant as the pion does not carry a net Baryon number, whereas the scalar diquark mass m qq decreases like m π − 2µ and the scalar anti-diquark mass mqq increases like m π + 2µ.
We compare our obtained masses of the (pseudo) Goldstone modes to the corresponding χPT predictions for twocolor QCD with staggered quarks [1], using µ c from the fit of the condensates in Section III (see Figure 11) and the lattice parameters aλ and am as input. The large error of the combined modeqq/f 0 mainly comes from the systematic error as the double-cosh fit is more sensitive to the fitting interval than the single-cosh fits of the other modes.
We find a general agreement of the scalar diquark mode qq and the π/ qq mode to their predictions. Deviations become notable at large chemical potential, where also the quark number density in Figure 2 deviates largely. The large discrepancy of theqq/f 0 mode for µ µ c might be due to omitting disconnected contributions, as at large chemical potential the scalar meson mode f 0 dominates this combined mode. Note that the scalar diquark mode qq has disconnected contributions of order O(λ 2 ), which we expect to have a small effect. Hence, we obtain similar results as in the previous study [1] for N f = 4. However, in Ref. [1] disconnected contributions have not been omitted and thus their obtainedqq/f 0 mode coincides better with the χPT prediction. In conclusion, we find that within the bulk phase fundamental staggered quarks resemble the chiral symmetry breaking pattern of adjoint QCD or G 2 -QCD in the continuum as most notably seen in the behavior of the pion branch above the onset at aµ c .

B. Goldstone spectrum outside bulk phase
Continuum two-color QCD with quarks in the fundamental representation obeys the pattern of symmetry breaking SU (2N f ) → Sp(2N f ). Applying χPT, it was found that the mass of the pion mode increases for µ > µ c [24] in this case, due to a swapping of the P S and P A branches when compared to the staggered action [1]. We now wish to test whether the pattern of symmetry breaking on the lattice will change to the continuum pattern in the limit a → 0. Thus, we calculate the spectrum again at a larger inverse gauge coupling and using an improved gauge action.
With the new parameters (β = 1.7, N s = 16 and N t = 32), where monopoles are strongly suppressed, we find a quite different behaviour of the (pseudo) Goldstone modes than in Sec. VI A. We observe that the combined modesqq/f 0 and π/ qq do not give any meaningful results. Instead, we therefore study the scalar diquark mode qq and the pion mode π individually, and additionally we extract the mass of the scalar anti-diquark qq from the operator χ T τ 2 χ −χτ 2χ T for µ < µ c . In the correlation function of the pion mode we find an additional contribution from an opposite parity state, which we filter out by applying a double cosh-fit of the form We first study the λ-dependence pion mass and the scalar diquark mass (see Figure 12). Here we neglected O(λ 2 ) contributions from the diagonal terms in Eq. (7) in the calculation of the correlation functions. We compare the results to χPT predictions for from continuum two-color QCD with quarks in the fundamental representation, for which we use the lattice parameters λ and m, and also µ c as obtained from measuring the pion mass at µ = 0, as input. For the smallest diquark source λ = 0.001 we find excellent agreement with χPT up to µ ∼ 1.7µ c where leading order χPT is not reliable anymore. We conclude that the diagonal terms are negligible for this choice of diquark source.
At β = 1.7 the pion mass stays constant for µ µ c , but it starts to increase for µ > µ c . This differs profoundly from the behaviour at β = 1.5 (cf. Figure 11). Again, the large error of the pion mode π for µ > µ c mainly comes from the sensitivity of the double-cosh fit to the fitting interval. We interpret the increasing pion mass for µ > µ c as strong evidence that the pattern of symmetry breaking changed to its continuum counterpart outside of the bulk phase. Figure 13 combines the results for all the considered meson channels, for λ = 0.001 where the agreement with χPT is nearly perfect. In this figure the O(λ 2 ) diagonal terms were in fact included, but their contribution is of similar magnitude as the statistical error.

VII. UNFOLDED LEVEL SPACINGS
The low energy features of QCD are dominated by its global symmetries. In this limit, QCD is well described with only the two lightest quarks. Low-energy QCD can be approximated with a random matrix theory very well, where all interactions among the degrees of freedoms of the theory are equally likely and completely determined by the global symmetries. In a random matrix theory (RMT), the matrix elements of the Dirac operator are replaced with uncorrelated random numbers, in a way such that the global symmetries of the operator prevail [53]. Any observable is then an average over random matrix elements and depends only on universal features of the theory, and not on microscopic interactions.
RMT has found wide application for studying the eigenvalue statistics in nuclear resonances [54]. The distribution of level spacings is cleared off of any microscopic interactions by rescaling the spacings such that their average is unity and removing the fluctuating part of the cumulative spectral function. This procedure is called unfolding. We wish to demonstrate here that the change of the Goldstone spectrum of twocolor QCD with staggered fermions when leaving the bulk phase is accompanied by a change of the unfolded level spacing distribution (ULSD) of the Dirac operator.
For different random matrix ensembles different forms of ULSD have been predicted [55,56]. For the Gaussian orthogonal ensemble (GOE), where the matrix elements are real and the Dyson index is β D = 1, the ULSD is

The ULSD of the GSE is given by
Using the Implicitly Restarted Arnoldi Algorithm, we have measured the ULSD both with the parameters from Ref.
The ULSD within the bulk phase are shown in Figure  14. We find that the entire spectral range is well described by Eq. (27), i.e. the distributions resembles the symplectic ensemble very well. Inside the bulk phase, the ULSD does not seem to contain any component distributed according to the orthogonal ensemble.
In the continuum limit the spectrum is expected to resemble the Gaussian orthogonal ensemble. With our improved parameters, we find however that neither Eq. (27) nor Eq. (26), fully describe the ULSD. In fact, our numerical data (seen in Figure 15) indicate some sort of intermediate state. When separating the low-lying (aλ = 0.000 . . . 0.036) from the high-lying eigenmodes (aλ = 0.040 . . . 0.160) it becomes clear that a large part of the low eigenmodes are now distributed according to the GOE, while the higher eigenmodes remain distributed according to the GSE (see Figure 16). Since the low eigen- modes govern the Goldstone modes, it is clear that our spectroscopic results should reflect the chiral-symmetry breaking pattern of the GOE. We suspect that in the continuum limit the symplectic part will vanish entirely.

VIII. CONCLUSION AND OUTLOOK
In this work we investigated the influence of bulk effects in QC 2 D with staggered fermions on the Goldstone spectrum and the unfolded level spacing distribution of the Dirac operator at finite density. We compared the Goldstone spectrum to predictions from leading order chiral perturbation theory for two sets of lattice parameters: • A 12 3 × 24 lattice with β = 1.5 and bare mass am = 0.025 using a standard Wilson gauge action. The density of Z 2 monopoles at µ = 0 is z ∼ 0.88 in this case.
• A 16 3 × 32 lattice with β = 1.7 and bare mass am = 0.01 using a treelevel improved Symanzik gauge action. The density of Z 2 monopoles at µ = 0 is z ∼ 0.27 in this case.
Our main result is that the Goldstone spectrum switches from that of any-color QCD with adjoint fermions in the bulk phase to that of twocolor QCD with fundamental quarks on the physical weak-coupling side of the bulk crossover, as most notably visible in a change of the pion branch. We show that this change is reflected in the unfolded level spacing distribution, which appears to obtain a larger and larger contribution from the Gaussian orthogonal random-matrix ensemble, starting with the low-lying eigenmodes, as one moves from strong to weak compling, while deeply in the bulk phase the distribution is completely dominated by the Gaussian symplectic ensemble. We conclude that a continuum limit leading to twocolor QCD with the correct chiral symmetry-breaking pattern is possible with rooted staggered quarks. We also observed that the standard connected susceptibility subtraction to obtain a renormalized chiral condensate cannot be used at finite µ, since the connected susceptibility contains a singular contribution at the diquark condensation transition and that a renormalization using a heavy quark condensate is also rendered unfea-sible. Developing a proper µ-dependent renormalization scheme might be possible using gradient flow techniques, but is left for future work.

ACKNOWLEDGMENTS
This work was supported by the Helmholtz International Center for FAIR within the LOEWE initiative of the State of Hesse.