Constraints on Strong Phase Transitions in Neutron Stars

We study current bounds on strong first-order phase transitions (PTs) along the equation of state (EOS) of dense strongly interacting matter in neutron stars, under the simplifying assumption that on either side of the PT, the EOS can be approximated by a simple polytropic form. We construct a large ensemble of possible EOSs of this form, anchor them to chiral effective field theory calculations at nuclear density and perturbative Quantum Chromodynamics at high densities, and subject them to astrophysical constraints from high-mass pulsars and gravitational-wave observations. Within this setup, we find that a PT permits neutron-star solutions with larger radii, but only if the transition begins below twice nuclear saturation density. We also identify a large parameter space of allowed PTs currently unexplored by numerical-relativity studies. Additionally, we locate a small region of parameter space allowing twin-star solutions, though we find them to only marginally pass the current astrophysical constraints. Finally, we find that sizeable cores of high-density matter beyond the PT may be located in the centers of some stable neutron stars, primarily those with larger masses.


INTRODUCTION
The detailed phase structure of Quantum Chromodynamics (QCD), the theory governing the strong nuclear interaction, is still largely unknown. A combination of lattice-field-theory calculations and two decades of ultrarelativistic heavy-ion-collision experiments carried out at RHIC and at the LHC has established the crossover nature of the deconfinement transition at vanishing and small baryon number chemical potentials (Aoki et al. 2006;Cheng et al. 2006). However, the phase structure of baryon-rich matter present in the cores of neutron stars (NSs), and probed in lower-energy heavy-ion collisions, remains uncharted.
Recent years have witnessed a rapid evolution in NS observations, which has for the first time permitted model-independent constraints of the equation of state (EOS) of the dense matter in NSs. This inference of the EOS has been based on the generation of large ensembles of model-agnostic EOSs that are then conditioned to be consistent with both NS observations and ab initio calculations Kurkela et al. 2014;Most et al. 2018;Annala et al. 2018;Tews et al. 2018;Landry & Essick 2019;Capano et al. 2020;Miller et al. 2020;Essick et al. 2020;Raaijmakers et al. 2020;Annala et al. 2020;Dietrich et al. 2020;Landry et al. 2020;Al-Mamun et al. 2021;Miller et al. 2021;Essick et al. 2021;Raaijmakers et al. 2021;Annala et al. 2022;Huth et al. 2022;Altiparmak et al. 2022;Lim & Holt 2022;Gorda et al. 2023). The observations that currently constrain the EOS the most include the existence of massive NSs with masses around and exceeding the two-solarmass limit (Demorest et al. 2010;Antoniadis et al. 2013;Cromartie et al. 2019;Fonseca et al. 2021), constraints on the tidal deformability obtained from the binary-NS merger GW170817 (Abbott et al. 2017(Abbott et al. , 2018(Abbott et al. , 2019, as well as the simultaneous mass-radius measurements performed with the NICER telescope (Miller et al. 2019;Riley et al. 2019;Miller et al. 2021;Riley et al. 2021). On the theoretical side, the EOS is constrained by ab initio calculations based on chiral effective field theory (EFT) at nuclear densities Tews et al. 2013;Lynn et al. 2016;Drischler et al. 2019Drischler et al. , 2020Keller et al. 2023) and perturbative QCD (pQCD) at very high densities (Kurkela et al. 2010;Gorda et al. 2021).
While the EOS is an interesting and fundamental quantity in itself, interest in the EOS also arises from the fact that its features can reflect the active degrees of freedom in the system and hence the physical phase of matter. In the past years, there have been several works that have suggested features of the EOS may be related to the onset of hyperonic matter [e.g., (Gal et al. 2016)], the onset of quark matter [e.g., (Annala et al. 2020)], the presence of a quarkyonic phase (McLerran & Reddy 2019) or of diquark pairing [e.g., (Leonhardt et al. 2020)].
The most dramatic connection between the EOS and the phase structure takes place if the EOS contains a strong first-order phase transition (PT) between a low-density hadronic phase and a new high-density phase.
Indeed, there have been multiple studies of EOSs featuring strong PTs, displaying, e.g., qualitatively different gravitational-wave signals arising from binary-NS mergers compared to a scenario with vanishing latent heat (Most et al. 2019;Bauswein et al. 2019;Fujimoto et al. 2023). These studies have, however, been restricted to individual models, so it is of great interest to systematically determine what bounds on first-order PTs follow from currently existing theoretical and observational constraints. Such results would be of great value for merger simulations, as they would map out the space of possible PTs and suggest what ranges of EOSs one should cover.
While very versatile, the model-agnostic EOS-inference studies performed so far have focused less on exploring the broad space of strong PTs. In this work, our objective is to perform a simple inference study of the NSmatter EOS including first-order PTs and setting no a priori limitations for the strength and onset density of the transition. We do so by generating a large ensemble of EOSs that continue from a chiral EFT EOS band using a two-segment polytrope with a PT inserted between them. This EOS is then constructed up to densities well beyond those reached in stable NSs, in practice ten times nuclear saturation densities n = 10n 0 (with n 0 = 0.16 fm −3 ). At this density two methods, detailed in Sec. 2, are used to restrict the range of allowed values for the energy density ϵ and pressure p to enforce a high-density constraint from pQCD. This simple setup, though less general than previous model-agnostic approaches without a first-order PT, allows us to clearly delineate the matter below and above the transition and remains flexible enough to allow a general exploration of PTs along the NS-matter EOS without resorting to a specific microscopic model at intermediate densities.
This article is structured as follows. Following the details of our EOS construction in Sec. 2, we divide the discussion into smaller sections, each focusing on specific findings. First, we investigate the effects of a PT on the set of possible masses and radii of stable NSs in Sec. 3, finding larger-radius stars than were found in previous model-agnostic studies. In Sec. 4, we then study the range of allowed PT parameters, namely the onset densities and latent heats consistent with current NS observations and ab initio calculations. Here we also compare these allowed ranges with specific models used in the merger literature and investigate the parameter space for twin-star (or third-family) solutions within our approach. The latter are discussed in detail in Sec. 5. Finally, we investigate the possibility of cores of high-density matter within stable massive NSs in Sec. 6, detailing the ranges of NS masses and radii consistent with them. Finally, we conclude with a summary in Sec. 7.

SETUP
We build our EOS ensemble using a continuous piecewise construction up to the baryon density of 10n 0 . For the density region n ≤ 0.57n 0 , we use the BPS crust EOS from Baym et al. (1971), followed by an EOS within the chiral EFT band spanned by the "stiff" or "soft" EOSs from Hebeler et al. (2013) up to a density n EFT ≡ 1.1n 0 . For simplicity, we generate intermediate EOSs between the stiff and soft EOSs with a linear interpolation of p, ϵ, and baryon chemical potentials at fixed n, which well approximates the results obtained from many-body calculations in that work. Beyond this density, we extend the EOS with a two-polytrope construction with a firstorder PT in between. The parameters of the PT, n PT and ∆n, specify the onset number density and jump in number density at the PT, respectively. Explicitly, this implies: 1. For n EFT < n ≤ n PT , the EOS is governed by a single polytropic form, namely, p(n) = p EFT (n EFT )(n/n EFT ) Γ1 , 2. For n PT < n ≤ n PT + ∆n, there is a first-order PT where the baryon density jumps by ∆n but the pressure remains constant, 3. For n PT + ∆n < n ≤ 10n 0 , we use another polytropic EOS, matched to the end of the PT: p(n) = p(n PT )[n/(n PT + ∆n)] Γ2 .
At the density of 10n 0 we use one of two possible highdensity constraints. First, we use the EOS ensemble generated in Annala et al. (2022) to discard all EOSs with inconsistent ϵ, p values at this density. This choice allows us to make a comparison to this work, whose ensemble did not feature explicit PTs, and therefore to directly assess the effect of first-order PTs on such EOS ensembles and the macroscopic properties of the corresponding NSs. Second, we match to the less restrictive 10n 0 region of Komoltsev & Kurkela (2022), which is obtained by demanding that the extrapolated EOSs can be connected to the perturbative QCD EOS at densities around 40n 0 (where this calculation is under quantitative control). We use this second matching when comparing the resulting allowed range of PTs with the models that have been used in recent binary-NS simulations. In all, our construction requires five parameters to fully specify a given EOS: 1. a continuous parameter s ∈ [0, 1] that linearly interpolates between the soft and stiff chiral EFT EOSs of Hebeler et al. (2013) at fixed n, 2. n PT ∈ [1.1, 10]n 0 , 5. p within the same region at 10n 0 .
These five parameters are sampled from uniform distributions and together uniquely fix Γ 1 and Γ 2 . Explicitly, s, Γ 1 , n PT , ∆n/n PT , and Γ 2 uniquely fix ϵ, p at 10n 0 , and this process can be inverted. Note that the upper bound on ∆n/n PT is a consequence of the construction: any EOS is constrained to have ∆n ≤ 10n 0 −n EFT . We additionally reject any EOS that becomes superluminal (with an acausal speed of sound c 2 s > 1) in the interpolated region. The five parameters are sampled uniformly, and for the present study, we sample in total about 1,500,000 possible causal EOSs.
Once we have generated our ensemble of EOSs, two astrophysical constraints are folded in as hard cuts: 1. Any viable EOS must be able to support a NS of 2M ⊙ mass, to be compatible with the observations of high-mass NSs. 1 2. Any viable EOS must be consistent with the LIGO/Virgo 90% credible interval for the tidal deformability Λ inferred from the GW170817 event.
Here, we use the fact that the chirp mass is known to very good accuracy M chirp = 1.186M ⊙ , while the mass ratio and tidal deformability are constrained by q = M 2 /M 1 > 0.73 andΛ < 720 (Abbott et al. 2019), with the latter defined bỹ Λ = 16 13 In this relation, M 1 and M 2 are the gravitational masses of the two NSs in the event (with M 2 < M 1 ), and Λ(M 1 ) and Λ(M 2 ) are their respective tidal deformabilities (assuming a common EOS) as defined in Hinderer (2008).
After implementing these two constraints, approximately 280,000 EOSs of our sample remain as viable candidate EOSs containing a first-order PT. This is the ensemble that we proceed to systematically investigate. EOSs within our ensemble, after imposing the two astrophysical constraints M TOV > 2M ⊙ andΛ GW170817 < 720. In the top panel, the cyan regions corresponds to our matching p-ϵ region at n = 10n 0 , while in the bottom panel, the black outline displays the allowed region of EOSs from Annala et al. (2022), which does not include explicit PTs. allowed regions (black lines) obtained in that work. This comparison has the advantage that we can directly see the impact of strong PTs. From the results of Fig. 1 we observe that PTs extend the range of allowed NS radii consistent with current astrophysical observations beyond the 13.5 km found in Annala et al. (2022). A closer inspection into those EOSs that feature the largest radii at masses around or above 1.4M ⊙ reveals that they originate from EOSs featuring strong PTs beginning at low densities. This is explicitly seen in Fig. 2 where we divide the results into EOSs with transition densities n PT < 2.0n 0 or n PT ≥ 2.0n 0 , demonstrating that the largest-radius solutions primarily stem from early PTs below 2n 0 NSs. Moreover, the early-PT EOSs are found  to primarily feature values of the parameter s ≲ 0.2, and hence correspond to matching to the softer part of the chiral EFT band of Hebeler et al. (2013). These EOSs then rapidly stiffen prior to the PT, which is necessary to fulfil the astrophysical constraints, in particular the existence of 2M ⊙ NSs. For the early PTs the majority of the EOSs with s > 0.2, arising from the harder part of the chiral EFT band, fail to pass the tidal-deformalibity constraint.
The occurrence of these larger-radius stars can also be traced to the broadening of the correlation between the radius and tidal deformability of a 1.4M ⊙ star. This is shown in the top panel of Fig. 3, again grouped into early PTs, n PT < 2.0n 0 , and ones occurring later, n PT ≥ 2.0n 0 . Phase transitions are seen to broaden the correlation by allowing stellar configurations built with early-PT EOSs to reach radii beyond 13.5 km, confirming the observations of Han & Steiner (2019). The bottom panel of Fig. 3 further demonstrates the linking of early PTs with larger radii. In this figure, we also observe an  Top panel: Correlation of R(1.4M ⊙ ) and Λ(1.4M ⊙ ), divided into EOSs with transition densities n PT < 2n 0 (blue) and n PT ≥ 2n 0 (red). The overlap of these regions is shown in blended color. Bottom panel: Correlation of Λ(1.4M ⊙ ) with the transition density n PT , with the points colored by their R(1.4M ⊙ ) values. Note that the largest-radius stars exhibit early PTs.
anticorrelation between R(1.4M ⊙ ) and n PT , as was previously found in Miao et al. (2020). Finally, we note that due to the necessary fine-tuning to fulfill all constraints, these early PTs form a very small fraction of our EOS ensemble, so that these large-radius stars are expected to obtain small posterior distributions in a corresponding Bayesian analysis.

MATCHING TO GENERAL n = 10n 0 REGION
We now turn to the more general EOS ensemble matched to the maximal n = 10n 0 region of Komoltsev & Kurkela (2022). In the top panel of Fig. 4, we show the allowed region of ϵ and p probed within stable NSs using this construction. For comparison, we also show the boundary from the analysis of Annala et al. (2022) (black lines) for the entire NS EOS. Similar to the M -R region discussed above, we observe that when PTs are included, the set of allowed EOS points extends beyond the previous p-ϵ region towards larger pressure values at low  points along the stable NS sequence with twin-star configurations removed, when using the more general n = 10n 0 matching region of Komoltsev & Kurkela (2022) (shown in salmon color), after imposing the two astrophysical constraints M TOV > 2M ⊙ andΛ GW170817 < 720. Bottom panel: Individual p-ϵ sequences corresponding to the twin-star solutions, with the unstable branches shown in light gray. and moderate energy densities. A similar increase in the p-ϵ region is also seen in Fig. 1. In particular, EOSs that are very stiff at low densities exit the allowed region from Annala et al. (2022) but are still permitted within our present construction. These are precisely the same types of EOSs that lead to the large-radius stars in Sec. 3. In addition, in our more general ensemble we find a small number of EOSs with twin-star (or third family) solutions (Gerlach 1968a,b;Schertler et al. 2000), in which the stable stellar sequence first becomes unstable just after the PT, but then becomes stable again at higher central densities following a stiffening of the EOS beyond the PT. These individual solutions are shown separately in the bottom panel of Fig. 4, and are discussed in detail in Sec. 5.  Fig. 1, while the salmon one is the more general n = 10n 0 matching region of Komoltsev & Kurkela (2022). The points denoted by "causality" correspond to our construction prior to applying astrophysical constraints.
To more closely examine the differences between our current construction and previous results, we study the p-ϵ matching region at n = 10n 0 in Fig. 5. We show the regions from Annala et al. (2022) and Komoltsev & Kurkela (2022), with the points from our current ensemble overlayed. The Annala et al. (2022) region indeed is a subset of the maximal region, as it must be. We also observe that before applying the astrophysical constraints, our current simple construction (denoted by "causality" in Fig. 5) covers about half the area of the two regions (either Annala et al. (2022) or Komoltsev & Kurkela (2022)), with a significant number of EOSs lying outside the Annala et al. (2022) region. More interestingly, even after applying the astrophysical constraints, we still find a significant number of EOSs that lie outside of that region, which emphasizes the effect of including PTs in a general EOS ensemble. We also note that the impact of the pQCD constraint from Komoltsev & Kurkela (2022) is visible in this figure as the cut to the upper right of the causal and astrophysically allowed regions.
In Fig. 6, we show the p-ϵ region with the individual PTs from our ensemble, distinguished by whether any point beyond the PT is stable (top panel: ϵ max > ϵ PT +∆ϵ) or whether the PT terminates the NS sequence (bottom panel: ϵ max = ϵ PT ), with ϵ max denoting the highest energy density reached along the stable stellar sequence. Practically, we define this ϵ max as the last central density whose forward finite-difference mass derivative satisfies dM/dn central > 0, using a high-resolution grid in number density. We note that our resulting definitions of stability or destabilization are not strictly equivalent  Komoltsev & Kurkela (2022) shown in salmon color. Top panel: EOSs for which the PT completes within the stable NS sequence, with twin solutions shown in green. Note that for most of the twin solutions, the central densities occur at much higher densities than the ends of the PTs. Bottom panel: EOSs for which the PT destabilizes the NS sequence.
to the microscopic stability analysis conducted in Seidov (1971); Schaeffer et al. (1983); Lindblom (1998). As discussed in Lindblom (1998), such microscopic stability does not necessarily lead to macroscopic cores of highdensity matter beyond the PT, and hence our definition will neglect some cores of negligible size. Also in this context, we note that the above classification into PTs which do or do not terminate the stellar sequence does not represent a clean delineation within the full ensemble since there are EOSs with arbitrarily small values of Γ 2 that effectively prolong the PT beyond the value of the ∆n parameter. Within these and future figures, when presenting results for EOSs that possess stable stellar configurations beyond the PT, we additionally remove any EOSs possessing Γ 2 < 0.15 from the ensemble. 2 We clearly observe a minimal ϵ PT ≈ 275 MeV/fm 3 for PTs that terminate the NS sequence, which is set by the astrophysical 2M ⊙ constraint. We also see that EOSs exhibiting twin-star solutions have a particular range of PT parameters, which will be further studied below. Moreover, Fig. 6 shows how the PTs in our general ensemble go beyond the region from Annala et al. (2022) without explicit PTs. Turning then to the M -R sequences in Fig. 7, we find the M -R region without twin solutions to be very similar to Fig. 1 above, following from matching to the highdensity region from Annala et al. (2022). This is not surprising, since the constraints at small and large R follow directly from the astrophysical 2M ⊙ and GW170817 constraints, respectively. The main effect of the more general 10n 0 matching region is to allow for twin solutions with PTs, which are plotted in the bottom panel of Fig. 7.
Next, we examine the range of allowed PT parameters ϵ PT and ∆ϵ, or n PT and ∆n, for the EOSs we have constructed. To make this discussion more concrete, we also compare our allowed ranges to values that have been studied in the literature, see Table 1. We show in Fig. 8 the allowed values of ϵ PT and ∆ϵ for PTs that either complete within stable NSs or destabilize the NS sequence, while the corresponding values of n PT and ∆n are shown in Fig. 9. In both figures, we separate the twin solutions from the remaining distribution of allowed ranges, and for comparison show the 99.5% confidence contour when matching to the Annala et al. (2022) region at n = 10n 0 (red dashed lines). We emphasize that, due to our choice of performing the high-density matching at n = 10n 0 , there is a range of parameters that cannot be reached in our simple model, indicated as the grey region in Fig. 9. Likewise there is a similar region in Fig. 8, but its boundary is not easy to parameterize.
Our first observation from these figures is that there is a substantial overlap in the PT parameters in the ∆ϵ- ϵ PT or ∆n-n PT planes that can either destabilize the sequence or complete within stable NSs, leading to highmass cores of matter beyond the phase transition. We will discuss these cores in more detail below in Sec. 6. Next, we compare our results to PTs studied in the literature, collected in Table 1. We mark these points in Figs. 8 and 9, placing them in panels that indicate how the corresponding EOSs were constructed: namely whether they were engineered to complete within the stable NS sequence or to destabilize it. We observe that past literature has mostly focused on parameters corresponding to lower transition densities, and that there remains a large range of allowed PT parameters for ϵ PT ≳ 800 MeV/fm 3 or n PT ≳ 4n 0 that has not been explored in merger simulations. Moreover, stronger (destabilizing) PTs with ∆ϵ ≈ 750 MeV/fm 3 and ∆n ≈ 3n 0 at intermediate densities have also not been probed. Additionally, we find a substantial range of destabilizing PTs that have not  Table 1). The dashed red line denotes the 99.5% confidence contour when matching to the Annala et al. (2022) region at n = 10n 0 .
been studied, but may lead to interesting dynamics in the post-merger phase of a binary-NS merger. In more detail, the PT scenarios explored in Bauswein et al. (2019) and Somasundaram et al. (2023) overlap nicely with our most likely range of PT parameters, while Most et al. (2019) studied a particularly strong PT 3 and Fujimoto et al. (2023) utilized one of intermediate strength. However, there remain even earlier and weaker PTs that destabilize the star and are allowed within our construction. Finally, Han & Steiner (2019) [see also Alford et al. (2013)] examined three different scenarios (labeled B, C, D) of EOSs with n max > n PT + ∆n, the parameters of which are shown in the top panel of Fig. 9. Their classes C and D correspond to stellar sequences that are connected or disconnected at the PT, respectively, meaning that the point just beyond the PT at n PT + ∆n is either stable or unstable. In the latter case (D), the sequence later includes another stable twin-star branch due to a stiffening of the EOS beyond the PT. Class B (for "both") describes EOSs where the point just after the PT is stable, but which later experience a destabilization and then a further stabilization with a twin-star branch. We find that classes B and C, where the point n PT +∆n is stable, lie within our region of nontwin-star solutions, while class D lies outside our region, somewhat in the vicinity of our twin-star solutions.

TWIN STARS
From the lower panel of Fig. 7 we see that two separate classes of twin-star solutions are generated in our ensemble. One class features a small PT that locally destabilizes the NS sequence, but the sequence is then quickly restabilized by a stiffening of the EOS. These solutions have high-mass twins all possessing M ≈ 2M ⊙ . For these solutions, the stable stellar sequence is almost connected, and as such, the PT parameters overlap with those of the EOSs without twin-star solutions. We note that these very few high-mass-twin EOSs lead toΛ ≈ 720, so that they only marginally pass the GW170817 tidaldeformability constraint and would obtain very small posterior distributions in a corresponding Bayesian analysis.
The second and more interesting class of twin-star solutions is more exotic, with a much larger phase transition that destabilizes the stellar sequence for a large range of central densities. These solutions have a low-mass twin branch beginning at a central density close to n = 10n 0 and remaining stable up to the matching density. The EOSs featuring low-mass twins have the peculiar feature that the low-density NS branch has large radii and does not on its own pass the GW170817 constraint, but it is the twin-star branch, with masses M ≈ 1.4M ⊙ , that allows the EOSs to pass the gravitational-wave constraint with their extremely small Λ values. For these EOSs, the GW170817 event would have been a merger between a twin star and a regular NS, or between two twin stars, as both of these scenarios lead to small binary tidal de-formabilitiesΛ ≪ 720. We finally remark that in this case it is not clear how these twin stars would form, as they have much smaller gravitational and baryonic masses than the highest-mass NSs along the same EOS. Hence, during their formation from a NS a large amount of material ∼ 0.5M ⊙ would need to be ejected. In addition, it is unclear whether such twin star solutions would be consistent with the electromagnetic counterparts of GW170817.
We conclude by observing that all of our twin-star solutions feature maximal masses only marginally exceeding 2M ⊙ . While within a more general EOS construction, the maximal masses may be somewhat increased, overall these twin-star solutions are not very probable, especially given possible observations of even higher-mass NSs (Linares et al. 2018;Romani et al. 2022).

CORES OF HIGH-DENSITY MATTER
In this section, we explore the cores of high-density matter in stable NSs, which we define as the part of the star described by the second polytrope after the phase transition. We do not speculate on the nature of the high-density phase. In Fig. 10 we show the maximum possible mass of the core for given PT parameters. Unsurprisingly, the largest cores are related to transitions with a very early onset, but cores with masses exceeding 0.5M ⊙ are possible for onsets up to 5n 0 (but small ∆n).
In Fig. 11, we display the mass fraction of the core, M core /M , in the M -R diagram, further dissecting the results by the minimal value of ∆n/n PT (corresponding . The mass fraction of the core M core /M for different cuts on ∆n/n. Note that the cuts shown correspond to the red lines in Fig. 10. to PTs above the corresponding red lines in Fig. 10). As one can readily observe, excluding the twin-star solutions the largest cores originate from weak phase transitions with ∆n/n PT < 0.5 and heavier NSs, while the largest PTs lead to small cores at masses M ≳ 2M ⊙ . This reflects the fact that in our ensemble all sizable PTs begin at fairly high densities, at which point the stars are already quite massive. Note also that the low-mass twinstar configurations possess sizeable cores of high-density matter.

SUMMARY
In this work, we have explored bounds for first-order PTs that are compatible with present EOS constraints from chiral EFT up to nuclear densities, current NS observations, and pQCD results at very high densities. To this end, we constructed a large ensemble of EOSs with a two-polytrope form, always featuring a first-order PT in between, demanding that the lower polytrope connects to the chiral EFT band from Hebeler et al. (2013), the higher polytrope connects to pQCD constraints at n = 10n 0 from Annala et al. (2022) or Komoltsev & Kurkela (2022), and that all physical EOSs are causal, support 2M ⊙ NSs, and are consistent with the tidaldeformability from GW170817. We find that within our PT setup all three constraints -chiral EFT, astrophysics, and pQCD -provide unique constraints to the set of allowed EOSs.
Our derived region of possible PTs covers nearly all individual EOSs considered in previous works on PTs and their inclusion in merger simulations, but is significantly larger and includes several scenarios that have not yet been explored in merger simulations. Moreover, we find that first-order PTs can extend the radius ranges of NSs to larger radii, but these are only compatible with astrophysics constraints if the PT begins early. Note that these solutions are also broadly consistent with the radius constraints from NICER (Miller et al. 2021;Riley et al. 2021). Our EOS ensemble also includes twin-star solutions, although their parameter space is very small, with the solutions just barely passing current astrophysical constraints. As expected, sizable cores, i.e., with a significant part of matter described by the second polytrope beyond the phase transition, mainly occur in stars with large masses.
Let us finally comment that the EOS ansatz used in this work is exploratory, so it is to be expected that the derived bounds for the EOSs and PT parameters can increase somewhat if one allows for more freedom in the EOS. Nevertheless, our systematic study with explicit PTs points to important qualitative features for NS properties possible with strong PTs.