Implicit correlations within phenomenological parametric models of the neutron star equation of state

The rapid increase in the number and precision of astrophysical probes of neutron stars in recent years allows for the inference of their equation of state. Observations target different macroscopic properties of neutron stars which vary from star to star, such as mass and radius, but the equation of state allows for a common description of all neutron stars. To connect these observations and infer the properties of dense matter and neutron stars simultaneously, models for the equation of state are introduced. Parametric models rely on carefully engineered functional forms that reproduce a large array of realistic equations of state. Such models benefit from their simplicity but are limited because any finite-parameter model cannot accurately approximate all possible equations of state. Nonparametric models overcome this by increasing model freedom at the cost of increased complexity. In this study, we compare common parametric and nonparametric models, quantify the limitations of the former, and study the impact of modeling on our current understanding of high-density physics. We show that parametric models impose strongly model-dependent, and sometimes opaque, correlations between density scales. Such interdensity correlations result in tighter constraints that are unsupported by data and can lead to biased inference of the equation of state and of individual neutron star properties.

A set of observations of different systems can be used to constrain a shared underlying property through a hierarchical inference scheme. The hierarchical formalism is derived in the context of combining data from different sources while faithfully incorporating their uncertainties and potential observational selection effects [20]; see, e.g., Refs. [7,21]. In the context of NS structure, the main objective is to obtain a posterior for the EoS as a shared variable among many astrophysical observations. The prior corresponding to this posterior is not necessarily straightforward to define because the space of potential EoS (i.e., the space of possible functions obeying basic physical constraints) that relate the pressure p and the baryon density ρ, p = p(ρ), is infinite dimensional. 1 The simplest way to define such a prior is through a parametrization of the EoS, which is a functional form of p(ρ) that typically depends on a few parameters. 1 We can equivalently use p(ε), with ε the internal energy density.
Common phenomenological models such as piecewisepolytrope [22], spectral [19,23] and speed-of-sound [24] parametrizations have been used to effectively sample candidate EoS for use in inference. The simplicity of a closed-form parametric expression comes at the cost, though, of being unable to faithfully represent many of the possible degrees of freedom in the true EoS. While many of these models can accurately represent most EoS derived from effective nuclear interactions [19,22], it is not always clear how to extend these parametrizations toward more general behavior in the EoS that may arise from phase transitions or new physics. This limitation of phenomenological parametric models has been recognized from the outset [22]. However, in this study we investigate another way that they may artificially restrict the inferred EoS. Parametric models use only a few parameters, which means that the values of p(ρ) at different densities are often correlated. These correlations represent a source of model dependence in the inference, which is undesirable insofar as it does not reflect true prior knowledge of the EoS at those densities. That is, the correlations induced by the choice of parametrization can constitute strong, unintentional prior beliefs about the EoS. This unwanted model dependence is a natural consequence of the phenomenological nature of the parametric models.
An alternative method for constructing a prior on the space of EoS, which we call nonparametric in what follows, targets more model flexibility by making use of Gaussian processes (GPs) [25]. This approach produces a multivariate Gaussian distribution for the function φ = log (c/c s ) 2 − 1 , where c s is the speed of sound and c is the speed of light. By conditioning the prior only weakly on existing nuclear-theory models, we generate a model-agnostic prior process for EoS. The chosen correlations between φ(p i ) and φ(p j ), or equivalently between the values of the EoS at different densities, are set by a kernel function, which is in turn described by a few parameters. Following Ref. [26], we consider a variety of possible kernel parameters to probe a range of different correlations and thus maximize model freedom. This approach allows us, in principle, to model any function p(ρ), and furthermore to probe a wide range of interdensity correlations and high-density EoS behavior.
Of course, completely unrestricted freedom in the EoS is neither desirable nor realistic, as certain physical constraints should be encoded into the EoS prior. For example, an EoS must be causal, and thermodynamically stable, dp dε = c 2 s > 0. (2) Imposing these constraints in the prior is desirable as it excludes unphysical models from the analysis. 2 In this paper, we examine common parametric and nonparametric EoS models to determine the extent to which each prior's assumptions impact inference of the EoS and NS properties. We find that the three parametric models we study (spectral, piecewise polytrope, and speed of sound) build additional interdensity correlations into the EoS beyond what can be attributed to causality and stability. These correlations between densities typically lead to more stringent constraints than are strictly supported by the data. On the other hand, the nonparametric model demonstrates the largest degree of model independence, restricted primarily only by causality and thermodynamic stability. We demonstrate that these strong, model-dependent interdensity correlations have already impacted inferred microscopic and macroscopic NS properties. Such effects are expected to become more severe as statistical uncertainties decrease with more data that probe different NS densities.
The remainder of the paper is organized as follows. In Sec. II, we describe our inference methods and our approach to investigating model dependence. In Sec. III, we examine the EoS and NS properties inferred with current data and show that they are influenced by correlations in the EoS prior. In Sec. IV, we illustrate the main limitations of parametric EoS inference with a toy model. In Sec. V, we quantify the implicit EoS correlations and demonstrate that the nonparametric model displays the largest degree of model independence. In Sec. VI, we study the correlations' potential impact on upcoming EoS inference using mock astrophysical measurements. In Sec. VII, we demonstrate that the limitations identified in parametric models cannot be resolved by making small modifications to the prior distributions. Finally, in Sec. VIII we discuss our conclusions.

II. METHODS AND MODELS
The posterior for the EoS depends on two elements: (i) the prior EoS process, and (ii) the data. Our goal in this study is to assess the effect of the prior as generated from different parametric and nonparametric models for the EoS. We therefore always employ the same data, which we briefly describe in Sec. II B. The hierarchical likelihood corresponding to this data is described in Refs. [21,27]. The EoS priors are described in detail in Sec. II A, where we discuss the different EoS models and parameter priors that generate each EoS prior process.

A. EoS prior
We wish to establish a prior process over candidate EoS. By this, we mean a probabilistic measure on the space of potential EoS. To do this, we use several models of the EoS. We distinguish parametric models, which provide a functional form for the EoS, from nonparametric models which do not impose such a functional form. We use three different phenomenological parametric models, a piecewise-polytrope [22] parametrization, a spectral parametrization [19], and a direct parametrization of the speed of sound [24]. The spectral and piecewisepolytrope parametrizations use a polytropic form for the EoS, so that In the piecewise-polytrope case, the polytropic index Γ is a piecewise-constant function of the pressure, while in the spectral case, log(Γ) is expanded as a polynomial in pressure. In the speed of sound parametrization, the speed of sound is expressed as a constant plus a Gaussian and a logistic curve which asymptotes to c 2 /3. Following past practice [28], we slightly relax the causality threshold and consider EoS with c 2 s < 1.1c 2 for all parametric models. See Appendix B for more details about each model and its implementation.
To establish a prior process, we must additionally supply a joint prior probability distribution on the parameters of each model from which a draw is a realization of the parameters and therefore a candidate EoS. For example, in the spectral model, the parameters are coefficients in the spectral expansion. In the piecewise polytrope, the parameters are the value of the polytropic index itself. For our headline results, we use standard priors for the parametric models [28,29], except for the speed-of-sound model, which we adapt to increase access to astrophysically relevant EoS; again, see Appendix B for details.
We compare the prior processes generated by the parametric models to a prior process from a nonparametric model [25]. While our nonparametric implementation does not assume a specific functional form for the EoS, it does parametrize the correlations between the sound speed at different densities. These correlations are described by a kernel function. In practice, we choose a large set of points, p i , and then the variable φ(p i ) is sampled from a multivariate Gaussian distribution. By changing the kernel's parameters and conditioning on different nuclear models, we can generate a range of GPs. We choose a model-agnostic prior, which is to say we average over multiple GPs with different correlations, each loosely informed by nuclear-theory models [25]. We do this to maximize the freedom of the model. See Appendix A for more details.
Due to its construction, the GP itself has parameters which control correlations. Such parameters have been termed hyperparameters [25], though we avoid this terminology here in order to avoid potential confusion with the term's use in hierarchical inference. In addition, the parametric models also have parameters which control the prior process; in general, such details are unique to each model. We instead focus primarily on the prior process induced by each EoS model with its chosen prior, returning to the subject of parameter distributions briefly in Sec. VII. For now, we simply note that our model implementations are typical of those used in the literature [19,24,28,30,31]. Lastly, we stress that our distinction between parametric and nonparametric models lies not in the existence of parameters but in the specification of a functional form for the EoS. In particular, we compare models with small, fixed numbers of parameters that are commonly used in the literature. For these models, the choice of functional form significantly impacts the range of EoS that can be represented.

B. Data and likelihood
Unless otherwise stated, all analyses in this paper make use of the same astronomical data as Ref. [27]. Specifically, we include two mass-tidal deformability measurements from gravitational wave (GW) detections of merging NSs [32,33], one heavy pulsar mass measurement with radio data [34], and two x-ray observations of NS masses and radii [10][11][12][13]. In the latter case, we also use the up-to-date radio mass measurement of the pulsar J0740+6620 [35]. Given these data d, the posterior probability density of a particular EoS ε is where I is any additional information we may have about the system, e.g., knowledge that the data originate from a NS, as is the case for the pulsar observations but not for the GWs. Here P (ε|d, I) is the posterior probability of the EoS given the data, P (d|ε, I) is the likelihood of the astrophysical data given the EoS, P (ε|I) is the prior probability of the EoS, and P (d|I) = P (d|ε, I)P (ε|I)Dε is the total probability of observing this data marginalized over all EoS in the prior, often called the evidence. For general astrophysical data, P (d|ε, I) must be computed by marginalizing over the astrophysical distribution of masses, spins, sky locations, and distances for individual events, which remains poorly constrained [36][37][38][39][40][41][42]. For the full expression, see Refs. [7,21]. The different datasets we use primarily inform the EoS at different densities. The heaviest pulsar mass measurements serve to downweight EoS which cannot support the observed NS masses; these constraints tend to most significantly impact inference near ∼ (4-6)ρ nuc and typically favor a stiffer EoS. The x-ray data provide constraints on the NS radius, and constraints so far have given information about the EoS mainly in the region ∼ (1-4)ρ nuc [12,21,27,43]. The GW observations provide constraints on the tidal deformabilities of the binary components, which are dominated by the loudest event observed so far, GW170817 [8,21]. In terms of densities, the relevant scale constrained by this measurement is ∼ (1-3)ρ nuc [21,27]. Future constraints with GWs are likely to lie in this density range, as the fractional uncertainty in Λ will be smallest for lower-mass NSs with less dense cores and larger tidal deformabilities. In principle, nuclear experiments or calculations could also be included in such an analysis, and would mainly constrain the EoS near or below ρ nuc [43][44][45][46][47]. However, we do not incorporate any in this work.

III. IMPACT OF EOS MODEL ON CURRENT EOS CONSTRAINTS
Following the above prescription, we analyze the existing data using the four different EoS priors and plot the resulting marginal posteriors for p(ρ) across a wide range of densities. Figure 1 compares the spectral and nonparametric models; similar plots for the other parametric models can be found in Appendix C. The posteriors differ in their predictions for the EoS. For instance, the spectral posterior is stiffer on average than the nonparametric one, especially above 4ρ nuc . Similar differences have also been pointed out in Refs. [12,24,43], where multiple EoS models were employed under identical analysis settings. Our goal here is to understand the origin of these discrepancies.
Since it is difficult to glean information about interdensity correlations from envelope plots like Fig. 1, we turn our attention to two macroscopic NS properties that roughly correspond to the EoS behavior at high and low densities: the maximum mass, M max , and the radius of a 1.4M NS, R 1.4 . In Fig. 2 Figure 1. Symmetric 90% credible region for the pressure p at each density ρ in units of the nuclear saturation density using the nonparametric and spectral prior processes. We show results including all astrophysical data (labeled "astro," solid lines) and restricting to the heavy pulsars only (labeled "psr," dashed lines). The latter choice ensures that prior choices on the Mmax supported by each model are irrelevant. Other parametric models are shown in Fig. 10. In all cases, we find that the p-ρ posterior depends on the EoS model even when identical data and inference schemes are employed.
plot shows that both the spectral prior and posterior seem to have more support for M max around 2.2−2.5M than the nonparametric case. However, this trend is reversed above 2.6 M . This observation suggests that the difference between the nonparametric and spectral posteriors cannot be trivially assigned to different marginal priors. To further demonstrate this, in the right panel we plot the same variables, but now reweighted to a flat marginal M max prior. As expected, the two posteriors differ, with the spectral model producing a narrower posterior.
To understand this discrepancy, we revisit the possible reasons M max is constrained on the high side. Though an upper limit on M max has been proposed based on the analysis of the counterpart of GW170817 [50], our analysis does not make use of it. Excluding an origin due to data, the upper limit on M max must be the result of the EoS prior. The two-dimensional M max -R 1.4 panel indeed shows that the upper limit on M max is related to the upper limit on R 1.4 [8,32]; the M max -R 1.4 prior does not cover the entire available region for either model, with larger M max requiring stiffer EoS and larger R 1. 4 .
The fact that larger M max requires large values of R 1.4 is not unexpected from causality considerations. Indeed, the causality condition and the pressure at twice saturation, p 2.0 , set an upper limit on the value of pressure at five times saturation, p 5.0 . Since p 2.0 and p 5.0 correlate with R 1.4 and M max , respectively [17], any causal EoS model should limit M max for certain low R 1.4 configurations [48]. To quantify this, we overplot the limiting M max -R 1.4 relation given by Ref. [49]: each point on the line represents a soft low-density EoS stitched to an EoS with c 2 s = 1 at different densities. This curve should be interpreted approximately, as the exact causality threshold depends on the details of the low-density EoS [51].
Nonetheless, the right panel shows that in the nonparametric case the prior fills more of the physically allowable M max -R 1.4 parameter space compared to the spectral prior. This indicates that the nonparametric prior has non-negligible support in the entire physically allowed region even if specific marginal priors might downweight some regions (left panel). The same is not true for the spectral model which cannot access certain regions of the M max -R 1.4 plane. 3 This demonstrates that the correlation between M max -R 1.4 that appears in the spectral model is not entirely due to causality considerations, but it is also affected by the specifics of the model. The nonparametric model, however, is able to produce EoS which fall near the causality limit, indicating physics rather than modeling artifacts are the primary limitation to model freedom. Figure 12 presents a qualitatively similar conclusion for the piecewise-polytrope and speed-ofsound models. The piecewise polytrope exhibits behavior similar to the spectral model, while the speed-of-sound parametrization exhibits the opposite problem: its prior does not include low M max for large R 1.4 values.
Crucially, the correlations we see in these corner plots are only those that are apparent from the twodimensional marginalized posteriors. They do not reveal the many hidden correlations within the parametric EoS models that are not as easily detected. It is possible for implicit correlations within the EoS prior to bias the inference in ways that are not obvious in low-dimensional projections.

IV. IMPACT OF INTERDENSITY CORRELATIONS: TOY MODEL
To better understand the effect of such implicit correlations, we first consider a simple toy model that demonstrates several of the issues with parametric models and introduce our techniques for diagnosing them.
We consider several simple linear parametrizations of the pressure as a function of energy density p(ε). This allows us to examine the prior processes induced by the assumption of linearity from various perspectives. We contrast this to a GP prior process in the same context, finding particularly striking differences in the effect of a precise measurement of the pressure at one density on our uncertainty in the pressure at other densities.
We begin with the simple parametric model of a linear relationship between the pressure and the energy density  . We show results with the nonparametric and spectral EoS models, and contours denote 90% credible regions. The black line in the two-dimensional plot represents a maximally stiff M -R curve [48] stitched to a fiducial low density EoS [49]. Due to the low-density stitching, this "causality" line should be interpreted as a fuzzy boundary and not a sharp line. Both panels demonstrate that the spectral and nonparametric EoS models produce different Mmax posteriors and that these differences cannot be attributed to the marginal priors. They are instead caused by correlations between low and high densities (equivalently, between Mmax and R1.4) imposed by the models. The correlations in the nonparametric case are due to causality, while the spectral case exhibits additional correlations and model dependence.
parametrized by the pressure at p a = p(ε a ) and the slope c 2 s . Ignoring causality constraints, we choose what appear to be uninformative priors and refer to this as the point+slope process. The top panel of Fig. 3 shows the envelope plot for this prior process, i.e., the marginal distributions of the pressure at each energy density.
The envelope plot appears reasonable. That is, the prior process assigns approximately equal uncertainty to each pressure. However, the envelope plot only shows the marginal distributions at each density. Figure 4 shows the correlations between pressures at different densities. From this we see that the prior process actually imposes strong correlations between all pressures. We quantify this correlation between p a and the pressure at some other density p b ≡ p( b ) with the mutual information [52], defined as where P (p a ) = dp b P (p a , p b ) is the marginal distribu-tion. For the point+slope parametrization, we compute The mutual information between p a and p b can be made arbitrarily small only in the limit σ a σ c 2 s |ε b −ε a |; however, this limit corresponds to vanishingly small marginal uncertainty for p a . We conclude that the assumption of a linear functional form can produce what seems to be a reasonable envelope plot in Fig. 3, but nevertheless induces model-dependent correlations between the pressure at different densities in Fig. 4.
In an attempt to remove the correlation between p a and p b , we consider the alternative parametrization described by the pressures at the two reference densities. We assume priors and refer to this as the two-point process. Figures 3 and 4 show envelope and marginal distributions, respectively. By construction, I(p a , p b ) = 0 for the two-point prior process, as also seen in Fig. 4. However, the envelope plot shows that the marginal prior actually tightens for pressures between the reference densities. That is, we are able to remove the correlation between two pressures only at the expense of asserting greater prior knowledge about other pressures. Additionally, Fig. 4 shows that there are still correlations between (p a , p b ) and other pressures. This hints at the fact that, when one assumes a specific functional form, it may be possible to remove the correlations between a small number of statistics, but it is generally difficult to make all correlations vanish simultaneously or to avoid making strong assumptions about specific values of the function.
In order to consider this effect more quantitatively, we introduce a generalization of the mutual information that considers three pressures [52] I(p a , p b , p c ) ≡ dp a dp b dp c P (p a , p b , p c ) ln P (p a , p b , p c ) P (p a )P (p b )P (p c ) = dp a dp b P (p a , p b ) dp c P (p c |p a , p b ) ln Even if one can choose parametrizations and priors such that I(p a , p b ) vanishes, there is another term when considering mutual information for three pressures. In fact, for both the point+slope and two-point prior processes, the integral over p c diverges as P (p c |p a , p b ) is a delta function (determined by the closed-form parametrization), while P (p c ) is a Gaussian with finite width. We conclude that the assumption of a linear relationship between the pressure and the density implies an infinite amount of information about the allowed relationships between variables. One cannot undo all these correlations at the same time by a clever choice of marginal prior distributions, although it may be possible to undo some of them. The failure of this reparametrization scheme anticipates the results of our investigation of alternative parametric priors in Sec. VII.
In general, the only way to undo all correlations simultaneously is to add more model freedom into the prior process. For example, one may add more reference densities to an existing model and generate a piecewise linear prior process. However, there will always be some densities between the (finite number of) reference densities, regardless of how many reference densities are chosen. In each of those regions, the piecewise-linear model is equivalent to our two-point prior process, and the strong correlations remain. One is then left with the question of how to extend the parametrization to remove all correla-tions in a scalable way. We show that a GP is a natural solution.
A GP, defined in terms of a mean function and a covariance kernel, describes our uncertainty in the infinitely many degrees of freedom in a function. With the assumption of Gaussianity, we can easily marginalize away uninteresting degrees of freedom, in our case retaining only the pressures on a dense grid of energy densities, and the GP reduces to a high-dimensional multivariate Gaussian distribution. Specifically, we consider the joint distribution induced over p a , p b , and p c by a GP with mean µ and covariance Σ, with matrix elements defined by a covariance kernel A common choice is the squared-exponential kernel although more complicated kernels are also used [26]. 4 Figures 3 and 4 show a GP assuming a squaredexponential kernel with parameters chosen to match the marginal distribution of the two-point prior process and l |ε a − ε b |. We also obtain which vanishes as exp[−2(ε a − ε b ) 2 /l 2 ] in the limit |ε a − ε b | l. Furthermore, P (p c |p a , p b ) is a normal distribution, and the generalization of the mutual information in Eq. (10) Therefore, if l |ε b − ε a | as well, I(p a , p b , p c ) → 0. We conclude, then, that the GP prior process can be made to simultaneously produce reasonable envelope plots (broad marginal distributions for all pressures) while retaining vanishingly small correlations between (reasonably separated) pressures. This is in stark contrast to the parametrized prior processes, where this is, in general, not possible.
We demonstrate one more useful diagnostic in this toy model through the conditioned distribution 4 It is worth noting that linear regression is a special case of a GP. That is, a GP can reproduce the linear model with an appropriate choice of covariance kernel: which shows how our knowledge of p i depends on p j . Figure 3 shows the envelope plots for the conditioned distributions corresponding to each of our prior processes when we condition on p c . We see that a constraint at ε c is broadcast to nearby densities in all cases, but the Gaussian process fills up the prior volume from the unconditioned marginal distributions the fastest. This is a visual manifestation of the correlations quantified by the mutual information. Indeed is just the Kullback-Leibler divergence D KL (P (b|a)||P (b)) from the unconditioned marginal to the conditioned marginal averaged over the possible a ∼ P (a).
In the case of realistic EoS inference, we have to consider even higher-dimensional spaces. A natural measure of correlations, then, is a generalization of the mutual information (sometimes called the total correlation, multivariate constraint, or multi-information [52]) (17) where H(x) = − dxP (x) ln P (x) is the entropy of the distribution P (x). Larger H imply broader distributions. We will consider this statistic in the context of real astrophysical constraints on the EoS in Sec. V. In general, one can make I small but still allow for very little model freedom (small H). We therefore seek prior processes with both large H(x 1 , · · · , x N ) and small I(x 1 , · · · , x N ). This is sometimes captured in the variation of information, defined as H − I, but we find it more useful to consider H and I separately.

V. IMPACT OF INTERDENSITY CORRELATIONS: IDEALIZED MEASUREMENT
Our toy model illustrates the potential impact of implicit correlations within EoS models on the results of EoS inference. In order to quantify the sensitivity of the parametric and nonparametric models to such modeldependent correlations between density scales, we now consider simulated NS observations. The macroscopic NS properties that astronomical observations target (masses, radii, and tides) are determined by a range of NS densities, it is therefore not straightforward to disentangle the effect of the data and the model dependence in the constraint that a single astronomical observation imposes on the EoS. Consequently, we begin with the same setup as Sec. IV: an idealized direct measurement of p(ρ) at a single density, while keeping in mind that a realistic astronomical measurement would correspond to a combination of many such constraints correlated across many  Figure 5. Similar to Fig. 1 but with a mock constraint injected directly into p(ρ = 2ρnuc) for each EoS prior process. The posterior after the simulated constraint is included ("astro+mock") is overplotted on the posterior with all current data ("astro"), and in a darker color. The constraint at a single density affects the parametric posteriors over a much wider range of density scales than the nonparametric one. The inset focuses around ρ = 2ρnuc. The two black straight lines provide an estimate of the constraints imposed by causality (c 2 s < 1) and thermodynamic stability (c 2 s > 0.1) around ρ = 2ρnuc, subject to the heavy pulsar measurements. The nonparametric posterior quickly "fills" more of the physically available region after satisfying the mock constraint, while the parametric posteriors do not. densities. 5 We consider a tight Gaussian constraint at p 2.0 with mean of 3.20 × 10 34 dyn/cm 2 based on a candidate EoS drawn from our GP prior that is consistent with all current parametric posteriors near 2ρ nuc . We arbitrarily choose the standard deviation, 2.61 × 10 33 dyn/cm 2 (∼ 8% relative uncertainty). We then plot the corresponding envelope for p(ρ) with this mock constraint and all other real astronomical data for each model in Fig. 5. In the nonparametric case, imposing this constraint pinches the p-ρ envelope around 2ρ nuc , but the uncertainty in the p(ρ) curve is unaffected beyond ≈ ±0.5ρ nuc . All the parametric models, though, change across several ρ nuc , indicating that the EoS at many scales is informed significantly by the EoS near 2ρ nuc . 5 An example of how one may obtain direct constraints on the pressure from nuclear experiments is demonstrated in Refs. [44,45].
In each panel, the inset zooms in around the ρ = 2ρ nuc region; to guide the eye the two black lines provide a rough estimate of the maximally causal (c 2 s = 1) and minimally stable (c 2 s = 0.1) EoS that can support the heavy pulsar observations (see Fig. 2 of Ref. [21]) around ρ = 2ρ nuc . The two lines were obtained by combining dp/dε = c 2 s and the first law of thermodynamics with the approximation that ε 2.0 = c 2 ρ 2.0 . The nonparametric prior process contains EoS draws that approach this limiting behavior near the constraint. Comparing the 4 panels, the nonparametric model fills more of the physically available space. The parametric models, on the other hand, are clearly subject to additional correlations between pressures besides causality and stability.
To quantify these correlations, we follow Sec. IV and compute the total correlation between the pressures at several reference densities. Table I shows the total correlation (I) and joint entropy (H) between ln p 1.0 , ln p 1.5 , ln p 2.0 , ln p 3.0 , and ln p 4.0 induced by the posterior process conditioned on the astrophysical data as well as the astrophysical data and the mock constraint on p 2.0 . 6 We consider these pressures as the central density of M max stars may be as low as 4ρ nuc [27], and therefore we focus on pressures that are confidently relevant for NSs. Although the precise values of I and H can be difficult to interpret, we notice some trends.
Overall, the nonparametric process consistently has the largest joint entropy and smallest total correlation, as desired. The parametric processes all have approximately equal I, which are larger than the nonparametric process by 1 nat. Additionally, the change in I when we additionally condition on a mock constraint on p 2.0 is much smaller for the nonparametric than for the parametric processes. This can be interpreted as the constraint on p 2.0 removing some correlations from the parametric processes by approximately fixing the value of p 2.0 .
What is more, the parametric processes have much smaller joint entropies than the nonparametric process in all cases. This is a manifestation of the reduced model freedom in the parametric processes as the nonparametric process explores more combinations of pressures than any parametric process. Although not exact, the exponential of the difference in entropies is an estimate of the ratio of the effective number of pressure combinations supported in each distribution: the nonparametric contains between 10 and 100 times as many possible pressure combinations as the parametric processes.
Additionally, the constraint on p 2.0 removes more entropy from each parametric processes than from the nonparametric process. While we expect the joint entropy to be smaller in all cases after measuring p 2.0 precisely, the additional entropy lost in the parametric processes is associated with the correlations between pressures. That is, knowledge of p 2.0 decreases our uncertainty in other pressures within the parametric processes, something that does not happen as strongly in the nonparametric process. This is apparent in Fig. 5 as well.
As a final note, Table I also reports the entropy of the marginal distributions over ln p 2.0 . We see smaller differences between these one-dimensional (1D) distributions, reinforcing the conclusion that the differences between the nonparametric and parametric processes arise mainly from correlations between multiple pressures.

VI. IMPACT OF INTERDENSITY CORRELATIONS: MOCK ASTROPHYSICAL OBSERVATIONS
Different astronomical probes provide information about different density scales, and therefore interden- 6 We estimate the entropies via Monte Carlo sums over kernel density estimates (KDEs) of the associated distributions. As such, the actual correlations may be smoothed by the KDE, which may act as upper limits on the estimates of the mutual information in some cases. sity correlations are likely to matter even more for realistic EoS inference than in the idealized case considered above. Implicit correlations in EoS models could artificially give the appearance of tension between observations of NSs or nuclear matter made via different channels. This has already been shown to be relevant in comparisons of nuclear experiments with astrophysical observations [44,45]. To investigate this possibility, we now repeat the previous sections' analysis for a simulated set of astronomical observations. We choose a candidate EoS with M max = 2.54 M and R 1.4 = 12.0 km. This EoS is relatively soft at low densities and stiff at high densities, but is consistent with the 90% 1D marginal pressure constraints for all of our models at all densities, except the speed-of-sound parametrization above 4ρ nuc . The combination of macroscopic parameters lies outside the spectral 90% credible region in Fig. 2, motivating its use in studying how tension appears in an analysis when such a mismatch arises. 7 We simulate three measurements of pulsar masses and radii with comparable uncertainty to the recent measurement for J0740+6620 [12,13]. This observation incorporated radio data to constrain the pulsar mass [35] and constrained the radius with x-ray data. We also simulate 20 GW detections of binary NS mergers at A+ detector sensitivity [53]. Note, however, that we do not impose prior knowledge of the NS nature of the components in our inference. The simulated pulsars are drawn from a uniform-in-central-density distribution, while the simulated binary NSs come from a uniform-in-mass distribution, under the condition that the NS masses lie below M max .
We analyze each dataset separately, folding the simulated measurements of each type onto all current astrophysical data. We plot the inferred posteriors in Fig. 6. We find that the nonparametric posterior for R 1.4 is centered on the correct value (12 km) with either x-ray or GW data. In the spectral case, though, while the GW measurements are consistent with the correct value, the x-ray posterior is in tension at 90% credibility. Moreover, the GW and x-ray posteriors are less consistent with each other, an observation that could lead to the erroneous conclusion of tension between different EoS probes.
We can understand this as follows. We expect GW measurements of high-mass NSs (≥ 1.7M ) to be less informative than lower-mass NSs, as the absolute impact of tidal parameters on the signal is weaker for more compact stars. In general, high-mass NSs will most likely be indistinguishable from black holes until the advent of next-generation detectors [27,54]. As such, high-mass systems offer little information for either M max or R 1.4 , and thus the GW data primarily probe only the lowmass/low-density part of the EoS. Indeed, Fig. 6 shows Table I. Total correlation (I) and entropy (H) of the joint distributions over ln p1.0, ln p1.5, ln p2.0, ln p3.0, and ln p4.0 induced by several processes as well as the entropy of the marginal distribution over only ln p2.0 (H(ln p2.0)). The nonparametric processes consistently have smaller I and (much) larger H than any parametric process, implying much more model freedom. This is the case even though the entropy of the marginal distributions for ln p2.0 can be comparable.  that each mock-GW M max posterior is similar to the respective posterior of Fig. 2, indicating that additional GW observations inform M max only weakly. On the other hand, x-ray measurements have already proven capable of bounding the radius of high-mass NSs [12,13]. Additionally, x-ray detection of pulsations in a compact object proves it is a NS, and thus its mass offers information about M max . Depending on the mass distribution of observed events, x-ray probes could thus probe the EoS at both low and high densities. Our mock x-ray dataset contains one such NS with mass 2.50 M . Fig-ure 6, then, shows that when we use a parametric model to fit all the x-ray data, biases can arise as no EoS in the prior process can simultaneously reproduce the correct values for both M max and R 1.4 . The bias is smaller in the GW-based results as the data there probe a narrower density range, resulting in the appearance of mild tension between the two datasets. By extension, a newly observed GW signal for a 1.4M NS would be in tension with the x-ray-based results, despite no real astrophysical inconsistency. Recent concerns of tensions between PREX-II [55] and astrophysical predictions may be influenced by a similar mechanism, as noted in Refs. [44,45].

VII. IMPACT OF PARAMETRIC PRIOR CHOICES
These investigations show that the parametric EoS prior processes include model-dependent interdensity correlations that influence the resulting inference. Such prior processes are constructed based on two ingredients: (i) a functional form for p(ρ) (or an equivalent quantity) and (ii) a prior for the parameters of the function. The former may be carefully engineered, while the latter can be changed more easily. As in Sec. IV, it is therefore reasonable to wonder if we can change the nature of the interdensity correlations by a trivial change in the parameter prior, or whether the correlations are inherent to the functional form. Below we argue for the latter, as also demonstrated in Sec. IV.
First, adding parameters does not necessarily always increase model freedom. Adding a parameter to any model we have shown so far will require choosing a distribution for that parameter, and a reasonable range will strongly depend on the functional form. For the piecewise-polytrope model, this process is somewhat easier, as additional adiabatic indices for new segments have clear physical meaning. Therefore reasonable ranges can be chosen. For the spectral model, with the addition of a new spectral component, there is no obvious mapping of parameter values to physics, and so tuning parameter ranges is much harder.
Second, of the existing parameters in the models, we typically find that only a few are meaningfully constrained. For example, the speed-of-sound model has only a single Gaussian bump, and thus current astrophysical data tightly constrain this bump to be at densities low enough to produce pulsars consistent with, e.g., Refs. [10,11,21]. This severely limits the flexibility of the model, as the logistic term is, by itself, not strong enough to support realistic NSs. In practice a 1 and a 2 are overconstrained in this model and a 4 , and a 5 are underconstrained. As a result, the speed-of-sound model is the least flexible (and leads to the most stringent constraints) even though it has the most parameters. We find similar behavior in the spectral and piecewise-polytrope models, suggesting that the effective number of parameters in the models is fewer than what is nominally stated.
Third, due to the fine-tuning of the parametric models, attempting to redefine priors on parameters is generally not an efficient way to expand model freedom. As an example, we consider the spectral EoS where we find that EoS candidates have a priori strong correlations between parameters in order to satisfy causality and stability. These correlations were noted in Ref. [29] and are also shown in Fig. 7. We find that the γ i are alternately strongly correlated or anticorrelated with each other. It is possible that other distributions, in particular distributions that upweight EoS further from the line of strongest correlation, reduce the strength of interdensity correlations.
We test this by upweighting more extreme γ i values. The reweighting procedure does indeed change the posterior distribution of γ i parameters, although the inferred distribution in M max -R 1.4 is effectively unchanged. Moreover, the M max -R 1.4 prior does not significantly extend into the previously excluded region closer to the stability threshold. We find similar results with a different reweighting of γ i that instead favors for central values. Figure 7 also shows that extending the prior range on γ i will not extend the reach of the spectral model as the parameter posteriors are not significantly affected by the prior cutoffs. We reach similar conclusions with the piecewise polytrope.
Overall, while it was possible to remove correlations between only a subset of pressures within our toy models in Sec. IV, it is generally difficult to do even that with real parametrized EoS models.

VIII. CONCLUSIONS
All models of the dense-matter EoS should contain some correlations between density scales due to causality and thermodynamic stability requirements. However, in this study we show that phenomenological parametric models such as the spectral, piecewise-polytrope, and speed-of-sound models impose even stronger correlations a priori. As a result, NS properties are constrained more tightly in parametric models than in nonparametric ones in ways that are not supported by the data. Regardless of whether these tighter constraints end up being compatible with the true EoS, their emergence is at-tributable to what are effectively model-dependent prior assumptions dictated by the phenomenological nature of the parametrizations. Viewed in this way, they deserve the same scrutiny as other prior choices imposed by the analyst.
The concerns about implicit correlations are alleviated by GP-based nonparametric models that enjoy extensive model freedom, restricted only by causality and thermodynamic stability. They allow us to generate, with no additional modeling effort, candidate EoS with complex phenomenology that could be associated with, e.g., a transition to quark matter in the cores of NSs. For example, Refs. [56,57] study EoS with complex speed-of-sound phenomenology, while Refs. [51,[58][59][60][61][62] consider strong first-order phase transitions that result in a discontinuity in the speed of sound and multiple stable branches. The GP prior process is able to recreate such behaviors generically.
The parametric models are relatively easier to implement. However this might come at the cost of finetuning which makes it harder to sample from the prior as many draws are unphysical. Extension to more complex phenomenology, such as phase transitions, is less straightforward and might need tailored parametric models [6,63,64] unless the parametrization supports such behavior inherently. The piecewise polytrope specifically, as implemented here, can lead to priors and posteriors with "kinks" [28,31], while Refs. [65,66] discuss its behavior in cases where the observed NSs do not reach high enough densities to probe all polytropic segments. More complicated parametric models exist (see Appendix B), but as we argue in Sec. VII, improving parametric models by adding parameters or extending the priors ranges is not always straightforward or efficient. However, extreme extensions to these models (for example a O(1000) parameter extension to the piecewise polytrope) could exhibit behavior that is closer to the nonparametric results than the few-parameter models they generalize.
In conclusion, commonly used parametric models of the EoS are hampered by built-in and often opaque correlations between density scales. These correlations already affect inferences based on these models, and these effects will only become more severe with additional astrophysical data. The impact of the EoS model on inference acts as an additional systematic error that must be addressed to achieve highly informative EoS constraints [67][68][69][70][71][72]. Our work shows that the nonparametric GP-based model addresses this EoS model systematic and restores model freedom by forgoing the use of specific functional forms for the EoS itself and instead parametrizing a wide range of possible correlations directly.

ACKNOWLEDGMENTS
We thank Les Wade for useful discussions on the implementation of the spectal model in LALsuite . R.E.  function. Each GP draw is a realization of a multivariate Gaussian distribution, which is loosely conditioned on nuclear models. The GP from which EoS candidates are sampled has a covariance which is governed by a kernel function through the parameters σ and that control the strength and length of the correlations respectively [26]. The EoS prior process includes EoS drawn from multiple underlying GPs with different parameters: log σ ∈ U (1, 10) and ∈ U (0.1, 0.9). However, our GPs' kernels contain additional terms as well. See Ref. [26] for more details. In total we use ∼ 2 × 10 6 draws, of which ∼ 3 × 10 5 contribute to the prior nontrivially.
In Ref. [12], a single GP is used to generate EoS realizations using the same method. This single GP is more tightly bound to the mean realization and nuclear models than corresponding piecewise-polytrope and spectral models due to the values of σ = 1 and = 1 chosen. Such values correspond to stronger correlations and over larger length scales than any GP we employ here. This demonstrates that although the GP model is flexible, it is not necessarily agnostic. This can be useful, for example, to examine the validity of a set of related nuclear models given astrophysical data [73]. In the piecewise-polytrope approach, consistent with Refs. [22,28], the polytropic exponent is a piecewise constant function, which changes value at two predetermined densities Here ρ 1 = 10 14.7 g/cm 3 and ρ 2 = 10 15 g/cm 3 are fixed via an optimization for a set of candidate EoS following from nuclear models [22]. The parameter K 1 is chosen to give some value p 1 ≡ p(ρ 1 ), and K 2 and K 3 are then fixed by continuity. Therefore {Γ 1 , Γ 2 , Γ 3 , p 1 } are the parameters in this model. Their corresponding priors are given in Table II. Extensions to this model with more polytropic segments or allowing the transition densities to vary are proposed in Refs. [74][75][76][77].
When {Γ 1 , Γ 2 , Γ 3 , p 1 } are sampled from a uniform distribution then the resulting total EoS will be neither necessarily causal or stable. Therefore, we have to enforce these constraints after the fact; specifically, we sample a set of parameters, compute the corresponding EoS, and save it only if it obeys causality and stability. 8 Overall, we retain ∼ 1.6 × 10 5 EoS. We verified that this number is enough to efficiently characterize the posterior by confirming that we get consistent results with half as many draws. For the computation of the p(ρ), and ε(p) relations, we used LALSimulation, a subsection of LALSuite [78]. For checks of NS properties such as causality, we used LALInference [78,79]. Our priors are slightly more restrictive than those used in Ref. [31] due to computational problems that arise for candidates with the highest Γ 2 , Γ 3 , which tend to represent acausal EoS candidates anyway.

Spectral parametrization
In the spectral approach, the polytropic exponent is expanded in a series of basis functions. Following the conventions of Ref. [19] which introduced the spectral parameterization, we take x ≡ p/p 0 where p 0 is the smallest pressure where the spectral parametrization will be used; the parametrization is matched to some other EoS at this density which serves as the low-density crust [80]. Then we set In most of the literature, and for our purposes n is set to 3. Note that the overall scaling of p(ρ) is fixed by γ 0 , and again we have four total parameters {γ 0 , γ 1 , γ 2 , γ 3 }.
In practice sampling individual parameters is impractical because generic combinations of parameters produce unphysical EoS, even if the parameter ranges are chosen carefully. Instead, following Ref. [29], we sample in a different parameter space r = (r 0 , r 1 , r 2 , r 3 ) and under an affine map construct samples in γ. The prior on r is given in Table II. Our analysis uses a total of ∼ 1.9 × 10 5 draws from the spectral model. We again use the LALSuite components LALSimulation and LALInference [78,79], with particular spectral components implemented by Ref. [28]. The spectral EoS is stitched to a model of the SLy EoS just below 0.5ρ nuc [28] (see Fig. 1).

Speed of sound parametrization
In this approach, the speed of sound is parametrized as a function of energy density. Taking z ≡ ε/(ρ nuc c 2 ), we write with a 1 , a 2 , a 3 , a 4 , a 5 real parameters, and a 6 fixed by matching to a low-density crust. In Ref. [24], the matching is done to a chiral effective field theory at ∼ ρ nuc with limits based on Fermi liquid theory enforced up to a density of 1.5ρ nuc . Since we do not wish to use more nuclear theory information for this model than others, we instead stitch to SLy at a density of 0.6ρ nuc , comparable to the stitching density of the spectral model. Because of this, the parameter ranges in our implementation must be adjusted to generate realistic EoS candidates. The prior on each parameter is given in Table II. Our analysis uses a total of ∼ 1.6 × 10 5 draws from this model. A similar model based on the speed of sound is presented in Ref. [81].

Causality in parametric models
Because of the relatively large uncertainties, we follow Ref. [28] in not excluding parametric EoS until they have a large violation of the speed of sound c s > 1.1c. This is the standard criteria used in LALinference for the piecewise-polytrope and spectral models as part of determining if an EoS is physical. For consistency we extend it to the speed-of-sound parametrization as well. The primary motivation for this is to allow a possibly acausal EoS to represent another, causal EoS which is not modeled effectively by the prior on EoS [28]. In addition, the LALinference implementation of the spectral and piecewise-polytrope models enforces this criterion only up to the central density of the maximum mass NS. In the speed-of-sound model, we require the EoS to be causal (or approximately causal) everywhere. The nonparametric model obeys exact causality (c s ≤ c) at all densities. speed-of-sound (astro-causal) speed-of-sound (prior-causal) speed-of-sound (prior) speed-of-sound (astro) causality Figure 9. Comparison between using strictly causal (c 2 s < c 2 ) parametric EoS with each model (gray) and the headline results allowing some violation of the causal limit (c 2 s < 1.1c 2 ). When restricting to only causal EoS, the issues of model dependence and insufficient coverage of the physically allowed Mmax-R1.4 space are more severe, especially for the piecewise-polytrope.
These choices were made for consistency with past work [28,32], but we still find that this extra model freedom does not enable to spectral and piecewise-polytropic models to fill in the physically available M max -R 1.4 space up to the causality threshold. Figure 9 shows that excluding the acausal models minimally affects the spectral and speed-of-sound results. However, the piecewise-polytrope results are noticeably tighter, and the M max − R 1.4 allowed parameter space is covered significantly less. This also explains why the nominal piecewise-polytrope prior supports larger pressures than the nonparametric prior in Fig. 10.  Figure 10. Symmetric 90% credible region for the pressure p at each density ρ in units of the nuclear saturation density (left) and the radius R as a function of the mass M . From top to bottom we show results with the spectral, piecewise-polytrope, and speed-of-sound parametric models. At each panel we overplot the corresponding nonparametric result for comparison. The spectral p-ρ panel is identical to Fig. 1, but we show it for completeness. As in Fig. 1 we show results with all astrophysical data (labeled "astro," solid lines) and restricting to the heavy pulsars only (labeled "psr," dashed lines).

Appendix C: Further results with the parametric models
Most results presented the main body of this paper were obtained using the spectral EoS model. In this Appendix we present similar results with the piecewisepolytropic and the speed-of-sound models. Figure 10 shows the pressure-density and mass-radius posteriors for all EoS models. As expected, we find that the posteri-ors differ due to the different models. Figure 11 shows the posteriors for the speed of sound as a function of the density where again the nonparametric case results in the less constrained results as an outcome of larger model flexibility. Figure 12 shows the equivalent of Fig. 2 for the piecewise-polytrope and the speed-of-sound models. We again find that both parametric models lead to tighter constraints on M max for high values.
The two- speed-of-sound (astro) Figure 11. The speed of sound squared as a function of baryon density in the spectral (top), piecewise-polytrope (middle) models, and speed-of-sound (bottom) models, compared to our nonparametric model. dimensional plots show that model-dependent correlations between the maximum mass and the radius (equivalently low and high densities) exclude certain regions of the M max -R 1.4 space in the parametric marginal priors. As with the spectral model, there is a gap between the piecewise-polytrope prior and the approximate causality threshold. The corresponding plots for the speed-ofsound parametrization show the opposite behavior: the prior reaches the causality threshold, but it fails to produce EoS with large R 1.4 and small M max . nonparametric (astro) nonparametric (prior) speed-of-sound (prior) speed-of-sound (astro) causality Figure 12. Similar to Fig. 2 but for the piecewise-polytrope (top panels) and the speed-of-sound (bottom panels) models. As with the spectral case, we find that the parametric model leads to tighter posteriors for Mmax compared to the nonparametric model. The two-dimensional plots show that this is again due to model-dependent correlations in Mmax-R1.4.