Bayesian nonparametric inference of neutron star equation of state via neural network
Abstract
We develop a new nonparametric method to reconstruct the Equation of State (EoS) of Neutron Star with multimessenger data. As an universal function approximator, the FeedForward Neural Network (FFNN) with one hidden layer and a sigmoidal activation function can approximately fit any continuous function. Thus we are able to implement the nonparametric FFNN representation of the EoSs. This new representation is validated by its capabilities of fitting the theoretical EoSs and recovering the injected parameters. Then we adopt this nonparametric method to analyze the real data, including masstidal deformability measurement from the Binary Neutron Star (BNS) merger Gravitational Wave (GW) event GW170817 and massradius measurement of PSR J0030+0451 by NICER. We take the publicly available samples to construct the likelihood and use the nested sampling to obtain the posteriors of the parameters of FFNN according to the Bayesian theorem, which in turn can be translated to the posteriors of EoS parameters. Combining all these data, for a canonical 1.4 neutron star, we get the radius km and the tidal deformability (90% confidence interval).Furthermore, we find that in the high density region (), the 90% lower limits of the ( is the sound speed and is the velocity of light in the vacuum) are above , which means that the socalled conformal limit (i.e., ) is not always valid in the neutron stars.
0000000190340866]MingZhe Han 0000000290787825]JinLiang Jiang 0000000191207733]ShaoPeng Tang 0000000289666911]YiZhong Fan
1 Introduction
Neutron stars (NSs) are natural laboratories to study the properties of the dense matter (see Refs.(Lattimer, 2012; Lattimer & Prakash, 2016; Özel & Freire, 2016; Oertel et al., 2017) for reviews), in which the extreme conditions (i.e., the conditions in the NS core) are hard to achieve in terrestrial experiments. Usually, the matter in NS (in most of its lifetime except for the birth and death/merger) can be well described by the zerotemperature pressuredensity relation, i.e., the equation of state (EoS). Up to now, the properties of the low density matter have been well constrained by the nuclear theories/experiments (Lattimer & Lim, 2013; Krastev & Li, 2019; Lim & Holt, 2018), and the behavior of super high density matter is robustly predicted by the perturbative quantum chromodynamics (pQCD) theory (Baym et al., 2018). However, the matter in the transition area (between the low and the super high density) is not clear yet. For instance, the presence of a hadronquark phase transition (Tang et al., 2021) or not is quite uncertain. The observations of NSs are essential to reveal the properties of the matter in the medium density region.
On August 17, 2017, the gravitational wave (GW) signal from the coalescence of binary neutron star (BNS) was detected by the LIGO/Virgo detectors (Abbott et al., 2017, 2019), which provides new avenues to constrain the EoS via the NS matter effects. Benefiting from this remarkable GW event (GW170817), various studies have been implemented to better understand the EoS of NSs (Sieniawska et al., 2018; Most et al., 2018; De et al., 2018; Abbott et al., 2018; Raaijmakers et al., 2019; Miller et al., 2019; Landry & Essick, 2019; Montaña et al., 2019; Miller et al., 2020; Tang et al., 2020; Raaijmakers et al., 2020; KanakisPegios et al., 2020; Essick et al., 2020; Landry et al., 2020; Zimmerman et al., 2020; Biswas et al., 2020; Huth et al., 2021; AlMamun et al., 2021). With the improvement on the sensitivities of GW detectors, the detectable volume will be significantly enhanced, and an increasing sample of BNS mergers will be available to constrain the EoS. Nevertheless, in this work we do not take into account GW190425 (Abbott et al., 2020) because of its weak tidal effect (due to the low signal to noise ratio and massive components). Moreover, its possible black holeneutron star merger nature (Han et al., 2020) may influence the inference. In addition to the GW signals, the massradius () measurements of the isolated NS PSR J0030+0451 were recently reported by the NICER collaboration (Gendreau et al., 2016). These measurements (Riley et al., 2019; Miller et al., 2019) were based on the pulse profile modeling of the Xray emission from the hot spots on the NS surface, which are more reliable and accurate than traditional/indirect spectroscopic measurements/estimates. Such data are valuable to constrain the NS EoS (Miller et al., 2019; Raaijmakers et al., 2019; Sieniawska et al., 2018). Tighter constraints on the EoS have been reported in the joint analysis of these multimessenger observations of NSs (Sieniawska et al., 2018; Raaijmakers et al., 2019; Miller et al., 2019; Jiang et al., 2020; Raaijmakers et al., 2020; Landry et al., 2020; Zimmerman et al., 2020; Biswas et al., 2020; Tang et al., 2021; Huth et al., 2021; AlMamun et al., 2021) under the assumption that all NSs share the same EoS.
The phenomenological parameterization methods, including for instance the spectral representation (Lindblom, 2010) and the piecewise polytropic expansion (Read et al., 2009; Özel et al., 2016; Raithel et al., 2017), have been widely adopted to study the EoS of NSs. These parametric models are described as , where is the probability density function (PDF) of the underlying probability distribution as a function of . The inference aims to get the credible intervals of . However, a specific parametric form may be modeldependent and the range and type of inferences are likely limited, too. In order to avoid model misspecification, the extension beyond the parametric assumptions may be needed to attain more robust models. In this situation, one can consider a nonparametric model where the functional space is so large that the approach will not depend on a specific parametric form. Such nonparametric approaches are considered as the models with infinite dimensional parameter. However, in practice we usually use finite dimensional parameter instead, since it can provide an approximation to the infinite sum. The nonparametric model via Gaussian process (GP) has been implemented to the studies of the EoS (Landry & Essick, 2019; Essick et al., 2020; Landry et al., 2020). These works use GP and Bayesian methods to infer the NS EoS from multimessager data. The posteriors are calculated by using the Monte Carlo integration, which is relatively inefficient in comparison to the Markovchain Monte Carlo (MCMC) method. The MCMC approach has not been adopted in these previous works because the associated jump proposals are thought to be nontrivial. In practice, one can use a finite dimensional EoS table to represent the real EoS, then the EoS model has a finite number of random variables (RVs) at the discrete set of densities drawn from the GP, whose covariance matrix is determined by a set of hyper parameters, and the hyper parameters is chosen based on the training set (i.e., the theoretical EoSs). However, to get a precise approximation one often take 100 dimensions, thus the model will have 100 parameters, and it is extremely expensive in computation to perform the Bayesian inference, while as one will see below the number of RVs of the model in this work is much fewer. Besides, the GP EoS prior relies on the choice of the training set, therefore, it might be limited by the current knowledge of EoS (see Landry & Essick, 2019; Essick et al., 2020, for more details). Fujimoto et al. (2018, 2020, 2021) have studied the methodology of the machine learning method for EoS. They take the use of the FeedForward Neural Network (FFNN) to map the neutron star data to the EoS parameters, while in their approach the training data are still generated by a specific parametric form. Different from their work, we do not consider the FFNN as an interpolation tool between the observations and the EoS, but as a representation of the EoS itself, then we can combine it with the Bayesian statistical framework, which is more explainable.
In this work we use the Bayesian nonparametric method via the FFNN to reconstruct the NS EoS with multimessenger data. We describe in details the nonparametric representation of NS EoS along with fitting the theoretical EoSs in Sec.2.1, and introduce the observed data used for the Bayesian inference in Sec.2.2. The main results of our work is presented in Sec.3, including the simulated/real data. Finally, Sec.4 is our summary and conclusion.
2 Method
2.1 Nonparametric representation of EoS
Consider a general regression problem of fitting an output as a function of input :
(1) 
where the is called the mean function and the is called the residual, usually they are independent. Based on which part of the model is nonparametric, the nonparametric models can be divided into three types: nonparametric mean function, nonparametric residual distribution, and fully nonparametric regression. Furthermore, the nonparametric mean function methods consist of basis expansion, bspline, Gaussian process, and the regression tree (for more details see Chapter 4 of Müller et al., 2015). In this work we use the nonparametric mean function method via the FFNN. An FFNN with one hidden layer can be written as
(2) 
where is the dependent variable, is the independent variable, is the number of neural nodes, and are the weights, is the bias, is the residual, and is the activation function (a nonlinear function). Cybenko (1989) and others pointed out that if we use a sigmoidal function as the activation function, the finite sums in the form of Eq.(2) are dense in the functional space of continuous functions defined in the ndimensional unit cube (for more details see Cybenko (1989)), which means that we can use Eq.(2) with a sigmoidal activation function and finite to fit any real continuous functions. This property exactly matches the feature of nonparametric models (i.e., the modeling should be flexible and robust) we have mentioned above. It thus gives us the possibility to take it as a nonparametric representation of the EoS. In this work we use the sigmoid function as the activation function
(3) 
which yields when . Both the input layer and the output layer have one node without activation function, as for the hidden layer, we have tried several settings (5, 10, 20, and more). We find that if the number of nodes is less than 10, the model could not represent all the EoSs with relatively small errors, and as the number of nodes increase, the efficiency of the sampling algorithm decreases. Therefore, we choose 10 nodes for this work. Besides, it should be noted that we call our model as nonparametric model for its ability of representing almost all theoretical EoSs (as shown below), even though it only has 31 parameters which is far away from infinity.
In the zero temperature case which we consider here, the NS EoS can be represented by the relation between the total energy density and the pressure . Thus in order to constrain the EoS, we need to reconstruct the relation from the data. However, not all the function can be a realistic EoS, unless it satisfies two conditions:

The microscopical stability condition, i.e.,

The causality condition, i.e., the sound speed , where is the speed of light in vacuum.
Fortunately, it has been shown that the auxiliary variable
(4) 
can automatically satisfy these two conditions (Lindblom, 2010), thus it is widely used in parameterizing the EoS (Landry & Essick, 2019; Essick et al., 2020; Landry et al., 2020). We therefore take it as the output variable, and the as the input variable of the FFNN. We standardize the inputs of the FFNN before the inference, i.e., , where and are the mean and standard deviation of , respectively. In the low density () range we match the constructed EoS with the known EoS SLy (Douchin & Haensel, 2001)^{1}^{1}1In practice, when integrating the macroscopic properties of neutron star, we divide the EoS constructed by the FFNN into a finite dimensional (100) EoS table (The pressures of these nodes are logarithmically uniform in , and by using Eq. 4 we can obtain the via fixed stepsize integration.), and between adjacent discrete density points we apply the method described in Appendix B of Lindblom & Indik (2014) to interpolate the EoS tables. The Eq. (4) gives an one to one map between the auxiliary variable and the first derivative of energy density with respect to the pressure. If the initial condition of pressure and energy density is further given, then this would yield a unique EoS. So at the matching point, we fix the first point of the constructed EoS by the corresponding pressure and energy density of the EoS SLy.. Additionally, according to the theoretical conjecture that the energy of the unitary gas is less than the energy of pure neutron matter and the various experimental constraints on the symmetry energy parameters(Tews et al., 2017; Lattimer & Steiner, 2014), we constrain the with (Jiang et al., 2019), where is the saturation density. And the pressure is limited to according to the theoretical conjecture that the two nucleon interaction pressure is an absolute lower bound because the threebody interactions in pure neutron matter are always repulsive (Özel et al., 2016). EoSs violating these conditions will be directly rejected in the sampling procedure.
To illustrate the accuracy and efficiency of this new nonparametric representation, we first implement the FFNN to fit the theoretical EoSs. We use the Python library Keras (Chollet & others, 2018) with TensorFlow (Abadi et al., 2016) as a backend. We take the mean square errors as the loss function, i.e.,
(5) 
where the is the true value and the is the predicted value. As for the optimization method we use Adam (Kingma & Ba, 2014). Three groups of EoSs that we fitted are characterized by the different compositions (see Fig.1):

Hadronic EoSs, i.e., BSK20, BSK21, BSK22, BSK23, BSK24, BSK25, BSK26 (Goriely et al., 2010, 2013); BSR2, BSR6 (Agrawal, 2010); DD2, DDHD, DDME2 (Banik et al., 2014; Gaitanos et al., 2004; Lalazissis et al., 2005); ENG (Engvik et al., 1996); GM1 (Glendenning & Moszkowski, 1991); KDE0V, KDE0V1 (Gulminelli & Raduta, 2015; Agrawal et al., 2005); MPA1 (Müther et al., 1987); NL3 (Lalazissis et al., 1997); RS (Friedrich & Reinhard, 1986); SK255, SK272, SKI2, SKI3, SKI4, SKI5, SKI6, SKMP, SKOP (Agrawal et al., 2003; Reinhard & Flocard, 1995; Nazarewicz et al., 1996; Bennour et al., 1989; Reinhard et al., 1999); SLY230A, SLY2, SLY9, SLY (Chabanat et al., 1997; Chabanat, 1995; Douchin & Haensel, 2001); TM1 (Sugahara & Toki, 1994).
The direct variable we fitted is the , we calculate the absolute errors in fitting , i.e., , and the values are as low as , while the absolute value of is . Furthermore, we use the fitted to calculate the corresponding , which is defined as , and the relative errors, i.e., , are as low as 0.1% 10%. Thus it is reasonable to conclude that the FFNN can well represent the various realistic EoSs.
2.2 Bayesian Inference
After constructing the nonparametric representation of EoS with FFNN, we can optimize the parameters (i.e., the weights, bias, and residual of the FFNN) with the observation data. The standard approach for neural networks is to approximate a minimal loss. From a statistician point of view, we can consider it as MAP (Maximum A Posteriori) or MLE (Maximum Likelihood Estimation) by constructing a proper likelihood function. These two point estimate methods are relatively easy to realize, but can not give us the credible intervals of the parameters that we are interested in. Then we use the Bayesian statistical framework to obtain the posterior distributions of the parameters.
Assuming that all neutron stars share the same EoS, the likelihood for Bayesian inference can be written as
(6) 
where , , and are mass, radius, and pseudo enthalpy of the star’s core, respectively. The and are respectively the parameters for GW and FFNN. is the interpolated GW likelihood that marginalized over the other nuisance parameters (Hernandez Vivanco et al., 2020), is the strain data of GW170817, and are primary/secondary mass and its tidal deformability. takes only if all the nuclear constraints are satisfied, otherwise it takes (In practice, we set the log(0) to be ). For the we follow the method in Sec. 3.1.2 of Miller et al. (2020), that is
(7) 
where the is the probability distribution for the mass of PSR J0740+6620 (Fonseca et al., 2021). While for constructing the , we use the Gaussian kernel density estimation (KDE) with the samples of PSR J0030+0451 from NICER (Riley et al., 2019; Miller et al., 2019).
We implement the enthalpybased formulas of Lindblom & Indik (2014) to solve the TolmanOppenhimerVolkoff and ReggeWheeler equations to get the mass, radius and tidal deformability of neutron star. The prior of is based on the method described in the Sec. 5.1 in Miller et al. (2021). We uniformly sample in (0, 1), the central enthalpy is , where the and are respectively the central enthalpy of a 1 and a maximummass nonrotating neutron star for the EOS under consideration. And the quadratic prior is because the central enthalpy changes rapidly near the maximum mass. Then the can be combined with to calculate a pair of macroscopic observables ( and ), and hence the . For the , we sample the chirp mass and the mass ratio () instead of the component masses with proper uniform priors. We can then map the transformed source frame masses to with the EoS parameters , and calculate the . As for the parameters of FFNN, we set all the 31 parameters to be uniformly distributed in , the bounds are chosen as the credible interval of the standard normal distribution, since the inputs have been standardized. We use the nested sampling algorithm PyMultiNest (Buchner et al., 2014) in Bilby (Ashton et al., 2019) to sample the posterior distributions of all those parameters, typically with 1000 live points, a sampling efficiency of 0.01 (without constant efficiency), and a tolerance of 0.5.
3 Results
3.1 Simulated data analysis
To illustrate the performance of our nonparametric method on recovering the EoS from observations, we simulate mock data calculated by theoretical EoSs and perform the same procedure applied to the analysis of the real data except for the nuclear constraints. We select three NSs’ masses based on the most probable masses in the posterior samples of GW170817 ^{2}^{2}2https://dcc.ligo.org/LIGOP1800061/public and PSR J0030+0451 (Riley et al., 2019; Miller et al., 2019), then use the theoretical EoS to calculate the corresponding tidal deformabilities and radius. The mock likelihood function is multidimensional Gaussian, with the means given by simulated values and the comparable variances based on current quality of the real data samples. We choose four tabulated EoSs as the theoretical model, including three hadronic EoSs (WFF1 (Wiringa et al., 1988), H4, MS1 (Müller & Serot, 1996)), and a hybrid EoS (ALF2), where the hadronic EoSs range from relatively soft to stiff.
The credible intervals of the recovered parameters are summarized in Tab.1. We can see that there are no egregious errors in the recovering process and all the injected values lie in the 90% credible intervals of the recovered values, which is in keeping with our expectations. Fig.2 shows the reconstructed EoSs () of all the injected theoretical EoSs, and the posterior regions comfortably include the curves of the corresponding injected EoSs. Therefore, it is rather appropriate to use this nonparametric model to reconstruct the EoS with the observational data.
EoS  

WFF1  
H4  
MS1  
ALF2 
3.2 Real data analysis
Our direct results are the posterior distributions of the weights and bias parameters of the hidden layer, but these parameters don’t have a strong correlation with the output variable (this is why we call it hidden layer). We thus transform them to the distribution of the output variable , and further calculate the parameters of the EoS that we are interested in. The results can be divided into five groups:

The microscopical stability condition, the causality condition, and the request of (i.e., Prior).

The constraints of 1) and the mass information of PSR J0740+6620 (i.e., PSR).

The constraints of 3) and the data of GW170817 (i.e., PSR+Nuclear+GW).

The constraints of 4) and the data of PSR J0030+0451 (i.e., PSR+Nuclear+GW+NICER).
The mass information of the PSR J0740+6620 raises the pressure at , while the possible region of the relation has been effectively narrowed down by the nuclear constraints (see upper left panel of Fig.3). The results yielded by adding the GW data indicate that stiffening at and softening at are not allowed. These characteristics can also be seen from the plot in Fig.3 where the relatively larger at and the relatively smaller at are not favored. Meanwhile, as also found in other researches (Raaijmakers et al., 2020; Jiang et al., 2020; Landry et al., 2020), the inclusion of the data of PSR J0030+0451 further improves the constraints.
For the bulk properties and , the PSR mass information significantly increase the radius at low mass region, while the nuclear constraints have ruled out a good fraction of high value regions (see lower panels in Fig.3). These larger and regions are further excluded by the GW data (see also the Tab.2). While the data of PSR J0030+0451 boost the lower bound of the and for and , respectively. This can be explained by the relatively larger radius of PSR J0030+0451 compared with the radii of NSs in GW170817 (Abbott et al., 2018), albeit they have similar mass.
It has been conjectured that, in NS the speed of sound is less than the speed of light in vacuum divided by , i.e., , which is also called the conformal limit. While the existence of the heavy NS () may call for the violation of conformal limit (Bedaque & Steiner, 2015). Our results show that the 90% lower limits of the are larger than 1/3 in the high density region, which means that the conformal limit is indeed violated at some densities.
Observable  Prior  PSR  PSR+Nuclear  PSR+Nuclear+GW  PSR+Nuclear+GW+NICER 

4 Summary
In this work, we introduce a new nonparametric representation of the NS EoS via the FFNN. Then we implement this nonparametric model to fit the theoretical EoSs including the hadronic EoSs, the hyperonic EoSs, and the quark matter (hybrid) EoSs. We find that our method performs well in reconstructing these theoretical EoSs. For the absolute value of of the order , the absolute fitting errors are on the order of . And the relative errors in calculating the corresponding are as low as 0.1% 10%. Then we make injections and recover them using the same procedure as what we applied on the real data analysis. We inject mock observations of three NSs considering the masses of GW170817 and PSR J0030+0451, and the tidal deformabilities and radius of these NSs are calculated by four theoretical EoSs, i.e., WFF1, H4, MS1, and ALF2. We find that there are no egregious errors on recovering both the macroscopic properties and microscopic relations with simulated data. After validating our model, we further adopt it to analyze the real observation data, i.e., massdeformability measurement of GW170817 and massradius measurement of PSR J0030+0451, and use Bayesian method to constrain the EoS parameters. We find that the radius and deformability of a canonical neutron star are correspondingly and (see the bottom panels of Fig.3), where the credible intervals are all 90%. These results are in agreement with those found in the previous investigations (Sieniawska et al., 2018; Raaijmakers et al., 2019; Miller et al., 2019; Jiang et al., 2020; Raaijmakers et al., 2020; Landry et al., 2020; Zimmerman et al., 2020; Biswas et al., 2020; Huth et al., 2021; AlMamun et al., 2021).
The speed of sound in NS is another significant property of the EoS. There is a conjecture that the sound speed in the NS satisfies the so called conformal limit, i.e., . Using our nonparametric method, we infer that the square of the sound speed in the core of the maximum mass configuration NS is larger than at 90% credible level (see the top right panel of Fig.3), that means the conformal limit is violated in the center of very massive NSs. Raaijmakers et al. (2019) have imposed a constraint that the speed of sound of each EoS converges to conformal limit at asymptotically high densities according to pQCD calculations (Fraga et al., 2014), in their speed of sound (CS) model. While our model does not guarantee this constraint, but notice that the central densities of the observed NSs are far away from pQCD density, this constraint will not significantly affect our results.
As mentioned in Sec.2.1, the finite sum of Eq.(2) can fit the continuous function at high resolution, and the larger the number of the nodes is the better the approximation gets. Therefore, we need to extend the width of the hidden layer(increase the number of the neural nodes) to get better approximation in the future. However, larger neural networks have more parameters, which may be extremely hard for the sampling algorithm to converge. Thus a ‘smarter’ sampling algorithm like PolyChord (Handley et al., 2015) that can support higher dimensional parameter space would be considered, and the balance between the performance of fitting and the computational cost would also be carefully considered. Besides, more complex model would lead to the overfitting problem (i.e., the model can better fit the current data, but worse predict the unknown data point.), and the posterior inference of more complex models may dramatically change as more data is accrued. Therefore, in future works we should pay more attention to the tradeoff between the variance and bias to get a better model.
5 Acknowledgments
We thank Mr. Lei Zu for discussing the neural networks. This work was supported in part by NSFC under grants of No. 11921003 and No. 11525313. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gwopenscience.org/), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the MaxPlanckSociety (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain.
References
 Abadi et al. (2016) Abadi, M., Barham, P., Chen, J., et al. 2016, arXiv eprints, arXiv:1605.08695. https://arxiv.org/abs/1605.08695
 Abbott et al. (2017) Abbott, B. P., et al. 2017, Phys. Rev. Lett. , 119, 161101, doi: 10.1103/PhysRevLett.119.161101
 Abbott et al. (2018) —. 2018, Phys. Rev. Lett. , 121, 161101, doi: 10.1103/PhysRevLett.121.161101
 Abbott et al. (2019) —. 2019, Physical Review X, 9, 011001, doi: 10.1103/PhysRevX.9.011001
 Abbott et al. (2020) —. 2020, Astrophys. J. Lett. , 892, L3, doi: 10.3847/20418213/ab75f5
 Agrawal (2010) Agrawal, B. K. 2010, Phys. Rev. C, 81, 034323, doi: 10.1103/PhysRevC.81.034323
 Agrawal et al. (2005) Agrawal, B. K., Shlomo, S., & Au, V. K. 2005, Phys. Rev. C, 72, 014310, doi: 10.1103/PhysRevC.72.014310
 Agrawal et al. (2003) Agrawal, B. K., Shlomo, S., & Kim Au, V. 2003, Phys. Rev. C, 68, 031304, doi: 10.1103/PhysRevC.68.031304
 AlMamun et al. (2021) AlMamun, M., Steiner, A. W., Nättilä, J., et al. 2021, Phys. Rev. Lett. , 126, 061101, doi: 10.1103/PhysRevLett.126.061101
 Alford et al. (2005) Alford, M., Braby, M., Paris, M., & Reddy, S. 2005, Astrophys. J. , 629, 969, doi: 10.1086/430902
 Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, Astrophys. J. Supp. , 241, 27, doi: 10.3847/15384365/ab06fc
 Banik et al. (2014) Banik, S., Hempel, M., & Bandyopadhyay, D. 2014, Astrophys. J. Supp. , 214, 22, doi: 10.1088/00670049/214/2/22
 Baym et al. (2019) Baym, G., Furusawa, S., Hatsuda, T., Kojo, T., & Togashi, H. 2019, Astrophys. J. , 885, 42, doi: 10.3847/15384357/ab441e
 Baym et al. (2018) Baym, G., Hatsuda, T., Kojo, T., et al. 2018, Reports on Progress in Physics, 81, 056902, doi: 10.1088/13616633/aaae14
 Bedaque & Steiner (2015) Bedaque, P., & Steiner, A. W. 2015, Phys. Rev. Lett. , 114, 031103, doi: 10.1103/PhysRevLett.114.031103
 Bennour et al. (1989) Bennour, L., Heenen, P. H., Bonche, P., Dobaczewski, J., & Flocard, H. 1989, Phys. Rev. C, 40, 2834, doi: 10.1103/PhysRevC.40.2834
 Biswas et al. (2020) Biswas, B., Char, P., Nandi, R., & Bose, S. 2020, arXiv eprints, arXiv:2008.01582. https://arxiv.org/abs/2008.01582
 Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, Astron. Astrophys. , 564, A125, doi: 10.1051/00046361/201322971
 Chabanat (1995) Chabanat, E. 1995, Effective interactions for extreme isospin conditions.
 Chabanat et al. (1997) Chabanat, E., Bonche, P., Haensel, P., Meyer, J., & Schaeffer, R. 1997, Nucl. Phys. A, 627, 710, doi: 10.1016/S03759474(97)005964
 Chollet & others (2018) Chollet, F., & others. 2018, Keras: The Python Deep Learning library. http://ascl.net/1806.022
 Cybenko (1989) Cybenko, G. 1989, Mathematics of Control, Signals and Systems, 2, 303, doi: 10.1007/BF02551274
 De et al. (2018) De, S., Finstad, D., Lattimer, J. M., et al. 2018, Phys. Rev. Lett. , 121, 091102, doi: 10.1103/PhysRevLett.121.091102
 Douchin & Haensel (2001) Douchin, F., & Haensel, P. 2001, Astron. Astrophys. , 380, 151, doi: 10.1051/00046361:20011402
 Engvik et al. (1996) Engvik, L., Osnes, E., HjorthJensen, M., Bao, G., & Ostgaard, E. 1996, Astrophys. J. , 469, 794, doi: 10.1086/177827
 Essick et al. (2020) Essick, R., Landry, P., & Holz, D. E. 2020, Phys. Rev. D , 101, 063007, doi: 10.1103/PhysRevD.101.063007
 Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, arXiv eprints, arXiv:2104.00880. https://arxiv.org/abs/2104.00880
 Fortin et al. (2016) Fortin, M., Providência, C., Raduta, A. R., et al. 2016, Phys. Rev. C, 94, 035804, doi: 10.1103/PhysRevC.94.035804
 Fraga et al. (2014) Fraga, E. S., Kurkela, A., & Vuorinen, A. 2014, Astrophys. J. Lett. , 781, L25, doi: 10.1088/20418205/781/2/L25
 Friedrich & Reinhard (1986) Friedrich, J., & Reinhard, P. G. 1986, Phys. Rev. C, 33, 335, doi: 10.1103/PhysRevC.33.335
 Fujimoto et al. (2018) Fujimoto, Y., Fukushima, K., & Murase, K. 2018, Phys. Rev. D , 98, 023019, doi: 10.1103/PhysRevD.98.023019
 Fujimoto et al. (2020) —. 2020, Phys. Rev. D , 101, 054016, doi: 10.1103/PhysRevD.101.054016
 Fujimoto et al. (2021) —. 2021, arXiv eprints, arXiv:2101.08156. https://arxiv.org/abs/2101.08156
 Gaitanos et al. (2004) Gaitanos, T., Di Toro, M., Typel, S., et al. 2004, Nucl. Phys. A, 732, 24, doi: 10.1016/j.nuclphysa.2003.12.001
 Gendreau et al. (2016) Gendreau, K. C., et al. 2016, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.W. A. den Herder, T. Takahashi, & M. Bautz, 99051H, doi: 10.1117/12.2231304
 Glendenning & Moszkowski (1991) Glendenning, N. K., & Moszkowski, S. A. 1991, Phys. Rev. Lett. , 67, 2414, doi: 10.1103/PhysRevLett.67.2414
 Goriely et al. (2010) Goriely, S., Chamel, N., & Pearson, J. M. 2010, Phys. Rev. C, 82, 035804, doi: 10.1103/PhysRevC.82.035804
 Goriely et al. (2013) —. 2013, Phys. Rev. C, 88, 024308, doi: 10.1103/PhysRevC.88.024308
 Gulminelli & Raduta (2015) Gulminelli, F., & Raduta, A. R. 2015, Phys. Rev. C, 92, 055803, doi: 10.1103/PhysRevC.92.055803
 Gusakov et al. (2014) Gusakov, M. E., Haensel, P., & Kantor, E. M. 2014, Mon. Not. Roy. Astron. Soc. , 439, 318, doi: 10.1093/mnras/stt2438
 Han et al. (2020) Han, M.Z., Tang, S.P., Hu, Y.M., et al. 2020, Astrophys. J. Lett. , 891, L5, doi: 10.3847/20418213/ab745a
 Handley et al. (2015) Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015, Mon. Not. Roy. Astron. Soc. , 453, 4384, doi: 10.1093/mnras/stv1911
 Hernandez Vivanco et al. (2020) Hernandez Vivanco, F., Smith, R., Thrane, E., & Lasky, P. D. 2020, Mon. Not. Roy. Astron. Soc. , 499, 5972, doi: 10.1093/mnras/staa3243
 Huth et al. (2021) Huth, S., Wellenhofer, C., & Schwenk, A. 2021, Phys. Rev. C, 103, 025803, doi: 10.1103/PhysRevC.103.025803
 Jiang et al. (2020) Jiang, J.L., Tang, S.P., Wang, Y.Z., Fan, Y.Z., & Wei, D.M. 2020, Astrophys. J. , 892, 55, doi: 10.3847/15384357/ab77cf
 Jiang et al. (2019) Jiang, J.L., Tang, S.P., Shao, D.S., et al. 2019, Astrophys. J. , 885, 39, doi: 10.3847/15384357/ab44b2
 KanakisPegios et al. (2020) KanakisPegios, A., Koliogiannis, P. S., & Moustakidis, C. C. 2020, Phys. Rev. C, 102, 055801, doi: 10.1103/PhysRevC.102.055801
 Kingma & Ba (2014) Kingma, D. P., & Ba, J. 2014, arXiv eprints, arXiv:1412.6980. https://arxiv.org/abs/1412.6980
 Krastev & Li (2019) Krastev, P. G., & Li, B.A. 2019, Journal of Physics G Nuclear Physics, 46, 074001, doi: 10.1088/13616471/ab1a7a
 Lackey et al. (2006) Lackey, B. D., Nayyar, M., & Owen, B. J. 2006, Phys. Rev. D , 73, 024021, doi: 10.1103/PhysRevD.73.024021
 Lalazissis et al. (1997) Lalazissis, G. A., König, J., & Ring, P. 1997, Phys. Rev. C, 55, 540, doi: 10.1103/PhysRevC.55.540
 Lalazissis et al. (2005) Lalazissis, G. A., Nikšić, T., Vretenar, D., & Ring, P. 2005, Phys. Rev. C, 71, 024312, doi: 10.1103/PhysRevC.71.024312
 Landry & Essick (2019) Landry, P., & Essick, R. 2019, Phys. Rev. D , 99, 084049, doi: 10.1103/PhysRevD.99.084049
 Landry et al. (2020) Landry, P., Essick, R., & Chatziioannou, K. 2020, Phys. Rev. D , 101, 123007, doi: 10.1103/PhysRevD.101.123007
 Lattimer (2012) Lattimer, J. M. 2012, Annual Review of Nuclear and Particle Science, 62, 485, doi: 10.1146/annurevnucl102711095018
 Lattimer & Lim (2013) Lattimer, J. M., & Lim, Y. 2013, Astrophys. J. , 771, 51, doi: 10.1088/0004637X/771/1/51
 Lattimer & Prakash (2016) Lattimer, J. M., & Prakash, M. 2016, Phys. Rept. , 621, 127, doi: 10.1016/j.physrep.2015.12.005
 Lattimer & Steiner (2014) Lattimer, J. M., & Steiner, A. W. 2014, European Physical Journal A, 50, 40, doi: 10.1140/epja/i201414040y
 Lim & Holt (2018) Lim, Y., & Holt, J. W. 2018, Phys. Rev. Lett. , 121, 062701, doi: 10.1103/PhysRevLett.121.062701
 Lindblom (2010) Lindblom, L. 2010, Phys. Rev. D , 82, 103011, doi: 10.1103/PhysRevD.82.103011
 Lindblom & Indik (2014) Lindblom, L., & Indik, N. M. 2014, Phys. Rev. D , 89, 064003, doi: 10.1103/PhysRevD.89.064003
 Miller et al. (2020) Miller, M. C., Chirenti, C., & Lamb, F. K. 2020, Astrophys. J. , 888, 12, doi: 10.3847/15384357/ab4ef9
 Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, Astrophys. J. Lett. , 887, L24, doi: 10.3847/20418213/ab50c5
 Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, NICER PSR J0030+0451 IllinoisMaryland MCMC Samples, 1.0.0, Zenodo, doi: 10.5281/zenodo.3473466
 Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, arXiv eprints, arXiv:2105.06979. https://arxiv.org/abs/2105.06979
 Montaña et al. (2019) Montaña, G., Tolós, L., Hanauske, M., & Rezzolla, L. 2019, Phys. Rev. D , 99, 103009, doi: 10.1103/PhysRevD.99.103009
 Most et al. (2018) Most, E. R., Weih, L. R., Rezzolla, L., & SchaffnerBielich, J. 2018, Phys. Rev. Lett. , 120, 261103, doi: 10.1103/PhysRevLett.120.261103
 Müller & Serot (1996) Müller, H., & Serot, B. D. 1996, Nucl. Phys. A, 606, 508, doi: 10.1016/03759474(96)00187X
 Müther et al. (1987) Müther, H., Prakash, M., & Ainsworth, T. L. 1987, Physics Letters B, 199, 469, doi: 10.1016/03702693(87)91611X
 Müller et al. (2015) Müller, P., Quintana, F. A., Jara, A., & Hanson, T. 2015, Bayesian Nonparametric Data Analysis, Springer Series in Statistics (Cham: Springer International Publishing), doi: 10.1007/9783319189680
 Nazarewicz et al. (1996) Nazarewicz, W., Dobaczewski, J., Werner, T. R., et al. 1996, Phys. Rev. C, 53, 740, doi: 10.1103/PhysRevC.53.740
 Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Reviews of Modern Physics, 89, 015007, doi: 10.1103/RevModPhys.89.015007
 Özel & Freire (2016) Özel, F., & Freire, P. 2016, Annu. Rev. Astron. Astrophys. , 54, 401, doi: 10.1146/annurevastro081915023322
 Özel et al. (2016) Özel, F., Psaltis, D., Güver, T., et al. 2016, Astrophys. J. , 820, 28, doi: 10.3847/0004637X/820/1/28
 Prakash et al. (1995) Prakash, M., Cooke, J. R., & Lattimer, J. M. 1995, Phys. Rev. D , 52, 661, doi: 10.1103/PhysRevD.52.661
 Raaijmakers et al. (2019) Raaijmakers, G., Riley, T. E., Watts, A. L., et al. 2019, Astrophys. J. Lett. , 887, L22, doi: 10.3847/20418213/ab451a
 Raaijmakers et al. (2020) Raaijmakers, G., Greif, S. K., Riley, T. E., et al. 2020, Astrophys. J. Lett. , 893, L21, doi: 10.3847/20418213/ab822f
 Raithel et al. (2017) Raithel, C. A., Özel, F., & Psaltis, D. 2017, Astrophys. J. , 844, 156, doi: 10.3847/15384357/aa7a5a
 Read et al. (2009) Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, Phys. Rev. D , 79, 124032, doi: 10.1103/PhysRevD.79.124032
 Reinhard et al. (1999) Reinhard, P. G., Dean, D. J., Nazarewicz, W., et al. 1999, Phys. Rev. C, 60, 014316, doi: 10.1103/PhysRevC.60.014316
 Reinhard & Flocard (1995) Reinhard, P. G., & Flocard, H. 1995, Nucl. Phys. A, 584, 467, doi: 10.1016/03759474(94)00770N
 Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, Astrophys. J. Lett. , 887, L21, doi: 10.3847/20418213/ab481c
 Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, A NICER View of PSR J0030+0451: Nested Samples for Millisecond Pulsar Parameter Estimation, v1.0.0, Zenodo, doi: 10.5281/zenodo.3386449
 Sieniawska et al. (2018) Sieniawska, M., Bejger, M., & Haskell, B. 2018, Astron. Astrophys. , 616, A105, doi: 10.1051/00046361/201833071
 Sugahara & Toki (1994) Sugahara, Y., & Toki, H. 1994, Nucl. Phys. A, 579, 557, doi: 10.1016/03759474(94)909237
 Tang et al. (2020) Tang, S.P., Jiang, J.L., Gao, W.H., Fan, Y.Z., & Wei, D.M. 2020, Astrophys. J. , 888, 45, doi: 10.3847/15384357/ab5959
 Tang et al. (2021) —. 2021, Phys. Rev. D , 103, 063026, doi: 10.1103/PhysRevD.103.063026
 Tews et al. (2017) Tews, I., Lattimer, J. M., Ohnishi, A., & Kolomeitsev, E. E. 2017, Astrophys. J. , 848, 105, doi: 10.3847/15384357/aa8db9
 Wiringa et al. (1988) Wiringa, R. B., Fiks, V., & Fabrocini, A. 1988, Phys. Rev. C, 38, 1010, doi: 10.1103/PhysRevC.38.1010
 Zimmerman et al. (2020) Zimmerman, J., Carson, Z., Schumacher, K., Steiner, A. W., & Yagi, K. 2020, arXiv eprints, arXiv:2002.03210. https://arxiv.org/abs/2002.03210