Equation of state and speed of sound of isospin-asymmetric QCD on the lattice

: We determine the QCD equation of state at nonzero temperature in the presence of an isospin asymmetry between the light quark chemical potentials on the lattice. Our simulations employ N f = 2 + 1 ﬂavors of dynamical staggered quarks at physical masses, using three diﬀerent lattice spacings. The main results are based on a two-dimensional spline interpolation of the isospin density, from which all relevant quantities can be obtained analytically. In particular, we present results for the pressure, the interaction measure, the energy and entropy densities, as well as the speed of sound. Remarkably, the latter is found to exceed its ideal gas limit deep in the pion condensed phase, the ﬁrst account of the violation of this limit in ﬁrst principles QCD. Finally, we also compute the phase diagram in the temperature – isospin density plane for the ﬁrst time. The data for all observables will be useful for the benchmarking of eﬀective theories and low-energy models of QCD and are provided in ancillary ﬁles for simple reuse.


Introduction
The theory of the strong interactions is Quantum Chromodynamics (QCD), featuring confinement of quarks and gluons at low energies, as well as asymptotic freedom at high scales.Albeit radically different in their properties, these two phases of strongly interacting matter are connected by a smooth crossover transition at zero net quark density according to lattice QCD simulations [1,2].How the dominant degrees of freedom transform from composite objects (hadrons) to colored quarks and gluons through this transition is described by the equation of state (EoS) of the system.In particular, the EoS gives a complete description of equilibrium QCD in terms of a relationship between thermodynamic observables including the pressure, the energy density or the entropy density.The phenomenological relevance of the EoS is manifold and ranges from heavy-ion physics to astrophysics and cosmology.The above observables control the evolution of the quark-gluon plasma in hydrodynamic models of heavy-ion collisions [3,4], the expansion of the early universe via the Friedmann equations [5] and also the mass and radius of stable neutron stars through the Tolman-Oppenheimer-Volkoff equation [6].In the latter context, a particularly relevant feature of the EoS is the speed of sound c s of QCD matter and the related polytropic index γ, which might serve as a proxy to distinguish stars with and without deconfined quark matter cores [7,8].
The above physical systems contain QCD matter in very different environments and therefore require the knowledge of the EoS as a function of different control parameters.These parameters include the temperature T , the chemical potentials µ conjugate to the conserved charges, as well as further variables like external electromagnetic fields.The relevant chemical potentials are charge µ Q , baryon µ B and strangeness µ S chemical potentials.While the baryon chemical potential is in most cases assumed to carry the dominant effect, in some cases µ Q can play the major role.This occurs for example for an early Universe featuring large lepton flavour asymmetries [9][10][11].Here the isentropic cosmological expansion leads to substantial charge chemical potentials, triggering the onset of pion condensation and producing characteristic signals for primordial black holes and gravitational wave spectra [12].Significant negative charge chemical potentials also arise for QCD matter in neutron stars due to an excess of down quarks over up quarks.
The charge chemical potential can be rewritten in terms of a nonzero isospin chemical potential µ I .While lattice simulations with a generic combination of quark chemical potentials suffer from the infamous complex action (or sign) problem, QCD at pure isospin chemical potential, i.e., at vanishing other chemical potential components, has a real action and is amenable to direct Monte-Carlo simulations [13].At low temperature, this setting exhibits a second-order phase transition to a phase with a Bose-Einstein condensate (BEC) of charged pions according to chiral perturbation theory [13].First lattice simulations with higher-than-physical quark masses have qualitatively confirmed this expectation [14][15][16][17] and also gave important insight to the structure of the phase diagram in the T -µ I plane [18,19] as well as pion dynamics at low T [20].In Refs.[21,22] we carried out a systematic investigation of this system with physical quark masses and determined the phase diagram in the continuum limit, revealing an interesting interplay of chiral symmetry breaking, deconfinement and Bose-Einstein condensation.The continuum phase diagram for the parameter space relevant for this study is shown in Fig. 1.
In this paper we determine the EoS throughout the phase diagram for a broad range of temperatures and isospin chemical potentials.Generalizing our approach at (approximately) zero temperature [23], we construct the pressure, the energy and entropy densities, the interaction measure and the speed of sound from the isospin density as primary observable.At low temperatures and high µ I , we find that the speed of sound increases above its conformal limit 1/ √ 3 (we use natural units, with the speed of light set to unity).This is the first evidence for the explicit violation in first principles QCD of this general bound expected from holography [24] and, together with the polytropic index, which we compute as well, might provide relevant information for the modeling of the EoS based on neutron star radii and masses [7,8].We also included a sketch of the region where c s exceeds this conformal bound in Fig. 1.Besides the EoS, we also use our results to draw the QCD phase diagram in the temperature -isospin density plane.This result, together with the com-  [21].The green band is the phase boundary to the phase with Bose-Einstein condensation (BEC), the blue band is the prolongation of the chiral crossover in the T -µ I plane and the red data point marks the pseudo-triple point, the meeting point of the crossover and the BEC phase.In yellow we also show a sketch of the region, where the speed of sound exceeds its conformal limit.
plete tabulated EoS as shown in the plots, is available in the ancillary files submitted to the arXiv along with the preprint of this paper.To facilitate the use of the EoS in phenomenological models, we provide an accompanying data publication [25], including the physical observables and the uncertainties for all temperatures and chemical potentials where results are available.The results can be used as benchmarks for low-energy models and effective theories of QCD as well as for comparing to functional approaches.First accounts of our findings have been given in Refs.[12,[26][27][28].

Simulation setup and main observables
In our simulations we use N f = 2 + 1 flavors of rooted staggered quarks with two steps of stout smearing at physical quark masses and the tree-level Symanzik improved gluon action.The line of constant physics for the bare quark masses is taken from Ref. [29].For the approach to the continuum limit we use lattices with temporal extents N t = 8, 10 and 12 and aspect ratios of N s /N t ≈ 3 (these ensembles have already been used for the phase diagram [21]), together with a number of additional ensembles at T = 0.More details concerning the run parameters are collected in appendix D. The isospin chemical potential µ I enters the light quark Dirac operator in an exponential form and is normalized such that pion condensation sets in at zero temperature at µ I = m π /2.As in our previous studies [12,[21][22][23] the simulations are performed including a pionic source parameter λ, in the light quark mass matrix (see also Refs.[14,15,17]), which serves as an infrared regulator and triggers pion condensation in a finite volume.Physical results are obtained by means of an extrapolation λ → 0, which is facilitated by improving the observables and reweighting the configurations.For details on this improvement, see Refs.[21][22][23].For computing uncertainties we use the bootstrap procedure with 1000 samples.As we will see below, the main observable is the isospin density, from which the full µ I -dependence of the EoS can be extracted.The strategy for the EoS computation will be outlined below.n I can be computed directly from the simulations as described in Ref. [23].To perform the λ-extrapolations we use the improvement program introduced in Ref. [21] with the application to n I as explained in Ref. [22,23].This improvement program results in fully controlled extrapolations and from now on we only discuss results which have already been extrapolated to λ = 0. We note that a well controlled λ-extrapolation is of particular importance to facilitate the following spline interpolations of n I in T and in µ I .

The EoS from an interpolation of the isospin density
Apart from the isospin density n I , the main task for the determination of the EoS is the computation of the pressure p and the interaction measure I.All other relevant quantities, apart from the speed of sound (to be discussed in Sec. 4) follow from these three quantities.
In particular, the energy and entropy densities are given by The EoS at vanishing chemical potential has been computed in the continuum limit in various setups and by different collaborations, see e.g.Refs.[30,31].It is thus convenient to separate the effects due to nonzero temperature from the modifications due to the presence of a nonzero isospin chemical potential.This is possible for the quantities of this section, for which the two types of contributions are added, but not for the speed of sound, defined via directional derivatives in the T -µ I plane.The pressure and the interaction measure can be written as p(T, µ I ) = p(T, 0) + ∆p(T, µ I ) , where the modifications of the EoS due to the isospin chemical potential, ∆p and ∆I, are the objects of interest in our study.Whenever we need to use results for p(T, 0) and I(T, 0), we use the results obtained from a reanalysis of the data of Ref. [30] with the parameterisation of I discussed in section 3 of Ref. [29].The correct inclusion of the correlations of the associated parameters are of particular relevance for the computation of the speed of sound.The details of the µ I = 0 data are discussed in appendix A.
A possible starting point for the computation of the EoS is the relation At vanishing chemical potential this equation is used to rewrite I as a derivative of the partition function with respect to the lattice scale [32][33][34][35].At nonzero chemical potential, one can follow a similar strategy to calculate the modifications of pressure and interaction measure using µ = 0 subtraction (see, e.g., Ref. [36] and Refs.[37,38]).The direct application of this method at nonzero µ I is discussed in appendix B, where we will see that it leads to large uncertainties for the interaction measure, due to the µ I = 0 subtractions.An alternative, which leads to more accurate results that we will present in the following, is to use a two-dimensional smooth interpolation of the results for n I to obtain the function n I (T, µ I ). 1 Since the isospin density is the derivative of the pressure with respect to µ I , see Eq. (2.1), the modification of the pressure can be computed from such an interpolation as [12,26,27] Inserting this into Eq.(2.4) gives for the interaction measure (see also Ref. [12]) The remaining task is to perform the interpolation and to obtain n I (T, µ I ) as a twodimensional function.Since the interpolation of an unknown function based on a discrete set of points subject to statistical uncertainties is an ill-posed inverse problem, the final interpolation will not be unique.This is already true for the interpolations used for the computation of the EoS at µ = 0 and will extend to the computation of the modifications of the EoS due to nonzero µ I discussed in the next section.The task at hand is to obtain an interpolation which is close to the actual physical solution while remaining as model-independent as possible.For the purpose of model-independence we average over all possible two-dimensional cubic spline interpolations with variable spline nodepoints (spline fits), weighted with the goodness of the description of the data.Here the goodness of the description is determined via the Akaike information criterion [39] and we have included a term to suppress unwanted (and unphysical) oscillatory solutions.Furthermore, the spline boundary conditions are chosen carefully to include the mandatory physical information on the interpolated function.The individual spline configurations are generated by a spline Monte-Carlo already introduced in Ref. [40] and discussed further in appendix C, where we also show a set of representative examples for this interpolation.
3 Results for the EoS Given the two-dimensional interpolation for the isospin density n I (T, µ I ), we can now proceed with the computation of the thermodynamic quantities we are interested in, namely the pressure p, the interaction measure I, as well as energy ϵ and entropy s densities using Eqs.(2.5), (2.6), (2.3) and (2.2).The modifications of the individual observables due to   µ I ≠ 0 are shown in Fig. 2. The pressure shows the strongest changes due to µ I for small temperatures. 2o show the magnitude of uncertainties and to have a more quantitative picture, we plot the pressure versus µ I for different temperatures in Fig. 3. Within the BEC phase the modification of the interaction measure shows an initial increase with µ I before it decreases towards larger µ I values.This phenomenon has already been observed at vanishing temperatures [23,28], where it leads to a negative interaction measure starting at around µ I /m π ≈ 0.84 [28], and is a clear signature of the presence of the BEC in the EoS, see also Ref. [12].We will discuss the full interaction measure at T ≠ 0 below.The decrease of ∆I is shifted to larger µ I values with increasing temperature, as can be seen from the right panel of Fig. 3, where we show the interaction measure versus µ I for a few different temperatures.The decrease is no longer visible at temperatures above around 150 MeVnote that this is just around the edge of the BEC phase for this particular lattice spacing.For these temperatures ∆I shows a peak in temperature direction if µ I /m π ≳ 0.5, indicating the strong influence of the phase transition on the interaction measure.The strong change of I in this region also translates to energy and entropy densities.The former generically shows a strong increase with µ I , which becomes less pronounced at larger temperatures outside of the BEC phase.The modifications of the entropy density mainly follow those of the interaction measure.
To get a feeling for the overall magnitudes, we plot the full results for p, I, ϵ and s in Fig. 4. To obtain the full results, we use the parameterization from Ref. [30] with the coefficients introduced in appendix A for the µ I = 0 quantities from Eq. (2.3). 3 We observe that the pressure generically rises with T and µ I .The exception is again the small T and large µ I region, where the increase seen in the plot, however, stems from the normalization with the temperature as explained above.A similar monotonic rise with T and µ I is generically observed for the energy density, even though the increase tends to become less pronounced for larger temperatures outside of the BEC phase.Note that the entropy density vanishes in the zero temperature limit for all values of µ I .
For the interaction measure we can clearly observe the decrease with µ I deep in the BEC phase for small temperatures.Furthermore, we also see a flattening and the onset of a peak-like structure for temperatures around the BEC phase boundary.This is an interesting observation.A peak of the interaction measure is present at µ I = 0 [30], see the gray curve in Fig. 5, at a temperature of around 190 MeV, i.e., above the phase transition.The observed flattening of I might thus be the onset of this peak structure, shifted towards smaller temperatures for increasing µ I .If this is the case, the peak approaches the thermal crossover at smaller values of µ I , until it becomes mostly consistent with the boundary of the BEC phase, see Fig. 1, around µ I /m π ≈ 0.8.This is further visualized in Fig. 5, where we show the interaction measure versus the temperature for different µ I /m π .Together with the shift of the onset of the plateau, we also observe the development of a narrower peak structure for larger µ I values.This might be an effect of the second order phase transition at the BEC phase boundary.A similar decrease of the temperature of the maximum with the chemical potential is also observed at nonzero baryon chemical potential, see [41,42] (as well as the large modifications of the interaction measure around T ≈ 155 MeV observed in Ref. [43]), for instance.

Lattice artifacts at N t = 8, 10 and 12
So far we have discussed results obtained for one value of N t , corresponding to one particular lattice spacing for each value of T .To obtain continuum results we have to increase N t while keeping the physical temperature and the aspect ratio N s /N t fixed.Unfortunately, for the isospin density a well controlled continuum extrapolation is not possible with the N t = 8, 10 and 12 lattices currently at our disposal.To show this, we plot the results for n I obtained from 8 × 24 3 , 10 × 28 3 and 12 × 36 3 lattices for two different temperatures in Fig. 6.For both temperatures, we observe that at small µ I the results of the N t = 12 lattices lie between those of the N t = 8 and 10 lattices.This indicates that higher order lattice Results for the interaction measure versus the temperature for different values of µ I /m π obtained on the 8 × 24 3 lattices from the two-dimensional spline interpolation for n I (T, µ I ) described in the text.We also compare to the µ I = 0 results from the parameterization of Ref. [30] described in appendix A.  artifacts are still present in this region and impede a proper continuum extrapolation.This behavior is in agreement with the large lattice artifacts observed for the leading order Taylor expansion coefficient in the direction of the isospin chemical potential, see, e.g.Ref. [44].The ordering of the results from the different N t values remains even for larger values of µ I at T = 149 MeV.At smaller temperatures the data rearranges in the vicinity of the phase transition to the BEC phase, such that the magnitude of n I increases from N t = 12 to N t = 10 and 8.
While a direct continuum extrapolation for the isospin density is not possible with the current dataset, we can still look at the magnitude of lattice artifacts for the other observables related to the EoS.We show the results for the pressure and the interaction measure from the lattices at different values of N t for three different temperatures in Fig. 7.The plot shows that the pressure, as the direct integral over isospin density, suffers from similar lattice artifacts as n I and again the results from the N t = 12 lattice are located between those of the N t = 8 and 10 lattices.For the interaction measure the situation is a bit different in the sense that all results typically overlap within the (comparably large) uncertainties.

The phase diagram in the T -n I plane
From the interpolation of n I we can also extract the phase diagram in the T -n I plane, of which a preliminary version has been presented in Ref. [28].To this end we use the phase boundaries in the T -µ I plane from Ref. [21] for the individual temporal extents, N t = 8 to 12, and determine the value of n I on these phase boundaries.The results for the chiral crossover and the BEC phase boundary for all N t are shown in Fig. 8.While the chiral crossover does not show any lattice spacing dependence, the BEC phase boundary exhibits slight lattice artifacts.
Due to the Silver blaze property at T = 0, the BEC phase transition starts at n I = T = 0.The BEC phase boundary at small temperatures and isospin denstities has been calculated via next-to-leading order chiral perturbation theory [45].We show those results up to T = 30 MeV (where the phase boundary in chiral perturbation theory starts to deviate significantly from the lattice result) together with the phase boundaries from the N t = 12 lattices in Fig. 9.
4 The speed of sound

Computation from the interpolation of the isospin density
Another interesting observable related to the EoS is the speed of sound c s .The isentropic speed of sound, i.e., the speed of a sound wave travelling isentropically through the plasma, is defined as where the subscript iso refers to the derivative taken in the direction of isentropic trajectories in phase space, i.e., for QCD at generic nonzero quark chemical potentials in the direction where and we have introduced the directional derivative ∂ ξ in this direction in parameter space.
For a pure isospin chemical potential the only relevant density is the isospin density, so that the condition (4.2) reduces to and we can write the directional derivative as where all quantities and derivatives can be obtained analytically from the spline interpolation for n I and the analytic form for the interaction measure at µ I = 0 from appendix A. Once α has been obtained, one can similarly analytically compute the directional derivatives of p and ϵ in Eq. (4.1).
Another quantity of interest related to the speed of sound, in particular for astrophysical and cosmological applications, is the polytropic index (see [8], as well as the lectures [46,47]) In the conformal limit, approached by QCD at asymptotically large densities or temperatures, it takes a value of γ = 1, while in the hadronic regime conformal symmetry is broken due to spontaneous chiral symmetry breaking, leading to large values of γ in the range of γ ≳ 2 (see the discussion in Ref. [8]).Consequently, γ can be seen as a measure in the EoS to distinguish between regions of hadron dominated matter (confined) or matter dominated by free quarks (quarkyonic/deconfined).In the study of Ref. [8] a value of γ = 1.75 has been introduced to distinguish between these two types of matter in neutron star cores.Finally, we also look at the normalized trace anomaly [48], which should be a number between −2/3 and 1/3 due to causality and thermodynamic stability.Furthermore, in Ref. [48] it has been argued that ∆ ≥ 0.

Speed of sound at vanishing temperature
Before discussing the results for the isentropic speed of sound in the parameter space of nonzero (T, µ I ), it is instructive to look at the limiting case of vanishing temperature.An initial study of the EoS at T = 0 on a coarse lattice with a ≈ 0.29 fm has already been presented in Refs.[23,28].Here we will present new results for the speed of sound at T = 0, obtained on 24 3 × 32 and 32 3 × 48 lattices at lattice spacings of a ≈ 0.22 fm and a ≈ 0.15 fm, respectively, including data up to µ I /m π ≈ 1.The results for these lattice spacings have already been presented partly in Ref. [49] where they also have been compared to the a ≈ 0.29 fm data.The starting point for the extraction of the EoS at zero temperature is again the isospin density, from which one can obtain the pressure and, consequently, all other thermodynamic quantities, using Eq.(2.5).Due to the Silver Blaze property, the isospin density vanishes outside of the BEC phase at T = 0.In practice, the simulations are performed at a small but non-vanishing temperature, so that residual temperature effects on n I need to be corrected in the vicinity of the transition.As already done in Ref. [23] we use chiral perturbation theory [13] to correct for these T ≠ 0 effects.In particular, we fit the results for n I for the two smallest values of µ I within the BEC phase, i.e., we include the data points up to µ I /m π ≈ 0.65, to the chiral perturbation theory expression for n I (see Eq. (B1) in appendix B of Ref. [23]).For our present ensembles this fit yields values of f π = 136(2) and 130(3) MeV, respectively, in good agreement with the physical value and the result obtained from the fit in Ref. [23] for the lattice with a ≈ 0.29 fm.To obtain a smooth interpolation for the isospin density, we again perform a model independent spline interpolation of the remaining data points via a spline Monte-Carlo as discussed above, where all splines are matched to chiral perturbation theory. 4The resulting interpolation of the isospin density is shown in the left panel of Fig. 10.
At T = 0 the condition of Eq. ( 4.3) is trivially fulfilled since s vanishes.Thus, the directional derivative is equivalent to the µ I -derviative, ∂ ξ = ∂ µ I .The resulting derivatives of p and ϵ in Eq. (4.1) can again be computed analytically.The results for the square of the isentropic speed of sound are shown in the right panel of Fig. 10, together with the conformal bound [24] as a gray dashed line.We observe that the squared speed of sound crosses the conformal bound at µ I /m π ≈ 0.643(4) and µ I /m π ≈ 0.646 (5), for a ≈ 0.22 and 0.15 fm, respectively, and reaches a peak at  We note that our finding of c 2 s > 1/3, as well as the development of a peak is in good agreement with recent results obtained in two-color QCD [37,38] (see also [50]).Furthermore, similar peaks in the speed of sound appear in quarkyonic models [51][52][53][54].At larger µ I , the speed of sound decreases and, on general grounds and according to perturbation theory [8], is expected to approach the conformal bound asymptotically from below.For our values of µ I , we currently do not see the decrease below the conformal bound.This would require simulations at yet higher isospin chemical potentials.
To allow for contact with recent studies on the EoS in neutron stars [8,[55][56][57][58][59], we plot c 2 s versus the energy density in the upper left panel of Fig. 11.Comparing to typical energy densities reached in the most massive neutron star cores (which are of the order of 10 3 MeV/fm 3 , see e.g.Refs.[55,59]), we see that the speed of sound reaches values which are larger than the conformal bound already for around one to two orders of magnitude smaller energy densities.In the upper right panel of Fig. 11, we show the polytropic index versus µ I .Just at the onset of pion condensation, γ assumes a value of γ ≥ 2, greater than the "quark matter bound" introduced in Ref. [8], in agreement with the prediction from chiral perturbation theory (yellow dashed curve).It then increases with µ I until it reaches its maximal value γ ≈ 2.4 around µ I /m π between 0.67 and 0.72, depending on the lattice spacing.Further increasing µ I , the polytropic index decreases below 1.75 and is seen to approach its conformal value of γ = 1 asymptotically.We note that the crossing of the "quark matter bound" might provide an alternative definition for the crossover [60,61], where effective degrees of freedom change from pions to Cooper pairs of u and d quarks.
Finally, the normalized trace anomaly ∆ is plotted in the bottom panel of Fig. 11.It starts at 1/3 at the onset of the BEC phase, in good agreement with chiral perturbation theory, and decreases towards larger µ I .Eventually it becomes negative between µ I /m π of 0.85 to 0.9 on the border of our parameter interval.Confirming the prediction of chiral perturbation theory, this shows a specific counter-example to the claim that ∆ would be strictly positive in QCD.

Speed of sound at nonzero temperatures
The results for the isentropic speed of sound at T ≠ 0 for the different lattices are shown in Fig. 12.For small temperatures, the speed of sound initially decreases slightly in the vicinity of the BEC phase boundary, before it starts to rise within the BEC phase.For N t = 8 and 10 it crosses the conformal bound around µ I /m π ≈ 0.75 to 0.8, i.e., at a somewhat larger µ I value than at T = 0.This is also visible in the left panel of Fig. 13, where we show the speed of sound obtained on the N t = 8 lattice for different temperatures including the  uncertainties.Contrary to the T = 0 case, we do not observe a maximum for c 2 s , likely due to the fact that it appears on the border or outside of our µ I range.
At this point we note that the speed of sound depends on the derivatives of n I , which are not well determined at the borders of our interpolation region at large µ I , as well as at our largest and smallest temperatures (at µ I = 0 it is determined by the parameterization from Ref. [30], see the lower left panel of Fig. 15).Consequently, the results for the speed of sound have to be taken with care beyond µ I /m π ≳ 0.85.To allow the reader to scrutinize the uncertainties at low T but large µ I , we show in the right panel of Fig. 13 c 2 s including its uncertainties for the different N t at a temperature of T = 120 MeV, close to the lower border of the temperature range.For N t = 12, the speed of sound does not reach the conformal bound, but we can still observe an increase of c s towards larger values of µ I .We interpret that as a shift of the maximum towards larger values of µ I when we approach the continuum.The question of the presence of the peak in the continuum limit can be answered once the results on finer lattices and at larger µ I become available.
Taking a closer look at the results from N t = 8, our most accurate results concerning the extraction of c s , albeit being furthest from the continuum, we show the squared speed of sound versus the energy density and the polytropic index γ versus µ I for a small temperature of T = 122 MeV in Fig. 14.Comparing the results for c s versus the energy density against those obtained at T = 0, shown in Fig. 11, we observe that the speed of sound crosses the conformal bound at about five times larger energy densities.Contrary to what is seen at T = 0, the polytropic index starts from a comparably small value γ = 0.82(4), and does not increase directly at the BEC phase boundary.The drastic increase towards larger values happens at µ I /m π ≈ 0.7 and it crosses the "quark matter bound" at around µ I /m π ≈ 0.85.Within our range of chemical potentials we see no sign of a maximum or the onset of a plateau.

Discussion and conclusions
In this paper we studied the QCD equation of state at nonzero temperatures and isospin chemical potentials using first-principles lattice simulations at physical quark masses.The results are shown in Figs. 2 and 4. We observe a strong rise of the pressure within the phase of condensing charged pions (BEC phase) at small temperature, which becomes milder when approaching the boundary of the BEC phase with T .We mention that our results for the pressure might also be useful to constrain the EoS at other chemical potentials.In our simulations the quark chemical potentials are set as µ u = −µ d = µ I and µ s = 0. Since flipping the sign of µ d merely amounts to a phase change in the fermion determinant, it is simple to see that for the partition functions Z(µ u = −µ d ) > Z(µ u = µ d ) holds.This implies that the pressure ∝ log Z at nonzero µ I provides an upper bound for the pressure in a setup with equal light quark chemical potentials.
The interaction measure was found to initially rise in the BEC phase, before it reaches a maximum and decreases.For small temperatures it eventually becomes negative deep in the BEC phase, providing an explicit counter-example to general positivity arguments [48].This effect diminishes as T grows, and around the high-T boundary of the BEC phase the impact of µ I is to shift the behavior of I(T ) to lower temperatures.We also determine the QCD phase diagram in the T -µ I plane, which for N t = 12 is shown in Fig. 9.
We have put a particular focus on the determination of the isentropic speed of sound c s , for which we show results both at zero and non-vanishing temperatures.At T = 0, c s initially increases, crosses the conformal bound of 1/ √ 3 around µ I /m π ≈ 0.64 and reaches a peak at µ I /m π ≈ 0.77 with a maximum of c 2 s ≈ 0.56, see Eq. (4.8) and Fig. 10.It then decreases again and is expected to approach the conformal bound asymptotically from below.To our knowledge, this is the first evidence for the explicit violation of this bound in first principles QCD.We note that similarly large c s values have also been observed at large µ B with the functional renormalization group approach in the point-like approximation for four-quark interactions [62].A remnant of the peak in c s remains visible at low temperatures (see Fig. 12), where it is shifted towards larger values of µ I .In Fig. 1 we also show a sketch of the region where the speed of sound exceeds the conformal bound in the phase diagram.
Besides c s , we also computed the polytropic index γ, which has been discussed as an indicator for the state of matter in the core of neutron stars.The results are shown in Figs.11 and 14 (right panels).At T = 0, the polytropic index starts from γ = 2, the value predicted by chiral perturbation theory [13], and then drops below γ = 1.75, the "quark matter bound" introduced in [8], around µ I /m π ≈ 0.8, eventually approaching the conformal value of γ = 1 asymptotically.At T ≠ 0 the behavior is quite different.The polytropic index starts from the µ I = 0 value, around γ = 0.8 to 1.0 -see Fig. 15 bottom right -from where it increases and approaches larger values of γ ≳ 1.75 only around µ I /m π ≳ 0.85.
Our tabulated results for the EoS will be useful for comparison to low-energy models and effective theories of QCD.To this end, our data for the full EoS, including a code to compute the observables, are available with the published version of this paper.Note that our calculations rely on three different lattice spacings, but owing to enhanced lattice artifacts at low temperatures for certain observables, we did not carry out a full continuum extrapolation here.
We finally comment on the consequences of our finding of the excess of the speed of sound over the conformal bound for the modeling of the EoS of neutron stars.With increasing amount of data on the masses and radii of the observed neutron stars in the Universe, several groups started to extract information on the QCD EoS at nonzero baryon density using this experimental data.In these approaches, the EoS is typically constrained at small and large densities from effective hadronic models (e.g.[63,64]) and perturbative QCD (e.g.[65]), respectively and then interpolated using a set of basis functions (e.g.[66][67][68][69] -see also [70,71] for non-parametric interpolations using neural networks).Newer studies use a large set of different types of basic functions and millions of different EoS interpolations [8,[55][56][57][58][59].While most recent studies indicate that experimental constraints favour a stiff EoS with a speed of sound that exceeds the conformal limit [55,[57][58][59]71], it has often been considered as extreme for the EoS to develop large speeds of sound of c 2 s > 0.5 to 0.6 or even to have an EoS which exceeds the conformal bound.In our study we provide direct evidence that an EoS with a speed of sound of this magnitude exists in QCD at small temperatures.Thus, such conditions for the EoS of cold dense QCD matter are certainly not unrealistic.
To compute the full EoS from the decompositions of Eq. (2.3), we need input at µ I = 0.Here we use the parameterization with t = T /(200 MeV), which has been employed in Refs.[29,30].Since we need to take the full correlations between the parameters into account for the correct computation of the uncertainties of derivatives and integrals, we use the parameters obtained from a reanalysis of the data from Ref. [30].The resulting parameters are listed in Tab. 1.We note the .Pressure, interaction measure, speed of sound squared and polytropic index in the temperature range relevant for this study, obtained from the parameterization at µ I = 0 of Ref. [30] discussed in the main text.
slight differences in the parameters compared to the ones obtained in Ref. [30].These differences in the reanalysis can be attributed to flat directions in parameter space, as already mentioned in Ref. [30], and do not lead to significant changes in the description of the data and the curves for the thermodynamic observables.For completeness, we show the results for the pressure, the interaction measure, the squared speed of sound and the polytropic index for this parameterization in the parameter region relevant for this study in Fig. 15.

B EoS from direct interaction measure calculations
An alternative to computing the EoS from the interpolation of the isospin density, is to rewrite the interaction measure via the derivative of the partition function with respect to the lattice scale, along similar lines as in Refs.[36][37][38].Using Eq. (2.6) and the derivatives of the lattice parameters with respect to the lattice scale, one obtains where β is the lattice coupling, S g the gauge action and am q the bare quark mass of flavour q in lattice units.We note that all the quantities appearing in Eq. (B.1) need to be renormalized properly, demanding, for instance, the knowledge of the quantities at T = 0, but nonzero µ I .Another way to ensure a proper renormalization is to make use of the decomposition of Eq. ( 2.3) and to compute only ∆I instead of I. Defining generically ∆O = ⟨O⟩ T,µ I − ⟨O⟩ T,0 , we obtain The remaining task is the computation of ∆S g and ∆ ψψ q .We show the results for ∆I(T, µ I ) versus µ I obtained for one particular temperature on a set of 24 3 ×6 lattices in the left panel of Fig. 16.For comparison we also show the result for ∆I obtained from the interpolation of the isospin density.We note that the qualitative behavior of the two data sets is similar, at least as far as comparison is possible due to the uncertainties for the direct computation, which are at least an order of magnitude larger compared to the uncertainties of the results obtained from the interpolation.The question is, where these large uncertainties orginate from.The problem is the combination of the λ-extrapolations at µ I ≠ 0 in combination with the subtraction of the µ I = 0 value (despite the fact that such an extrapolation is not necessary there).We show the λ-extrapolation of the Symanzik improved gauge action in the right panel of Fig. 16 in comparison to the result at µ I = 0 at the same temperature, which when subtracted give ∆S g .As is evident from the plot, the large uncertainties come from the subtraction of two quantities of similar magnitude, so that relative uncertainties are enhanced by orders of magnitude compared to the uncertainties of the individual quantities.This is particularly pronounced for the gauge action, but a similar behavior is also seen for the quark condensates.We note, that such a subtraction is absent for the isospin density, leading to way more accurate results for the quantity to interpolate and, consequently, for the EoS.

C Isospin density at λ = 0 and model-independent spline interpolations
The basic data for the isospin density, Eq. (2.1), is obtained at non-vanishing pion source, λ ≠ 0. For the extrapolations to λ = 0 we use the improvement program from Refs.[21][22][23].The resulting λ-extrapolations are basically flat and can be done using either a linear function in λ 2 or a constant.As an estimate for the systematic uncertainty associated with the extrapolation, we use the maximal deviation of the final result with either the extrapolation using the alternative (linear or constant), functional form or any of the two data-points at the lowest λ values.Note, that this is particularly important for the extraction of the equation of state, since an underestimation of uncertainties might lead to unphysical fluctuations which significantly affect the interpolation using spline fits.
For the extraction of the equation of state we use an average over all possible cubic spline interpolations of the isospin density in the two-dimensional parameter space (T, µ I ) with less grid-than data points (i.e., spline fits), weighted with an estimator for the "goodness" of the spline fit.Note, that for our spline fits the positions of the spline nodepoints generically do not coincide with the positions of the data points.As already discussed in Ref. [40], this average for an observable A (for instance the isospin density n I (T, µ I ) for given values T and µ I ) can be written as Here N G is the total number of spline nodepoints and N x (N G ) the number of nodepoints which can be varied in the particular spline setup.Note that N x and N G do not need to be equivalent (but always N x ≤ N G ), since some of the nodepoint positions can be fixed.This is the case for the nodepoints at the lower end of the n I -splines in the µ I direction, for instance, which are kept at µ I = 0.In Eq. (C.1), ⃗ x is the vector of (two-dimensional) nodepoint positions for the N x variable nodepoints.Note that typically the allowed range for the nodepoint values is restricted, as outlined below.The action S spl (⃗ x, N G ) represents the estimate for the "goodness" of the spline fit.Possible choices have already been discussed in Ref. [21].As the basic action we use the Akaike information criterion [39] (see also Ref. [72]), where N P is the number of parameters of the fit.One of the major problems for any spline interpolation or spline fit is the possible appearance of oscillatory solutions, i.e., solutions with additional minima and maxima as the spline attempts to capture all of the datapoints.These solutions can in particular be triggered by statistical fluctuations of data points and are particularly problematic for the equation of state, since additional unphysical minima and maxima might have strong effects on quantities like the speed of sound.To suppress those solutions we include another term in the action, following the spirit of Ref. [73].The term signifies the stability of the spline solution under small variations of the nodepoints.The parameters of the spline f n (with n = 1, . . ., N P ) are given either by the value of the spline on one of the spline nodepoints or by the derivatives on the nodepoints of the spline boundaries, depending on the particular spline setup.If we vary the nodepoints slightly and have a stable, non-oscillatory spline solution, we expect those values to not change significantly.I.e., given a variation α of one of the nodepoints, here with nodepoint index k and the variation in direction i, where ϵ is a small (not necessarily positive) number compared to the typical distance between two datapoints, we expect the spline parameters f α n to differ only slightly from the previous parameters f n .The parameter variation with respect to the typical statistical uncertainties for the parameters can be estimated by Here σ(f n ) is the statistical uncertainty of parameter f n , obtained from applying the same spline fit to the individual bootstrap samples for the data points.For stable fits we expect D α to be a number not much larger than one, of course depending on the typical order of magnitude of the statistical uncertainties and the typical change of the value of two consecutive data points.To suppress unwanted oscillatory solutions in the sum of Eq. (C.1), we add to the action the average of D α over all possible spline variations α, so that the total action is given by with the tunable parameter γ.
The tunable parameters of the spline average outlined above concern the possible numbers of nodepoints in each direction, possible constraints on the nodepoint locations, the boundary conditions of the spline, i.e., for each spline boundary one derivative for a cubic spline, as used here, and the parameter γ as well as the size of ϵ.For the interpolation of the isospin density, we use three to five nodepoints in each direction and demand that always two data points reside between the outermost and the consecutive nodepoint on each border of the grid and at least one data point lies between two consecutive nodepoints in each direction.The application of these constraints to the spline grids is straightforward if the data points themselves form a rectangular grid.The outer nodepoints in µ I direction are fixed at µ I = 0 and on this whole boundary we impose n I = 0. To account for the isospin density being an uneven function in µ I , we also impose ∂ 2 n I /(∂µ I ) 2 | µ I =0 = 0.The second derivatives in all other directions have been kept as free parameters for the spline fit and the positions of the outer gridpoints are allowed to vary.
To efficiently perform the sum from Eq. (C.1), we use Monte-Carlo methods as proposed in Ref. [21].In particular, we employ a Metropolis algorithm with a symmetric proposal probability for changes in the spline nodepoints.For ϵ we have chosen a random number between a tenth of the distance between the two nearest data points in positive or negative direction in such a way, that the nodepoint remains between the two data points after variation.To tune the parameter γ, we have performed several runs starting from small values of γ and monitored the resulting splines in the Markov chain.We stopped increasing γ when we found that no significant oscillations leading to local minima/maxima structures between two nodepoints showed up in the final average over spline configurations.As the final value we choose γ = 10.For our final results we have first tuned the number of nodepoints to the optimal value by performing 100 independent thermalisations allowing changes in the number of nodepoints and the nodepoint locations with respect to the data points.We then search for the spline with the lowest action and restrict ourselves to these number of nodepoints and the intervals in which the spline nodepoints reside with respect to the grid of data points.The final results are then obtained from 100 splines obtained by 20 independent runs where we vary the nodepoints in this constrained setup with 20000 thermalisation updates and with 5 spline configurations each, separated by 10000 spline updates.G , as well as the intervals in which the nodepoint positions have been varied are given in Tab. 2. We show the location of the 100 final nodepoint sets used in the analysis in Fig. 17.Note that the outer spline nodepoints can lie anywhere outside of the data point interval and that the outer nodepoint of the lower boundary in µ I direction has been held fixed at µ I = 0.The resulting interpolation for N t = 8 lattices is shown for a set of temperatures in Fig. 18.The uncertainties include the uncertainty due to the individual data points (computed using the bootstrap procedure with 1000 samples) and from the Monte-Carlo over spline interpolations.As a comment concerning the T = 0 splines, in this case we are speaking about a one-dimensional spline interpolation, which is rather well behaved and uncritical concerning the spline Monte-Carlo.In this case, the number of nodepoints and their location has been allowed to vary freely up to the maximal possible number of nodepoints given the number of available data points with the matching to chiral perturbation theory as discussed in section 4.2.

D Simulation Details
We provide the run parameters for the simulations at T ≠ 0 and T = 0 in Tab. 3.For each of the mentioned parameter values we have simulated at up to five different values of the pion source parameter λ with 1.4 > λ/m ud ≳ 0.05.At T = 0 we have used a fixed number of three different λ values.For each of these parameter sets we generated between 500 and 2000 trajectories, measuring observables (here the isospin density n I ) on every fifth configuration.Typically we use two to three independent chains to acquire the full statistics.We estimate autocorrelations using the integrated autocorrelation time of the plaquette expectation value.The values we obtain on the different ensembles are mostly between 5 and 10 in molecular dynamics units, with the tendency to large values for smaller temperatures and λ-values, as well as for larger values of µ I and N t .While this indicates that typically two consecutive configurations are correlated, the autocorrelation times obtained for n I , the main observable of our study, are much smaller.In rare cases we observed autocorrelation times of the order of 20 in molecular dynamics units, so that 4 consecutive measurements are correlated.Generically, we have checked the error analysis by using binning prior to the bootstrap procedure, but did not observe a significant dependence of the uncertainties on the binsize.Furthermore, in the improved λ-extrapolation, values at different λ are combined using a linear function to extract the λ = 0 result, averaging out the fluctuations of the individual simulation points.

Figure 1 .
Figure1.Phase diagram of isospin asymmetric QCD determined in Ref.[21].The green band is the phase boundary to the phase with Bose-Einstein condensation (BEC), the blue band is the prolongation of the chiral crossover in the T -µ I plane and the red data point marks the pseudo-triple point, the meeting point of the crossover and the BEC phase.In yellow we also show a sketch of the region, where the speed of sound exceeds its conformal limit.

3 Figure 2 .
Figure 2. Results for the modifications of the pressure (top left), the interaction measure (top right), the energy density (bottom left) and the entropy density (bottom right) due to nonzero isospin chemical potential on the 8 × 24 3 lattices.Uncertainties are not shown for better visibility.

Figure 3 .
Figure 3. Results for the modifications of the pressure (left) and the interaction measure (right) versus µ I for different temperatures obtained on the 8 × 24 3 lattices from the two-dimensional spline interpolation for n I (T, µ I ) described in the text.

3 Figure 4 .
Figure 4. Results for the pressure (top left), the interaction measure (top right), the energy density (bottom left) and the entropy density (bottom right) in units of the temperature obtained on the 8 × 24 3 lattices from the two-dimensional spline interpolation for n I (T, µ I ) described in the text.Uncertainties are not shown for better visibility.

4 µ
Figure 5.Results for the interaction measure versus the temperature for different values of µ I /m π obtained on the 8 × 24 3 lattices from the two-dimensional spline interpolation for n I (T, µ I ) described in the text.We also compare to the µ I = 0 results from the parameterization of Ref.[30] described in appendix A.

Figure 7 .
Figure 7. Results for the pressure (left) and the interaction measure (right) for temperatures T = 120, 145 and 165 MeV (top to bottom) from N t = 8, 10 and 12 lattices.

n I /m 3 πTN t = 8 N t = 10 N t = 12 Figure 8 .Figure 9 .
Figure 8. Results for the chiral crossover temperature (left) and the BEC phase boundary (right) in the T -n I plane for different values of N t , based on our results from Ref. [21].

Figure 10 .
Figure 10.Left: Results for the isospin density together with the spline interpolation at T = 0 obtained on 24 3 × 32 and 32 3 × 48 lattices and lattice spacings of a ≈ 0.22 and 0.15 fm.The yellow part of the curve is obtained directly from the chiral perturbation theory expression for n I .Right: Results for the isentropic speed of sound at T = 0, obtained from the spline interpolation of the left panel.Also shown are the chiral perturbation theory result (dashed yellow line) with the pion decay constant obtained from the fit discussed in the text as well as the conformal bound[24] (dashed gray line).

µFigure 11 .
Figure 11.Left: Results for the isentropic speed of sound at T = 0 versus the energy density in the BEC phase.Also shown is the conformal bound[24] (gray dashed line).Right: Results for the polytropic index γ at T = 0. Also shown is the result from chiral perturbation theory (yellow dashed line) and the "quark matter bound" introduced in Ref.[8] (gray dashed line).Bottom: Results for the normalized conformal anomaly ∆ at T = 0, together with the result from chiral perturbation theory (yellow dashed line).

/m π c 2 sFigure 12 .
Figure 12. Results for the isentropic speed of sound at T ≠ 0, obtained on the N t = 8, 10 and 12 lattices (from top left to bottom left) in the T -µ I plane, excluding uncertainties.

Figure 13 .
Figure 13.Results for the isentropic speed of sound at T ≠ 0 at different temperatures on the 8 × 24 3 lattices (left panel) and at a temperature of 120 MeV for N t = 8, 10 and 12 (right panel).

Figure 14 .
Figure 14.Left: Results for the isentropic speed of sound versus the energy density at T = 122 MeV obtained on the 8 × 24 3 lattices.Also shown is the conformal bound[24] (gray dashed line).Right: Results for the polytropic index γ at T = 122 MeV.Also shown is the "quark matter bound" introduced in Ref.[8] (gray dashed line).

Figure 15
Figure15.Pressure, interaction measure, speed of sound squared and polytropic index in the temperature range relevant for this study, obtained from the parameterization at µ I = 0 of Ref.[30] discussed in the main text.

Figure 16 .
Figure 16.Left: Results for the interaction measure obtained from a direct computation explained in the text (blue points) and obtained from the two-dimensional interpolation of the density (red band) on 24 × 6 lattices at a temperature of 136 MeV.Right: λ-extrapolation (blue band) of the Symanzik improved gauge action (blue open circles) on a 24 3 × 6 lattice with T = 136 MeV and µ I /m π = 0.25.We also show the result at µ I = 0 for the same temperature (black filled box and black band) which needs to be subtracted from the λ = 0 result to obtain the contribution to ∆I.

Figure 17 .
Figure 17.Summary of the location of spline nodepoints in the final run.Note that the spline nodepoints can only be generated in certain regions of the parameter space, indicated by the accumulation of dots.The final numbers of nodepoints in the different directionsN G = N (T ) G × N (µ I )

3 µFigure 18 .
Figure 18.Comparison of the results for the isospin density n I versus µ I for different temperatures (left) and T for different µ I (right) obtained from the simulations on 8 × 24 3 lattices and the associated spline interpolation described in the text.

Table 2 .
Number of spline nodepoints in the T and µ I directions and the associated intervals in which the Monte-Carlo-generated nodepoints reside.