New Covariant Density Functionals of Nuclear Matter for Compact Star Simulations

We generate three families of extended covariant density functionals of nuclear matter that have varying slope of symmetry energy and skewness at nuclear saturation density, but otherwise share the same basic parameters (symmetry energy, compressibility, saturation parameters, etc.) with the standard DDME2, DD2, and MPE functionals. Tables of the parameters of these new density functionals are given, which can be straightforwardly used in DDME2, DD2, and MPE parameterization-based codes. Furthermore, we provide tables of a large number of equations of state (81 for each family) that can be used in astrophysical simulations to assess the impact of variations of the not-well-known slope of symmetry energy and skewness of nuclear systems on the astrophysics of compact objects. We also provide tables of computed integral parameters (mass, radius, and tidal deformability) that can be used, e.g., for modeling gravitational waveforms. Finally, for the extended DDME2-based parameterization, we implement a first-order phase transition to quark matter to obtain a family of equations of state that accommodates a phase transition to quark matter. Analogous tables of the equations of state and integral parameters are provided for this case as well.


INTRODUCTION
The detection of gravitational waves from binary compact sources by first-and second-generation gravitational-wave detectors and the prospect of higher-sensitivity observations by the third-generation gravitational observatories, such as the Einstein Telescope, provide a strong motivation for modeling the astrophysical sources of gravitational waves.Binary neutron star mergers offer significant information via the inference of the tidal deformation in the premerger phase, as already observed in the GW170817 event (Abbott et al. 2017a(Abbott et al. , 2018(Abbott et al. , 2019a)), as well as through observations of oscillations of the postmerger objects (Abbott et al. 2017b(Abbott et al. , 2019b)).
An important input in binary neutron star modeling is the equation of state (EoS) of dense matter in compact stars.To date, a large number of EoSs have been proposed which can be classified in the following groups: (a) the microscopic EoS which are based on first-principle computations, such as those based on potentials fitted to the nucleon-nucleon scattering data or those based on QCD, either in the perturbative or low-energy limits (for recent reviews, see Oertel et al. 2017;Burgio et al. 2021;Drischler et al. 2021); (b) density functional-based EoSs, which start with some form of energy density functional with adjustable parameters that are then fitted to bulk properties of nuclei and/or nuclear matter (for recent reviews, see Oertel et al. 2017;Blaschke & Chamel 2018;Sedrakian et al. 2023); (c) matter-compositionagnostic models, which are based on some form of a parameterization of the EoS, well-known examples being polytropes (Ozel & Psaltis 2009;Read et al. 2009), constant speed of sound (CSS) parameterizations (Alford et al. 2013;Zdunik & Haensel 2013), metamodels based on Taylor expansions of energy density close to saturation density (Margueron et al. 2018;Ferreira et al. 2022;Zhang & Li 2023), and nonparametric EoSs (Landry & Essick 2019;Landry et al. 2020;Essick et al. 2021), etc.
The covariant density functionals (CDFs) provide a fast and reliable method to incorporate physical constraints from bulk properties of nuclear systems and astrophysical information starting from baryon-meson Lagrangians (Serot & Walecka 1997;Bender et al. 2003;Vretenar et al. 2005;Meng et al. 2006;Oertel et al. 2017).They provide access to microscopic information, such as self-energies (decomposed according to their Lorentz structure), matter composition, chemical potentials, and effective masses.At the same time, they allow for flexible adjustments of model parameters with the inflow of new information and constraints on the EoS or other aspects of micro-and macro-physics.A good example is how CDF-based theories confront the so-called "hyperon puzzle"; see Sedrakian et al. (2023) and references therein.
To obtain a CDF, the partition function of a many-body system is often evaluated in the lowest-order Hartree approximation, which leads to relativistic mean-field models (Serot & Walecka 1997;Bender et al. 2003;Vretenar et al. 2005;Meng et al. 2006;Oertel et al. 2017).However, those CDFs that reproduce the available data to high accuracy do not strictly represent mean-field realizations of any given Lagrangian, as the correlations beyond the mean-field are built into the fitting procedure via additional elements, such as (a) nonlinear mesonic fields (e.g., Boguta & Bodmer 1977;Lalazissis et al. 1997;Todd-Rutel & Piekarewicz 2005;Chen & Piekarewicz 2014) and (b) the density dependence of the coupling constants (e.g., Typel & Wolter 1999;Long et al. 2004;Lalazissis et al. 2005;Typel et al. 2010;Typel 2018).As a result, the self-energies are not linear functions of density, as is the case in the mean-field theory.For this reason, we will avoid the term "relativistic mean-field model" and use CDF theory instead.
In this work, we provide a large set of EoSs that are based on CDF descriptions of nuclear matter.Our point of depar-ture will be three standard CDFs with density-dependent couplings given by Lalazissis et al. (2005), Typel et al. (2010) and Typel (2018), which were fitted to the laboratory data (commonly, the binding energies, charge radii, spin-orbit splittings, the neutron skin thickness, etc. for a set of magic and semimagic nuclei) using a least-squares minimization procedure.Note that the fitting procedures are not identical, but are based on the calibration of parameters on a particular set of the physical observables of finite nuclei and/or nuclear matter.Next, we generate a large number of CDFs by varying the density dependence of the coupling constants.These variations are then mapped onto the variations of the coefficients of the expansion of the energy density of isospin asymmetric nuclear matter in Taylor series around the nuclear saturation and isospin symmetrical limits.As a proof of principle, we have previously provided a smaller sample of CDF EoSs (Li & Sedrakian 2019), which were mapped on the EoS models following from purely phenomenological expansions of the nuclear matter EoS near the saturation density that is most prominently used in the metamodels of EoSs (Margueron et al. 2018), but also in the Bayesian analysis of EoS parameters (Char et al. 2020;Traversi et al. 2020;Malik et al. 2023;Zhu et al. 2023) and for providing sensible constraints on parameters of nuclear and neutron star matter (Lattimer & Prakash 2004;Lattimer 2023).Thus, the first goal of this work is the construction of three new and large families of nucleonic EoSs that at saturation reproduce a given standard CDF, such as DDME2 (Lalazissis et al. 2005), DD2 (Typel et al. 2010), and MPE (Typel 2018).They, however, possess an alternative density dependence of the couplings and, thus, the resulting coefficients of Taylor expansion of energy density.These families are constructed to cover a substantial range of the mass-radius diagram of compact stars that is compatible with current astrophysical constraints.This will allow us to carry out systematic studies of the dependence of astrophysical simulations on the important not-well-determined parameters of the nuclear EoS, such as the slope of symmetry energy and skewness of symmetric matter at nuclear saturation density.The second goal extends the DDME2 functional-based analysis to allow for a phase transition to quark matter.The quark matter EoS is modeled using simple CSS parameterization (Alford et al. 2013;Zdunik & Haensel 2013).We assume a strong, firstorder phase transition to quark matter in this case.Thus, this particular set of EoSs allows astrophysical implementations of QCD phase transition in dense matter.Finally, all the generated data for the EoSs and the integral parameters of compact stars, such as mass, radius, and tidal deformability, are provided in the form of tables, that can be downloaded from the https://github.com/asedrakian/DDCDFs repository, while static versions are available in Zenodo at doi:10.5281/zenodo.8350680.This work is structured as follows.In Section 2 we briefly review the formalism of CDF theory at the saturation density and its standard realizations in the case of density-dependent couplings.In Section 3 we propose new CDFs with modified couplings and present their predictions for the EoS, the mass-radius diagram, and the tidal deformabilities of static, spherically symmetrical stars.Section 4 describes the case where a first-order phase transition to quark matter is allowed for.Here we present the same parameters of static hybrid stars as in the case of nucleonic models.Our conclusions are collected in Section 5.

COVARIANT DENSITY FUNCTIONALS
In this work, we use the CDF approach based on the Lagrangian of stellar matter with nucleonic degrees of freedom L = L N + L m + L l + L em , where the nucleonic Lagrangian is given by where ψ N are the nucleonic Dirac fields with masses m N , and σ, ω µ , and ρ µ are the mesonic fields that mediate the interaction among nucleon fields.The remaining pieces of the Lagrangian correspond to the mesonic, leptonic, and electromagnetic contributions, respectively (Sedrakian et al. 2023).
The meson-nucleon couplings which are density-dependent and are given by (Typel & Wolter 1999;Typel 2018) with a constant value g m (ρ ref ) at a reference density ρ ref , and a function f m (r) that depends on the ratio r = ρ/ρ ref .
For the isoscalar channel, one has a rational function with conditions which reduce the number of free parameters.The density dependence for the isovector channel is taken in an exponential form: The energy density of isospin asymmetric matter is customarily split into an isoscalar and an isocector term: where ρ = ρ n + ρ p is the baryonic density, with ρ n(p) denoting the neutron (proton) density, δ = (ρ n − ρ p )/ρ is the isospin asymmetry, and E 0 (ρ) and E sym (ρ) are, respectively, the energy of symmetric matter and the symmetry energy.At densities close to the saturation ρ sat Equation ( 6) can be further Taylor expanded as where χ = (ρ − ρ sat )/3ρ sat .The coefficients of the expansion are known as nuclear matter characteristics at saturation density, namely, incompressibility K sat , the skewness Q sat , the CDF symmetry energy J sym and its slope parameter L sym .The coefficients of the leading-order terms in the expansion, namely K sat and J sym are well determined, whereas the coefficients of higher-order terms, namely Q sat and L sym are not well known.By definition, these coefficients allow one to assess the properties of nuclear matter predicted from various models.
In the most widely-used CDF parameterizations, e.g., DDME2 (Lalazissis et al. 2005) and DD2 (Typel et al. 2010), the reference density is chosen as the vector density at saturation ρ ref = ρ v sat .If we fix in the Lagrangian (1) the nucleon and meson masses, the properties of infinite nuclear matter (the seven characteristics listed in Table 1) can be computed uniquely in terms of seven adjustable parameters (Li & Sedrakian 2019).These are three coupling constants at saturation density (g σ , g ω , g ρ ) and four parameters (a σ , b σ , a ω , a ρ ) that control their density dependence.Besides these choices, we also use the newly proposed MPE parameterization (Typel 2018), for which dependence on the scalar density is considered for the coupling of scalar meson σ and the reference density ρ ref = ρ s sat .In this case, the third condition in Eq. ( 4) is dropped.The predicted nuclear matter parameters with these three CDFs are summarized in Table 1.It is evident that the choice of the reference density is part of the definition of the CDF and thus characterizes it.

NUCLEONIC MODELS
In this work, we use the DDME2, DD2, and MPE CDF parameterizations as reference CDFs on which to build our extension.Because the nuclear matter coefficients J sym and L sym are strongly correlated while the value of symmetry energy E sym at density ρ c = 0.11 fm −3 is almost identical for a variety of nuclear models (e.g., Trippa et al. 2008;Ducoin et al. 2011;Li et al. 2018), we hold the value of E sym (ρ c ) constant for each family of CDF when L sym is being varied.Consequently, J sym varies as well.There are also correlations between other pairs of coefficients.For example, K sat -Q sat pair is correlated but to a less noticeable degree (e.g., Typel 2018;Malik et al. 2020).A correlation among the pair L sym -K sym may also exist (e.g., Ducoin et al. 2011;Providência et al. 2014), such correlations being a property of CDF itself.Having fixed the values of ρ sat , E sat , K sat , and E sym (ρ c ), we proceed to quantify the uncertainties in the CDFs by covering a range of terms of not-well-known higher-order characteristics of nuclear matter, specifically Q sat and L sym .To this end, for isoscalar σand ω-mesons we fix their values of couplings at saturation density, but modify their density dependences; while for isovector ρ-meson we fix its coupling value at ρ c but modify its density dependence, so the value of coupling at saturation density is thus changed accordingly.The parameters entering each set of nuclear matter coefficients are obtained by a numerical search.
We then compute the EoS of cold stellar matter by implementing the conditions of β-equilibrium and charge neutrality, where both electrons and muons are included.We further match smoothly our core EoS to that of inner crust EoS given in Baym et al. (1971b) at the crust-core transition density ∼ ρ sat /2.For the outer crust, we use the EoS given in Baym et al. (1971a).With this input, the integral parameters of compact stars, in particular, the mass M and the radius R, are then computed from the Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman 1939;Oppenheimer & Volkoff 1939).The tidal deformability λ critically characterizes how easily an object can be deformed under an external tidal field (Hinderer 2008;Binnington & Poisson 2009).It is given as λ = 2/3 k 2 R 5 , 0 2 0 0 4 0 0 6 0 0 8 0 0 1 0 0 0 1 2 0 0 0 2 0 0 4 0 0 6 0 0 where k 2 is the dimensionless Love number that can be computed simultaneously with the TOV equations.We work with the dimensionless Λ = λM −5 for convenience.
Figure 1 shows the EoS of stellar matter generated by DDME2 family of CDF models [panel (a)] and for three families of CDF models, DDME2, DD2 and MPE, for pairs of (Q sat , L sym ) values as indicated, thus showing the impact of the form of density functional on the EoS [panel (b)].
Our collection in Figure 1 (a) contains in total 81 EoSs for each fixed set of low-order characteristics of DDME2, which are generated by varying −600 ≤ Q sat ≤ 1000 MeV and 30 ≤ L sym ≤ 110 MeV with steps ∆Q sat = 200 MeV and ∆L sym = 10 MeV.Tables 2 and 3 in Appendix provide the corresponding parameter sets for each CDF.The range of these parameters is chosen as follows.For Q sat = −600 MeV and independent of the value L sym the maximum mass is about 2.0 M ⊙ , which matches the mass measurement of PSR J0740+6620 (Cromartie et al. 2019;Fonseca et al. 2021).For Q sat ≥ 600 MeV, the maximum mass of a compact star reaches the value 2.59 +0.08 −0.09 M ⊙ deduced for the light companion of GW190814 (Abbott et al. 2020) if interpreted as a neutron star.The lower and upper values of L sym are those predicted by the majority of the density functionals and are within the range currently discussed (see, e.g., Adhikari et al. 2021Adhikari et al. , 2022;;Lattimer 2023;Sedrakian et al. 2023).
Next, let us assess the influences of the low-order characteristics and the different functional forms of the couplings that we have fixed within each family of CDFs on the EOS and the gross properties of neutron stars.From Figure 1   we observe that: (a) the differences in the low-order characteristics of nuclear matter among the three families of CDFs, which control the properties of finite nuclei and nuclear matter at subsaturation densities, have negligible impact on the EoS, especially at high densities, which are more relevant for neutron stars.This is clearly seen from the comparison of the DDME2 and DD2 CDFs which become identical for the stiff EoS in each family; and (b) the high-density behavior of the EoS depends strongly on the specific chosen family of CDFs (which differ by their functional form) and/or the choice of density dependencies of couplings on scalar or vector density.As seen from the comparison of the DDME2 and MPE CDFs, sizable differences emerge in the high-density domain, especially for those parameters that lead to soft EoSs.
Figure 2 shows the maximum mass M max [panel (a)] and the radius R 1.4 of a canonical 1.4 M ⊙ mass star [panel (b)] for static, spherically symmetrical stellar configurations predicted by the three families of EoSs generated from the three different CDFs we used.Here we vary continuously Q sat for some representative values of L sym .In agreement with the observed difference among EoSs, it is clearly seen that the choice of the specific form of the CDF and corresponding low-order characteristics has a small impact (< 2%) on the radius of canonical 1.4 M ⊙ mass star.The variations of such a magnitude expected for the low-mass domain of compact stars are significantly smaller than current measurement accuracy of ∼ 10%.At the same time, the differences in the maximum masses predicted by DDME2/DD2 and MPE CDFs can reach up to 5%, which correspond to ∼ 0.1 M ⊙ , and G W 1 7 0 8 1 7 P S R J 0 0 3 0 + 0 4 5 1 L s y m = 3 0 1 1 0 Q s a t = 1 0 0 0 -6 0 0 P S R J 0 7 4 0 + 6 6 2 0 G W 1 9 0 8 1 4 ' s s e c o n d a r y In panel (a) the color regions show the 90% confidence interval (CI) ellipses from each of the two NICER modeling groups for PSR J0030+0451 and J0740+6620 (Riley et al. 2019(Riley et al. , 2021;;Miller et al. 2019Miller et al. , 2021)), the 90% CI regions for each of the two compact stars that merged in the gravitational-wave event GW170817 (Abbott et al. 2019a), and finally the 90% CI for the mass of the secondary component of GW190814 (Abbott et al. 2020).In panel (b), the constraint for a 1.36 M ⊙ star deduced from the analysis of GW170817 event (Abbott et al. 2019a) is shown too.
are comparable to the error bars in the mass measurements of very massive (close or beyond two solar mass) compact stars.Therefore, the choice of the functional form of the CDF (among many other factors) influences the predictions for the properties of compact stars and justifies considering different families of CDFs in astrophysical studies of compact stars.
Figure 3 shows the mass-radius curves [panel (a)] and mass-tidal deformability relations [panel (b)] for static stellar configurations based on the 81 EoS of the DDME2 family.Constraints from multimessenger astronomy are highlighted for comparison.These concern the mass and radius measurements for PSR J0030+0451 (Riley et al. 2019;Miller et al. 2019) and J0740+6620 (Riley et al. 2021;Miller et al. 2021) by NICER, the compactness and tidal deformability constraints extracted from the binary compact star merger GW170817 (Abbott et al. 2019a), and the mass of the secondary component of GW190814 event (Abbott et al. 2020).For these observations, all the regions/limits are given at 90% credible intervals.As seen from Figure 3, the softness of the EoSs at the intermediate density implied by the GW170817 event and the stiffness at high density suggested by the massive pulsars/objects can be reconciled by an appropriate choice of the parameters of L sym and Q sat .
Similarly, we generate the collection from the MPE family, for which the parameters are provided in Tables 4 and 5 in the Appendix.Note that with this functional form, a slightly higher value Q sat ≥ −500 MeV is required to reach

HYBRID MODELS
About half of our models from the previous section have tidal deformabilities that are outside the inference for GW170817.However, these EoSs can be reconciled with this inference if a phase transition to quark matter occurs.In this section, we demonstrate the utility of such models by supplementing them with a quark EoS using as a model the CSS parameterization (Zdunik & Haensel 2013;Alford et al. 2013) where p tran is the transitional pressure with energy density ε tran at nucleonic density ρ tran , ∆ε is the energy discontinuity, and c s is the sound speed in the quark phase.Figure 4 illustrates how such a transition to quark matter allows a nuclear EoS with a large L sym to be consistent with astrophysical constraints.We show two scenarios in which the phase transition occurs at low or high enough densities, respectively, considering different values of c s .For stiff quark EoS (c s ≈ 1.0), which allows a large enough jump in energy density ∆ε under consideration, there can be a second stable branch of compact stars.This leads to the emergence of twin configurations -two stable compact stars have an identical mass, but different radii and, consequently, tidal deformabilities (e.g., Alford et al. 2013;Benic et al. 2015;Paschalidis et al. 2018;Alvarez-Castillo et al. 2019;Christian et al. 2019;Montana et al. 2019).In the first scenario, the phase transition occurs at ρ tran /ρ sat ≲ 2.0, causing the nucleonic branch to end at a lower mass ≲ 1.4 M ⊙ , and the hybrid branch extends up to M max > 2 M ⊙ (e.g., Paschalidis et al. 2018;Montana et al. 2019;Blaschke et al. 2020;Li et al. 2021;Christian & Schaffner-Bielich 2022;Li et al. 2023a).In the second scenario, in contrast, the phase transition occurs at ρ tran /ρ sat ≳ 3.0, so that the heaviest star on the nucleonic branch has M max > 2 M ⊙ , and the hybrid branches lie at lower masses, which yield ultracompact hybrid stars (Li et al. 2023b).
Note that ellipses/limits referring to GW170817 have been obtained under the assumption of nucleonic stars, therefore they are relevant, strictly speaking, only for constraining the nucleonic EoS models.The ultracompact hybrid star models shown in Figure 4, however, could satisfy the mutual dependence of tidal deformabilities Λ 1 − Λ 2 of members of a binary for GW170817 (Li et al. 2023b).In this scenario, the multimessenger data can be accounted for if PSR J0030+0451 and J0740+6620 are purely nuclear stars with large radii, whereas one of the components of GW170817 binary is ultracompact hybrid star (Li et al. 2023b).Ultracompact star arise in scenarios where the conformal limit of QCD is reached at quite high densities.If this density is low, then the mass-radius curves are limited to the regions with larger radii, around 10-12 km (Sedrakian 2023).

CONCLUSIONS
In this work, we constructed three families of CDFs starting from DDME2, DD2, and MPE parameterization allowing for variations of the slope of the symmetry energy and skewness of symmetric matter at saturation density.The new CDFs allow for systematic studies of the sensitivity of various astrophysical scenarios toward variations of these parameters.One way to utilize the results is to modify the codes that are based on the original CDFs with density-dependent couplings to generate the EoS and other observables.Tables 3 and 5 in the Appendix provide the parameter values of extended density functionals for such purpose.Another way to utilize our results is to use directly the provided tables of EoSs in numerical simulations.In addition, we provided tables of integral parameters, such as the mass, radius, and tidal deformability, which can be used in astrophysical analysis, for example, to generate gravitational-wave templates.The effects of first-order phase transition in neutron stars were studied using the simple CSS parameterization of the EoS of quark matter within two scenarios of low-and high-density phase transition to quark matter.The EoS and integral parameter tables are also provided in this case.
So far, we have concerned only been concerned with purely nucleonic models.Of course, hyperonization and the onset of resonances are a possibility in compact stars, in which case both the EoS and integral parameters are modified (e.g., Chatterjee & Vidaña 2016;Tolos & Fabbietti 2020;Sedrakian et al. 2023).Clearly, our alternative parameterizations can be used to study hypernuclear matter, eventually admixed with ∆-resonances, to study the sensitivity of the EoS on the highorder characteristics of nuclear matter.2 and 4, respectively.As we already discussed in the text, these models are well constrained with respect to the lower-order saturation characteristics.In Tables 3 and 5 we further present sets of alternative parameterizations that produce different values of Q sat and/or L sym , but preserve these lower-order values predicted by the original DDME2 and MPE parameterizations.To this end, one needs to modify only the parameters in functions f σ and f ω [see Equation (3)] which control the density dependence of the couplings in the isoscalar sector, and f ρ [see Equation ( 5)] in the isovector sector.Note.
-The masses are in units of MeV.
Table 3 Alternative parameterization of the density dependence of the couplings in the isoscalar and isovector channels for the indicated values of Q sat (MeV) and/or L sym (MeV).The values of g σ and g ω are the same as in the DDME2 parametrization that given at the saturation density, while that for g ρ here refers to the value at the cross density ρ c = 0.11 fm −3 .-For ρ meson coupling, the value of g ρ (ρ sat ) can be simply obtained through formula g ρ (ρ sat ) = g ρ (ρ c )e −aρ(1−ρc/ρsat) , and a larger value of a ρ corresponds to a lower value of the slope parameter L sym .Note.
-The masses are in units of MeV, and density in fm −3 .
Table 5 Alternative parametrization of the density dependence of the couplings in the isoscalar and isovector channels for the indicated values of Q sat (MeV) and/or L sym (MeV).The values of g σ and g ω are the same as in the MPE parametrization that given at the saturation density, while that for g ρ here refers to the value at the cross density ρ c = 0.11 fm −3 .
Q -For ρ meson coupling, the value of g ρ (ρ sat ) can be simply obtained through formula g ρ (ρ sat ) = g ρ (ρ c )e −aρ(1−ρc/ρsat) , and a larger value of a ρ corresponds to a lower value of the slope parameter L sym .

Figure 1 .
Figure 1.EoS for purely nucleonic stellar matter.In panel (a) the models are generated with DDME2 family of CDF models by varying the parameters Q sat ∈ [−600, 1000] MeV and L sym ∈ [30, 110] MeV.The effects of parameter L sym on the low-density region of EoS are shown in the inset for illustration.In panel (b) the same is shown for three families of CDF models with fixed values of pairs (Q sat , L sym ) (in MeV) as indicated in the plot.
(b), where we vary the form of the CDF by fixing pairs of values (Q sat , L sym ),

Figure 2 .
Figure 2. The maximum mass M max [panel (a)] and the radius R 1.4 of a canonical 1.4 M ⊙ mass star [panel (b)] that are predicted by the three families of EoSs generated from each CDF used for various values of L sym as a function of Q sat (in MeV).In panel (a) the shadings show the masses of PSR J0470+6620(Cromartie et al. 2019;Fonseca et al. 2021) and the secondary component of GW190814(Abbott et al. 2020), while in panel (b), the vertical line indicates the upper limit on the radius set by the analysis of GW170817(Abbott et al. 2018).

Figure 3 .
Figure 3. Mass-radius [panel (a)] and mass-tidal deformability [panel (b)]relations for nucleonic EoS models with different pairs of values of Q sat and L sym (in MeV).In panel (a) the color regions show the 90% confidence interval (CI) ellipses from each of the two NICER modeling groups for PSR J0030+0451 and J0740+6620(Riley et al. 2019(Riley et al. , 2021;;Miller et al. 2019 Miller  et al.  , 2021)), the 90% CI regions for each of the two compact stars that merged in the gravitational-wave event GW170817(Abbott et al.  2019a), and finally the 90% CI for the mass of the secondary component of GW190814(Abbott et al. 2020).In panel (b), the constraint for a 1.36 M ⊙ star deduced from the analysis of GW170817 event(Abbott et al. 2019a) is shown too.

Figure 4 .
Figure 4. Mass-radius [panel (a)] and mass-tidal deformability [panel (b)] relations for hybrid EoS models constructed from a stiff nucleonic EoS (Q sat = 600, L sym = 90 MeV) with phase transitions at low and high densities.For all the hybrid branches, the values of sound speed squared c 2 s are as indicated in the plot.The dotted thin lines indicate unstable configurations.Constraints from multimessenger astronomy are shown too.

Support
Program for Chongqing Overseas Returnees (grant No. CX2021007). A. S. acknowledges the support from the Deutsche Forschungsgemeinschaft grant No. SE 1836/5-2, the Polish NCN grant No. 2020/37/B/ST9/01937 at Wrocław University, and the Institute for Nuclear Theory at the University of Washington, for its kind hospitality, stimulating research environment, and partial support from the INT's U.S. Department of Energy grant No. DE-FG02-00ER41132.APPENDIX PARAMETER VALUES OF THE DDME2 AND MPE FAMILIES In this appendix we present the values of the parameters determining the CDFs.The parameters of the original DDME2 and MPE parameterizations are shown in Tables

Table 1
Nuclear matter characteristics at saturation density of the reference CDF parameterizations.

Table 2
Nucleon, meson masses and meson-nucleon coupling constants in the DDME2 parameteriztation, where g m (m = σ, ω, and ρ) refer to the values at the saturation density.

Table 4
Nucleon, meson masses and meson-nucleon coupling constants in the MPE parametriztation, whereby g m (m = σ, ω, and ρ)refer to the values at the saturation density, and ρ s sat the corresponding scalar density.