Hybrid and quark star matter based on a non-perturbative equation of state

With the recent dawn of the multi-messenger astronomy era a new window has opened to explore the constituents of matter and their interactions under extreme conditions. One of the pending challenges of modern physics is to probe the microscopic equation of state (EoS) of cold and dense matter via macroscopic neutron star observations such as their masses and radii. Still unanswered issues concern the detailed composition of matter in the core of neutron stars at high pressure and the possible presence of e.g. hyperons or quarks. By means of a non-perturbative functional renormalization group approach the influence of quantum and density fluctuations on the quark matter EoS in $\beta$-equilibrium is investigated within two- and three-flavor quark-meson model truncations and compared to results obtained with common mean-field approximations where important fluctuations are usually ignored. We find that they strongly impact the quark matter EoS.


I. INTRODUCTION
The recent clean mass determinations for two pulsars, PSR J1614-2230 and PSR J0348+0432 [1], confirm the existence of neutron stars (NS) with a mass of about -and possibly even beyond, see [2] -two solar masses. Since the NS mass and radius depend strongly on the underlying equation of state (EoS), see e.g. the reviews [3][4][5], any model for the NS interior should produce an EoS leading to a maximum NS mass at least compatible with the above observations. The maximum mass is thereby most sensitive to the EoS at the highest densities. The central densities of the maximum mass configuration can reach values of about 1 fm −3 , well above the nuclear saturation density, ρ 0 ∼ 0.16 fm −3 , and other degrees of freedom than neutrons, protons and electrons might appear. Although additional degrees of freedom soften the EoS and thus lower the maximum mass, it has been demonstrated that the appearance of hyperons [6], mesons, or ∆-baryons [7], is not excluded by the observation of the two massive pulsars. Numerous studies demonstrate in addition the possibility of a transition to quark matter at high density in such massive stars and the formation of so-called hybrid stars [8,9]. Under the hypothesis of absolutely stable strange quark matter [10], even pure quark stars, also referred to as strange stars, might exist [11].
The onset of a new degree of freedom causes not only a reduction of the maximum mass but in general leads to smaller radii of the stars, too. For strange quark stars, being composed of self-bound matter, the mass-radius relation is qualitatively different from neutron and hybrid stars. In particular it follows M ∝ R 3 for small masses, such that observing a pulsar with very small radius would be a strong indication for a strange quark star. Radii have presently been essentially determined from x-ray observations, see [12] for a recent review, and ex- * E-Mail:konstantin.otto@physik.uni-giessen.de † E-Mail:micaela.oertel@obspm.fr ‡ E-Mail:bernd-jochen.schaefer@theo.physik.uni-giessen.de tracted from the tidal deformability of GW170817 measured by the LIGO-Virgo collaboration [13][14][15][16][17]. The obtained values lie in the range of 10-14 km for a fiducial 1.4 M star and are perfectly compatible with neutron or hybrid stars. In the near future further precise radius determinations are to be expected from the NICER mission and additional binary NS merger data from LIGO-Virgo detectors.
From a theoretical point of view, it is not possible up to now to derive the dense matter EoS from first principles over the entire range necessary for describing compact stars, neither on the hadronic nor on the quark side. Therefore simplified EoS at finite densities need to be constructed. Currently two main strategies are pursued. The first one consists in parameterizing the EoS in the unknown density domains, putting the least possible number of model assumptions. The parameters are then adjusted to existing constraints from nuclear experiments, observations and/or theoretical calculations, including attempts to extract the density dependence of the EoS directly from NS mass and radius data. Examples are the EoS by [18,19]. The second, more traditional strategy is based on modelling dense matter. It is less flexible than the aforementioned but has the advantage of allowing to track among others the matter's particle content. Although some points remain open, for instance on inhomogeneous matter in the crust, decades of considerable effort together with constraints from experimental data and theoretical calculations have led to sophisticated and reliable models for nuclear matter up to roughly twice the saturation density. Above this density, not only the models are less under control, but non-nucleonic degrees of freedom might appear and the situation becomes more complicated.
Hybrid EoS thereby include both hadronic as well as quark matter and are usually obtained by a combination of the corresponding EoS for both sides. The quark matter part is subject to even larger uncertainties than the hadronic one, being per se a non-perturbative QCD issue. It is still under development in particular at high baryonic densities and small temperatures appropriate for NS. Examples are based on perturbative QCD [20], density-dependent quark-mass models [21], NJL-models [22], quark-meson model investigations [23], Dyson-Schwinger type approaches [24], holographic [25], or quasiparticle models [26]. Lattice QCD calculations, which have allowed to considerably improve our understanding of the EoS at low baryon chemical potentials, are afflicted by the sign problem at finite density and can thus not directly be applied to neutron stars. For the moment some models include constraints from lattice results at low chemical potentials, see e.g. [27] or use QCD-like theories not subject to the sign problem, for instance G2-QCD [28], but much effort is still needed.
In this work an EoS for quark matter within a model based on the underlying chiral symmetry breaking of QCD is derived and the consequences for the structure of a non-rotating star are analyzed. In contrast to many previous works using the mean-field approximation, the functional renormalization group (FRG) approach is employed here to incorporate nonperturbative effects. An accessory benefit of the FRG is that quantum and thermal fluctuations are taken into account that are of particular importance in the vicinity of phase transitions and are usually ignored in mean-field approximations. These might be some reasons for the recent growing interest in the application of the FRG method to neutron star matter [29,30].
A phase transition from hadronic matter to quark matter can be modeled which might be of relevance for hybrid stars. In that context, strange quarks might be suppressed due to their relatively large effective mass [22,[31][32][33]. In addition, large effective strange quark masses often render hybrid stars containing a strange quark matter core unstable against gravitational collapse [22,32], such that we consider the two-and three-flavor cases separately here.
With our approach it becomes feasible to investigate the influence of quantum as well as density fluctuations on observables for neutron stars such as the mass-radius relation in a systematic manner. We will consider here non-magnetized matter, since magnetic effects on the EoS are expected to play a minor role for pulsars [34].
The paper is organized as follows: after a brief setup of the used effective quark-meson model for quark matter, in Sec. II, three different approximations of the grand potential are presented which incorporate various contributions of certain quantum, thermal and density fluctuations. The numerical results of the phase structure are summarized in Sec. III. As input for the analysis of the mass-radius relation of a neutron star the EoS is needed. In Sec. IV various EoS for symmetric as well as for β-equilibrated and charge neutral matter are calculated and compared with a non-perturbative EoS obtained with the FRG with and without strange quarks. In addition, a phase transition from hadronic to quark matter is implemented and the parameter independent consequences are studied. The speed of sound in quark matter for various approximations is analyzed in Sec. IV A. The resulting mass-radius relations for all the used approximations are presented in Sec. IV B. We conclude in Sec. V and summarize our findings. Parameter choices and numerical details can be found in the Appendices A and B.

II. CONSTRUCTING A NON-PERTURBATIVE EOS
As mentioned above, we will be mainly concerned with determining the EoS of the quark core of a hybrid star or that of a strange star. To that end, we will consider a quark part and a leptonic part. The latter will be treated as a relativistic non-interacting Fermi gas. In the following, we focus on the quark part which will be calculated within an effective twoand three-quark flavor quark-meson model framework and we will investigate different approximations. The quark-meson model has been widely used as an effective model of lowenergy QCD since it successfully incorporates dynamical chiral symmetry breaking and the generation of constituent quark masses. In contrast to NJL-type models often employed in this context, the quark-meson model used here does not implement a repulsive vector interaction channel which would stiffen the EoS [35]. We confront two different mean-field approximations of the effective potential with the results obtained within the functional renormalization group. This comparison enables a systematic and parameter independent analysis of the influence of certain quantum and density fluctuations on the EoS for vanishing temperature.

A. Quark-meson Models
The quark-meson model consists of N f flavors of constituent quarks q and dynamical (pseudo-)scalar meson fields φ a encoded in the meson matrix with φ a := σ a + iπ a and the generators T a of the U (N f ) group transformations. The model features besides the usual kinetic terms for all involved fields a Yukawa-type interaction between the quarks and mesons as well as mesonic selfinteractions through the chirally invariant potential which is in general a function of N f independent chiral invariants For two light and one heavy quark flavors, N f = 2 + 1, the generators of the corresponding U (3) transformations in flavor space can be chosen as the usual Gell-Mann matricesλ a , i.e. T a =λ a /2. Omitting in the effective potential the highest chiral invariant ρ 3 whose associated coupling constant is of negative mass dimension, the three-quark flavor Euclidean Lagrangian reads wherein additionally the lowest order axial U (1) A symmetry breaking term and two explicit chiral symmetry breaking terms −c l σ l −c s σ s have been added. We assume here a perfect SU (2) isospin symmetry, such that the two lightest quark flavors up and down can be replaced by one index l = u = d while the strange quark flavor is denoted by the index s. The relation of the singlet-octet basis (σ 0 , σ 8 ) and the nonstrange-strange basis (σ l , σ s ) in the scalar meson sector is governed by a rotation such that the isospin-symmetric vacuum condensate evaluates to For only two quark flavors, N f = 2, the model simplifies to in the mesonic chiral effective potential U with one scalar field σ and three pseudoscalar pions π summarized in ϕ T := (σ, π). Implicitly a maximal U (1) A -symmetry breaking is assumed in this representation since the remaining chiral (pseudo)scalar multiplets, the η and a fields, are neglected. For more details see e.g. [36]. The grand partition function Z in thermal equilibrium is defined by a path integral over the quark/antiquark and meson fields wherein the temperature is introduced via the Matsubara formalism. The Euclidean Lagrangian generally contains N f independent quark chemical potentials µ f which are added to the quark-meson Lagrangian in standard thermodynamic manner The three quark flavor chemical potentials are not independent. Since we are interested in cold neutron star matter, we assume weak equilibrium including β-equilibrium with neutrinos freely leaving the star, where µ denotes the quark chemical potential related to baryon number, µ = µ B /3, and µ e the electron chemical potential which in the present case is the negative charge chemical potential. In addition, electrical charge neutrality has to be fulfilled, such that only one independent chemical potential remains. We choose µ as such. Note that even though the chemical potentials introduced in Eq. (11) break isospin symmetry, we still assume only one light condensate σ l as an approximation. This leads to equal masses m u = m d = m l even in isospin asymmetric matter. Finally, the logarithm of the grand partition function yields the total grand potential which in general incorporates the thermal, density and quantum fluctuations of the quark and meson fields

B. Mean-field Approximation
In mean-field approximation (MFA) the grand potential for the quark-meson model basically splits into a quark loop and a static meson contribution For N f = 2 + 1 quark flavors the meson contribution Ω m , ignoring the ρ 3 invariant, reads with the chirally symmetric meson potential U χ (ρ 1 , ρ 2 ) = m 2 ρ 1 + λ 1 ρ 2 1 + λ 2 ρ 2 introducing in this manner a mass parameter m 2 and two quartic couplings λ i .
For only two quark flavors the chiral potential simplifies to The employed input parameters can be found in Appendix A. The integration over the quark loop yields a UV finite and explicitly temperature dependent contribution with the usual Fermi-Dirac distribution functions , the quark energies E f = p 2 + m 2 f and corresponding light, m l = gσ l /2, and strange quark masses m s = gσ s / √ 2. The light and the strange condensates σ l and σ s are the temperature and quark chemical potential dependent minima of the full thermodynamic potential. Again, for only two quark flavors the strange quark contribution in this expression is simply dropped and the light quark mass is replaced by m l = gσ/2 with one chiral order parameter σ.
In general, the UV divergent vacuum contribution of the quark loop to the potential can be completely absorbed in the model parameters because the quark-meson model is renormalizable. Ignoring this divergent zero point part yields the standard MFA (sMFA). However, this vacuum contribution is automatically included in the fully renormalized quark flow equation, In this way important vacuum fluctuations in addition to the thermal and density fluctuations from the quark loops are included in the grand potential [37]. In the following, we denote this approach by renormalized MFA (rMFA). In both meanfield approximations the fluctuations and back-reactions on the mesonic sector are completely left out and the same static (tree-level) meson potential is used. However, all quark and meson fluctuations can finally be considered by applying the functional renormalization group method.

C. Functional Renormalization Group
As mentioned in the introduction a suitable framework to incorporate quantum fluctuations in a consistent way is based on the non-perturbative functional renormalization group in terms of the Wetterich equation [38] with the RG time t = ln(k/Λ). The effective average action Γ k interpolates between a microscopic or bare UV action S Λ = Γ k→Λ and the full quantum effective action Γ = Γ k→0 in the infrared and governs the dynamics of the field expectation values after the integration of quantum fluctuations from the UV scale Λ down to the infrared scale k IR . The infrared regulator R k specifies the regularization of quantum fluctuations near an infrared momentum shell with momentum k and Γ (2) k denotes the second functional derivative of the effective average action with respect to the fields of the given theory. The trace represents a summation and integration over all discrete and continuous indices. The highly non-linear Wetterich equation has a simple one-loop structure but includes higher loop contributions in perturbation theory since the full nonperturbative propagator enters the loop diagram. One advantage of this approach is that it does not rely on the existence of a small expansion parameter and thus is applicable in the nonperturbative regime. QCD-related reviews on the functional RG approach can be found in [39].
In order to solve the functional equation numerically some truncations are required that turn it into a finite-dimensional partial differential equation. This truncation might induce a certain dependence of physical observables on the regulator, but this can be minimized by choosing optimized regulators or by implementation of RG consistency [40]. In this work, a modified three-dimensional flat regulator has been used [41].
The flow equation for Γ k must be supplemented with an initial condition at k = Λ which according to Eq. (4) reads for three flavors withρ 2 := ρ 2 − ρ 2 1 /3. This truncation for the effective action corresponds to a leading-order derivative expansion with standard kinetic terms for the meson fields. Note that in this local potential approximation (LPA) no scalar wave function renormalizations and no scale-dependence in the Yukawa couplings between quarks and mesons are taken into account. However, in this truncation the important dynamical back-reaction of the mesons on the quark sector of the model is already included. For more details see [36].
Finally, this yields for three quark flavors the IR and UV finite flow equation for the effective potential where now the flow of the mesonic degrees of freedom is taken into account. The meson energies include the RG scale dependent meson masses m b which are obtained by diagonalizing the mass entries of the matrix Details and the lengthy explicit expressions of the eigenvalues can be found in [36]. Evolving the system towards the infrared yields the full thermodynamic potential evaluated at the solution of the gap equation, i.e., the minimum of the grand potential.
As initial UV condition for the flow Eq. (21) the meson potential is parameterized as with the scale-dependent chiral potential which differs only in the second argumentρ 2 from the corresponding chiral potential in mean-field approximation Eq. (15). Note that only the expansion coefficients a ij in the potential are scale dependent while all remaining parameters are kept constant. See App. A for the initial parameter fixing. For two quark flavors the flow, Eq. (21), simplifies to with only one scalar σ and three pion degrees of freedom. Moreover, for T = 0 and N f generic quark flavors the flow reduces to (26) In this limit the Fermi-Dirac distributions of the fermionic threshold functions in Eq. (21) or in Eq. (25) become a sharp Heaviside function such that the Silver Blaze property is satisfied [42]. Hence, for µ 2 f > m 2 f , only scales above the Fermi sea f contribute to the corresponding quark loop and are integrated out yielding a finite quark density. Hence, increasing the chemical potential suppresses more and more the quark dynamics of the model.
Note that at vanishing temperature the quark-meson model is equivalent to its Polyakov-loop extended version, the Polyakov-quark-meson (PQM) model [43]. In such PQM models the deconfinement phase transition is captured statistically by including an effective potential for the gluon background field in terms of the order parameter fields for deconfinement, the Polyakov loops. There are basically two major modifications by the Polyakov loops which vanish in the zero temperature limit exactly: all known variants of the effective Polyakov loop potential [44][45][46][47] are proportional to the temperature and the implicit dependence of the Polyakov-loop variables on the quark loop dynamics degenerate to a standard Fermi-Dirac distribution without a Polyakov-loop contribution. Hence, exactly at T = 0 the Polyakov-loop in this context is irrelevant for the thermodynamics at µ > 0. However, a phenomenological finite density generalization of the Polyakov loop potential at zero temperature can stiffen the EoS [48]. For a recent review see e.g. [49] and for recent (P)QM phase structure investigations with the FRG see e.g. [50].

III. PHASE STRUCTURE
The bulk thermodynamics of the system is determined by the effective potential which depends in the case of two flavors on one (non-strange) condensate and for three or, more precisely, (2+1) flavors on two (one non-strange and one strange) condensates. The condensates are determined by minimizing the total effective potential with respect to the corresponding fields, here generically denoted by Φ, This yields the chiral non-strange σ l and strange σ s condensates or expectation values as functions of the external parameters, i.e. the temperature and the quark chemical potentials. The thermodynamic pressure is just the negative of the effective potential evaluated at the minimum and normalized in vacuum In the FRG case, the infrared-evolved potential Ω k=0 − Ω 0 k=0 has to be used. By varying the external parameters the phase structure of the chiral transition can now be analyzed. For simplicity, we consider here a common chemical potential for all quark flavors, i.e., we set µ e = 0 in Eq. (11). Our numerical findings are collected in Fig. 1 where the phase diagrams for two quark flavors (solid lines) and for three flavors (dashed lines) are shown. The numerical input parameters and contingent cutoff dependency are collected in App. A; details on the numerical implementation can be found in App. B.
Typical for the phase structure obtained with the FRG is the back-bending of the chiral first-order phase transition line (dark blue lines) for small temperatures below the critical endpoints (CEPs) denoted as black dots in the figure. For the chosen vacuum input parameters and in the LPA truncation of the FRG equation, see App. A, the CEP is located at very small temperatures, around T ∼ 10 MeV for two and three quark flavors. In mean-field approximations the thermodynamical behavior at small temperatures is different and the first-order transition line hits the chemical potential axis perpendicularly [51].
Furthermore, since the inclusion of fluctuations generically washes out the chiral phase transition, the crossover transition line is shifted to higher temperatures if more fluctuations in the thermodynamic potential are taken into account. This is nicely demonstrated in Fig. 2(a) where both the light (solid lines) and the strange (dashed lines) chiral condensates are shown as a function of temperature for vanishing quark chemical potential.
In sMFA where only the thermal quark loop contribution is considered, the pseudocritical crossover temperature T c ≈ 140 MeV at µ = 0 is smallest. Already the inclusion of the vacuum quantum fluctuations of the quarks, labeled as rMFA in the figures, lifts the pseudocritical temperature by about 30 MeV. Interestingly, the whole chiral phase transition is shifted constantly towards higher temperatures by roughly this amount for all chemical potentials, cf. Fig. 1. As a consequence, the first-order transition at T = 0 is also pushed to higher chemical potentials as visible in Fig. 2(b). This trend is continued at least for moderate chemical potentials when additionally meson fluctuations with the FRG are taken into account. However, for smaller temperatures and due to the back-bending of the transition line, cf. Fig. 1, the critical chemical potential is pushed to smaller values which will be of relevance for the EoS later.
All condensates exhibit for T = 0 a first-order phase transition close to µ ≈ 300 MeV corresponding to the light quark mass in the vacuum. In sMFA, the first-order transition is strongest, in FRG weakest. Hence, the gap in the FRG light condensate is quite small and melts only moderately after the transition, still signaling a chirally broken phase in this density regime of the phase diagram [52]. In rMFA the light condensate is constant until µ = 300 MeV and melts down before the first-order jump which is consistent with the Silver Blaze property. It is likely that for a sigma mass below 560 MeV the rMFA condensates immediately jump at the light quark masses as well.
Just below a quark chemical potential of about 430 MeV, the value that coincides with the strange quark mass in the vacuum, a further decrease is seen in all three strange condensates and a smooth chiral phase transition takes place.
However, when strange quarks are added to the system an opposite behavior is found for vanishing and moderate chemical potentials and the three-flavor crossover line is pushed down again, see the dashed line in Fig. 1. The difference to the two-flavor phase structure shrinks for decreasing temperatures. Below T < 50 MeV almost no influence of the strange quark on the transition line is observed where the dashed line merges with the solid two-quark flavor line.

IV. EOS FOR QUARK AND HYBRID STARS
The equations of state (EoS) for symmetric quark matter, i.e., for equal chemical potentials, obtained in MFAs and with the FRG, are compared to each other in Fig. 3(a). Solid lines are the two-quark flavor findings and the dashed lines the corresponding three-flavor calculations. The numerical results are almost insensitive to the strange quark before the onset of the strange chiral phase transition around energy densities ε ≈ 550 MeV/fm 3 but start to deviate thereafter, cf. the inlay in Fig. 3(a). In MFAs the transition is more gradually realized and the deviation is less pronounced than within the FRG. Furthermore, it is obvious that vacuum fluctuations reduce the slope of the EoS, i.e. the sound speed, and over most of the shown density range the EoS obtained in FRG has still smaller slope, see section IV A.
With the inclusion of a free relativistic electron gas and the conditions for weak equilibrium and charge neutrality, cf. Eqs. (11) and (12), we obtain slightly modified EoS. The results for β-stable and neutral matter are presented in Fig. 3(b). In comparison to symmetric quark matter that does not satisfy the above conditions, pressure is reduced at given energy density. The effect is, however, negligible in mean-field approximation for N f = 2. In the three-flavor case, the transition is pushed to smaller energy densities. Note that by using only one chiral condensate σ l for both up and down quark flavors a certain mismatch is caused.
Due to the splitting of the chemical potentials, the isospin symmetry is broken, but our approximation yields in all cases degenerated up-and down-quark masses. The mismatch is largest in the vicinity of the chiral phase transition. For large µ, the restoration of chiral symmetry in the light quark sector suppresses both quark masses such that only a mass differences which is induced by the explicit symmetry breaking terms is left. However, for a more complete analysis including isospin symmetry breaking the introduction of an additional third chiral condensate is necessary which is beyond the scope of the present work. In order to allow for a description of hybrid stars with a phase transition from hadronic to quark matter in the interior of the star we combine the quark matter EoS with a nuclear one. The transition is achieved with a standard Maxwell construction that maximizes the pressure for a given chemical potential. For the nuclear EoS we consider some representative models compatible with several nuclear physics constraints as well as the maximum neutron star mass and the GW170817 tidal deformability. Three of them are energy density functional models, one is based on a non-relativistic Skyrme parameterization, RG(SLy4) [53], and two are relativistic mean field models, HS(DD2) [54] and SFHo [55]. The BL EoS [56] is formulated in the framework of the Brueckner-Bethe-Goldstone many-body theory with chiral nuclear forces.
In Fig. 4 a comparison of different nuclear EoS (dash-dotted lines) with the N f = 2 (solid) and N f = 2 + 1 (dashed) EoSs evaluated with the FRG respecting β-equilibrium and charge neutrality is given. Obviously, all nuclear EoS except the HS(DD2) EoS produce higher pressure at a given baryon chemical potential (µ B = 3µ) for the entire range of interest for compact stars. Hence, no hybrid stars could exist with these model combinations. The pressure of the HS(DD2) EoS intersects the two-flavor FRG pressure curve twice. The intersection for a higher pressure around µ B /m n ≈ 1.45 displays the appropriate physical transition from nuclear to quark matter while the lower intersection point has been dropped because it is unexpected to find a quark matter phase at such low densities. Nevertheless, its existence can be explained by the attractive meson interactions in the QM model which, due to their binding energy, lead to a first-order transition at a baryon chemical potential smaller than the neutron mass.
By construction, the combination of the two-flavor QM EoS from the FRG with the HS(DD2) EoS leads to a first-order phase transition between the confined nuclear matter and the deconfined quark matter which is characterized by a discontinuity in the energy density. This can be observed in Fig. 5 which depicts the constructed hybrid N f = 2 (DD2+QM2) and N f = 2 + 1 (DD2+QM2+1) EoSs in comparison to the hybrid EoS QHC19 [57] and a parameterized EoS [58]. 1 The QHC19 EoS features a smooth crossover quark-hadron tran- sition. For the parameterized EoS the HS(DD2) EoS is used for the hadronic regime and the quark matter side is parameterized as which describes a first-order transition at {ε c , p c } with energy gap ∆ε and a constant slope s = ∂p/∂ε, i.e. a constant quark matter sound speed squared, thereafter. Following Ref. [58], we choose p c = 1.89 × 10 35 dyn cm −2 and ε c = 9.02 × 10 14 g cm −3 which corresponds to n c ≈ 3n 0 and ∆ε/ε c = 0.6. This ensures a phase transition at similar densities to those found in the employed QM model construction together with a large gap in energy density which is a necessary criterion for the possible occurrence of twin stars, i.e. an additional branch of hybrid stars with the same mass but different radius than their nuclear counterpart [60]. A comparison to the FRG calculation in Fig. 5 reveals that the gaps in our construction are much too small to produce a disconnected second branch but might favor a single connected hybrid-nuclear branch. As explained in [60], if ∆ε becomes too large, there might also not be any stable hybrid stars at all.
For the slope, we consider two extreme parameterizations, one with s = 1/3 corresponding to the asymptotical QCD value, and one with s = 1 corresponding to the maximally allowed sound speed by causality.

A. Sound speed
More information on dense matter can be gained through a detailed investigation of the speed of sound [61]. It measures the stiffness of the EoS for a one fluid flow by the thermodynamic derivative of the pressure with respect to the energy density at constant entropy and particle numbers and can be identified as the speed of propagation of sound waves. Causality implies an upper bound c 2 s ≤ 1 and thermodynamic stability a lower bound c 2 s > 0. For an ideal gas composed of point-like ultrarelativistic (massless) components the squared speed of sound is equal to one third, c 2 s = 1/3. This is common to all systems with conformal symmetry of which an ideal massless gas is just an example. Even for any strongly interacting system the vanishing of the trace of the energymomentum tensor, a feature of conformal theories, implies that the energy density is connected to the pressure by ε = 3p hence yielding c 2 s = 1/3 independently of density, temperature, or interactions. The speed of sound is decreased such that c 2 s < 1/3 when a mass for the components is included or when (perturbative) interactions among the components take place. In the case of QCD at asymptotically high densities or temperatures, far exceeding the densities in the core of compact stars, a weak coupling expansion is valid (pQCD) such that c 2 s is expected to reach the conformal limit with increasing density from below [62]. This behavior is confirmed in QCD lattice calculations at finite temperature as well as at zero and small baryon chemical potentials [63].
The speed of sound has also been investigated in alternative theories for which for example the AdS/CFT correspondence holds where calculations in the strong coupling limit are feasible, see e.g. [64]. It has been conjectured that c 2 s is always bounded from above in such classes of strongly coupled field theories by the conformal value of 1/3 [65] although recently counterexamples have been presented [66]. For more details of this conjecture see [62].
The speed of sound of the QM model in both mean-field approximations and the FRG calculation is found generally to be always smaller than c 2 s = 1/3. An alternative scenario could be the presence of a bump in c 2 s at intermediate densities before approaching the upper bound from below asymptotically and thus implying the existence of a maximum and a local minimum of c 2 s as a function of the chemical potential. This scenario is supported by another recent FRG analysis including diquark condensation [29] where a maximum in c 2 s above 1/3 is found. The additional inclusion of vector interactions in the quark-meson model [67] is also expected to stiffen the EoS.
Our result for the speed of sound of quark matter with a flavor-symmetric chemical potential is shown in Fig. 6(a). In the N f = 2 mean-field approximation (solid lines) the speed of sound converges to the limit c 2 s = 1/3 while the addition of strange quarks (dashed lines) leads to a reduction of c 2 s around scales of the strange chiral phase transition. In the FRG solution, the speed of sound is generally smaller than the asymptotic mean-field values beyond the transition.
Furthermore, the strength of the first-order chiral phase transition, i.e. the gap in the order parameter, is found to correlate with the size of the jump in the speed of sound. Hence, the strong first-order transition in sMFA leads to a jump of c 2 s close to its asymptotic Stefan-Boltzmann value c 2 s = 1/3 which leads to the almost linear behavior of the EoS even at low pressure, see Fig. 3(a). The more washed-out transition in rMFA induces an initially smaller slope of the EoS. This becomes more significant in the FRG calculation: due to an even smoother transition a comparably small gap ∆ε is found in Fig. 3(a). In agreement with Fig. 6(a) the slope is consistently smaller than that of the mean-field calculations.
For N f = 2 + 1, c 2 s is found to be sensitive to the numerical error caused by the employed solution method of the flow equation which leads to visible fluctuations at high µ. Therefore, Fig. 6(a) depicts for µ > 350 MeV averaged values in conjunction with error intervals displayed as blue band. For more technical details see App. B.
The speed of sound for β-stable and charge-neutral quark matter is displayed in Fig. 6(b). Note that due to the usage of only one light condensate, the first-order transition cannot be resolved exactly in this approximation and hence the drop of the speed of sound to zero at low chemical potential is not shown. Due to the numerical uncertainties mentioned above, we postpone a careful N f = 2 + 1 analysis to a future work. Qualitatively, we observe the same behavior as for a flavor-symmetric chemical potential. However, in mean-field approximation we find that the reduction of the speed of sound due to the onset of strangeness takes place at lower chemical potentials and more gradually than in symmetric quark matter. In Fig. 6(b) we also show the speed of sound for the HS(DD2) nuclear EoS and indicate the transition point to quark matter in the DD2+QM2 EoS by a vertical line. As suggested in [68], this discontinuity in the speed of sound can be related to a δfunction singularity in the fundamental derivative, leading to possibly non-convex thermodynamics.

B. Neutron star models
In order to determine the influence of fluctuations on the mass-radius relation of a neutron star, we employ hydrostatic equilibrium solutions for a relativistic spherically symmetric compact star composed of a perfect fluid, which have been derived from Einstein's equations by Tolman, Oppenheimer and Volkoff (TOV) [69]. In Schwarzschild coordinates the corresponding equations determining the pressure p and the enclosed gravitational mass M of the star as function of the radius r read 2 dp(r) dr = − G r p(r) + ε(r) M (r) + 4πr 3 p(r) r − 2GM (r) wherein ε denotes the energy density and G the gravitational constant. For a given EoS in terms of p(ε) as input these equations can be integrated from the origin for a given choice of a central pressure p 0 , i.e. p(r = 0) = p 0 and M (r = 0) = 0.
The value of the radius where the pressure vanishes, i.e. p(r = R) = 0, defines the surface and thus mass M (R) and radius R of the star. Varying the unknown central pressure p 0 yields the mass-radius relation. The mass-radius relations for the pure β-stable and chargeneutral quark matter EoS are shown in Fig. 7(a) for three different approximations. Please keep in mind that quark matter is not absolutely stable within our setup and that thus such pure quark stars could not exist. We, nevertheless, show the In general, all three-flavor calculations yield a smaller maximum mass than the corresponding two-flavor results, which can be understood by the softening of the EoS due to the additional strange degrees of freedom. Only the two-flavor sMFA and both the two-and three-flavor FRG results yield a maximum mass above 2M . Furthermore, the inclusion of the renormalized vacuum fluctuations in the rMFA in contrast to the sMFA leads to smaller masses and radii. The additional consideration of mesonic fluctuations via the full FRG computation increases the maximum mass even slightly beyond the sMFA result but also leads to significantly larger radii. For the sake of completeness the causality constraint R ≤ 2.87GM [4] is also displayed in the figure.
The mass-radius relations from the combined EoSs for a hybrid star with the N f = 2 or N f = 2 + 1 quark-meson matter side employing the FRG, respectively, and a hadronic phase parameterized by the HS(DD2) EoS are presented in Fig. 7(b), labeled again as DD2+QM2 and DD2+QM2+1. For N f = 2, the onset of quark matter leads to the visible kink in the DD2+QM2 curve slightly below M = 2M corresponding to a central baryon number density of approximately 0.47 fm −3 . Below this value, the hybrid star mass-radiusrelation coincides with the nuclear HS(DD2) one as it should. For N f = 2 + 1, the DD2+QM2+1 curve exhibits a similar behavior, but the onset of quark matter occurs at a smaller baryon number density of approximately 0.43 fm −3 . Thus, since the transition occurs well above nuclear saturation density, both hybrid EoSs satisfy constraints from nuclear physics as implemented in the HS(DD2) EoS. The maximum hybrid star mass of about 2.1M for N f = 2 complies well with current observations, whereas the N f = 2 + 1 curve does not satisfy the 2M limit. Since the quark matter onset occurs only for masses slightly below 1.8M and higher, the value of the GW170817 tidal deformability obtained from both hybrid EoS does not change significantly with respect to the HS(DD2) valueΛ ≈ 795 for a mass ratio of 0.8 of the two coalescing stars. It is in slight tension with recent LIGO/Virgo data [13,14] but in agreement with the observation that the HS(DD2) EoS leads to a relatively large radius for intermediate mass stars. The QHC19 model leads to a radius smaller by almost 2 km. For comparison, the M -R relations for four different pure nuclear RMF models (dash-dotted lines) are also shown in Fig. 7(b). As mentioned above, since the stiffness of the quark matter EoS obtained with the FRG is not large enough a hybrid star EoS with these nuclear EoS is not possible within our approximation.
As expected from Fig. 5, the parameterized EoS leads to a kink in the mass-radius relation at a mass slightly above the DD2+QM2 curve. Contrary to the connected hadronic and hybrid branches in the latter, the large energy gap of the s = 1 parameterized EoS leads to a disconnected hybrid branch and therefore twin stars at masses of about 2M . For s = 1/3, the pressure in the quark matter phase is not sufficient to counteract the strong gravitational pull due to the large energy density of the quark core, cf. Ref. [58], and thus does not support a stable hybrid star branch. Hence, we can rule out the occurrence of twin stars in our model due to the small energy gap at the phase transition from nuclear matter to quark matter and due to the small stiffness of the quark matter EoS. In case of a Maxwell construction with parameters that feature a larger energy gap, there might not even be any stable hybrid stars with a QM model quark matter description.
In summary we found that it is feasible to construct a hybrid nuclear-quark EoS where the quark matter part is obtained from an effective theory within the functional renormalization group approach. In particular, we would expect that the inclusion of repulsive interaction channels suspected to play a significant role might stiffen the quark matter EoS to a degree that allows for a realistic description of combined hadronic and quark matter with other hadronic models. Furthermore, it might deepen our understanding of the role and possible abundance of strange matter in compact stars.

V. SUMMARY AND CONCLUSIONS
The core of neutron stars contains strongly interacting matter at extreme densities, reaching several times nuclear matter saturation density for the most massive ones. Observations, in particular precise high mass determinations [1] and the GW170817 tidal deformability [13,14], thereby put constraints on the equation of state for conditions not accessible to experiments. Although much recent progress has been achieved, many questions remain open, in particular on the composition of the inner core: does it contain nonnucleonic degrees of freedom? In addition to hyperons, nuclear resonances or mesons, a quark matter core might appear. In this paper we have performed a non-perturbative study of the quark matter EoS. We have considered a twoand three flavor quark-meson model, employing different approximations from mean-field to the functional renormalization group. Apart from including fluctuations in a consistent non-perturbative way, the FRG has the additional benefit of reducing the parameter dependence of the usually applied EoS models.
The quark-meson model fully incorporates chiral symmetry breaking and we have performed a first study of the impact of quantum and density fluctuations on the EoS for vanishing temperature. Since the different approximations of the grand potential were fixed to the same input parameters, our numerical findings are solely attributed to the impact of the fluctuations. As anticipated from studies of the phase diagram and confirmed by our investigations, fluctuations tend to wash out the chiral phase transition. Within the EoS, the softening due to the appearance of strange quarks is therefore pushed to higher densities. Quark stars obtained from the FRG, including fluctuations in the quark matter EoS, have higher maximum masses and radii compared with their mean field counterparts. Furthermore, we have constructed a hybrid star EoS, combining our FRG quark matter EoS with a nuclear one via a Maxwell construction. The results for hybrid stars with a two-flavor quark matter core are in reasonable agreement with existing constraints. In contrast to many studies within the mean-field NJL model, see e.g. [32], our FRG EoS allows for gravitationally stable hybrid stars with a three-flavor core, which however leads to a maximum mass below the highest observed pulsar masses. We notice that the inclusion of a repulsive vector interaction in the quark-meson model is susceptible to additionally stiffen the quark matter EoS and allow for constructing hybrid stars with nuclear EoS with smaller radii and tidal deformabilities.
We have presented here one of the first constructions of a non-perturbative EoS for high densities within the FRG that is a very promising non-perturbative method for computing the EoS directly from quark-gluon degrees of freedom, cf. [29,30]. Please note that for small temperatures, subject of the aforementioned analysis, dynamical baryonic correlations, i.e. the dynamical interrelation of three-quark and baryonic degrees of freedom, become an increasingly important issue. In particular, for densities relevant for the neutron star interior such correlations are likely to be of utmost importance and certainly modify the physical findings. Generally, the improvement of the FRG calculations towards a broader density regime can be achieved in a systematic manner. The elaboration of the clustering of quarks into baryons as well as the emergence of long-range correlations between baryons makes the QCD-based FRG approach towards lower densities increasingly auspicious and also sophisticated. First successes could in parts already be achieved in a similar QCD related context [70]. Due to the expected rich phase structure in this area of the QCD phase diagram a related and important aspect addresses the inclusion of further additional interaction channels like, for example, vector-axial-vector and/or diquark-quark channels which allow for further phases such as e.g. 2SC or CFL phases [58]. Work along these lines is ongoing.
MeV. For N f = 2, all strange quantities are omitted. g, c l and the remaining two input parameters for the chiral potential are set to reproduce the subset m l , m π , m σ and σ k,0 in the same fashion as above.
The parameter set for the FRG and rMFA flow equations have been optimized by a global differential evolution algorithm [71] with an initial UV cutoff of Λ = 1 GeV. In the full FRG case we stopped the evolution around k IR = 80 MeV where the condensates are already frozen and in the rMFA case we stop at k IR = 1 MeV. Note that all obtained numerical results are insensitive to IR values when chosen in this region while a UV cutoff dependence for the rMFA results can still be seen. However, when choosing a UV cutoff larger than Λ > 2 GeV for the rMFA results any cutoff dependence disappears [72].
The used input parameters for the mean-field potentials can be found in Tab. I. The UV coefficients a ij for the chiral potentials in the FRG calculations are listed in Tab. II. Further details on the input parameters can be found in Refs. [36].  68   TABLE II. FRG input parameters similar to Tab. I. Additionally, the utilized Yukawa coupling g and sigma mass mσ are also quoted. The parameter a01 corresponds to the modified invariantρ2 by analogy with the mean-field parameter λ2 in Tab. I.

Appendix B: Numerical Solution Techniques
The numerical methods employed for the solution of the FRG flow equations and for the TOV equations are outlined in this appendix.
The FRG flow equations (21) and (25) are partial differential equations (PDEs) including a partial derivative w.r.t. the RG scale k and field derivatives encoded in the mass terms. Setting up a two-dimensional grid for the chiral potential in the variables x = σ 2 l and y = 2σ 2 s − σ 2 l , cf. [36], and interpolating the derivatives from the discrete grid points, it is possible to reduce the PDE to a set of coupled ordinary differential equations (ODEs). For the interpolation, we use cubic splines in each of the two grid directions, respectively. The ODE solution is obtained from an explicit Runge-Kutta type step algorithm. For N f = 2, the grid reduces to one dimension.
In the context of this work, the T = 0 limit of the flow equation, Eq. (26), is utilized wherein the fermionic threshold function reduces to a Heaviside function. This also encodes the Silver Blaze property of the theory because E l ≥ m l = gσ l,0 /2 (B2) implies that θ(E f − µ) = 1 for all µ < m f and hence the flow at the IR minimum σ l,0 does not change with respect to the vacuum flow. Here we have assumed for simplicity a flavor-independent chemical potential. Of course, this property only holds for µ < µ c where µ c signifies the chiral firstorder transition. In both the sMFA and FRG solutions, we observe µ c < m l , cf. Fig. 2(b). Unfortunately, due to the utilized grid method the Silver Blaze property is subject to a numerical error. For any chemical potential, all grid points located at σ l < 2µ/g display a different running than in vacuum. Since the bosonic energies incorporate field derivatives that are approximated from an interpolation of all grid points, the flow experiences small modifications at the IR minimum even if µ < m l . This effect aggravates when µ approaches m l from below, µ m l . It leads to fluctuations of the chiral condensate around the vacuum IR value. Those fluctuations are small and hardly visible in Fig. 2(b). However, they lead to an unphysical phase of very small but non-zero pressure. Furthermore, for a flavor-dependent chemical potential the first-order transition is additionally distorted by the error of our approximation, see Sec. IV. Thus, in the numerical treatment of the FRG EoS, data points close to the phase transition are omitted and the EoS from the physical phase with restored chiral symmetry is polynomially extrapolated down to p = 0. This procedure only affects the low-pressure outer region of the calculated pure quark stars. The dependence of the star radius on the extrapolation error has been checked and found to be negligible.
For N f = 2+1, similar numerical fluctuations as discussed above are also found when µ s approaches the order of the strange quark mass m s at the current IR minimum. They are most prevalent in the determination of the speed of sound. A strong dependence of these fluctuations on the exact grid point configuration is observed and gives strong evidence for the claim as a numerical artifact. Therefore, in Fig. 6(a) for µ > 350 MeV only the average derivativec 2 s := ∆p/∆ε is shown as dots, uniformly spaced at a distance of 15 MeV. For each dot,c 2 s has been calculated from the (p, ε) tupels at the dot to its left, itself, and the dot to its right. Furthermore, the highest deviation of the microscopically calculated, fluctuating speed of sound c 2 s fromc 2 s in the respective interval is indicated by the border points of the shaded region such that all c 2 s data points lie within this region.
In the pure quark matter star calculations, the TOV equation is solved with an explicit Runge-Kutta algorithm. The evolu-tion is stopped when the radial pressure p(r) reaches a value of 10 −5 relative to the central pressure. The EoS data points are interpolated with cubic splines, utilizing two separate splines in case of a discontinuity due to a first-order transition such as observed in the rMFA EoS.