Phase transition phenomenology with nonparametric representations of the neutron star equation of state

Astrophysical observations of neutron stars probe the structure of dense nuclear matter and have the potential to reveal phase transitions at high densities. Most recent analyses are based on parametrized models of the equation of state with a finite number of parameters and occasionally include extra parameters intended to capture phase-transition phenomenology. However, such models restrict the types of behavior allowed and may not match the true equation of state. We introduce a complementary approach that extracts phase transitions directly from the equation of state without relying on, and thus being restricted by, an underlying parametrization. We then constrain the presence of phase transitions in neutron stars with astrophysical data. Current pulsar mass, tidal deformability, and mass-radius measurements disfavor only the strongest of possible phase transitions (latent energy per particle $\gtrsim 100\mathrm{MeV}$). Weaker phase transitions are consistent with observations. We further investigate the prospects for measuring phase transitions with future gravitational-wave observations and find that catalogs of $O(100)$ events will (at best) yield Bayes factors of $\sim 10:1$ in favor of phase transitions even when the true equation of state contains very strong phase transitions. Our results reinforce the idea that neutron star observations will primarily constrain trends in macroscopic properties rather than detailed microscopic behavior. Fine-tuned equation of state models will likely remain unconstrained in the near future.

Figure 1

(Left) one-dimensional 90% symmetric marginal posterior credible regions for the radius as a function of mass conditioned on current data. We show results with only pulsar masses (denoted PSR) and pulsar masses, GW observations, and NICER x-ray pulse profiling (denoted PGX). We additionally show maximum-likelihood EOSs from subsets of the prior conditioned on the size of the latent energy per particle $\mathrm{\Delta }(E/N)$ of phase transitions that overlap with the central densities of NSs between $1.1–2.3M_{\odot }$ [small: $\mathrm{\Delta }(E/N)\le 10\mathrm{MeV}$ and large: $\mathrm{\Delta }(E/N)\ge 100\mathrm{MeV}$]. (Right) correlations between the radius at two reference masses, $M=1.4$ and $2.0M_{\odot }$. While the one-dimensional marginal distributions are similar, EOSs with small $\mathrm{\Delta }(E/N)$ show stronger correlations between $R_{1.4}$ and $R_{2.0}$ than EOSs with large $\mathrm{\Delta }(E/N)$. This is because the radius can change rapidly when $\mathrm{\Delta }(E/N)$ is large, as is evident in the maximum-likelihood EOS.

Figure 2

Examples of CSS EOSs based on DBHF [66] with a causal extension ($c_s=c$) beyond the end of the phase transition. We show examples with (top) weak and (bottom) strong phase transitions, defined by whether there are multiple stable branches. For each EOS, we show (top left) the pressure and (bottom left) the sound-speed as a function of baryon density, (top center) the moment of inertia and (bottom center) the novel feature introduced in Sec. 2b [Eq. (2)] as a function of gravitational mass, and (top right) the $M\text{−}\mathrm{\Lambda }$ and the (bottom right) $M\text{−}R$ relations. Stable (unstable) branches are shown with dark solid (light dashed) lines. Each curve is labeled with connections between macroscopic phenomenology and microphysical features. (black annotations) The maximum mass of cold, nonrotating stars ($M_{\mathrm{TOV}}$) and, where relevant, the beginning and end of stable branches. (red annotations) The beginning and end of features as identified by the procedure in Sec. 2b. (red shading) The extent of the identified features.

Figure 3

Analogous to Fig. 2 but for more complicated phase transition phenomenology associated with mixed phases (Gibbs construction) from Han et al. [67], obtained by implementing specific hadronic and quark models. Again, the features introduced in Sec. 2b correctly identify the beginning and end of the phase transition even though there is no discontinuity in $c_s$ at the onset and the phase transition corresponds to a wide range of masses. The broad extent of the phase transition is not readily apparent from the macroscopic properties alone, which show a sharp feature only at the end of the phase transition.

Figure 4

The feature extraction algorithm: (left) the sound-speed as a function of baryon density and (right) $\arctan ({\mathcal{D}}_M^I)$ [Eq. (2)] as a function of the gravitational mass. The algorithm progresses from top to bottom, first with the identification of local minima in $\arctan ({\mathcal{D}}_M^I)$ and then pairing each with a corresponding running local maximum in $c_s$. The number of features reported corresponds to the number of unique running local maxima in $c_s$ selected; in this case 1. The multiplicity of each feature corresponds to the number of local minima in $\arctan ({\mathcal{D}}_M^I)$ that are paired with the same running local max in $c_s$; in this case 3. For demonstration purposes, we show how the algorithm would progress if we had $R_{c_s^2}>1.7$. If the threshold on the drop in the sound-speed $R_{c_s^2}$ was $\le 1.7$, the algorithm would accept the first pairing (second row) and instead report two features: one at lower densities with multiplicity two and one at higher densities with multiplicity one. This would be the case for the main results presented in Secs. 3 and 4, which use a threshold $R_{c_s^2}=1.1$.

Figure 5

Correlations between the divergence between macroscopic properties caused by a phase transition $\mathrm{\Delta }\ln I-⟨\mathrm{\Delta }\ln I⟩$ and the latent energy per particle of the associated phase transition $\mathrm{\Delta }(E/N)$ for all transitions that begin at masses greater than $0.7M_{\odot }$. Color indicates the proximity of the phase transition’s end to $M_{\mathrm{TOV}}$. Large divergences in macroscopic properties can only be caused by phase transitions with large $\mathrm{\Delta }(E/N)$, but not all phase transitions with large $\mathrm{\Delta }(E/N)$ cause large divergences in macroscopic properties.

Figure 6

Marginalized (unshaded) priors and (shaded) posteriors for parameters that characterize phase transitions based on current astrophysical data from pulsar masses, GWs, and x-ray mass-radius measurements. For each EOS we report the properties of the transition with the largest $\mathrm{\Delta }(E/N)$ that overlaps with each mass interval. We report (left to right), the latent energy ($\mathrm{\Delta }(E/N)$), the onset energy density (${ϵ}_t$), the onset pressure ($p_t$), the energy density at the end of the transition (${ϵ}_e$), and the onset mass scale ($M_t$) for three mass-overlap regions: $0.8–1.1M_{\odot }$, $1.1–1.6M_{\odot }$, and $1.6–2.3M_{\odot }$.

Figure 7

Ratios of probabilities conditioned on different numbers of features. Compare to Table 1; see Eqs. (3) and (4) for an explicit definitions of our notation. (Left) distributions over the number of stable branches and (right) distributions over the number of ${\mathcal{D}}_M^I$ features for EOSs with $\mathrm{\Delta }(E/N)\ge 10$, 50, and 100 MeV, respectively for different mass-overlap regions: (top) $0.8–1.1M_{\odot }$, (middle) $1.1–1.6M_{\odot }$, and (bottom) $1.6–2.3M_{\odot }$. We show the ratio of maximum likelihoods (black dots) and the posterior divided by the prior (circles and x’s). As in Table 1, we consider (PGX, red circles) the ratio of the posterior conditioned on PSR masses, GW coalescences, and x-ray timing and compare it to our nonparametric prior as well as (blue x’s) the posterior conditioned on only PSR masses. Error bars approximate $1\sigma$ uncertainties from the finite size of our prior sample. In general, a single stable branch without strong ${\mathcal{D}}_M^I$ features is preferred.

Figure 8

Bayes factors vs catalog size comparing (left-most column) multiple stable branches vs a single stable branch and (right three columns) at least one ${\mathcal{D}}_M^I$ feature vs no ${\mathcal{D}}_M^I$ features. We consider features that overlap with three mass ranges: (top row) $0.8–1.1M_{\odot }$, (middle row) $1.1–1.6M_{\odot }$, and (bottom row) $1.6–2.3M_{\odot }$. We also show three different injected EOSs: (blue, no phase transition) DBHF, (orange, weak phase transition at $\sim 1.9M_{\odot }$) DBHF_3504, and (green, strong phase transition at $\sim 1.5M_{\odot }$) DBHF_2507. Shaded regions denote $1\sigma$ uncertainties from the finite size of our Monte Carlo sample sets. Different realizations of catalogs will also produce different trajectories; these should only be taken as representative.

Figure 9

Sequences of one-dimensional marginal posteriors for $\mathrm{\Lambda }(M)$ at (left to right) 1.2, 1.4, 1.6, 1.8, and $2.0M_{\odot }$ for different simulated EOSs: (top, blue) DBHF, (middle, orange) DBHF_3504 (phase transition at $\sim 1.9M_{\odot }$) and (bottom, green) DBHF_2507 (phase transition at $\sim 1.5M_{\odot }$). These posteriors show the distributions of $\mathrm{\Lambda }(M)>0$ (i.e., they only consider EOSs with $M_{\mathrm{TOV}}\ge M$). These posteriors are conditioned only on simulated GW events (no real observations), and a line’s color denotes the number of simulated GW events within the catalog (light to dark: fewer to more events) along with the true injected values (vertical black lines). The prior is shown for reference (gray shaded distributions). For very small $\mathrm{\Lambda }$, primarily associated with DBHF_2507 at high masses, the true value falls near the lower bound in the prior. The primary effect of additional observations is to reduce support for larger values of $\mathrm{\Lambda }$. While significant uncertainty in $\mathrm{\Lambda }(M)$ remains after 100 events, the nonparametric prior is able to correctly infer $\mathrm{\Lambda }(M)$ at all $M$ simultaneously, including sharp changes in $\mathrm{\Lambda }(M)$ over relatively small mass ranges.

Figure 10

Joint posteriors for $\mathrm{\Delta }(E/N)$ and transition onset mass ($M_t$) inferred from simulated GW catalogs for (left, blue) DBHF and (right, green) DBHF_2507. Gray curves denote the (reweighed) prior, color denotes the size of the catalog, and contours in the joint distribution are 50% highest-probability-density credible regions. Solid lines denote the true parameters for DBHF_2507; there are no such lines for DBHF because it does not contain a phase transition. As in Fig. 6, extracted parameters correspond to the feature with the largest $\mathrm{\Delta }(E/N)$, but here we only require features to overlap the broad range $0.8–2.3M_{\odot }$.

Figure 11

Stellar sequences for incompressible two-phase Newtonian stars with ${\rho }_L=2{\rho }_{\mathrm{sat}}=5.6\times {10}^{14}\mathrm{g}/{\mathrm{cm}}^3$, $p_T=5\times {10}^{34}\mathrm{dyne}/{\mathrm{cm}}^2$, and various values of ${\rho }_H$. We plot (top) the $M\text{−}I$ relation and (bottom) $\arctan ({\mathcal{D}}_M^I)$ as a function of the stellar mass. Stable branches are shown with solid lines, and unstable branches are shown with dotted lines. The bottom panel inset focuses near the discontinuity for curves with; ticks on the y-axis correspond to the values in Eq. (a11).

Figure 12

An additional example of the impact of thresholds within the feature extraction algorithm with an EOS realization with a relatively short correlation length. (Top) trivial thresholds; (middle) threshold on the size of $\mathrm{\Delta }\arctan ({\mathcal{D}}_M^I)$; (bottom) threshold on the amount $c_s^2$ must decrease (analogous to Fig. 4). The rapid oscillations in $c_s^2$ are identified when selecting based on ${\mathcal{R}}_{c_s^2}$ but they are rejected when selecting based on $\mathrm{\Delta }\arctan ({\mathcal{D}}_M^I)$; their relatively small $\mathrm{\Delta }(E/N)$ do not produce significant changes in the $M\text{−}I$ relation.

Figure 13

The effective number of EOS samples from the posterior process as a function of catalog size for (solid) catalogs comprised of only mock GW observations and (dashed) catalogs that include real pulsar mass measurements in addition to mock GW observations. For each of the three true EOS considered in Sec. 4, we find an approximately exponential decrease of the number of effective samples with the catalog size.

Figure 14

Distributions of radii and tidal deformabilities at reference masses as well as $M_{\mathrm{TOV}}$ conditioned on current data. These distributions de facto exclude EOSs with $M_{\mathrm{TOV}}<2M_{\odot }$ by requiring ${\mathrm{\Lambda }}_{2.0}>0$ (enforced through the logarithmic scale). As in Fig. 1, there are much weaker correlations between low-mass and high-mass observables.

Figure 15

An additional example of an EOS with mixed phases (Gibbs construction) from Han et al. [67], analogous to Fig. 3.

Figure 16

Several realizations from our nonparametric prior, each with a single stable branch but with different numbers of phase transitions.

Figure 17

Additional realizations from our nonparametric prior, each with multiple stable branches. Typically, we always identify a phase transition associated with the loss of stability between stable branches, even if the stable branches are small (bottom row).