QCD phase diagram for nonzero isospin-asymmetry

The QCD phase diagram is studied in the presence of an isospin asymmetry using continuum extrapolated staggered quarks with physical masses. In particular, we investigate the phase boundary between the normal and the pion condensation phases and the chiral/deconfinement transition. The simulations are performed with a small explicit breaking parameter in order to avoid the accumulation of zero modes and thereby stabilize the algorithm. The limit of vanishing explicit breaking is obtained by means of an extrapolation, which is facilitated by a novel improvement program employing the singular value representation of the Dirac operator. Our findings indicate that no pion condensation takes place above $T\approx 160$ MeV and also suggest that the deconfinement crossover continuously connects to the BEC-BCS crossover at high isospin asymmetries. The results may be directly compared to effective theories and model approaches to QCD.


I. INTRODUCTION
The thermodynamics of strongly interacting matter, as described by quantum chromodynamics (QCD), plays a characteristic role for the phenomenology of heavy-ion collisions, the structure of compact stars and the evolution of the early Universe. The relevant parameters that impact QCD physics in these settings include the temperature and the net densities of the individual quark flavors . In the light quark sector ( = , ), the characteristic combinations are the baryon density = ( + )⇑3 and the isospin density = − . While the former measures the overall excess of strongly interacting matter over antimatter, the isospin density describes the asymmetry between the up and down quarks or, equivalently, between protons and neutrons in the system.
The above-mentioned physical settings are all thought to exhibit a strong isospin asymmetry. The initial state of typical heavy-ion collisions has around twice as many neutrons as protons, which has implications, for example, for the imbalance between the generated charged pions [1]. Cold neutron stars are characterized by even lower proton fractions [2]. A high isospin asymmetry is carried by charged pion degrees of freedom that might be relevant for compact stars or for nuclear matter in general [3]. Although the early Universe is typically assumed to undergo an evolution with almost vanishing densities, a large lepton asymmetry might propagate into the baryon sector and shift equilibrium conditions around the time of the QCD transition [4].
In the grand canonical approach to QCD, the densities are traded for the corresponding chemical potentials and the isospin chemical potential is defined as = ( − )⇑2. The phase diagram in the temperature -isospin chemical potential plane is known to exhibit at least four different phases that can be relevant for the physical systems mentioned above. We discuss these phases qualitatively in Fig. 1. In the vacuum (for low values of and of ) QCD is confining and the effective degrees of freedom are hadrons. On the one hand, if the temperature is raised at low isospin chemical potential, deconfinement sets in and the quark-gluon plasma phase is realized. Lattice QCD simulations have demonstrated that the transition to deconfinement is an analytic crossover [5] and occurs near the chiral pseudocritical temperature of ≈ 155 MeV [6]. On the other hand, if is increased at low temperature, another phase transition line is encountered, which can be best understood using an effective description of QCD based on pionic degrees of freedom. Charged pions are the lightest hadrons that couple to the isospin chemical potential, and their effective dynamics is described by chiral perturbation theory ( PT) [7]. At = 0, if the isospin chemical potential exceeds the critical value , = ⇑2, sufficient energy is pumped into the system so that charged pions can be created. Due to the bosonic nature of pions, a Bose-Einstein condensate (BEC) is formed at this point. PT also predicts that the transition between the vacuum and the BEC state is of second order with the universality class O(2) [7]. As is usually the case, the temperature tends to destroy the condensate so that this phase transition line is expected to bend towards higher values of as grows. While PT predicts a monotonous rise of , ( ), our preliminary lattice simulations reveal a flattening of this curve around the zero-density deconfinement crossover [8,9]. This suggests that deconfinement is in some sense stronger than the condensation mechanism, and in the quark-gluon plasma phase no pions can be created by increasing .
Yet another transition is expected to occur at even higher isospin chemical potential. Perturbation theory predicts that the attractive gluon interaction forms Cooper-pairs (BCS superconductivity) of and¯quarks in the pseudoscalar channel. Since the resulting pair has the same quantum numbers as the pion condensate, the deconfinement-type transition between the BEC and BCS states is expected to be an analytic crossover [7].
The particular structure of the transition lines separating these four phases (hadronic, quark-gluon plasma, BEC and BCS) in the phase diagram is not yet known. The results we present below favor the scenario depicted in Fig. 1, where the = 0 deconfinement transition continuously connects to the BEC-BCS crossover at high isospin chemical potentials. We note that for asymptotically large values of , perturbative arguments suggest [7,10] a decoupling of the gluonic sector and the emergence of a first-order deconfinement phase transition. This region is not included in Fig. 1.
The QCD phase diagram at nonzero temperature and isospin chemical potential has been studied using a multitude of approaches. Besides the already mentioned chiral perturbation theory as effective description [7,11], hard thermal loop perturbation theory [12] and different low-energy models of QCD have been employed. The latter include the Nambu-Jona-Lasinio (NJL) model in four [13] and in two dimensions [14], the linear sigma or quark meson model [15], scalar field theory [16], random matrix models [17], the thermodynamic bag model [18], the hadron resonance gas model [19] and by means of the functional renormalization group or Dyson-Schwinger equations [20]. The impact of the isospin chemical potential was also investigated in the limit of large number of colors [21] and employing the holographic principle [22]. For large , the phase diagram was argued to be related to the analogous phase diagram with baryon chemical potential via the orbifold equivalence [23]. In particular, this equivalence is expected to hold outside the pion condensation phase only, so that the precise knowledge of the location of the BEC transition line is relevant from this point of view as well.
A significant advantage of the > 0 system is that -unlike for nonzero baryon chemical potential -there is no sign problem and lattice QCD simulations can be employed to investigate the phase diagram directly. Following the pioneering work of Refs. [24,25], the phase diagram was studied both in the grand canonical [26][27][28] and in the canonical ensemble [29,30]. The low isospin density region was also discussed by means of a Taylor expansion in [31,32], while the strong coupling regime was investigated using the strong coupling expansion [33]. Finally, we mention that QCD with > 0 has deep relations to two-color QCD with baryon chemical potentials. The latter theory has also been the subject of intensive research, both on the lattice [34] as well as in model and functional approaches [35].
We mention that, besides the broad range of applications for this system, an additional motivation to study the > 0 theory is its similarity with nonvanishing baryon chemical potential, > 0, both on the conceptual and on the technical level. At = 0, both systems exhibit the so-called silver blaze phenomenon [36] and undergo a phase transition from the vacuum to a phase where colorless composite particles are created. This transition is accompanied by the accumulation of nearzero eigenvalues of the fermion matrix in both cases, leading to an ill-conditioned inversion problem and a breakdown of naive simulation algorithms. This necessitates the use of an infrared regulator that we denote by below. Understanding these concepts and facing these technical challenges in the (sign-problem-free) ≠ 0 theory may give us insight on how to assess the ≠ 0 system in the future.
In this paper we follow the grand canonical approach and simulate QCD with 2 + 1 flavors of dynamical staggered quarks at various values of and . Besides the inclusion of the strange quark, we improve over currently existing studies in the literature by using physical quark masses and employing an improved lattice action in the simulations. The latter enables a faster scaling towards the continuum limit. Furthermore, we present a novel method to perform the extrapolation of the infrared regulator to zero -the main technical achievement in this project. This method involves the singular values of the massive Dirac operator, which we discuss for the first time on the lattice in this context. Preliminary results of our simulations were presented in Refs. [8,9] and applied in Ref. [37].

II.1. Symmetries
First of all, it is instructive to discuss the pattern of the symmetry breaking that leads to Bose-Einstein condensation at high isospin chemical potentials. In Euclidean spacetime, our choice for the continuum action =¯ℳ for the light quarks = ( , ) ⊺ includes Here is the gluon field and denote the Pauli matrices. Besides the isospin chemical potential and the quark mass , we also included an explicit symmetry breaking term in ℳ that couples to the charged pion field ± , (2) The role of the parameter -referred to as a pionic source in the following -will be elucidated below. At = = 0 (but nonzero ), the action is symmetric under the chiral group SU (2) × U (1). The U (1) symmetry, corresponding to baryon number conservation, is not affected by the chemical potential nor by the pionic source and is not discussed in the following. The SU (2) symmetry group is broken down to U 3 (1) by the chemical potential. The subscript indicates that the generator of the remaining symmetry is 3 . Pion condensation is signaled by the spontaneous breaking of this U 3 (1) symmetry, by any of the expectation values ∐︁¯5 1̃︁ , ∐︁¯5 2̃︁ .
The spontaneous breaking of the continuous U 3 (1) symmetry implies the presence of a Goldstone mode. The introduction of the pionic source in (1) corresponds to an additional explicit breaking that selects the 2 direction for the ground state and makes the would-be massless mode a pseudo-Goldstone boson. In fact, in a finite volume, such a trigger is necessary for the spontaneous breaking to occur. The physical limit → 0 is obtained subsequently by means of an extrapolation. Notice that for the U 3 (1) symmetry, its spontaneous (by ∐︀ ±̃︀ ) and explicit (by ) breaking is completely analogous to the one of the standard chiral symmetry at = 0, together with its spontaneous (by ∐︁¯̃︁) and explicit (by ) breakings. One very important difference is that while in nature > 0, the parameter is unphysical and the limit → 0 needs to be taken.
For vanishing isospin chemical potential, the pionic source can be rotated into the mass parameter. In particular, we can perform the chiral rotation (3) This rotation brings the action into a flavor-diagonal form Thus, at = 0 the mass and the pionic source parameter are indistinguishable and only their squared sum plays a physical role. This will have important consequences for the ultraviolet structure of the theory, see Sec. II.3 below. Notice that for = 0 the pionic source has the same interpretation as the twisted mass parameter used in the context of Wilson fermions [38].

II.2. Lattice setup
To investigate pion condensation and the associated spontaneous symmetry breaking, we study 2+1-flavor QCD with > 0 and > 0 nonperturbatively on the lattice. We work on 3 × lattices with spacing and sites indexed by the coordinates ( , , , ). The temperature and the spatial volume read = 1⇑( ) and = ( ) 3 , respectively. To discretize the Dirac operator we use the staggered formulation; here the equivalent of 5 is the local spin-flavor structure 5 = S 5 ⊗ F 5 = (−1) + + + . At = = 0, this operator couples to the Goldstone boson in the chiral limit.
The partition function of this system is given in terms of the path integral over the gluon links = exp( ), where we employed the rooting procedure (for a discussion on the theoretical issues related to rooting, see Ref. [39]). In Eq. (5), = 6⇑ 2 denotes the inverse gauge coupling, is the tree-level Symanzik improved gluon action, ℳ is the light quark matrix in the basis of the up and down quarks and ℳ is the strange quark matrix, Here, the argument of ⇑ indicates the chemical potential , introduced on the lattice by multiplying the forward (backward) timelike links by ( − ). The pionic source parameter induces an explicit symmetry breaking and serves to trigger pion condensation as we discussed above in Sec. II.1.
The integrand of needs to be positive to allow the importance sampling of the gluon configurations. The strange quark determinant is real and positive due to the standard 5 -hermiticity relation 5 ℳ 5 = ℳ † . To show the positivity for the light sector, we need to discuss the symmetry properties of the staggered Dirac operator at ≠ 0. The staggered equivalent of chiral symmetry implies that holds. In addition, the Dirac operator satisfies so that the light fermion matrix is 1 5 -hermitian Thus, taking the determinant of both sides shows that det ℳ is real.
One can also show that the determinant is positive by considering where = diag(1, 5 ) and we used Eq. (8). Indeed, since has unit determinant, we have (11) Thus, both determinants in the measure of the path integral (5) are positive.
The fourth root of the determinants in (5) is approximated via rational functions. The simulation setup with > 0 was first introduced for the = 2 theory in the pioneering work of Ref. [24] for the quenched case and in Ref. [25] for dynamical QCD. The same technique was also used in Ref. [28]. Here we extend this setup by including the strange quark as well. In addition, we improve the lattice action by using the tree-level Symanzik gauge action and by employing two steps of stout smearing in the Dirac operator. The quark masses are tuned to their physical values along the line of constant physics (LCP) ( ), as determined in Ref. [40], with the pion mass ≈ 135 MeV. Our simulation algorithm is based on Ref. [41]. In addition we implement a Hasenbusch-type improvement scheme that is typically used in the context of mass preconditioning [42]. In our setup this amounts to the replacement det ℳ ( ) = det ℳ ( ) ⋅ det ℳ ( )⇑ det ℳ ( ) with > , which allows to use a larger step size (and a lower precision) in the simulation algorithm for the second (more expensive) factor.

II.3. Observables and renormalization
Our primary observables are the pion condensate and the quark condensate. Both are obtained from the partition function via differentiation, Inserting the partition function (5) and rewriting the light quark determinant using Eq. (11), we obtain for the condensate operators, The relation (13) allows for direct measurements using noisy estimators. The observables of Eq. (12) are subject to additive renormalization. This is necessary, since log contains ultraviolet divergences in the inverse lattice spacing. The structure of these divergences can be determined based on dimensional arguments (see Ref. [43] for a discussion at where we suppressed further divergences that contain . Note that the divergences are independent of , since the chemical potential couples to a conserved charge [44]. Thus, it suffices to consider the case of = 0. Above in Eq. (4) we have seen that for vanishing isospin chemical potential, the mass and the pionic source may be rotated into each other, so that the two parameters can only appear in Eq. (14) in the form 2 + 2 . The quark condensate and the pion condensate inherit the quadratic and the logarithmic divergences from log . However, for ∐︀ ±̃︀ these vanish at = 0 -the point of interest for the physical theory. Since the light quark mass is nonzero, the additive divergences remain in ∐︁¯̃︁ and need to be subtracted even in the limit → 0. The standard choice is to consider the difference to ∐︁¯̃︁ measured in the vacuum, i.e. at = = 0. Finally we need to address the multiplicative renormalization of our observables. The quark condensate and the pion condensate have nontrivial renormalization constants,¯= −1 and = −1 . The equivalence of and at = 0 and the -independence of the renormalization constants, however, imply that = , so that a multiplication of the condensates by cancels the multiplicative divergence. 1 Altogether, the renormalized observables read where we also included a normalization factor involving the pion mass = 135 MeV and the chiral limit of the pion decay constant = 86 MeV and added unity to the quark condensate for convenience. In this normalization (which follows Ref. [45]), Σ¯= 1 at = = 0 due to the Gell-Mann-Oakes-Renner relation. In addition, zero-temperature leading-order PT [7] predicts a gradual rotation of the condensates so that Σ 2 + Σ 2 = 1 holds irrespective of , which can also be observed to some extent in the full theory.
In addition to the fermionic observables we also consider the Polyakov loop, as a measure for deconfinement. The multiplicative renormalization of amounts to with an arbitrary choice for ⋆ = ( ⋆ , 0) and ⋆ . Different choices correspond to different renormalization schemes; we choose ⋆ = 162 MeV, where the bare Polyakov loops are already significantly nonzero, and set ⋆ = 1. This renormalization prescription for was developed at = 0 [46] and also put into practice at nonzero background magnetic fields [47].
The pionic source parameter > 0 is introduced in order to trigger pion condensation and to stabilize the simulation algorithm. However, in order to obtain physical results needs to be extrapolated to zero. Most of the observables exhibit a pronounced dependence on , making this extrapolation cumbersome. As a typical example, in the top panel of Fig. 2 we show the pion condensate Σ of Eq. (15) as measured using various values of on our 24 3 × 6 ensembles at = 113 MeV. Especially below and around the critical chemical potential , ≈ ⇑2, a direct fitting of the data appears to be hopeless. This necessitates various improvements, both in the valence sector (the definition of the operators) and in the sea sector (reweighting of the configurations), to eliminate the -dependence in the observables. The result of the improvement, which we will describe below, is also included in the top panel of Fig. 2 -clearly revealing the transition between the vacuum and the BEC phase, which was hidden for the direct data at high .
The simulations at low values of suffer from a numerical problem. Since the pionic source acts as an infrared regulator, it has a substantial impact on the condition number of the fermion matrix ℳ . In particular, the average iteration count cg for the convergence of the conjugate gradient algorithm used for inverting ℳ grows significantly as decreases. This is especially the case in the pion condensed phase (see the bottom panel of Fig. 2), where the inversion at our smallest is an order of magnitude slower than at the largest pionic source. 2 Moreover, in order to maintain a reasonable acceptance in the Metropolis step, the simulations at low need to be performed with a reduced step size along the Monte-Carlo trajectory due to increasing fluctuations in the fermion force. This leads to a slowing down by a further factor of 3 − 4 for the lowest pionic source.
Thus, improving our observables so that they lie closer to their → 0 limit is necessary -both for controlling the extrapolation and for sparing simulation time. It turns out that the improvement for the pion condensate is drastically different from the improvement of the other 2 We also mention that simulating with = 0 is not feasible even for our smallest isospin chemical potentials due to the occasional appearance of configurations with fermion matrices of very high condition numbers, for which the numerical inversion breaks down. observables due to the fact that ∐︀ ±̃︀ plays the role of the order parameter for Bose-Einstein condensation and thus has a highly nontrivial finite volume scaling. There is, nevertheless, a common element in the improvement program, which is the singular value representation of the massive Dirac operator.

III.1. Singular values
To get acquainted with the notion of singular values, let us begin with the eigenvalue equation of the massive Dirac operator. For the up quark, the eigensystem reads where the eigenvalues are complex numbers. Using chiral symmetry (7) and the hermiticity relation (8)  Note that ⇑ ( ) is not a normal operator thus (︀ ⇑ ( ), ⇑ † ( )⌋︀ ≠ 0 and its left and right eigenvectors do not coincide. Nevertheless, Eqs. (18) and (19) reveal that for each eigenvalue in the up quark sector there is a complex conjugate pair in the down quark sector, which is why the determinant of the total light quark matrix is real and positive as we have seen in Sec. II.2. We can create a hermitian operator by taking the modulus squared of the Dirac operator [48]. The square roots of the eigenvalues of this operator are referred to as the singular values 3 of the Dirac operator. The eigenvalue equation reads Note that due to the non-normality of ⇑ ( ), there is no relation between the singular values and the eigenvalues . We note that the singular values of ⇑ ( ) + and of ⇑ (− ) + coincide due to the hermiticity relation (8).

III.2. Banks-Casher type relation for the pion condensate
An important characteristic of chiral symmetry breaking at = 0 is reflected by the Banks-Casher relation [49], which states that the chiral limit of the quark condensate is related to the density of the Dirac eigenvalues around zero. Below we demonstrate that the → 0 limit of the pion condensate can be similarly obtained by looking at the density of the singular values around zero. The derivation follows Ref. [48], where massless quarks at > 0 were considered. Here we generalize the discussion to arbitrary . Throughout this section we neglect exact zero modes, whose contribution is discussed in detail in Ref. [48].
Using Eq. (20), we can write down the singular value representation of the pion condensate. Writing the trace of Eq. (13) in the basis of the modes we obtain, Here in the first step we considered the volume to be large enough so that the singular values become suffi- 3 The name refers to the role of these in the singular value decomposition of non-hermitian matrices.
ciently dense and the sum can be replaced by an integral introducing the density of the singular values, In the second step of Eq. (21) we performed the → 0 limit, which lead to a representation of the -function and resulted in the density (0) around zero. Eq. (21) tells us that having a nonzero pion condensate is equivalent to an accumulation of the near-zero singular values of the massive Dirac operator. Notice that the ordering of the two limits in Eq. (21) -just like in the usual Banks-Casher relation -is essential, and setting = 0 in the operator directly would give ± = 0. Equation (21) also connects the slowing down of the operator inversion and the high condition number of the fermion matrix in the BEC phase (as we demonstrated in the bottom panel of Fig. 2) with a physical notion: the emergence of a nonzero pion condensate.

III.3. Improvement of the quark condensate
For the quark condensate operator of Eq. (13), the spectral representations reads In this case the → 0 limit can be taken explicitly, without approaching the thermodynamic limit first. However, unlike the pion condensate, the quark condensate necessitates the calculation of various matrix elements and not only that of the singular values. Besides the values of the operators at = 0, it will be advantageous to have access to the dependence of the operator as well. In particular, we define the change in the operator between and = 0, We will use¯to improve the dependence of theō perator calculated using noisy estimators.

III.4. Leading-order reweighting
The above improvements eliminated the explicit dependence of the operators (i.e., of the valence quarks) on the pionic source. The remaining -dependence originates from sea quarks, i.e. from the nonzero value of in the path integral measure used to define expectation values ∐︀̃︀ . In particular, this involves the -dependence of the light quark determinant in Eq. (5). To get rid of this contribution, we need to manipulate the distribution of configurations by introducing the reweighting factors, where we used Eq. (11). This way the > 0 determinant is canceled in the expectation values so that we mimic the distribution that would have been obtained via a simulation directly at = 0. Note that Eq. (25) only holds if the distributions at > 0 and at = 0 have sufficiently large overlap.
The full determinant is a very expensive object as in principle it necessitates the calculation of all singular values. However, since we only need ( ) for small pionic sources, we can expand the reweighting factor in . Rewriting the logarithm of the reweighting factor in this manner, we obtain More specifically, here we replaced the logarithm of the determinant ratio with its Taylor-expansion in the pionic source (the expansion variable is denoted ) and evaluated the result at = . We then exploited the fact that odd terms in the expansion vanish. On comparison with Eq. (13), the leading order term in the expansion is found to be proportional to the pion condensate. The resulting leading-order reweighting factor is denoted by LO ( ). Thus we conclude that the reweighting of an observable, to leading order in , involves the exponential of the pion condensate (measured at > 0) times the fourvolume. Since the pion condensate is anyway measured on the individual configurations to compute ∐︀ ±̃︀ , this improvement comes with no extra costs. We test this improvement on small lattices, where a calculation of the complete spectrum of singular values is feasible, enabling a direct comparison between and LO . Specifically, we employ the low-temperature 8 4 ensembles generated in Ref. [28] and plot the two reweighting factors against each other in Fig. 3 for two values of around the critical isospin chemical potential. The scatter plot clearly shows the strong correlation between and LO and that the two factors become identical in the limit → 0.

IV. RESULTS: IMPROVEMENT PROGRAM
In the following we always perform the leading order reweighting with LO = exp(︀− ± ⇑(2 )⌋︀, where ± is 4 Notice furthermore that LO tends to overestimate if the reweighting factors are below average (and vice versa). This can be understood using the singular value representation of Eqs. (25) and (26), which can be used to show that the difference LO − is a positive and monotonously decreasing function of at fixed . Thus, a fluctuation that reduces the singular values (i.e. that pushes below average), inevitably increases the deviation LO − , thereby explaining the tendency visible in Fig. 3. calculated from (13) employing noisy estimators to evaluate the traces. The measurements are carried out on the reweighted ensembles.

IV.1. Improved pion condensate
We describe our strategy to determine (0) in more detail. The singular values of the massive Dirac operator (i.e., the eigenvalues of ⋃︀ ⇑ ( ) + ⋃︀ 2 , cf. Eq. (20)) are calculated using the Krylov-Schur algorithm. We find it sufficient to work with the lowest 50-150 singular values for each configuration. Using this set of singular values we build a histogram for the integrated spectral density where we discard bins that lie below the average lowest singular value. The statistical error of ( ) in each bin is estimated via the jackknife procedure. In Fig. 4 we plot ( )⇑ for two ensembles just around and slightly beyond the transition to the pion condensed phase.
To obtain (0), we need to extrapolate ( )⇑ down to zero. This is performed via polynomial fits of the histogram. To take into account the correlation between the individual bins, we minimize the correlated 2 corr , which involves the inverse of the correlation matrix of the data. To avoid numerical problems during this inversion, we smear the lowest eigenvalues of following the strategy of Ref. [50]. We perform the fits with polynomials of various degrees and over various fit ranges. The resulting values for (0) are weighted by exp(− 2 corr ⇑ dof ) with dof the number of degrees of freedom in the fit. The weighted results are used to build yet another histogram (red area in Fig. 4), similar to the analysis conducted in Ref. [51]. The central value for (0) is obtained by the median, while the statistical plus systematical error of the fit is estimated by the range that contains the middle 68% of the histogram (orange area in Fig. 4). The analysis gives drastically different results for the two chemical potentials considered in the figure. While the intersect of the fits is large for ⇑ = 0.53, it is consistent with zero for ⇑ = 0.48. (In the latter case (0) and its error is taken as the positive part of the orange histogram.) This demonstrates -via the relation (21) -the emergence of a nonzero pion condensate in the BEC phase. The so obtained results for ∐︀ ±̃︀ are included in the top panel of Fig. 2 above.

IV.2. Improved quark condensate
In principle, for the improvement scheme of the quark condensate, outlined in Sec. III.3, we need all singular values . This is in contrast to the improvement of the pion condensate, where only the spectral density of singular values around zero is of importance. Nevertheless, computing all singular values for each configuration is unfeasible in practice and, in fact, not necessary. The important simplification follows from the observation that the difference¯from Eq. (24) is dominated by the dependence of the terms involving the lowest singular values. It might thus be sufficient to perform the explicit → 0 extrapolation in the operator for only the first low modes and to leave the contribution from the remaining singular values untouched. To this end we introduce the truncated operator and the associated truncated differencē to rewrite the expectation value as The second term vanishes when we perform theextrapolation, so that we obtain This type of improvement will be efficient if is large enough to ensure that¯≈¯so that the subtraction effectively removes most of the dependence of the observable.
The optimal value for to achieve a good balance between the improvement for the dependence of the observable and the computational cost will, in general, depend on the operator, the lattice size and the lattice spacing, as well as on the magnitude of the -values in use. The results for¯versus for different values of are shown in the top panel of Fig. 5. The figure indicates that moderate values of suffice to ensure that most of the dependence in the valence sector is absorbed by the improvement. In the bottom panel of Fig. 5 we plot the combination̂︀¯−¯̃︂, revealing no significant dependence for ≳ 100. For comparison, the dependence of the unimproved observable (marked by the short horizontal lines in the figure) is much more pronounced.
We find that ∼ 100 does suffice to warrant well controlled -extrapolations for all ensembles and observables considered so far. Thus, the same set of singular values can be used for the improvement of the pion condensate as well as for the quark condensate.

IV.3. Final -extrapolations
Having performed the improvements described above, we can carry out the final extrapolation of the results to = 0. Since the improvements do not achieve a perfect elimination of the pionic source [the reweighting is done only at leading order, the singular value sums are truncated at a finite , and the finiteness of the volume distorts the relation (21)], a mild dependence of the results on is still expected. In Fig. 6 we show a few representative examples for the extrapolations of Σ and of Σ¯, where we also include the unimproved observables [i.e. the quantities obtained via the full traces, Eq. (13)]. Contrary to the pronounced (and thus uncontrolled) dependence of the latter on , we find that the improved observables can be reliably fitted by either a constant or a linear function (on general grounds, Σ is an odd function of , whereas Σ¯is even in the pionic source). We also attempted to perform a fit of the unimproved data by means of a spline with Monte-Carlo-generated nodepoints (for the details of this fit procedure, see Ref. [8]). For comparison, the so obtained extrapolations are also included in Fig. 6. Unlike the improved extrapolations, these fits are always dominated by the points at low , making them unreliable in several cases.

V.1. Phase diagram
In the following we work with the → 0 extrapolated observables to determine the phase diagram of the system in the − plane. We work with four lattice ensembles  Fig. 7 shows our = 6 results for Σ and for Σ¯from Eq. (15). The onset of pion condensation is characterized by the abrupt change at the boundary between the vacuum and the BEC phase, mapping a critical line , ( ). The quark condensate exhibits a similarly pronounced change at the pion condensation phase boundary. In addition, it also features a smooth dependence on the temperature around the chiral crossover transition ( ). In the following, we refer to the point, where these two transition lines meet as the pseudo-triple point at = , and = . The notation refers to the fact that the chiral transition is pseudocritical and does not correspond to a true phase transition (in which case this point would be a true triple point).
The boundary of the pion condensation phase, , ( ), is determined by the points where Σ acquires a nonzero expectation value. We estimate its location by interpo- The pion condensate (top panel) and the quark condensate (bottom panel) as functions of the temperature and the isospin chemical potential as measured on our = 6 lattices. The color coding reflects the magnitude of the observables: blue (Σ = 0) towards yellow (Σ = 1). Notice that the orientation of the -axis is different in the two panels. lating 5 the results for Σ using a two-dimensional cubic spline fit with Monte-Carlo-generated nodepoints [8]. From PT it is known that , (0) = ⇑2 [7]. We observe that , ( ) is independent of up to a temperature of about 130 to 140 MeV -see also the results shown in [8,9]. For a temperature of about 155 to 160 MeV the phase boundary starts to flatten out, and above ≈ 161 MeV the pion condensate is found to vanish within errors for all ensembles. Here we follow the phase boundary up to = 120 MeV, but this flat behavior is found to persist at least up to ≈ 190 MeV on our = 6 ensemble, see below. These findings contradict 5 A technical issue regarding this fit is the non-analyticity of Σ at the critical chemical potential. Although it is regulated by the finite volume, the sharp behavior around , cannot be captured by the smooth splines. To avoid this issue, here we fit the results from Sec. IV.1, allowing for negative intersects for (0) (i.e. the full orange histograms of Fig. 4). This enables a more stable determination of the Σ = 0 contour. the expectations from PT, which predicts a monotonous rise in , ( ) [7].
To perform the continuum extrapolation, we parameterize , ( , ) by a polynomial in ( − 0 ) with lattice spacing dependent coefficients. We set 0 = 140 MeV, below which , ( , ) is found to be independent of . To capture both the flatness of the critical chemical potential for low temperatures and its abrupt rise around ≈ 155 MeV, we found it necessary to include terms proportional to ( − 0 ) 2 , ( − 0 ) 3 and ( − 0 ) 4 in the fit and to fix , ( 0 , 0) = ⇑2. We find that for 140 MeV ≲ ≲ 155 MeV, our = 6 results are outside of the scaling region for lattice artefacts of ( 2 ) so that we excluded these from the continuum extrapolation. The final continuum curve, extrapolated including terms proportional to 2 , is shown in Fig. 8 (top panel), together with the data for the phase boundary which has been included in the fit. ≈ 161 MeV serves as an upper bound for the critical temperature, as Σ is found to vanish within the uncertainties for all lattice spacings and for all values of considered in this region, as mentioned above.
We proceed with the chiral crossover transition temperature ( ), which we define as the location of the inflection point of Σ¯with respect to . 6 The data shows an initial reduction in ( ) up to about = 70 MeV, followed by a slight increase. Beyond the pseudo-triple point, where the crossover line and the pion condensation boundary meet, we observe the two transition lines to coincide within errors. To be more quantitative, we once again determine ( ) using a twodimensional spline fit where the nodepoints have been generated via Monte-Carlo. In the spline fit 7 we include the constraint that Σ¯is an even function of . To perform the continuum extrapolation, we parameterize ( , ) by a polynomial which is even in , including data up to , (0) = ⇑2. The behavior of ( ) beyond this value is difficult to capture with polynomials. For the range of chemical potentials included in the fit, the data is already well described by a quadratic polynomial in and ( 2 ) lattice artefacts, see Fig. 8 (bottom panel). Adding a 4 -term or higher-order lattice discretization errors was found to only increase the uncertainties and not lead to significant differences.
The precise location of the pseudo-triple point is also of interest. We define , ( ) to be the value for 6 Note that in Refs. [8,9] we employed a different procedure, defining as the point where Σ¯( , ( )) = Σ¯(0, (0)) holds. This definition is only valid for < , (0), since beyond that the condensate is affected by already at = 0. On the contrary, the inflection point can be used to define ( ) at any value of . 7 Here we only take into account the fit for > 135 MeV, where the spline captures the behavior of Σ¯( , ) sufficiently accurately. Below this temperature the condensate changes abruptly at the BEC phase transition, which is not captured by our smooth interpolation. where the two transition lines become consistent within errors -i.e. where the uncertainties overlap. This provides a conservative lower bound for the pseudo-triple point. The results for , ( ) are shown in Fig. 9. As for the pion condensation phase boundary, the continuum extrapolation is performed including lattice artefacts of ( 2 ) and excluding the = 6 results. The resulting continuum extrapolation is indicated by the shaded orange band in Fig. 9. We found this extrapolation to be more stable than the fit to all available lattice spacings, including lattice artefacts of ( 4 ). Both extrapolations lead to similar results, see the comparison in Fig. 9.
We are now in the position to draw the continuum phase diagram, which we display in Fig. 10. The chiral crossover transition starts from a temperature of (0) = 159(4) MeV, which is consistent with the crossover temperatures from [6] within uncertainties. The results ex- Continuum extrapolations for the isospin chemical potential , corresponding to the pseudo-triple point, where the chiral crossover line meets the pion condensation boundary. The orange curve corresponds to the continuum extrapolation for = 8, 10 and 12 including lattice artefacts of ( 2 ) and the gray curve is the continuum extrapolation for all points including an additional ( 4 ) term. ≳ 160 MeV, we do not observe pion condensation up to = 120 MeV. The two transition lines meet at the pseudo-triple point, for which we obtain , = 70(5) MeV in the continuum limit, indicated by the yellow point in Fig. 10. The corresponding temperature is determined conservatively by taking into account the upper bound for  (7) MeV. From what we observe at finite lattice spacings, we expect that chiral symmetry restoration and the pion condensation phase boundary coincide from the pseudo-triple point on. To demonstrate this, in Fig. 11 we plot the pion condensate together with the quark condensate for > , . The figure indicates that pion condensation (defined by the point where Σ = 0) occurs together with chiral symmetry restoration (the inflection point of the condensate). The initial rise of the chiral condensate in Fig. 11 is an interesting feature in the pion condensation phase. A similar tendency has been observed in a study of the phase diagram of a related two-color NJL model [52]. We interpret it as a remnant of the relation between pion and chiral condensate Σ 2 + Σ 2 = 1, discussed in Sec. II.3, which follows from PT to leading order [7].

V.2. Polyakov loop
Next, we elaborate on the properties of the deconfinement transition in terms of the renormalized Polyakov loop . In contrast to the quark condensate, the Polyakov loop exhibits no pronounced inflection point. To capture how deconfinement depends on the isospin chemical potential, we consider the curves in the − plane, where = const. is satisfied. Considering our definition (17) for the renormalization, the contour with = 1 is a possible choice for the transition temperature. In addition, the distance between the various contour lines is related to the slope of around the transition point. We show the contours of constant in the phase diagram for = 6 in Fig. 12. The contour lines are observed to be insensitive to the BEC phase transition, as they enter the pion condensation phase without any sign of critical behavior. In this plot we also included results obtained at high isospin chemical potentials, up to ≈ 200 MeV for the pion condensation boundary and up to ≈ 240 MeV for a few Polyakov loop contours. After crossing the phase boundary, the contours continue to decrease, in agreement with the sketch of the phase diagram in Fig. 1. We postpone a more detailed discussion of the deconfinement transition to a future publication with an extended range in .

V.3. Order of the BEC transition
Finally, we investigate the nature of the transition between the vacuum and the pion condensed phase via a finite size scaling analysis. To this end three = 6 ensembles with = 16, 24 and 32 are employed. The results for Σ around the critical isospin chemical potential , are shown in Fig. 13. Both at = 113 MeV and at = 136 MeV a sharpening of Σ ( ) is visible, confirming that a real phase transition takes place in the → ∞ limit. To be more quantitative, we proceed with the assumption that the transition is of second order. In this case, our results at > 0 and ≈ , must reflect the critical behavior that emerges as → 0 and → , in the vicinity of the transition point. For the following analysis we use the data obtained on the 24 3 ×6 ensemble at the lowest temperature = 113 MeV.
First we compare our (unimproved) > 0 results for Σ to the prediction of PT. This involves two free parameters: the critical isospin chemical potential and a normalization factor , which corresponds to the magnitude of the renormalized quark condensate Σ¯in the vacuum, where the vacuum angle is determined implicitly by the second equation (see Refs. [28,53]). The data is very well described by the predicted behavior, see the top panel of Fig. 14. Taking into account data points for ⇑ < 0.63 and the two smallest pionic sources ⇑ < 0.2, we obtain a reasonable fit with 2 ⇑dof ≈ 1.6. For higher isospin chemical potentials PT tends to underestimate Σ , as has already been noted in Ref. [28]. The critical chemical potential is , = 71(2) MeV and the normalization factor = 1.02(1), both very close to the values expected in the zero-temperature and thermodynamic limits ( ⇑2 and unity, respectively).
We also consider the critical behavior of the O(2) uni-versality class to fit the order parameter as a function of and . In particular, this dependence involves the universal scaling function , with the critical exponents and and the free parameters , , 0 , and 0 . For the construction of , we follow Ref. [54]. The result of the fit, using the two lowest pionic sources and isospin chemical potentials 0.4 < ⇑ < 0.63, is shown in the bottom panel of Fig. 14, demonstrating the collapse of the data on the universal curve. The critical chemical potential is found to be , = 74(2) MeV. However, since the fit quality is rather low ( 2 ⇑dof ≈ 3), we extended the fit function to include scaling violations in the spirit of Ref. [54]. We found it necessary to work with which allowed us to achieve reasonable fit qualities ( 2 ⇑dof ≈ 1.5) for broader fit intervals ⇑ < 0.5 and 0.2 < ⇑ < 0.7. The resulting fit is included in the top panel of Fig. 14. We find , = 70(2) MeV, again lying very close to the expected value of ⇑2. Summarizing, our results strongly indicate that the transition into the pion condensed phase is of second order in the infinite volume limit. In addition, the behavior of the order parameter at nonzero values of the explicit breaking parameter is consistent with the critical scaling in the O(2) universality class, as expected due to the spontaneous symmetry breaking pattern discussed in Sec. II.1.

VI. SUMMARY
In this paper we have presented the first study of the phase diagram of QCD at nonzero isospin chemical potentials in the continuum limit with dynamical , and quarks at physical quark masses. Our main result is the continuum extrapolated phase diagram, as shown in Fig. 10. The key features of the phase diagram are the chiral crossover transition, the pion condensation boundary and their meeting point, the pseudo-triple point. We observe that the critical chemical potential for pion condensation remains constant up to ≈ 140 MeV, followed by a drastic change in behavior and a saturation at around ≈ 160 MeV. Above this temperature, we no longer observe pion condensation, at least up to = 120 MeV. This finding might be particularly relevant for the orbifold equivalence in the large limit, as the pion condensation region is smaller than previously expected.
The chiral crossover temperature decreases slightly until it reaches the vicinity of the pion condensation phase boundary. The two transition lines meet at , = 70(5) MeV and are found to be on top of each other within uncertainties for higher isospin chemical potentials. Our scaling analysis of the pion condensation boundary shows consistency with a second order phase transition in the O(2) universality class.
Using our = 6 lattice ensembles, we find first indications that the deconfinement transition temperature -defined in terms of the renormalized Polyakov loopdecreases and smoothly penetrates into the pion condensation phase. This behavior is depicted schematically in Fig. 1. If this tendency persists in the continuum limit and for larger values of , the scenario where the deconfinement transition connects continuously to the BEC-BCS crossover will be favored. Our results for the phase diagram can be compared directly to effective theories and models of QCD and serve as a highly nontrivial check of the validity of these approaches.
The main technical novelty of this study is the development of an improvement program for the extrapolations in the infrared regulator , discussed in detail in Secs. III and IV. Without the use of these improved extrapola-tions the reliable extraction of the phase diagram would have hardly been possible. We would like to emphasize that a similar infrared regulator will likely be necessary to enable simulations at finite baryon chemical potential -granted that the sign problem has been solved. In this case, a generalization of the methods presented here might be helpful to give insight to the investigations at finite .