Mapping neutron star data to the equation of state using the deep neural network

TL;DR

This paper introduces a deep neural network approach to estimate the neutron star equation of state from observational data, providing insights consistent with nuclear models and gravitational wave constraints.

Contribution

The study presents a novel supervised learning method using neural networks to map neutron star observational data to the equation of state, advancing the analysis of dense matter.

Findings

Results align with traditional nuclear model extrapolations

Consistent with gravitational wave tidal deformability bounds

Demonstrates the feasibility of neural networks in astrophysical inference

Abstract

The densest state of matter in the universe is uniquely realized inside central cores of the neutron star. While first-principles evaluation of the equation of state of such matter remains as one of the longstanding problems in nuclear theory, evaluation in light of neutron star phenomenology is feasible. Here we show results from a novel theoretical technique to utilize deep neural network with supervised learning. We input up-to-date observational data from neutron star X-ray radiations into the trained neural network and estimate a relation between the pressure and the mass density. Our results are consistent with extrapolation from the conventional nuclear models and the experimental bound on the tidal deformability inferred from gravitational wave observation.

Tables

Table 1. Table 1: Our neural network architecture in this work. In the input layer 56 neurons correspond to parameters of 14 points of the mass, the radius, and their variances. In the output layer 5 neurons correspond to 5 parameters of the EoS.

Layer	Number of neurons	Activation function
0 (Input)	56	N/A
1	60	ReLU
2, 3	40	ReLU
4 (Output)	5	$\tanh$

Full text

I Introduction

Neutron stars provide us with a natural laboratory to study the densest state of matter in our universe (see Refs. Haensel et al. (2007); Lattimer (2012); Özel and Freire (2016); Baym et al. (2018); Blaschke and Chamel (2018) for recent reviews). The essential ingredient for neutron star structures is the equation of state (EoS) $p=p(\rho)$, i.e., a relation between the pressure $p$ and the mass density $\rho$. It is a longstanding challenge to evaluate the EoS from the first-principles theory.

In the cores of neutron stars the baryon density may reach $\gtrsim 5\rho_{0}$, where $\rho_{0}$ is the normal nuclear density $\rho_{0}\simeq 2.7\times 10^{17}$ kg/m3. At such high density properly dealing with quantum chromodynamics (QCD) is indispensable. Symmetries of QCD imply a speculative duality at high density between hadronic and quark states, called quark–hadron continuity Schäfer and Wilczek (1999). The duality at high density has been confirmed also in a particular limit of large colors of quarks, and the dual state was named quarkyonic matter McLerran and Pisarski (2007). The EoS construction founded on quarkyonic matter has been proposed Fukushima and Kojo (2016); McLerran and Reddy (2018), which is consonant with the phenomenological EoS constructions Masuda et al. (2013); Alvarez-Castillo et al. (2014); Baym et al. (2018); Oter et al. (2019).

Although QCD is the established theory, the first-principles calculations of the EoS have serious problems. Among various theoretical approaches the most powerful method is the lattice-QCD simulation; however, the standard Monte-Carlo algorithm breaks down at finite density, dubbed the sign problem (see Ref. Aarts (2016) for a review). The perturbative QCD (pQCD) calculation is also feasible Kurkela et al. (2010), but it is valid only at asymptotically high density.

Thanks to the recent advances in observations, the quality of the neutron star observables is being improved steadily (see, e.g., Refs. Gendreau et al. (2016); Abbott et al. (2017) for NICER and GW170817). To circumvent the above-mentioned difficulties, current theoretical efforts are directed toward the EoS inference from these stellar observables, especially masses $M$ and radii $R$ (pairs of which are called $M$-$R$ relation). This is mediated by the Tolman–Oppenheimer–Volkoff (TOV) equation, and the mapping from $M$-$R$ to the EoS is in one-to-one correspondence Lindblom (1992). Now, the Bayesian analysis is one standard method to implement such inference Özel et al. (2010); Steiner et al. (2010, 2013); Raithel et al. (2016); *Raithel:2017ity; Alvarez-Castillo et al. (2016). If the number of available observational data is sufficiently large, the likelihood would be well localized such that a choice of the prior distribution scarcely affects the result. In reality, however, the number of data is limited, as tabulated in Refs. Ozel et al. (2016); Özel and Freire (2016); Bogdanov et al. (2016) and plotted in Fig. 2 (Left), and we may not exclude such factors. Hence, it would be complementary to develop an independent analysis based on a different principle than the Bayesian analysis.

Here, we propose a new method to utilize the neural network in the deep learning machinery to estimate the EoS from real observational $M$-$R$ data, as an extension from Ref. Fujimoto et al. (2018). Deep learning provides us with a way to find a regression function for complex nonlinear systems, and there are many physics applications, which include QCD physics Pang et al. (2018); Mori et al. (2018), nuclear physics Niu and Liang (2018), and gravitational waves George and Huerta (2018) (see also Ref. Allen et al. (2019) and references therein). As we explicate below, an advantage to employ the deep learning method lies in the fact that the numerical implementation is straightforward, so we are relatively free from implicit biases.

II Methods

II.1 Compilation of observational data

Ideally, with sufficient computational resources, machine learning would be capable of directly dealing with full multidimensional data from the observation. Figure 2 (Left) shows only a single contour for each neutron star, but the full data is available in the form of the probability distribution as exemplified in Fig. 2 (Right) for (arbitrarily chosen) two representatives out of 14 observations.

In the present work we simplify our analysis by approximately characterizing one probability distribution with four parameters. We project the two-dimensional distribution onto the one-dimensional $M$-axis (and $R$-axis) by integrating over $R$ (and $M$, respectively); in other words, we make marginal distributions with respect to $M$ and $R$. Such marginal distributions are represented by blue shaded shapes outside the frame on Fig. 2 (Right). Then, these two distributions along the $M$-axis and the $R$-axis are fitted by Gaussians as overlaid by red curves. Since each Gaussian has two parameters, namely, the mean and the variance, we sample $2\times 2\times 14=56$ parameters out from the raw $M$-$R$ data of 14 neutron stars. Now, our task is to find a mapping from these 56 observational parameters onto the most likely EoS, i.e., $p=p(\rho)$.

II.2 Training and validation data with fluctuations

We will utilize the neural network to represent such a mapping, and for the optimization, we generate training dataset; many sets of randomly generated EoS and the corresponding observational data. It is important to note that this mapping is not necessarily invertible; even for the same EoS the observational data points may fluctuate according to the probability distributions originating from observational errors. We need to train the neural network to “recognize” that the observational data points could depart from the $M$-$R$ relation.

Here we outline how we generated the training and validation data for our supervised learning. The first step is the EoS generation; we divide a sufficiently wide density range, $[\rho_{0},\,8\rho_{0}]$ in this work, into five segments equally separated in the logarithmic scale, that is, $[\rho_{i-1},\,\rho_{i}]$ with $i=1,2,\dots 5$ and $\rho_{5}=8\rho_{0}$. We then randomly choose an average sound velocity in each segment, $c_{s,i}^{2}$, with a uniform distribution in the causal range $\varepsilon<c_{s,i}^{2}<1-\varepsilon$ (in the natural units $c=1$), where we introduced a regulator, $\varepsilon=0.01$, to avoid singular behavior in solving the TOV equations. Note that the uniform distribution is chosen to cover wide parameter regions efficiently. Now we have 5 pressure values as $p_{i}=p_{i-1}+c_{s,i}^{2}(\rho_{i}-\rho_{i-1})$ for $\rho=\rho_{i}$.

Up to $\rho=\rho_{0}$ we adopt a conventional nuclear EoS, for which we chose SLy4 Douchin and Haensel (2001), one of the standard EoSs for nuclear matter (meaning that $p_{0}=p(\rho_{0})$ is fixed by SLy4), and for $\rho>\rho_{0}$ the pressure is interpolated with a polytrope function, i.e., $p=p(\rho)=K_{i}\,\rho^{\Gamma_{i}}$ for $\rho_{i-1}<\rho<\rho_{i}$, where two parameters, $K_{i}$ and $\Gamma_{i}$, are solved with two boundary conditions, $p_{i}=p(\rho_{i})$ and $p_{i-1}=p(\rho_{i-1})$.

For a given EoS, the $M$-$R$ relation follows from the TOV equations, which we call the genuine $M$-$R$ curve. We randomly sample 14 data points along the genuine $M$-$R$ curve in a region $M>M_{\odot}$ (whose lower bound $M_{\odot}$ is chosen loosely so that the region is large enough to cover masses from the actual observations). Then, the variances of the Gaussian distribution, denoted by $\sigma_{M_{i}}$ and $\sigma_{R_{i}}$, are randomly sampled from the uniform distribution on the ranges $[0,\,M_{\odot})$ and $[0,\,5\,\text{km})$, respectively. These ranges are sufficient for our purpose in view of Fig. 2. The real data are not necessarily centered on the bare data point ($R_{i}^{(0)}$, $M_{i}^{(0)}$), and we shall shift each distribution by $\Delta M_{i}$ and $\Delta R_{i}$ that we chose randomly from the Gaussian distributions with $\sigma_{M_{i}}$ and $\sigma_{R_{i}}$. To summarize the above, one observation for the training data consists of 14 probability distributions of the Gaussian shape whose center is $(R_{i}^{(0)}+\Delta R_{i},\,M_{i}^{(0)}+\Delta M_{i})$ and variances are $\sigma_{R_{i}}$ and $\sigma_{M_{i}}$ along the $R$-axis and the $M$-axis, respectively. For the neural network to learn the correlation between the variances ($\sigma_{R_{i}}$, $\sigma_{M_{i}}$) and how far the actual data is off from the genuine $M$-$R$ curve, we prepare 100 ensembles of different variances for each EoS and sampled 14 data points, and then prepare 100 ensembles of shifts, $\Delta M_{i}$ and $\Delta R_{i}$, for each generated set of variances. This means that we prepared $100\times 100$ ensembles of data for each EoS and sampled 14 data points. For the training dataset we repeated the above process 500 times to cover a wide variety of EoSs; the total training dataset is thus $100\times 100\times 500=5,000,000$ sets of the EoS and the 14 data points. For the validation dataset we generate $1\times 1\times 100$ sets to monitor the convergence and avoid the overfitting; for each step of learning process, we calculated the loss functions for the training and the validation data (see Ref. Fujimoto et al. (2018) for technicalities).

II.3 Neural network design

We specify the setups for the actual calculation. For numerics we employ a Python library, Keras Chollet (2015) using TensorFlow Abadi et al. (2016) backend. The design of our feedforward neural network is summarized in Tab. 1. Our objective is to construct a network that can convert the neutron star data to the EoS parameters, so the input and the output layers have 56 and 5 neurons, respectively. These correspond to 56 parameters of observed 14 neutron stars, ($M_{i}$, $R_{i}$, $\sigma_{M,i}$, $\sigma_{R,i}$) ($i=1,2,\dots,14$), and 5 parameters of the EoS, $c_{s,i}^{2}$ ($i=1,2,\dots,5$). We chose the activation function at the output layer as $\sigma^{(4)}(x)=\tanh(x)$, so that the sound velocity automatically satisfies the causal bound. For hidden layers the activation function is the ReLU, i.e., $\sigma^{(k)}(x)=\max\{0,x\}$ ($k=1,2,3$), which is known to evade the vanishing gradient problem and a standard choice in deep learning LeCun et al. (2015). We implement the loss function by the mean square logarithmic errors (msle). The optimization method of our choice is Adam Kingma and Ba (2014) with the batch size 1000. We initialized neural network parameters with the Glorot uniform distribution Glorot and Bengio (2010).

II.4 Uncertainty estimate from credibility of reproducibility

In our strategy we took care of the probability distribution in the observational side only, but the deduced EoS also has such a probability distribution around the most likely curve. To implement that, instead of randomly generating EoSs, we could have generated some distributions on the $\rho$-$p$ plane and sample fluctuating EoSs according to the generated distribution, which would, however, increase the size of the training dataset tens of thousands larger and require gigantic computational resources.

Here, we employ an alternative practical way to quantify the credibility of the deduced EoS with less efforts. We generate 10 independent training datasets to prepare 10 independent neural network models. For the same real experimental data, those 10 neural network models output 10 deduced EoSs. If a part of the EoS is insensitive to the $M$-$R$ observation, different neural network models would lead to different EoSs in such an unconstrained region. From the dispersion over 10 deduced EoSs, therefore, we can estimate the credibility of our results. Strictly speaking, this dispersion is not the probability distribution of the likely EoS but a measure to quantify how much the same deduced EoS is reproduced with the same analysis. In other words, this measure is to be regarded as the credibility of reproducibility within the present setup of machine learning. If the physical error bar is large, the credibility band would be large, but a small credibility band does not always guarantee small physical error bar. In this sense our uncertainty estimate gives a lower bound. Here we note that the uncertainty estimated in this way accounts for the statistical part (see the band labeled by “10 NNs” in Figs. 2 and 3). Uncertainty including systematics can be quantified by the root-mean-square deviation between the guessed and true values using the validation data (see the band labeled by “validation” in Figs. 2 and 3), as addressed in Ref. Fujimoto et al. (2018). This leads to an uncertainty width of 1.7 km for $R$ at $M=1.4M_{\odot}$ in the $M$-$R$ plane, which is comparable to our inferred width of 1.3 km (68% CL).

III Results and discussions

III.1 Deduced equation of state

In Fig. 2 we present the deduced EoS by the blue line and its credibility by the light blue shade (labeled by “10 NNs”). Uncertainty quantified in a different way is also overlaid by the light red shade in Fig. 2 (labeled by “validation”). Our results are in favor of standard EoSs calculated within the nuclear many-body model, such as APR4 Akmal et al. (1998), BSk20 Goriely et al. (2010), ENG (Dirac–Brueckner–Hartree–Fock method) Engvik et al. (1996), and SLy4 (non-relativistic potential) Douchin and Haensel (2001), some of which are overlaid on Fig. 2. Our results indicate that the constraints from currently observed neutron stars do not have enough resolution to probe a possibility of the first-order phase transition as encoded in QHC18 (hybrid phenomenological construction) Baym et al. (2018). The gray band represents an estimate from the chiral effective theory ($\chi$EFT) Hebeler et al. (2013), and our results lie within this band. In Fig. 2, for reference, we show MS1b (relativistic mean-field) Mueller and Serot (1996), WFF1 (variational) Wiringa et al. (1988), and several other phenomenological EoSs. The Bayesian analyses Steiner et al. (2013); Özel and Freire (2016) are also overlaid in Fig. 2. Note that while Özel & Freire Özel and Freire (2016) and our present analysis use the same astrophysical data, Steiner et al. Steiner et al. (2013) employs eight X-ray sources.

It is an interesting question how the corresponding $M$-$R$ curve looks like because even the knowledge of the existence of the $M$-$R$ curve is not provided to the neural network during the supervised learning by the $M$-$R$ points and the EoS parameters. Figure 3 shows the $M$-$R$ curves corresponding to the EoSs in Fig. 2. We see that our deduced EoS (blue curve) certainly supports massive neutron stars above two solar mass Demorest et al. (2010); *Fonseca:2016tux; Antoniadis et al. (2013); Cromartie et al. (2019).

III.2 Discussions

One may want to know why the uncertainty band of our deduced EoS looks such narrow. A part of the reason lies in the boundary condition in the low density side; we assumed SLy4 for $\rho\leq\rho_{0}$ because up to this density the EoS is well constrained by nuclear properties accessible by terrestrial experiments. So our results should be more precisely regarded as the most likely extrapolation from SLy4 with help of the observational data of 14 neutron stars. It shall be a future work to inspect possible bias effect induced by such a choice of the EoS up to $\rho_{0}$. Also we can in principle remove such an assumption by extending the neural network architecture including data from nuclear physics experiments (e.g., symmetry energy; see discussions in Ref. Gandolfi et al. (2012)) on top of neutron star data. Such a global analysis over all available data from astrophysics and nuclear physics experiments would be an ambitious future challenge.

At the same time, we can argue from a different point of view. The light blue band in Fig. 2 may look small at first glance, but the resolution is not yet good enough to justify/falsify a first-order phase transition. In view of the light blue band in Fig. 3, the corresponding uncertainty for the $M$-$R$ relation is $\sim 1\,\text{km}$.

Another important physical quantity derived from the EoS is the sound velocity, $c_{s}$, which is plotted in Fig. 5. Interestingly, the deduced sound velocity is smaller than $1/\sqrt{3}$ (the conformal limit value, viz., a naive upper bound for massless ultrarelativistic systems) for $\rho\lesssim 2\rho_{0}$. With further increasing $\rho>2\rho_{0}$, the sound velocity becomes significantly greater than $1/\sqrt{3}$, and eventually the increasing behavior is saturated beyond $\sim 4\rho_{0}$. Such a sharp increase of the sound velocity around $2\rho_{0}$ appears in accordance with the recent studies Bedaque and Steiner (2015); Tews et al. (2018). At even higher densities $>4\rho_{0}$ it is likely that the sound velocity starts decreasing and approaches the conformal limit of asymptotically free quarks and gluons. This in turn implies that the saturation seen around $4\rho_{0}$ hints a transition to weakly interacting quark matter.

Finally, we shall confirm that our deduced EoS is consistent with the recent gravitational wave experiment, specifically the tidal deformability $\Lambda$. Once the EoS is given, the tidal deformability can be calculated following the method outlined in the Ref. Hinderer et al. (2010). The experimentally determined bound, $\Lambda(1.4M_{\odot})=190^{+390}_{-120}$ Abbott et al. (2018), is indicated by a red bar in Fig. 5. Our deduced EoS leads to $\Lambda(1.4M_{\odot})=320\pm 120$ which is entirely consistent with the GW170817 measurement within the error bar as it should be. For the moment we utilize the tidal deformability as a benchmark test, but in the future the neural network should be better designed to implement what is called the multi-messenger observation, inclusive of gravitational waves as well as electromagnetic waves.

IV Summary

In this work we successfully utilized a new method based on the machine learning to infer neutron star EoS in a way independent of the existing methods. In our method the deep neural network can deal with nonlinear mapping from masses $M$ and radii $R$ of neutron stars to the EoS parameters. The neural network model was optimized with training datasets of size 5,000,000, and the convergence was monitored with an independent validation dataset. In this way, from available $M$-$R$ data from 14 neutron stars, we deduced an EoS to find it compatible with the conventional nuclear EoS and the currently existing constraints. Dealing with two-dimensional $M$-$R$ distribution for the neural network input would be an important extension for the future. Still, our successful results would be a first step toward further refinements to incorporate the gravitational wave measurements and nuclear physics experiments. Machine learning’s advantage lies in handling such a large set of complex data, and this direction deserves investigations.

Acknowledgements.

We thank Toru Kojo, Andrew Steiner and Wolfram Weise for encouraging discussions. K. F. was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. 18H01211.