Equation of state for hot QCD and compact stars from a mean field approach

TL;DR

This paper develops a mean-field model for hot QCD matter that accurately describes lattice QCD data, nuclear matter, and neutron star properties, revealing multiple phase transitions and their effects on matter behavior.

Contribution

The paper introduces a comprehensive CMF model that unifies lattice QCD results with nuclear and neutron star matter properties, capturing multiple phase transitions in a single framework.

Findings

Predicts three consecutive cross-over transitions in QCD matter.

Accurately describes neutron star properties and merger observations.

Reproduces lattice QCD thermodynamics at zero chemical potential.

Abstract

The thermodynamic properties of high temperature and high density QCD-matter are explored within the Chiral SU(3)-flavor parity-doublet Polyakov-loop quark-hadron mean-field model, CMF. The quark sector of the CMF model is tuned to describe the $\mu_B=0$ thermodynamics data of lattice QCD. The resulting lines of constant physical variables as well as the baryon number susceptibilities are studied in some detail in the temperature/chemical potential plane. The CMF model predicts three consecutive transitions, the nuclear first-order liquid-vapor phase transition, chiral symmetry restoration, and the cross-over transition to a quark-dominated phase. All three phenomena are cross-over, for most of the $T-\mu_B$-plane. The deviations from the free ideal hadron gas baseline at $\mu_B=0$ and $T\approx 100-200$ MeV can be attributed to remnants of the liquid-vapor first order phase transition in nuclear matter. The chiral crossing transition determines the baryon fluctuations at much higher $\mu_B\approx1.5$ GeV, and at even higher baryon densities $\mu_B\approx2.4$ GeV, the behavior of fluctuations is controlled by the deconfinement cross-over. The CMF model also describe well the static properties of high $\mu_B$ neutron stars as well as the new neutron star merger observations. The effective EoS presented here describes simultaneously lattice QCD results at $\mu_B=0$, as well as observed physical phenomena (nuclear matter and neutron star matter) at $T\cong0$ and high densities, $\mu_B>1$ GeV.

Figures

Tables

Table 1. Table 1: Best fit values of parameters extracted from a scan over the parameter space.

$T_0$ (MeV)	$a_1$	$a_2$	$b_3$	$g_{q\sigma }=g_{s\zeta }$
180.0	-11.67	9.33	-0.53	-1.0

Table 2. Table 2: Initial state properties obtained from the one-dimensional stationary case of the central heavy ion collision – Taub adiabat  Taub ( 1948 ); Thorne ( 1973 ) . The entropy per baryon S / A 𝑆 𝐴 S/A , temperature T 𝑇 T , the initial baryon density n B / n 0 subscript 𝑛 𝐵 subscript 𝑛 0 n_{B}/n_{0} and the quark fraction 1 3 ​ n q / n B 1 3 subscript 𝑛 𝑞 subscript 𝑛 𝐵 \frac{1}{3}n_{q}/n_{B} are presented for various colliding energies.

$\sqrt{s_{NN}}$ (GeV)	$E_{\mathrm{lab}}$ (GeV)	$S/A$	$T$ (MeV)	$n_B/n_0$	$\frac{1}{3}{n}_q/n_B$
2.2	1.6	2.9	63.3	2.6	0.0
2.4	2.1	3.5	71.7	3.2	0.1
2.6	2.7	4.0	79.5	3.7	0.2
3.0	3.9	5.1	96.6	4.9	0.3
3.5	5.6	6.4	120.0	6.0	0.4
4.5	9.9	8.3	158.9	8.6	0.7
6.2	19.6	10.5	205.2	14.9	0.9
7.7	30.7	13.0	245.7	18.9	1.0
9.1	43.2	15.0	278.9	21.3	1.0

Equations

$$m_{\text{i}\pm }^{\ast }$$

$$V=V_0+\frac{1}{2}{k}_0{I}_2-k_1{I}_2^2-k_2{I}_4+k_6{I}_6,$$

$$I_2=({\sigma }^2+{\zeta }^2)$$

$$I_6$$

$${\Omega }_q=-T\sum _{i\in Q}\frac{d_i}{(2\pi {)}^3}\int d^{3}k\ln \left(1+\Phi \exp \frac{-\left(E_i^{\ast }-{\mu }_i^{\ast }\right)}{T}\right)$$

$${\Omega }_{\overline{q}}=-T\sum _{i\in Q}\frac{d_i}{(2\pi {)}^3}\int d^{3}k\ln \left(1+{\Phi }^{\ast }\exp \frac{-\left(E_i^{\ast }+{\mu }_i^{\ast }\right)}{T}\right).$$

$$m_q^{\ast }$$

$$m_s^{\ast }$$

$$U$$

$$+b(T)\log [1-6\Phi {\Phi }^{\ast }+4({\Phi }^3+{\Phi }^{\ast 3})-3(\Phi {\Phi }^{\ast }{)}^2],$$

$$a(T)$$

$${\rho }_i=\frac{{\rho }_i^{\mathrm{id}}(T,{\mu }_i^{\ast }-v_{i}p)}{1+\sum _j{v}_j{\rho }_j^{\mathrm{id}}(T,{\mu }_j^{\ast }-v_{j}p)},$$

$$\frac{I}{T^4}=\frac{\epsilon -3P}{T^4}.$$

$${\chi }_n^B=\frac{{\partial }^n(P/T^4)}{(\partial {\mu }_B/T{)}^n},$$

$$T_f({\mu }_B)=T_0\left(1-{\kappa }_2^f{\left(\frac{{\mu }_B}{T_0}\right)}^2-{\kappa }_4^f{\left(\frac{{\mu }_B}{T_0}\right)}^4\right).$$

$${\kappa }_2^f$$

$$f_{2n}(T)=\frac{1}{(2n)!}{\frac{{\partial }^{2n}f(T,{\mu }_B)}{\partial {\left({\mu }_B/T\right)}^{2n}}|}_{(T,0)}.$$

$$c_{s/n_B}^2=\frac{n_B^2{\partial }_{T}s-2n_{B}s{\partial }_{{\mu }_B}s+s^2{\partial }_{{\mu }_B}{n}_B}{(ϵ+P)\left({\partial }_{{\mu }_B}{n}_B{\partial }_{T}s-{\left({\partial }_T{n}_B\right)}^2\right)}$$

$$(P_0+{\epsilon }_0)(P+{\epsilon }_0)n^2=(P_0+\epsilon )(P+\epsilon )n_0^2,$$

$${\gamma }^{\mathrm{CM}}=\frac{\epsilon n_0}{{\epsilon }_{0}n},{\gamma }^{\mathrm{CM}}=\sqrt{\frac{1}{2}\left(1+\frac{E_{\mathrm{lab}}}{m_N}\right)}.$$

$$Q_{\text{ij}}=-\lambda {\mathcal{E}}_{ij}.$$

$$\lambda =\frac{2}{3}{k}_2{R}^5.$$

$$\Lambda =\frac{\lambda }{M^5}=\frac{2}{3}{k}_2{\left(\frac{R}{M}\right)}^5.$$

Full text

I Introduction

At low temperatures and high baryonic densities a transition from hadronic matter to a cold deconfined state with quark degrees of freedom is expected McLerran and Reddy (2019), that is supported by perturbative QCD calculations Kurkela et al. (2010). This transition is a particular case of appearance of a deconfined state of quark-gluon plasma (QGP) that expected at large densities and large temperatures Shuryak (1980). The study of the unknown QCD phase diagram is one of the main motivations for the state-of-the-art research in the nuclear and particle physics community. A rich phase structure is conjectured for finite temperatures and chemical potentials Rajagopal (1999); Alford (2003); Buballa (2005); Schäfer (2003); Fukushima and Hatsuda (2011). The QCD phase structure is a most important ingredient for the understanding of the Early Universe, ultra relativistic heavy ion collisions, and of the evolution, structure and inspirals of neutron stars Takami et al. (2015). Even though QCD is a well-established theory with only a few parameters, perturbative calculations are inappropriate in the crossover regions of $T$ and $\mu_{B}$ discussed here due to the large values of the QCD coupling constant at scales relevant for most of these applications Braaten and Nieto (1996). The infamous sign problem de Forcrand (2009) prohibits lattice QCD (LQCD) calculations at finite densities.

The running experiments at the Large Hadron Collider (LHC), at the Relativistic Heavy Ion Collider (RHIC), at the Super Proton Synchrotron (SPS), and at the Heavy Ion Synchrotron (SIS) provide state-of-the-art data of the measurements of properties of matter produced in heavy-ion collisions (HIC). The detailed information of the particle production as measured in these experiments allow to extract both thermodynamic and kinetic characteristics of the system that is created.

Astrophysical observations of compact stars, along with data from the recent gravitational-wave detection by LIGO provide an additional tool to probe the equation of state of dense nuclear and possible quark matter Annala et al. (2018); Most et al. (2018); Abbott et al. (2018); Most et al. (2019); Bauswein et al. (2019); Hanauske et al. (2017, 2019) in the region of moderate temperatures and high baryon densities, close to those as created in HIC.

First principle LQCD calculations suggest a smooth crossover transition at vanishing baryochemical potential $\mu_{\rm B}=0$ from hadronic to partonic degrees of freedom Aoki et al. (2006). Although there is no indication of a first or a second order phase transition in the energy density, pressure, entropy density and speed of sound, there are other observable in LQCD – chiral susceptibilities which seem to indicate a chiral cross-over at a pseudocritical temperature $T_{\rm pc}\approx 155$ MeV Borsanyi et al. (2014); Bazavov et al. (2014). However, the extension of LQCD calculations to finite $\mu_{\rm B}$ is a difficult problem. There are approximate lattice methods, such as extrapolations by Taylor expansion and the analytic continuation from purely imaginary to real $\mu_{B}$ are reasonable only for small baryon densities. The exploration of higher baryon densities require effective QCD models, which respect the known symmetries of QCD and describe appropriately the known phenomenology of strong interactions.

The current knowledge of the properties of strongly interacting matter suggests a number of features that ought to be incorporated in any reasonable effective QCD model:

• •

First, nuclear matter and the nuclear liquid-vapor phase transition at moderate temperatures, close to the nuclear saturation density Pochodzalla et al. (1995);

• •

Second, the chiral symmetry restoration should lead to the Stefan-Boltzmann limit for the thermodynamic properties at high temperature and at high chemical potential Borsanyi et al. (2010);

• •

The transition from hadronic- to quark-gluon degrees of freedom, at high temperatures and/or chemical potentials, is a crucial ingredient for the consistent description of QCD matter.

Often, these different aspects of QCD are modeled within separate frameworks, which are then merged through various constructions.

The present work formulates a single combined framework to describing the QCD thermodynamics, which simultaneously satisfies all the constraints, from lattice QCD, and known nuclear matter properties, as well as neutron star observations.

The resulting equation of state is then used to estimate various properties of systems created in both heavy-ion collisions and neutron star physics. Sec. II presents the description of the Chiral SU(3)-flavor parity-doublet Polyakov-loop quark-hadron mean field model, CMF. Section III describes the fine-tuning of model parameters to the $\mu_{\rm B}=0$ LQCD data on the trace anomaly, and presents a comparison with the CMF model predictions for various conserved charge number fluctuations with the corresponding lattice data. The QCD phase diagram deduced from the CMF model is studied in Sec. IV. The creation of hot and dense QCD matter as created in heavy-ion collisions at various collision energies is studied using the 1-dimensional hydrodynamics in Sec. V, the respective trajectories along the QCD phase diagram are explored as well. Section VI presents the CMF model predictions for the observed neutron star properties.

II Chiral SU(3)-flavor parity-doublet Polyakov-loop quark-hadron mean-field model, CMF

The Chiral SU(3)-flavor parity-doublet Polyakov-loop quark-hadron mean-field model, CMF, is an extension of the previously proposed $\sigma$-$\omega$ model with parity doubling for nuclear and hadron matter Detar and Kunihiro (1989); Hatsuda and Prakash (1989); Papazoglou et al. (1997, 1998, 1999); Sasaki and Mishustin (2010). The CMF model was extended to include quark degrees of freedom Steinheimer et al. (2011a, b); Dexheimer et al. (2013); Mukherjee et al. (2017); Motornenko et al. (2019). This model is a phenomenological effective unified approach to describe interacting hadron-quark matter. The Lagrangian includes essential symmetries and features of QCD. These include:

• •

Chiral symmetry restoration in the hadronic sector, in particular the baryon parity doubling so an explicit mass term for baryons is possible, even when chiral symmetry is restored. This leads to a restoration of mass degeneracy among baryons and their respective parity partners Aarts et al. (2017, 2018).

• •

Eigenvolume corrections for hadrons, which allow for an effective modeling of their repulsive interactions. This suppresses hadronic densities and ensures a transition to a parton-dominated matter at large densities, when quark and gluon d.o.f. appear.

• •

Chiral symmetry restoration for quarks and a dynamical generation of their masses.

• •

The Polyakov loop via a QCD-motivated potential incorporates the deconfinement transition.

A detailed description of the hadronic part of the CMF model can be found in the literature Dexheimer et al. (2008); Steinheimer et al. (2011a, b); Mukherjee et al. (2017). It is based on a realization of a $\sigma$-$\omega$ model in mean-field description. Here the relevant fermionic degrees of freedom are baryons that interact through mesonic mean-fields. The version of this model used here includes all states in the SU(3)f baryon octet, together with their parity partners, i.e. states with the same quantum numbers but opposite parity. The LQCD data suggests that the same mechanism should be implemented for higher baryonic states – the baryon octet Aarts et al. (2018), this is a plan for future studies. In the limit of chiral symmetry restoration these parity partner states should be degenerate and their masses equal, which then serves as a signal for chiral symmetry restoration. To allow for such a behavior, the baryon masses are dynamically generated by their couplings to the scalar $\sigma$-field and the scalar strange $\zeta$-field, which serve as the order parameters for the chiral transition:

[TABLE]

Here $+$ stand for positive and $-$ for negative parity states, $g^{\text{(j)}}_{\text{i}}$ are the coupling constants of baryons to the two scalar fields, $m_{0}=759$ MeV is the baryon mass at the restored phase. In addition, there is an SU(3)f symmetry-breaking mass term proportional to the strangeness content of the baryons, there $n_{s}$ is the number of strange quarks in the baryon, and $m_{s}=130$ MeV is the mass of the strange quark. The couplings $g^{\text{(j)}}_{\text{i}}$ are tuned to reproduce the vacuum masses of baryons.

The mean-field values of the chiral fields are driven by the thermal contribution from baryons and quarks, and controlled by the scalar meson interaction, driving the spontaneous breaking of the chiral symmetry:

[TABLE]

with

[TABLE]

where $V_{0}$ is fixed by demanding that the potential vanishes in the vacuum.

The parameters of the scalar and the vector interactions are fitted to describe nuclear matter properties Steinheimer and Schramm (2014); Mukherjee et al. (2017). Contributions of all established hadronic resonances are included here with their vacuum masses Tanabashi et al. (2018). These states can be coupled to meson fields as parity doublets as well. However, this is not done in the current implementation, they only interact with the other particles via their excluded volume only.

The quark degrees of freedom are incorporated similarly with PNJL approach Fukushima (2004). The appearance of quarks is controlled by the value of Polyakov loop $\Phi$, which plays the role of the order parameter for the deconfinement transition. The coupling of the quarks to the Polyakov loop is introduced through the thermal energy of the quarks. Their thermal contribution to the grand canonical potential $\Omega$ is given by

[TABLE]

and

[TABLE]

The sums run over all light quark flavors (u, d, and s), $d_{i}$ is the corresponding degeneracy factor, $E^{*}_{i}=\sqrt{m^{*2}_{i}+p^{2}}$ is the energy, and $\mu^{*}$ is the chemical potential of the quark. Note, two- and three- quark contributions to $\Omega$ are omitted in the CMF model since hadronic excitations are explicitly included.

The effective masses of the light quarks are generated by the $\sigma$ field (non-strange chiral condensate) as well, the mass of the strange quark is generated by the $\zeta$ field (strange quark-antiquark state). The small explicit mass terms $\delta m_{q}=5$ MeV, and $\delta m_{s}=150$ MeV for the strange quark, and $m_{0}q=253$ MeV which corresponds to an explicit mass term which does not originate from chiral symmetry breaking:

[TABLE]

Dynamics of the Polyakov-loop is controlled by the effective Polyakov-loop potential $U(\Phi,\Phi^{*},T)$ Ratti et al. (2006):

[TABLE]

The parameters of this potential can be fixed to the lattice QCD data in the pure gauge sector Ratti et al. (2006). However, this yields an unsatisfactory description of the (2+1)-flavor QCD thermodynamics when the hadrons are explicitly included in the model. Therefore the parameters of the Polyakov Loop potential are adjusted in the present work to describe properly the (2+1)-flavor lattice data.

The CMF model incorporates excluded-volume effects, in order to suppress the hadronic degrees of freedom in the regions of the phase diagram where physically quarks and gluons dominate Steinheimer et al. (2011b). Consequently, all the thermodynamics densities, including the quark contribution, are reduced as parts of the system are occupied by EV-hadrons:

[TABLE]

the $v_{j}$ are the eigenvolume parameters for the different species. $p$ is a system pressure without contribution of mean fields, $\mu^{*}$ is the chemical potential of the hadron. The $v$ is assumed to be $v_{B}=1$ fm3 for (anti-) baryons, $v_{M}=1/8$ fm3 for mesons, and is set to zero $v_{q}=0$ for quarks.

III Constraining the CMF model to the lattice data

To introduce the constraints on the CMF model parameters from lattice QCD at high temperature and zero net-baryon density we use the QCD trace anomaly $I$, “the interaction measure”, as a reference:

[TABLE]

The free parameters of the present model are tuned to reproduce this quantity. The LQCD trace anomaly permits to calculate all observable thermodynamic quantities. Analysis of other Lattice data (chiral susceptibility) seem to show that for chirally related observables there is a crossover transition with a “pseudo”-critical temperature at $T\approx 156$ MeV. The analysis of this data by a phenomenological model suggests a half-hadron, half-quark composition in that region Albright et al. (2014).

The parameters of the CMF model’s quark sector needed to reproduce the trace anomaly data from LQCD are found by a least mean squares fitting procedure for the parameters of the Polyakov loop potential $U(\Phi,\Phi^{*},T)$ and for the coupling constants $g_{q\sigma}$ and $g_{s\zeta}$ of the quarks to the chiral condensates $\sigma$ and $\zeta$, respectively. All in all this fixes 5 model parameters, $T_{0},a_{1},a_{2}$, $b_{3},~{}g_{q\sigma}=g_{s\zeta}$ (we set $g_{q\sigma}$ and $g_{s\zeta}$ to the same value). The quark parameter fitting is performed through a scan over the parameter space on a $8\times 6\times 7\times 6\times 6$ sized grid, minimizing the root mean square deviation of the CMF model data on $I/T^{4}$ and computed on the lattice results. The resulting parameter values are presented in Table 1. The comparison of the CMF model to the lattice data is shown in Fig. 1.

The values in Table 1 suggest couplings $g_{q\sigma},~{}g_{s\zeta}$ of quarks to the chiral fields of about $1/3$ of the baryons, this one may expect from the additive quark model. Larger values of quark couplings would significantly influence the size of the peak in interaction measure $I/T^{4}$, as studied in Ref. Mukherjee et al. (2017). In case of large values of $g_{q\sigma},g_{s\zeta}$ the strong interplay between the chiral symmetry restoration and the deconfinement transition would result in too large values of the interaction measure. The large peaks in the baryon number susceptibilities are in contrast to the lattice data.

Higher order baryon number susceptibilities $\chi_{n}^{B}$ that are a LQCD measure of the particle number fluctuations

[TABLE]

as well as the curvatures of various lines of constant physical quantities are also interesting in scope of LQCD data.

The behavior of the $\chi_{4}^{B}/\chi_{2}^{B}$ at $\mu_{B}=0$ is presented in this section, the study for finite values of $\mu_{B}$ is presented in Sec. IV. The comparison of the CMF model with the available LQCD data for the $\chi_{4}^{B}/\chi_{2}^{B}$ is shown in Fig. 2, indicating a fair agreement of the CMF model with the lattice data.

Lattice QCD studies often explore regions of finite $\mu_{B}$ by using the Taylor series expansion. The Taylor expansion in series of $T$ and $\mu_{B}$ up to $\mathcal{O}(\mu_{B}^{4})$ was used in Bazavov et al. (2017a) to calculate “lines of constant physics”: i.e. lines in the $T-\mu_{B}$ plane where certain thermodynamic quantities like pressure, energy density, and entropy density $P,\varepsilon,s$, are constant. The coefficients $\kappa^{f}_{2}$ and $\kappa^{f}_{4}$ ($f\equiv P,\varepsilon,s$) represent these contour lines in the $T-\mu_{B}$ plane using the following parametrization Bazavov et al. (2017a):

[TABLE]

Here the coefficients $\kappa^{f}_{2}$ and $\kappa^{f}_{4}$ are calculated from Eqs. 12 and 13, see Ref. Bazavov et al. (2017a) for details:

[TABLE]

Here

[TABLE]

The coefficients $\kappa^{f}_{2}$ and $\kappa^{f}_{4}$ are calculated in the CMF model for the pressure $P$, the energy density $\varepsilon$, and the entropy density $s$ as functions of the temperature $T$. The CMF model predictions are in a reasonable agreement with recent LQCD calculations Bazavov et al. (2017a). The rather low values of $\kappa_{2}^{f}$ and $\kappa_{4}^{f}$ suggest also small curvatures of lines of constant physical observables in the temperature region studied here. Effects of the finite chemical potential are small, therefore these lines are almost horizontal in the $T$-$\mu_{B}$ plane. The coefficients for the entropy and the energy density indicate that $\kappa_{2}^{s}<\kappa_{2}^{\varepsilon}$, meaning a decrease of the entropy density along the lines of constant energy density.

IV The CMF model phase diagram

Two order parameters, the chiral condensate $\sigma$ and the Polyakov loop $\Phi$, plus the interacting baryon octet within the SU(3)-flavor $\sigma$-model permit four different phases within the CMF model. These phases are characterized as:

• •

A dilute gas of interacting hadrons;

• •

A hadronic liquid – a dense hadronic phase, the transition from the hadron gas to the hadronic liquid is the nuclear liquid-vapor phase transition. Quarks start to appear in the hadronic liquid, but they are negligible;

• •

A chirally restored phase, where the mass symmetry between the parity partners is restored. Here the quark masses are decreased, hence quarks give a sizable contribution to the thermodynamics;

• •

A quark dominated phase, where most of the energy density is carried by quark and gluon degrees of freedom. The gluon contribution is modeled by the Polyakov loop potential Ratti et al. (2006) and $\epsilon_{q}/\epsilon_{\rm tot}\approx 1$.

The baryon number susceptibilities $\chi_{n}^{B}$, that can be calculated using Eq. 10, are proportional to the respective cumulants of the baryon number distribution. Higher-order baryon number susceptibilities do increase in proportion to the increasing powers of the correlation length Stephanov (2009). Such an increase in correlation length would be reflected in large values of 2nd and higher-order susceptibilities in the vicinity of a critical point and in the region of phase transition. Hence, these quantities are useful indicators of a critical behavior in the CMF model. Deviations of $\chi_{n}^{B}$ from the corresponding baselines indicate a transformation between different phases, which is reflected usually in a non-monotonic behavior of these observables, e.g. skewness $\chi_{3}^{B}/\chi_{2}^{B}$ and kurtosis $\chi_{4}^{B}/\chi_{2}^{B}$.

The calculated skewness ($\chi_{3}^{B}/\chi_{2}^{B}$) and kurtosis ($\chi_{4}^{B}/\chi_{2}^{B}$) (Fig. 4) in the CMF model exhibit non-trivial structures in the $T$-$\mu_{B}$ phase diagram. The regions of deviations from baseline separate regions with quantitatively different properties which are often dubbed – “phases”. Note the sharp phase boundaries indicate first order phase transitions, FOPT, these are only observable at quite moderate temperatures $T<50$ MeV. The hadron phase located at both low temperatures $T$ and baryon chemical potential $\mu_{B}$ represents a dilute gas of interacting hadrons. There, the fluctuation measures $\chi_{2}^{B}/\chi_{1}^{B}$, $\chi_{3}^{B}/\chi_{2}^{B}$, $\chi_{4}^{B}/\chi_{2}^{B}$ are quite close to unity, consistent with the Skellam distribution baseline. The system exhibit a FOPT to a dense hadronic liquid phase with rising chemical potential $\mu_{B}\approx 1$ GeV. Here, fluctuations are reduced due to the repulsive interactions. Quarks start to appear in moderation. The liquid phase exhibits additional FOPT at $\mu_{B}\approx 1.5$ GeV, and second order transition at $\mu_{B}\approx 2.4$ GeV, which show up in the structure of the baryon number susceptibilities at these high $\mu_{B}$. The transition at $\mu_{B}\approx 1.5$ GeV is due to the restoration of chiral symmetry. The transition at $\mu_{B}\approx 2.4$ GeV is due to the quark dominated phase where baryonic density is mainly contributed by quarks. The non-monotonic behavior of the fluctuation measures $\chi_{3}^{B}/\chi_{2}^{B}$ and $\chi_{4}^{B}/\chi_{2}^{B}$ reflect the transitions. In contrast to the liquid-vapor transition, those two transitions do not change the Skellam baseline. Hence, the fluctuation measures are rather small $\lesssim 1$ before and after the “transition” .

The chiral critical point of the CMF model is located at a rather low temperature $T_{\rm chiral}^{\rm CP}\approx 17$ MeV. This value is close to the critical temperature of the nuclear liquid-gas transition in the same model, the critical $\mu_{B}^{CP}$ is remarkably different, though: the appearance of the parity partners controls the dynamics of the chiral fields: as the parity partners - in the CMF model - obey the same repulsive interaction strength as the nucleons, the critical point appears at that low temperature. This phenomenon has been observed in various mean field models before.

The different phases shown in Fig. 5 in the $T-\mu_{B}$ plane are related to the chiral field $\sigma$ and the quark fraction. The chiral field is close to its vacuum value, $\sigma=\sigma_{0}$, at the hadron gas region, here the quark fraction is close to zero, as expected. Both observables deviate from their vacuum values at higher densities and temperatures only.

The chiral field drops off more slowly or $\mu_{B}=0$ than seen in lattice QCD calculations, where chiral field rapidly drops around $T=160$ MeV. The reason for this discrepancy is due to the fact that in the present CMF model the thermodynamics at these temperatures are strongly influenced by baryonic resonances which are not coupled to the chiral fields. Baryon resonances like the $\Delta$ ought to be coupled to chiral fields including their chiral partners. This brings down the chiral condensate of lower temperature as can be seen in Zschiesche et al. (2005).

The speed of sound is another important signature for the simulations of the dynamics of heavy ion collisions and neutron star mergers. $c_{s/n_{B}}^{2}$ presents a derivative $c_{s/n_{B}}^{2}=\frac{\partial P}{\partial\varepsilon}$ at constant $S/A=s/n_{B}$ entropy per baryon, that allows to estimate a speed of propagation of sound-like excitations in non-dissipative hydrodynamic evolution. The isentropic speed of sound can be calculated as Floerchinger and Martinez (2015):

[TABLE]

The partial derivatives with respect to the chemical potential and to the temperature are performed at constant temperature and at constant chemical potential, respectively. The calculated speed of sound shows three local minima which correspond to three locally softest points of the EoS. These three minima correspond to phase boundaries, where the baryon number susceptibilities present a non-monotonic behavior. Note, that the speed of sound reaches quite large values, $c_{s}^{2}\approx 0.7$, in the higher density region of nuclear matter. This high speed of sound results due to the strong repulsion between the baryons, before the onset of deconfinement. Thereafter, the vector repulsion and baryon excluded volume cease, as such terms have not been predicted for the quarks.

V Application to heavy-ion collisions

The presented EoS is used as input for hydrodynamical simulations of both heavy-ion collisions and neutron star mergers. To illustrate which regions of the phase diagram can be reached in collisions at low and moderate collision energies the stationary 1-dimensional Taub adiabat model is used Taub (1948); Stoecker et al. (1978); Thorne (1973). The expansion is described at lines of constant entropy per baryon $S/A=\,$const (isentropes). These lines depict the isentropic matter evolution of ideal fluid dynamics at different collision energies.

The entropy is produced in the earliest stage of a heavy ion collision by the shock violent compression Stoecker and Greiner (1986). During the system’s expansion there is only a moderate increase of entropy due to the rather small viscosity Csernai et al. (2006); Romatschke and Romatschke (2007), hence, an isentropic expansion scenario is a reasonable approximation Steinheimer et al. (2008).

The expansion of the equilibrated matter then continues until the system becomes so dilute that the chemical as well as the kinetic freeze-out occur and the chemical composition is fixed.

The entropy per baryon ($S/A$) is calculated in the 1-dimensional stationary scenario of central heavy ion collisions – the two colliding slabs of cold nuclear matter Baumgardt et al. (1975); Stoecker et al. (1978, 1980, 1981a, 1981b); Hahn and Stoecker (1988); Stoecker and Greiner (1986) conserve the baryon number, energy and momentum across the shock front in accord with the relativistic Rankine-Hugoniot equation (Taub adiabat), RRHT, Taub (1948); Thorne (1973). Thus, the produced entropy is directly associated to the collision energy. The thermodynamic properties across the shock front are described by the RRHT-equation

[TABLE]

where $P_{0}$, $\varepsilon_{0}$ and $n_{0}$ correspond to the initial pressure, energy density, and baryon density in the local rest frame of each of the two slabs. The two symmetric slabs consist of the nuclear matter in the ground state, $P_{0}=0,~{}\varepsilon_{0}/n_{0}-m_{N}=-16$ MeV and $n_{0}=0.16~{}{\rm fm^{-3}}$. With any known relation $P=P(\varepsilon,n)$, Eq. 15 can be solved. Furthermore, the collision energy is related to the created density as follows:

[TABLE]

Here $\gamma^{\rm CM}$ is the Lorentz gamma factor in the center of mass frame of the heavy ion collisions and $E_{\rm lab}$ is the beam energy per nucleon in the laboratory frame of a fixed target collision. This relation can be obtained from the full stopping condition Stoecker et al. (1978, 1980, 1981a, 1981b); Hahn and Stoecker (1988); Stoecker and Greiner (1986); Merdeev et al. (2011). The initial state thermodynamics (density, temperature and entropy) of the hot, dense participant matter is obtained from Eqs. 15 and 16 as a function of the collision energy. The known initial entropy yields the lines of constant entropy which give the trajectories of the heavy ion collisions in the phase diagram.

The predicted isentropic expansion trajectories are shown in the $T-\mu_{B}$ plane phase diagram in Fig. 7.

Note that 1-dimensional stationary RRHT-adiabat scenario predicts a very strong compression and heating already at intermediate lab (fixed target) bombarding energies. The heavy ion participant system crosses the weak chiral transition predicted by the present CMF model already at $E_{\rm lab}\approx 2$ A GeV, i.e. at GSI’s SIS18 accelerator facility. Here the specific total entropy is predicted to reach $S/A\approx 3$, in accord with previous RMF-calculations Hahn and Stoecker (1988) which also used the 1-D RRHT-scenario. The $T-\mu_{B}$ values, $T\approx 70$ MeV, $\mu_{B}\approx 1.2$ GeV, with net baryon densities $n_{B}/n_{0}\approx 3$, reached here in heavy ion collisions, coincide with the $T-\mu_{B}$ values reached in binary neutron star collisions, as recent general relativistic fully 3+1-dimensional hydrodynamical calculations have confirmed Hanauske et al. (2017); Most et al. (2019) for the gravitational wave event GW170817. At these temperatures and densities, $T\approx 70$ MeV and $n_{B}/n_{0}\approx 3$, the RRHT model predicts that there are abut 20% of the dense matter are already transformed to quarks.

At $E_{\rm lab}=5.6$ A GeV, $\sqrt{s_{\rm NN}}=3.5$ A GeV roughly 40% of the CMF-matter is in the quark state in the RRHT model – a prerequisite for hot quarkyonic matter. Hence, this energy, which is presumably only reachable by the BMN detector at the Nuclotron at JINR Dubna, and the FXT fixed target setup at the STAR-detector at RHIC, is of great interest: here the matter starts to be dominated by quarks, rather then by in-medium baryons, at $T>100$ MeV and $\mu_{B}\gtrsim 1.5$ GeV. This is predicted by the present CMF model, when using the 1-D RRHT ideal hydrodynamics. This model predicts that the quarkyonic transition is crossed also at higher energies, using the isentropic expansion of the matter at specific total entropy $S/A>6$. In fact, non-equilibrium viscous effects may increase the specific entropy of the system. However, pre-freezeout radiation, e.g. of Kaons and other hadrons with small scattering cross-sections which can escape early from the semi-equilibrated baryon-rich, dense system can considerably lower the specific entropy during the expansion. So an answer to the question whether the local entropy per baryon increases or decreases during the time evolution awaits more detailed microscopic/macroscopic modeling.

Hence, heavy ions fixed target experiments of SIS at FAiR and SPS at CERN as well as STAR BES program at RHIC probe temperatures from $50<T<280$ MeV and chemical potentials from $500<\mu_{B}<1700$ MeV for the collision energy range $\sqrt{s_{\rm NN}}<10$ GeV considered here. In this region the CMF model shows not an additional phase transition, but the remnants of the nuclear liquid-vapor transition at $T\approx 20$ MeV. The chiral transition at larger chemical potentials can influence the dynamical evolution, too. The present results suggest that heavy-ion collisions mostly probe regions where the nuclear matter liquid-vapor critical point dominates – hence, the observed baryon fluctuations are largely due to remnants of the nuclear liquid-vapor phase transition. This had been suggested also in previous works Fukushima (2015); Mukherjee et al. (2017); Vovchenko et al. (2017, 2018). The CP associated with the chiral symmetry restoration in the CMF model lies at $\mu_{B}\approx 1.5$ GeV and $T\approx 20$ MeV. This high density region is to the best of our knowledge is reachable in the interior of neutron stars, NS, and in binary general relativistic NS mergers Dietrich et al. (2015); Radice et al. (2017); Most et al. (2019); Bauswein et al. (2019); Hanauske et al. (2019).

Fig. 8 presents the square of the isentropic speed of sound, at fixed specific entropy, $c_{s}^{2}$ and the quark fraction as function of the temperature for the studied collision energies, i.e. it shows how respective observable quantities evolve during the cooling of the system while it expands. Those isentropic lines which belong to $\sqrt{s_{NN}}>4.5$ GeV probe the softest point of the EoS, which is attributed to the chiral symmetry restoration. At this energy region there are strong local maxima and minima of the speed of sound squared after which $c_{s}^{2}$ rapidly increases during the expansion due to the decrease of the quark fraction, as a result of the rapid appearance of baryons, the EoS stiffens quickly doe to the hard-core baryon-baryon repulsion. For collision energies $\sqrt{s_{NN}}<4.5$ GeV the initial state is not dominated by quarks, hence, the systems starts to expand at rather high values of $c_{s}^{2}$, which then monotonously decrease during the expansion, as a result of the diminishing repulsion between the baryons.

VI Application to neutron stars

The densities in the neutron star interiors also can exceed by several times the nuclear matter saturation density. At these high densities, the lack of a detailed knowledge of the equation of state and the appropriate microscopic degrees of freedom is similarly disturbing as in the relativistic heavy ion collisions discussed above. The discussion of the role of hyperonic, quarkyonic and strange quark degrees of freedom at these NS densities is ongoing.

The CMF model can be employed directly to describe the neutron star matter. Here we work out without any changes to the coupling constants and parameters used to describe the $\mu_{B}=0$ LQCD results. The temperatures isolated in neutron star interiors are negligibly small, in comparison to what we observe in heavy ions collisions and hot QCD scales. The calculations here are done in the limit $T=0$. In contrast to ordinary isospin symmetric nuclear matter, the neutron star matter is in $\beta$-equilibrium which preserves the total electric neutrality of the NS matter and locally ensures stability with respect to $\beta$-decay. As a consequence, strangeness and hypercharges assume finite non-zero values. These constraints require the presence of leptons. In addition, one must allow for quark-hadron degrees of freedom.

Fig. 9 depicts the CMF model predictions of the relative particle species number content of all different particle species present inside a CMF neutron star, at $T=0$, as function of baryo-chemical potential. One feature of the present CMF model is the absence of baryon resonances (deltas etc. and of hyperons) even though they are included in the CMF model. Their total absence in the present CMF model calculation at $T=0$ is due to the very strong hard-core repulsion by the excluded volume corrections.

The calculated EoS at $T=0$ can be used as the input for the Tolman-Oppenheimer-Volkoff (TOV) equation, which allows to relate the mass and the radius of any static, spherical, gravitationally bound object Tolman (1939); Oppenheimer and Volkoff (1939), i.e. here a static neutron star, NS. The outer layer of neutron stars presumably consist of mostly neutron rich nuclei and clusters in chemical and $\beta$- equilibrium. Those nuclei are not yet included in the CMF model. Hence, another input for the EoS of the NS crust is needed. Here, we use the classical crust-EoS Baym et al. (1971) matched to the CMF-EoS at $n_{B}\approx 0.05~{}{\rm fm^{-3}}$.

Figure 10 presents the results on the NS mass-radius relations obtained by solving the TOV equation with the present CMF EoS, matched just to that crust-EoS. The total fraction of the star’s mass which consists of light and strange quarks is presented in color-code. The most massive stable solution of the TOV equation contains only $<30\%$ deconfined quarks, i.e. for lighter NS only a small fraction of the star’s mass originates from deconfined quark matter. If the quark fraction is increased above $30\%$ the stars become unstable. The central density of the stable stars can never exceed $n_{B}=6\,n_{0}$, as shown in the lower part of Fig. 10. Here again the maximum mass indicates the “last stable star”. The continuous slow transition from NS matter to a sizable deconfined quark phase implies a smooth appearance of quarks in the star structure and prevents a “second family” of stable solutions to appear. This does prohibit a strict separation between a quark core and the hadronic interior of the star. This is a CMF result due to the Polyakov loop implementation of the deconfinement mechanism. Similar results are obtained in the Quarkyonic Matter-model, where the deconfinement is realized by the appearance of the quarks from inside of the Fermi sea while the hadrons there reside exclusively on the surface shell in momentum space McLerran and Reddy (2019). A similar approach to deconfinement was suggested in Marczenko et al. (2018), however, there the produced mass-radius diagram differs from the CMF-model Motornenko et al. (2019) due to the different realization of the chiral symmetry restoration.

The response of a neutron star to non-spherical gravitational fields is reflected in the tidal deformability coefficient $\lambda$ Hinderer et al. (2010), which depends strongly on the EoS. During the inspiral phase of binary neutron star merger, both neutron stars experience tidal deformations induced by the other respective accompanying neutron star partner. The tidal deformability $\lambda$ is a measure of the induced quadruple moment $Q_{ij}$ in a response to the external tidal field $\mathcal{E}_{ij}$:

[TABLE]

$\lambda$ is directly proportional to the second Love number $k_{2}$:

[TABLE]

For convenience, usually the dimensionless tidal deformability $\Lambda$ is presented as:

[TABLE]

Here, $M$ and $R$ are the mass and radius of the neutron star. A proper value of $\Lambda$ is important for the description of the inspiral stage during the merger of two neutron stars.

Various estimates of $\Lambda$ emerged after the detection of GW170817 by the LIGO collaboration Abbott et al. (2017). Ref Annala et al. (2018) argued that for a $1.4M_{\odot}$ neutron star the tidal deformability and star radius are constrained to $\Lambda_{1.4M_{\odot}}>120$ and $R_{1.4M_{\odot}}<13.6$ km. It was concluded by means of a Bayesian analysis that for a $1.4M_{\odot}$ star the deformability should be $375.5<\Lambda_{1.4M_{\odot}}<800$ and the radius is at $12.00<R_{1.4M_{\odot}}<13.45$ km, with respective $2\,\sigma$ confidence levels, see Ref. Most et al. (2018). A recent analysis by the LIGO and Virgo collaborations Abbott et al. (2019, 2018) provides detailed constraints, by using a Bayesian analysis based on the reproducing of the details of the gravitational wave signal.

VII Summary

A unified and consistent approach to QCD thermodynamics is presented which is appropriate for essentially all temperatures and densities relevant for both heavy ion collisions and for neutron star matter. The Chiral SU(3)-flavor parity-doublet Polyakov-loop quark-hadron mean-field model, CMF model, includes the main features of QCD hadron phenomenology, as well as a very good description of known QCD thermodynamics. The CMF model allows for a simultaneous description of many nuclear (astro-) physical data, consistent with astrophysical observations as well as heavy ion collisions and of compact stars. The CMF model is improved by the fixing of the relevant CMF parameters for the quark sector to the state-of-the-art LQCD data on the interaction measure. Here, the parameters of the Polyakov-loop potential and of the quark couplings to the chiral fields have been fixed. A good agreement is found between the CMF model predictions for LQCD data on both the baryon number susceptibilities and the “lines of constant physics”. The CMF model is used to explore the phase diagram of strongly-interacting matter at a wide range of $T$ and $\mu_{B}$. Three critical regions are found, which are connected to the nuclear liquid-vapor phase transition, to the chiral symmetry restoration, and to the deconfinement. The region of phase diagram accessible to experiments of high energy heavy-ions collisions is dominated by remnants of the nuclear liquid-vapor phase transition. Other critical regions may be probed by neutron star structure and in binary neutron star mergers. The calculated properties of neutron stars, like the mass-radius relation, the chemical composition of the stars and the tidal deformabilities are in a good agreement with recent experimental observations. The applicability of the improved CMF model to such a wide range of strongly interacting systems is impressive. For the first time, a QCD-motivated EoS is presented which precisely describes the thermodynamic observables for the whole QCD phase diagram.

Acknowledgements

The authors thank HIC for FAIR, HGS-HIRe for FAIR, BMBF, and DFG for support. J.S. appreciates the support of the SAMSON AG, WGG-Forderverein, and the C.W. Fück-Stiftungs Prize 2018. H.St. acknowledges the support through the Judah M. Eisenberg Laureatus Chair at Goethe University, and the Walter Greiner Gesellschaft, Frankfurt. Computational resources have been provided by the Center for Scientific Computing (CSC) at the J. W. Goethe-University, Frankfurt.