Nonparametric model for the equations of state of neutron star from deep
  neural network
It is of great interest to understand the equation of state (EOS) of the
neutron star (NS), whose core includes highly dense matter. However, there are
large uncertainties in the theoretical predictions for the EOS of NS. It is
useful to develop a new framework, which is flexible enough to consider the
systematic error in theoretical predictions and to use them as a best guess at
the same time. We employ a deep neural network to perform a non-parametric fit
of the EOS of NS using currently available data. In this framework, the
Gaussian process is applied to represent the EOSs and the training set data
required to close physical solutions. Our model is constructed under the
assumption that the true EOS of NS is a perturbation of the relativistic
mean-field model prediction. We fit the EOSs of NS using two different example
datasets, which can satisfy the latest constraints from the massive neutron
stars, NICER, and the gravitational wave of the binary neutron stars. Given our
assumptions, we find that a maximum neutron star mass is $2.38^{+0.15}_{-0.13}
M_\odot$ or $2.41^{+0.15}_{-0.14}$ at $95\%$ confidence level from two
different example datasets. It implies that the $1.4 M_\odot$ radius is
$12.31^{+0.29}_{-0.31}$ km or $12.30^{+0.35}_{-0.37}$ km. These results are
consistent with results from previous studies using similar priors. It has
demonstrated the recovery of the EOS of NS using a nonparametric model.
Wenjie Zhou, Jinniu Hu, Ying Zhang, Hong Shen
# Nonparametric model for the equations of state of neutron star from deep neural network

###### Abstract

It is of great interest to understand the equation of state (EOS) of the neutron star (NS), whose core includes highly dense matter. However, there are large uncertainties in the theoretical predictions for the EOS of NS. It is useful to develop a new framework, which is flexible enough to consider the systematic error in theoretical predictions and to use them as a best guess at the same time. We employ a deep neural network to perform a non-parametric fit of the EOS of NS using currently available data. In this framework, the Gaussian process is applied to represent the EOSs and the training set data required to close physical solutions. Our model is constructed under the assumption that the true EOS of NS is a perturbation of the relativistic mean-field model prediction. We fit the EOSs of NS using two different example datasets, which can satisfy the latest constraints from the massive neutron stars, NICER, and the gravitational wave of the binary neutron stars. Given our assumptions, we find that a maximum neutron star mass is \(2.38^{+0.15}_{-0.13}M_{\odot}\) or \(2.41^{+0.15}_{-0.14}\) at 95% confidence level from two different example datasets. It implies that the \(1.4M_{\odot}\) radius is \(12.31^{+0.29}_{-0.31}\) km or \(12.30^{+0.35}_{-0.37}\) km. These results are consistent with results from previous studies using similar priors. It has demonstrated the recovery of the EOS of NS using a nonparametric model.

Neutron Star, Deep neural network, Gaussian process regression

## 1 Introduction

Neutron stars, remnants of very massive stars at the end of their lifecycle, are one of the most compact objects in the universe, attracting a lot of attention within the fields of astrophysics and nuclear physics (Oertel et al., 2017). Rapid developments in space observation technologies and gravitation-wave detection have proven advantageous to the measurement of neutron star properties, such as mass, radius, and tidal deformability. Three massive neutron stars with masses of around \(2M_{\odot}\), PSR J1614-2230 (Demorest et al., 2010; Fonseca et al., 2016; Arzoumanian et al., 2018), PSR J0348+0432 (Antoniadis et al., 2013), and PSR J0740+6620 (Cromartie et al., 2020) have been discovered in the past decade. The gravitational wave from a binary neutron star merger, the GW170817 event, was first detected by the LIGO and Virgo collaborations in 2017, providing a constraint on the tidal deformability of a neutron star at \(1.4M_{\odot}\)(Abbott et al., 2017, 2018, 2019). Furthermore, simultaneous measurements of the mass and radius of PSR J0030+0451 and PSR J0740+6620 were recently analyzed by the Neutron Star Interior Composition Explorer (NICER) (Riley et al., 2019; Miller et al., 2019; Riley et al., 2021; Miller et al., 2021).

These studies have improved our knowledge of neutron stars, providing insight into their interior structure and components. A neutron star can be divided into the atmosphere, outer crust, inner crust, outer core, and inner core regions. Its properties are strongly dependent on the equation of state (EOS) of dense nuclear matter (Lattimer & Prakash, 2000; Glendenning, 2001; Weber, 2005; Lattimer & Prakash, 2007; Baym et al., 2018). In the core region, the density approaches 5-10 times the nuclear saturation density. Therefore, in this high-density region, the EOS plays an essential role in investigations, yet it cannot be well determined by current terrestrial methodologies.

Conventionally, the EOS of neutron stars can be extrapolated by the nuclear many-body approaches, such as the density functional theory (Ring, 1996; Bender et al., 2003; Meng et al., 2006; Stone & Reinhard, 2007;
Niksic et al., 2011; Dutra et al., 2012, 2014) and _ab initio_ method (Akmal et al., 1998; Van Dalen et al., 2004; Sammarruca, 2010; Sammarruca et al., 2012; Wei et al., 2019; Wang et al., 2020), which can describe the ground-state properties of finite nuclei and nuclear saturation properties very well. However, there are large uncertainties, when these methods are extended to calculate the high-density EOS. They generate many kinds of neutron star mass-radius relations. Furthermore, the isospin dependence of EOS, i.e., the symmetry energy effect, is strongly correlated to the radii of low-mass neutron stars (Li et al., 2008; Bao et al., 2014; Bao and Shen, 2015; Sun, 2016; Ji et al., 2019; Li et al., 2019; Hu et al., 2020). With present observations of neutron stars, a smaller slope of symmetry energy, \(L\) is preferred. Meanwhile, several exotic hadronic degrees of freedom and/or hadron-quark transitions may appear in the core region of a neutron star because of the phase diagram of strong interaction (Yang and Shen, 2008; Xu et al., 2010; Chen et al., 2013; Orsaria et al., 2014; Wu and Shen, 2017; Ju et al., 2021; Huang et al., 2022). Hence, it is very difficult to generate a unified EOS in a self-consistent theoretical framework.

Recently, the data-driven methodologies, such as Bayesian inference (Ozel et al., 2010; Raithel et al., 2017; Steiner et al., 2010; Alvarez-Castillo et al., 2016; Miao et al., 2021), deep neural network (DNN) (Fujimoto et al., 2018, 2020, 2021; Farrell et al., 2022; Ferreira et al., 2022), nonparametric EOS representation (Landry and Essick, 2019; Essick et al., 2020, 2020), support vector machines (Murarka et al., 2022; Ferreira and Providencia, 2021), and so on, have been introduced to generate the possible EOSs using the latest observables of neutron stars. In Bayesian inference, the EOS is parameterized and the corresponding parameters are obtained with a marginal likelihood estimation on the posterior probability in terms of model parameters (Ozel et al., 2010). Fujimoto et al. proposed a scheme that can map the finite observation data of neutron stars onto the EOS with a feed-forward DNN. They present the EOS as a polytropic function with different speeds of sound at distinct density segments (Fujimoto et al., 2018). To avoid the limitations of a parametric EOS, Landry _et al._ developed a nonparametric method to generate the EOS from the observables of gravitation waves by combining the Gaussian process and Bayesian inference methods, where the EOS of the neutron star is represented by the Gaussian process with finite points (Landry and Essick, 2019). The matching between the EOS and neutron star observations was carried out by the Bayesian inference. Han et al. also reconstructed the EOS of a neutron star using another Bayesian nonparametric inference method where the EOS was produced by the neural network with a sigmoid type as the activation function (Han et al., 2021).

In this work, a new machine learning methodology is proposed to reconstruct a nonparametric model for the EOSs of neutron stars based on the scheme proposed by Fujimoto et al., where the complete EOS is generated by the Gaussian process regression method with finite data points about the pressure-energy relation. A DNN is trained with the constraints of neutron star mass-radius relations, the masses of the heavy neutron stars, and the measurements of NICER. In Section 2, the framework of the Gaussian process regression method and the construction of the DNN is given in detail. The nonparametric EOS model of neutron stars generated by the DNN is shown in Section 3. A summary is presented in Section 4.

## 2 Gaussian Process Regression and Neural Network

### Gaussian Process Regression

The Gaussian process (GP) (Huang et al., 2022; Williams and Rasmussen, 2006), a random process, is a series of normal distributions of random variables in an index set combination. If the set of random variables \(\{f(x):x\in\chi\}\) is taken from the GP with the mean function \(m(x)\) and the covariance function \(k(x_{1},x_{2})\), the corresponding random variables \(f(x_{i})\) satisfy the multivariate Gaussian distribution for any finite set, \([x_{1},\cdots,x_{m}]\in\chi\),

\[\left[\begin{array}{c}f(x_{1})\\ \vdots\\ f(x_{m})\end{array}\right]\sim\ \ \mathcal{N}\left(\left[\begin{array}{c}m(x_{1})\\ \vdots\\ m(x_{m})\end{array}\right],\left[\begin{array}{ccc}k(x_{1},x_{1})&\cdots&k (x_{1},x_{m})\\ \vdots&\ddots&\vdots\\ k(x_{m},x_{1})&\cdots&k(x_{m},x_{m})\end{array}\right]\right), \tag{1}\]

which can be simply expressed as

\[f(\cdot)\sim GP(m(\cdot),k(\cdot,\cdot)). \tag{2}\]

All linear combinations of random variables in GP obey the normal distribution. For each finite-dimensional set, its probability density function on the continuous exponential set is the Gaussian measurement of all random variables. Therefore, it is regarded that the infinite-dimensional set can be generalized by the extension of the multivariate Gaussian distribution.

Hence, the GP can be applied to solve a normal regression problem,

\[y^{(i)}=f(x^{(i)})+\epsilon^{(i)}, \tag{3}\]

where \(X\) is defined as the training set and its components \((x^{(1)},...,x^{(m)})\) are independently and identically
distributed with unknown distribution. \(\epsilon^{(i)}\) is an independent noise variable, which is also given by a normal distribution with variance \(\sigma^{2}\), \(N(0,\sigma^{2})\). This scheme is called the Gaussian process regression (GPR) method. Usually, it is assumed that \(f\) follows the GP with a mean value of zero for notation simplicity,

\[f(\cdot)\sim GP(0,k(\cdot,\cdot)). \tag{4}\]

The test set \(X^{*}=(x^{(1*)},...,x^{(m*)})\), has the same independent co-distribution as \(X\), marked as \(X\to X^{*}\). Therefore, the posterior distribution \(p(y^{*}|X,\;X^{*})\) is predicted in GPR as the Gaussian distribution of the results, which is different from the general linear regression.

According to the properties of the GP, a joint distribution of the training and test sets is obtained,

\[\left[\begin{array}{c}\vec{f}\\ \vec{f}^{*}\end{array}\right]\Bigg{|}\,X,X^{*}\sim\mathcal{N}\left(\vec{0}, \left[\begin{array}{cc}K(X,X)&K(X,X^{*})\\ K(X^{*},X)&K(X^{*},X^{*})\end{array}\right]\right) \tag{5}\]

where the matrix elements \(K(X^{A},X^{B})_{i,j}=k(x_{i}^{A},x_{j}^{B})\). In the GP, the covariance function \(k_{ij}\) is also called the kernel function. The standard choice is the squared-exponential kernel,

\[k_{se}(x_{1},x_{2})=\sigma^{2}\exp\left(-\frac{||x_{1}-x_{2}||^{2}}{2l^{2}} \right). \tag{6}\]

Meanwhile, their noises obey similar distributions,

\[\left[\begin{array}{c}\vec{\epsilon}\\ \vec{\epsilon^{*}}\end{array}\right]\sim\mathcal{N}\left(\vec{0},\left[ \begin{array}{cc}\sigma_{wn}^{2}I&\vec{0}\\ \vec{0}^{T}&\sigma_{wn}^{2}I\end{array}\right]\right). \tag{7}\]

Here, \(\sigma_{wn}^{2}\) is the hyper-parameter corresponding to white noise, which is different from the signal variance parameter, \(\sigma\) in \(k_{se}\). The summation of two independent multivariate Gaussian variables is still a multivariate Gaussian variable,

\[\left[\begin{array}{c}\vec{y}\\ \vec{y^{*}}\end{array}\right]\Bigg{|}\,X,X^{*}=\left[\begin{array}{cc}\vec{ \epsilon}\\ \vec{\epsilon^{*}}\end{array}\right]+\left[\begin{array}{c}\vec{\epsilon} \\ \vec{\epsilon^{*}}\end{array}\right]\sim \tag{8}\] \[\mathcal{N}\left(\vec{0},\left[\begin{array}{cc}K(X,X)+\sigma_ {wn}^{2}I&K(X,X^{*})\\ K(X^{*},X)&K(X^{*},X^{*})+\sigma_{wn}^{2}I\end{array}\right]\right)\]

Based on the properties of multivariate Gaussian distribution, the conditional distribution over the unknown \(y^{*}\) is,

\[y^{*}|y,\;X,\;X^{*}\sim\mathcal{N}\left(\mu^{*},\Sigma^{*}\right), \tag{9}\]

where,

\[\mu^{*} = K(X^{*},X)(K(X,X)+\sigma_{wn}^{2}I)^{-1}\vec{y},\] \[\Sigma^{*} = K(X^{*},X^{*})-K(X^{*},X)\] \[(K(X,X)+\sigma_{wn}^{2}I)^{-1}K(X,X^{*}).\]

\(\mu^{*}\) and \(\Sigma^{*}\) are the mean and covariance functions of the probability distribution for our prediction results, respectively. Therefore, given the hyper-parameters \(\sigma\) and \(l\) in the kernel function, a probability distribution describing the whole test set by the GPR method can be obtained. In principle, the mean function should be selected as the "actual data curve". However, it is strongly dependent on the hyper-parameters, \(\sigma\) and \(l\) that are determined by maximizing the marginal log-likelihood, defined as,

\[\log p(\mathbf{y}|\sigma,l) = \log\mathcal{N}(0,K_{yy}(\sigma,l))\] \[= -\frac{1}{2}\mathbf{y}^{T}K_{yy}^{-1}\mathbf{y}-\frac{1}{2}\log|K_{yy}|- \frac{N}{2}\log(2\pi),\]

where \(K_{yy}=K(X^{*},X^{*})\). Therefore, with a small number of data points, a relatively reasonable EOS curve and its confidence range can be predicted in the framework of the GPR method.

The direct matching between the EOS of a neutron star, i.e., the pressure-energy relation, and the observables of a neutron star may generate nonphysical solutions, such as the speed of sound of neutron star matter being less than zero or larger than the speed of light, \(c_{s}<0\) or \(c_{s}>c\), or the energy density becoming less than zero in some extreme conditions.

Recently, a new intermediate variable \(\phi\) was proposed to construct the EOS of a neutron star (Lindblom, 2010; Landry & Essick, 2019). \(\phi\) is defined as,

\[\phi=\textbf{log}\left(c^{2}\frac{d\epsilon}{dp}-1\right). \tag{12}\]

It avoids the aforementioned weird behaviors, as when \(\phi\in\textbf{R}\), the speed of sound obeys \(0\leq c_{s}^{2}=dp/d\epsilon\leq c^{2}\), which automatically satisfies the physical requirements. When \(p>0\), the \(\epsilon>0\) can be kept. Due to the large pressure magnitude of pressure, the \(\phi\) is regarded as a function of \(\log p\) so that it is easier to determine the hyper-parameters. Therefore, Eq. (12) will be expressed as,

\[\phi=\textbf{log}\left(\partial\textbf{log}\epsilon\frac{e^{\textbf{log} \epsilon}}{p}c^{2}-1\right), \tag{13}\]

where \(\partial\textbf{log}\epsilon=\left.\frac{\partial\log\epsilon}{\partial\log p }\right|_{p=p_{i}}\).

In the training set, \(n\) data points (\(\phi_{i},\log p_{i}\)) are randomly chosen. Once the optimal hyper-parameters are obtained by GPR, the continuum \(\phi-\log p\) curve can be generated. The corresponding EOS of the neutron star, \(\epsilon(p)\) is provided by numerically integrating

\[\frac{\partial\epsilon}{\partial p}=\frac{1+e^{\phi}}{c^{2}}. \tag{14}\]
### DNN method

In the available investigations of the structure of neutron stars, the EOS of neutron star matter was first calculated by either the nuclear many-body method or the parameterization function under the conditions of \(\beta\)-equilibrium and charge neutrality. The EOS was then input to the Tolman-Oppenheimer-Volkoff (TOV) equation (Tolman, 1939; Oppenheimer and Volkoff, 1939), which describes a spherically symmetric and isotropic star in a static gravitational field with general relativity.

\[\begin{split}\frac{dp}{dr}=&-\frac{G\epsilon(r)m(r )}{c^{2}r^{2}}\left[1+\frac{p(r)}{\epsilon(r)}\right]\\ &\times\left[1+\frac{4\pi r^{3}p(r)}{m(r)c^{2}}\right]\left[1- \frac{2Gm(r)}{c^{2}r}\right]^{-1}\\ \frac{dm}{dr}=&\frac{4\pi r^{2}\epsilon(r)}{c^{2}}, \end{split} \tag{15}\]

where \(r\) is the radial coordinate, representing the distance to the center of the star. The functions \(p(r)\) and \(\epsilon(r)\) are pressure and energy density (i.e., mass density), respectively. We can easily integrate these differential equations starting at \(r=0\), with the initial condition \(p(r=0)=p_{c}\). When it is integrated into the surface of the neutron star, i.e., the radius \(R\) and \(p(r=R)=0\), then \(M=m(R)\) corresponds to the total mass of the neutron star. Therefore, a continuum mass-radius (\(M\)-\(R\)) relation of a neutron star can be generated by the TOV equation.

A functional mapping between the EOS space and \(M\)-\(R\) space is constructed through the above framework, in a process called "TOV mapping". In principle, such mapping is invertible; thus, there should be a relevant inverse mapping (Lindblom, 1992), where the EOS can be uniquely reconstructed from the observed \(M\)-\(R\) relationship of the neutron star. However, in actuality, the complete \(M\)-\(R\) curve cannot be directly obtained from the observed data due to the discontinuities and uncertainties inherent in neutron star observations (Fujimoto et al., 2021). Therefore, a more likely EOS can be inferred from the neutron star observations with uncertainties.

The DNN is a powerful machine learning method to connect the EOS with observed data, following the idea of Fujimoto et al. (Fujimoto et al., 2021). The neural network (NN) is a representation of the fitting parameters of a function. Deep learning, e.g., the machine learning method using a DNN, is a process of optimizing the parameters contained in the function represented by an NN. Deep learning can be divided into supervised learning and unsupervised learning. The supervised learning that we adopted needs to have specific inputs and outputs before it can complete the fitting process with the training data (i.e., regression).

Compared with general fitting methods, the advantage of deep learning lies in the generalization properties of NNs. It does not need to rely on any prior knowledge about the proper form of the fitting function. Due to a large number of neurons (and neuron layers) and fitting parameters, an NN with a sufficient number of neurons can generate any continuous function (Cybenko, 1989; Hornik, 1991).

The model function of a feed-forward NN can be expressed as,

\[\begin{split}\mathbf{y}=f(\mathbf{x}|&\{W^{(1)},b^{(1)}, \cdots,W^{(l)},b^{(l)},\cdots,\\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad
of \(\mathcal{B}\), and the approximate derivative is,

\[\frac{\partial\mathcal{L}(W^{(l)})}{\partial W^{(l)}}\approx\frac{1}{|\mathcal{B }|}\sum_{n=1}^{|\mathcal{B}|}\frac{\partial\ell\left(\boldsymbol{y}_{n},f\left( \boldsymbol{x}_{n}|\left\{W^{(l)},b^{(l)}\right\}_{l}\right)\right)}{\partial W ^{(l)}}, \tag{18}\]

where, the batch size \(|\mathcal{B}|\) represents the number of sample points in \(\mathcal{B}\). Since each optimal choice varies from case to case, its error will be shown later as a part of our estimations on the EOS confidence. The epoch denotes the number of scans of the entire training data set \(\mathcal{D}\). Parameters are updated with each small batch, so an epoch is equivalent to iterating \(|\mathcal{D}|/|\mathcal{B}|\) small batches of data until all iterations are completed. In addition, the derivative \(\frac{\partial\ell}{\partial W}\) that appeared in Eq. (18) was calculated by the back-propagation method.

In this training process, the mean square logarithmic error (msle) is regarded as the loss \(\ell(\boldsymbol{y},\boldsymbol{y}^{\prime})\) in Eq. (17),

\[\ell_{\text{msle}}(\boldsymbol{y},\boldsymbol{y}^{\prime})\equiv|\log \boldsymbol{y}-\log\boldsymbol{y}^{\prime}|^{2}. \tag{19}\]

With a loss function, our NN can begin the basic training. The parameter initialization of NN will be discussed later, in detail.

Therefore, it is useful to compare our method with other methods proposed to generate the EOS of neutron star. In the present framework, the fitted EOSs are obtained by the DNN. The neutron star observation data is chosen as the input layer, while the constraint EOS is set up as the output layer. The training process is finished with the observations' likelihoods and the EOS priors generated by the theoretical model. The EOSs in the priors and the output layer are presented by several discretized points in \(\phi\)-function to satisfy the constraint of the speed of sound and are smoothly connected by GP. On the other hand, the EOSs in the work of Fujimoto et al. were parameterized as a polytrope function dependent on the speeds of the sound of neutron star matter. Furthermore, the fitted EOSs in the work of Landry and Essick were produced by Bayesian inference with a set of nuclear-theoretic models.

## 3 The Numerical Details and Results

To prepare the training data set, the EOSs from relativistic mean-field (RMF) models were used to obtain the generation interval of GPR fitting data points. Nine RMF parameterizations were selected: BigApple, DD2, DDLZ1, DDME1, DDME2, DDMEX, NL3, PKDD, and TW99 (Fattoyev et al., 2020; Typel et al., 2010; Wei et al., 2020; Niksic et al., 2002; Lalazissis et al., 2005; Taninah et al., 2020; Lalazissis et al., 1997; Long et al., 2004). All of these RMF parameter sets can provide neutron stars, whose maximum masses are larger than \(2.0M_{\odot}\)(Huang et al., 2020). The EOS from the NL3 set generated a maximum mass of neutrons star around \(2.78M_{\odot}\).

The \(\epsilon\)-\(p\) relation in the EOS was transferred into the \(\phi\)-\(\ln p\) function, where \(\ln p\) is the natural logarithm of pressure. After calculating the means and variances of the \(\phi\)-\(\ln p\) relations from the above nine EOSs, it was found that their mean value is very close to the EOS from the DDME1 set (Niksic et al., 2002). To investigate the stability of initial values in the present framework, two schemes were adopted to generate the fitting interval with the GPR method:

1. _Scheme 1_ - After obtaining the mean and variance of \(\phi\)-\(\ln p\) functions from nine RMF parameter sets, the 95% confidence interval of the variance was selected as the generation range of \(\phi_{i}\). As shown in panel (a) of Fig. 2, this interval encloses all EOSs from the RMF model.
2. _Scheme 2_ - The \(\phi\)-\(\ln p\) function provided by DDME1 set was regarded as the standard, and \(\phi\pm 0.3\phi\) are chosen as the upper and lower bounds of the generation range of \(\phi_{i}\). Such an interval is consistent with the one obtained by scheme 1, to a large extent.

In Fig. 3, the corresponding \(\epsilon\)-\(p\) relations of scheme 1 and 2 are compared to the model-informed and model-agnostic priors in the Bayesian inference method by Landry and Essick (Landry and Essick, 2019). The \(\epsilon\)-\(p\) relations from scheme 1 and scheme 2 in the present work are almost identical, which are also consistent with the model-informed prior. Since all of them are more strictly constrained by the theoretical EOSs. On the contrary,

Figure 1: The NN flow chart of present framework.
the model-agnostic prior has a loose boundary. It may consider more range of plausible EOSs.

To produce an EOS of neutron stars (including the high-density region) with the GPR method and aforementioned schemes, seven pressure points \(\ln p_{i}\) (\(i=1,~{}2,~{}\cdots,~{}7\)) were selected, with the same interval, in the range \(\ln p\in[1,7]\). \(\phi_{i}\) was randomly generated in the training interval at each \(\ln p_{i}\) point as an initial data set (\(\phi_{i},~{}\ln p_{i}\)). The EOS below nuclear saturation density was chosen as the one from the SLy4 set. A smooth and continuous \(\phi(\ln p)\) function is fitted by the GPR method, where the hyper-parameters, \(l\) and \(\sigma\) are obtained by maximizing the marginal log-likelihood, as shown in Eq. (11). Furthermore, the star point, \(\phi_{1}=\phi(\ln p=1)\) was fixed as the magnitude from the DDME1 parameter set.

The \(M\)-\(R\) relation of a neutron star can be calculated using the EOS from the GPR method by solving the TOV equation. In the present framework, the training data set of the DNN should assemble the points on the \(M\)-\(R\) curve, which correspond to the observables. The method proposed by Fujimoto et al. (Fujimoto et al., 2021) is used in this work to generate training data. Firstly, the maximum masses of neutron stars less than \(2.2M_{\odot}\) and the \(M\)-\(R\) relations that did not satisfy the radii constraints of PSR J0740+6620 and PSR J0030+0451 (Miller et al., 2019, 2021) were excluded from the training data. Then, 14 points in the mass regions, \([M_{\odot},M_{max}]\) on the \(M\)-\(R\) curve were randomly chosen as "the original data points" \((M_{i},R_{i})\) to simulate the real observations of the 14 available neutron stars. To consider the errors in the observations, the variances of the Gaussian distributions about the mass and radius, \(\sigma_{M_{i}}\) and \(\sigma_{R_{i}}\), were randomly taken from the uniform distribution in the ranges, \([0,M_{\odot}]\) and \([0,5\mathrm{km}]\). The deviations of mass and radius \((\Delta M_{i},\Delta R_{i})\) were calculated by the Gaussian distribution with the variances of \(\sigma_{M_{i}}\) and \(\sigma_{R_{i}}\). Finally the "real data point" \((M_{i}+\Delta M_{i},R_{i}+\Delta R_{i})\) was obtained. The set \((M_{i}+\Delta M_{i},~{}R_{i}+\Delta R_{i},~{}\sigma_{M_{i}},~{}\sigma_{R_{i}})\) can be compared to the observational data of neutron stars.

A group of \(i=14\) data points \((M_{i},~{}R_{i})\) was selected from the \(M\)-\(R\) curve generated by each EOS, and \(j=100\) groups of different variances \((\sigma_{M_{ij}},\sigma_{R_{ij}})\) were randomly sampled for each \(M_{i}\)-\(R_{i}\) data point. Later, \(k=100\) groups of deviations, \(\Delta M_{ijk}\) and \(\Delta R_{ijk}\) were provided by each variance set, \((\sigma_{M_{ij}},\sigma_{R_{ij}})\). In this way, \(100\times 100\) sets of data for each EOS were prepared and 14 data points were sampled. The above process was repeated by 500 times to include as wide a range as possible, resulting in \(500\times 100\times 100=5,000,000\) sets, where one set includes 14 data points.

Figure 3: The corresponding \(\epsilon-p\) relations of scheme 1 and 2 in Figure 2 and the model-informed and model-agnostic priors in the Bayesian inference method by Landry and Essick (Landry and Essick, 2019).

Figure 2: The generation range of \(\phi\)-\(\ln p\). We will randomly select points within this range and then use GPR method to generate EOS. In panel (a ) the nine EOSs are treated to obtain the mean \(\mu\) and variance \(\sigma\), whose 95% confidence interval is taken to obtain the fitting range. In panel (b), the generation range is based on DDME1 curve, with a fluctuation of 0.3.
For the architecture of the NN, the Python library, Keras (Chollet et al., 2015) was employed, with TensorFlow (Abadi et al., 2016) as the backend. The number of NN layers, their corresponding neurons, and the activation functions are shown in Table 1. The hyperbolic tangent function of the output layer makes the results fall between \((-1,1)\), speeding up the training. The msle is chosen as the loss function, given in Eq. (19). The optimization method was Adam (Kingma and Ba, 2014) by taking the batch size as \(1000\). The default initialization NN argument was the Glorot Uniform distribution (Glorot and Bengio, 2010).

The DNN models for a full training set of \(5,000,000\) data were compare with a random sampling of \(1,000,000\) data in the training set, giving similar results, but with the latter greatly improving the training efficiency. In addition, for all models, the changes in loss functions for the training of epoch were almost identical. The loss functions estimated for the validation data and training data are shown as an example in Fig. 4. When the epoch \(>10\), the verification loss is consistent with the training loss, whereas when the epoch \(>100\), the verification loss is stable. Therefore, each DNN model was trained with \(1,000,000\) data. The validation set was taken as the \(10,000\) sets from the rest \(4,000,000\) sets to check the convergence. Once the epoch \(=100\), the model was considered finished.

Due to the differences in initial input and training data, there was some uncertainty about the output results of the DNN. Therefore, the process was repeated \(100\) times to generate \(100\) independent DNN models. The uncertainties in the training results were estimated from the fitted \(100\) EOSs. In Fig. 5, \(200\) relations about \(\phi\)-\(\ln p\) from scheme 1 in panel (a) and scheme 2 in panel (b) are reconstructed through the training data of the DNN. Each curve is smoothly connected with seven output points by the GPR method, as shown in the inserts. It was found that most of these curves have similar pressure-dependence behaviors. Their differences increase in the high-density region due to the observation discrepancies associated with the \(14\) neutron stars.

The \(\phi\)-\(\ln p\) relations must be converted to the \(\epsilon\)-\(p\) function by integrating the Eq. (14) to obtain the EOS of the neutron star. In Fig. 6, the neutron star EOSs with the \(68\%\) and \(95\%\) confidence levels from the DNN with scheme 1 in panel (a) and scheme 2 in panel (b) are shown and compared to those joint constraints from the GW170817 and GW190814 events (Abbott et al., 2020) and the EOS from DDME1. In the inserts, the original \(200\) EOSs from the DNN training are plotted. To analyze the uncertainties of the EOSs, it was assumed that the pressures at each energy density from the machine learning model satisfy the Gaussian distribution. Therefore, the mean EOS was obtained as the dashed curve with the dark blue shadow representing the \(68\%\) confidence level and the light blue shadow, the \(95\%\), respectively. In the low-density region, our estimations are consistent with the joint constraints on the EOS from the GW170817 and GW190814 events. With density increasing, present EOSs are softer than the joint constraints, since the maximum masses of the \(14\) neutron stars are just around \(2M_{\odot}\). Furthermore, the fitted EOS differs slightly from the EOS of DDME1 in scheme 2, despite this being regarded as the mean value of the training data. In the mediate region of energy density, the EOS generated by the DDME1 is harder than the fitted one, since the radius of the neutron star from DDME1

\begin{table}
\begin{tabular}{c|c|c} \hline \hline Layer & Number of neurons & Activation function \\ \hline \(1\)(Input) & \(56\) & N/A \\ \hline \(2\) & \(60\) & ReLU \\ \hline \(3\) & \(40\) & ReLU \\ \hline \(4\) & \(40\) & ReLU \\ \hline \(5\)(Output) & \(6\) & tanh \\ \hline \hline \end{tabular}
\end{table}
Table 1: The setup of present DNN. The number of input and output neurons can be modified according to different network conditions. Here, the number of neurons at output layer is \(6\), because \(\phi(\ln p=1)\) has been fixed as the value obtained from DDME1 set.

Figure 4: The Loss probabilities as functions of epoch with the training data and validation data.
is a little larger when compared with the observations of the 14 neutron stars, as shown later. These results demonstrate that the EOS of the present framework is independent of the initial input of the training set.

Here, it must be emphasized that the inconsistencies in EOSs fitted by LIGO-Virgo-KAGRA (LVK) collaborations from GW170817 and GW190814 events, and present work are generated by the different theoretical frameworks and priors. In the LVK analysis, the EOSs in the priors were given by the spectral representation and are determined by the adiabatic index \(\Gamma\) as shown in Refs. (Read et al., 2009) and (Lindblom, 2010). The EOS parameters of the prior ranges in LVK were choices from the 34-neutron star matter EOSs, including the PAL6, APR1-4, WFF1-3, MS1-2, and so on (Read et al., 2009). The maxim masses of the neutron star from these EOSs are in the range of \(1.47\sim 2.78M_{\odot}\) and the radii at \(1.4M_{\odot}\) are \(9.36\sim 15.47\) km. Correspondingly, the prior of EOSs space in the present framework is taken from the 9 RMF parameter sets, which only can generate the maximum masses of the neutron stars from \(2.0\sim 2.4M_{\odot}\). Therefore, the harder EOSs were fitted by LVK at high-density regions.

Once the EOS of the neutron star were determined, its \(M\)-\(R\) relation was obtained by solving the TOV equation. The \(M\)-\(R\) relations from our deduced EOSs are plotted in Fig. 7, with 68% (dark blue) and 95% (light blue) confidence levels. The corresponding \(M\)-\(R\) distributions of the observed 14 neutron stars are given as contour plots. The masses of massive neutron stars, PSR J0348+0432, PSR J0740+6620, and PSR J1614-2230; the secondary compact object of the GW190814 event; and the radii of PSR J0030+0451 and PSR J0740+6620 from the NICER are given and compared. The fitted EOSs from schemes 1 and 2 nicely reproduce the neutron star observations and are able to generate massive neutron stars. Their radii are consistent with the results of the 14 observed neutron stars and the mass-radius simultaneous measurements from NICER. Furthermore, the \(M\)-\(R\) relation from the DDME1 set is shown as a solid line, which was chosen as the mean value to generate the training data set in scheme 2. Its

Figure 5: The 200 DNN models about \(\phi\)-\(\ln p\) from schemes 1 and 2.

Figure 6: The EOSs from the nonparametric machine learning methods with scheme 1 and 2 and comparing to those from the joint constraints from GW170817 and GW190814 events, and from the DDME1 set.
radius at the mediate mass region is a little larger when compared with the 14 observed neutron stars. The output EOSs of the DNN from scheme 1 provide smaller radii, which coincide with the distribution of observables. This shows that the final results of present framework is independent of the generating scheme for the training data.

In a binary neutron star merger, one neutron star will be deformed by the external gravitational field of another star. The magnitude of deformation is denoted as the tidal deformability, which is dependent on the EOS of the neutron star and can be extracted from the gravitational wave provided by the binary neutron star. In the GW170817 event, the dimensionless tidal deformability at \(1.4M_{\odot}\) was inferred as \(\Lambda_{1.4}=190^{+390}_{-120}\)(Abbott et al., 2018). In Fig. 8, the dimensionless tidal deformabilities as functions of neutron star masses from schemes 1 and 2, with 68% and 95% confidence levels, are plotted and compared to the constraint from the GW170817 event and the results from the DDME1 set. The \(\Lambda\) decreases with the neutron star mass since it is proportional to \(R^{5}/M^{5}\) of the neutron star. Therefore, the \(\Lambda\) from the DDME1 is relatively larger. The \(\Lambda_{1.4}\) from the reported machine learning framework completely satisfies the measurements from the gravitational wave detection.

Table 2 lists the properties of neutron stars fitted by the DNN with nonparametric training data: namely, the maximum masses of neutrons stars, the corresponding radii, the radii at \(1.4M_{\odot}\) and \(2.08M_{\odot}\), and the dimensionless tidal deformability at \(1.4M_{\odot}\) with 68% and 95% confidence levels in schemes 1 and 2. These variables were compared to the results from the DDME1 parameter set. Both of these two schemes can generate the massive neutron star with a mass close to \(2.55M_{\odot}\). The radius of the \(1.4M_{\odot}\) neutron star is around 12.30 km, which is consistent with the value extracted from the GW170817 of \(R_{1.4}=11.9\pm 1.4\) km (Abbott et al., 2019). The radius of \(2.08M_{\odot}\) neutron star is fitted around 12.0 km now. The radius and mass of PSR J0740+6620 were analyzed as \(12.39^{+1.30}_{-0.98}\) km and \(2.072^{+0.067}_{-0.066}M_{\odot}\), from

Figure 8: \(\Lambda\)-\(M\) relation, generated by the fitted EOSs and compared to that from DDME1 and the values extracted from GW170817 events.

Figure 7: The mass-radius relation of neutron star from the nonparametric machine learning method, the observation distributions from 14 neutron stars, the masses of massive neutron stars, and the radii constraints from the NICER.
NICER, by Riley et al. (Riley et al., 2021). The results from the two schemes are similar, with differences are less than 2%.

It can be found that present fits about the properties of the neutron are comparable with those generated by model-informed priors in the works of Landry and Essick (Landry and Essick, 2019), while they are much more constrained than the ones from model-agnostic prior. It is because our training data is just prepared to reproduce the theoretical EOSs, while the possibility that the EOS might be quite different from current theoretical fits was considered in model-agnostic prior.

Finally, the \(M\)-\(R\) relations from the two schemes to generate the training set, were compared and given in Fig. 9. Their behaviors are quite similar. The only difference is that the radii of the neutron stars and the uncertainties from scheme 2 are a little larger than those of scheme 1 because of the influence of the DDME1 set. This demonstrates that the fitted EOSs in the present framework is strongly independent of the choice of initial training data values using the GPR method.

## 4 Summaries and perspectives

A nonparametric methodology has been proposed to infer the EOSs of neutron star matter from recent observations of the neutron stars. A DNN was designed to map the mass-radius observables to the energy-pressure relation of dense matter. The GPR method was applied to construct the EOSs, and this method was completely independent of any apparent function form.

To generate the training data set, two schemes of the example data were adopted to provide the initial EOS. The mean values and variances of EOSs from nine successful relativistic mean-field model parameter sets were considered in the first scheme; whereas in the second, the mean value was chosen from the DDME1 set and the derivation was fixed as 0.3. A 5-million training data set was constructed by including the uncertainties in the mass and radius of neutron stars. Furthermore, in the training set, the constraints of the massive neutron star and the mass-radius simultaneous measurements were also taken into account in the training set.

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline  & C. L. & \(M_{max}[M_{\odot}]\) & \(R_{max}\) [km] & \(R_{1.4}\) [km] & \(R_{2.08}\) [km] & \(\Lambda_{1.4}\) \\ \hline DDME1 & & 2.45 & 11.83 & 12.99 & 12.98 & 692 \\ \hline \multirow{4}{*}{scheme 1} & 68\% & \(2.38^{+0.07}_{-0.07}\) & \(11.07^{+0.16}_{-0.17}\) & \(12.31^{+0.15}_{-0.16}\) & \(11.95^{+0.23}_{-0.23}\) & \(459^{+37}_{-46}\) \\ \cline{2-7}  & 95\% & \(2.38^{+0.15}_{-0.13}\) & \(11.07^{+0.34}_{-0.32}\) & \(12.31^{+0.29}_{-0.31}\) & \(11.95^{+0.44}_{-0.47}\) & \(459^{+82}_{-81}\) \\ \hline \multirow{4}{*}{scheme 2} & 68\% & \(2.41^{+0.08}_{-0.07}\) & \(11.15^{+0.21}_{-0.20}\) & \(12.30^{+0.17}_{-0.19}\) & \(12.03^{+0.27}_{-0.27}\) & \(448^{+55}_{-43}\) \\ \cline{2-7}  & 95\% & \(2.41^{+0.15}_{-0.14}\) & \(11.15^{+0.41}_{-0.39}\) & \(12.30^{+0.35}_{-0.37}\) & \(12.03^{+0.53}_{-0.54}\) & \(448^{+110}_{-86}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: The maximum masses of neutrons star, the corresponding radii at \(1.4M_{\odot}\) and \(2.08M_{\odot}\), and the dimensionless tidal deformability at \(1.4M_{\odot}\) from the nonparametric EOS models with 68% and 95% confidence levels in scheme 1 and 2 and compared to those from DDME1.

Figure 9: The \(M\)-\(R\) relation comparisons between two schemes with 95% confidence interval and the constraints from the massive neutron star and NICER.
One hundred independent NN models were generated with different training data sets, producing one hundred EOSs of a neutron star. These were analyzed with the standard statistical method and EOSs with the 68% and 95% confidence levels were obtained. They were softer when compared with the join constraints from the GW170817 and GW190814 events. The mass-radius relations from our fitted EOSs fully satisfy the present various astronomical observations of neutron stars. The dimensionless tidal deformability at \(1.4M_{\odot}\) was also consistent with the data extracted from the GW170817. Finally, concerning the creation of training data, the results from both schemes were almost identical. This shows that the present fitted EOSs are strongly independent of the initial set of training data set.

Our nonparametric NN framework can be naturally extended to other supervised learning fields to avoid the limitations of specific function forms. In the future, the original data on the gravitational wave from the binary neutron star will be included in the input layer to simulate the observations more realistically. The hadron-quark phase transition was excluded in the present training data set, and this too will be considered in future work.

## 5 Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (Grant Nos. 11775119 and 12175109), and the Natural Science Foundation of Tianjin (Grant No: 19JCYBJC30800). We are grateful to the referee for his constructive comments and suggestions.
