A Deep Learning Approach to Extracting Nuclear Matter Properties from Neutron Star Observations

Understanding the equation of state of dense QCD matter remains a major challenge in both nuclear physics and astrophysics. Neutron star observations from electromagnetic and gravitational wave spectra provide critical insights into the behavior of dense neutron-rich matter. The next generation of telescopes and gravitational wave observatories will offer even more detailed observations of neutron stars. Utilizing deep learning techniques to map neutron star mass and radius observations to the equation of state allows for its accurate and reliable determination. This work demonstrates the feasibility of using deep learning to extract the equation of state directly from neutron star observational data, and to also obtain related nuclear matter properties such as the slope, curvature, and skewness of the nuclear symmetry energy at saturation density. Most importantly, we show that this deep learning approach is able to reconstruct \textit{realistic} equations of state, and deduce \textit{realistic} nuclear matter properties. This highlights the potential of artificial neural networks in providing a reliable and efficient means to extract crucial information about the equation of state and related properties of dense neutron-rich matter in the era of multi-messenger astrophysics.


I. INTRODUCTION
The quest to determine the equation of state (EOS) of dense neutron-rich matter is a paramount challenge facing modern physics and astrophysics, representing one of the most pressing and critical unanswered questions [1][2][3].The EOS has significant implications for a broad range of phenomena, including heavy ion collision dynamics, binary neutron star mergers, supernovae, and gravitational waves.Both nuclear physics (see, for example, Refs.[4][5][6][7][8][9][10][11][12][13]) and astrophysics (see, for example, Refs. ) communities have made it a priority to investigate this fundamental problem and have established a wide range of research facilities, including telescopes, observatories, and gravitationalwave detectors, in order to advance our understanding of the EOS [1,2].
The nucleonic component of the EOS of cold neutron star matter can be expressed in terms of the energy per nucleon E(ρ, δ) as [36] where E SN M (ρ) is the energy per nucleon of symmetric nuclear matter (SNM), E sym (ρ) is the nuclear symmetry energy, and δ = (ρ n − ρ p )/ρ is the isospin asymmetry.In the above equation, ρ n , ρ p , and ρ represent the neutron, proton, and total density, respectively.Currently, the EOS of cold nuclear matter under extreme conditions is still uncertain and controversial, especially at supra-saturation densities, mainly because of the unknown high-density behavior of the nuclear symmetry energy E sym (ρ) [4,5].
In order to obtain the equation of state (EOS) of nuclear matter from first principles, one must solve for quantum chromodynamics (QCD), the fundamental theory of strong interactions.However, current model-independent results are limited to a narrow density range.At low densities around 1 − 2ρ 0 (where ρ 0 = 0.16 fm −1 is the saturation density of symmetric nuclear matter), ab initio methods can be combined with nuclear interactions derived from Chiral Effective Theory (χEFT) with controlled uncertainty estimates [37][38][39][40][41][42][43][44][45][46].For densities at ρ > ∼ 50ρ 0 , perturbative QCD calculations provide reliable results [47][48][49][50][51][52][53].However, for intermediate densities between 2 − 10ρ 0 , no reliable QCD predictions exist [54].To determine the EOS in this region, non-perturbative methods such as Monte Carlo simulation of QCD on a lattice (lattice QCD) must be developed, which faces significant challenges such as the sign problem in finite-density systems [55].Consequently, construction of the EOS at intermediate densities still relies on phenomenological approaches using many-body methods and effective interactions such as the relativistic mean field (RMF) theory and density functionals based on Skyrme, Gogny, or Similarity Renormalization Group (SRG) evolved interactions.
In recent years, there has been significant progress in determining the EOS at high densities from both nuclear laboratory experiments and multi-messenger astrophysical (MMA) observations of neutron stars (NSs).The experimental data from heavy-ion reactions collected from intermediate to relativistic energies, specifically related to nucleon collective flow and kaon production, has already significantly constrained the EOS of symmetric nuclear matter up to around 4.5 ρ 0 [6].The cooperation between the nuclear physics and astrophysics communities has resulted in substantial advancements in constraining the symmetry energy around and below the saturation density of nuclear matter using a combination of terrestrial nuclear experiments and astrophysical observations [5,[9][10][11][56][57][58][59].However, the density dependence of the nuclear symmetry energy E sym (ρ) at supra-saturation densities and the possible hadron-toquark phase transition remain the most uncertain aspects of the high-density EOS [4,5,21,22,24].The presence of new particles, such as hyperons and resonances, is also highly dependent on the high-density trend of E sym (ρ) [60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75].
The recent MMA observations of NSs have offered a unique opportunity to explore the high-density EOS.This represents an alternative way to independently extract the EOS by means of statistical approaches (as highlighted in Refs.[35,[76][77][78][79][80][81]).These observations encompass a wide range of methods, such as the Shapiro delay measurements of massive ∼2M pulsars [82][83][84], the radius measurement of quiescent low-mass X-ray binaries and thermonuclear bursters [76-78, 85, 86], X-ray timing measurements from the NICER mission [87][88][89], and the detection and inference of gravitational waves from compact binary mergers involving NSs [90][91][92] (by the LIGO/VIRGO/KAGRA [93][94][95] collaboration).Common observables for NSs include mass M , radius R, moment of inertia I, quadrupole moment Q, dimensionless tidal deformability Λ (and its derivatives such as Love number k 2 and tidal deformability λ), and compactness M/R.The NICER mission, for instance, targets the compactness M/R of neutron stars by measuring the gravitational lensing effect of the thermal emission from the star's surface.Meanwhile, gravitational-wave (GW) observations of binary neutron star (BNS) and neutron star-black hole (NSBH) mergers provide information about the tidal disruption of the star in the presence of its companion, which is quantified through the tidal deformability parameter λ.
There exist various statistical approaches to determine the most likely EOS from neutron star observational data.Of these, the use of Bayesian inference is widespread [35,[76][77][78][79][80].Gaussian processes also provide a non-parametric representation of the EOS [81].However, the uncertainty in Bayesian analyses raises questions regarding the true nature of dense matter EOS [54].In light of this, alternative model-independent methods are being sought.Deep neural networks (DNNs) [96,97] have garnered attention in the research community, where deep learning (DL) algorithms have displayed exceptional proficiency in tasks such as image recognition [98] and natural language processing [99].Furthermore, these techniques have been applied to various physics and astrophysics domains, including the analysis of GW data for detection [100][101][102][103][104][105][106][107], parameter estimation [108,109], and denoising [110].In previous works we employed Convolutional Neural Network (CNN) [111] algorithms to detect and infer GW signals from BNS [112,113] and, very recently, from NSBH [114] mergers.Additionally, the use of DNNs as a tool to extract the dense matter EOS from neutron star observations has also been explored in a growing number of studies [54,[115][116][117][118][119][120].
In a recent investigation [121], we presented an innovative approach to determine the nuclear symmetry energy, E sym (ρ), by utilizing DL techniques in conjunction with astronomical observations of neutron stars.Our results demonstrate that deep neural networks have the capacity to accurately extract E sym (ρ) from a set of M − R or M − Λ NS observations.This approach offers a promising avenue for exploring the high-density behavior of E sym (ρ), which remains a challenging task in nuclear physics.
In this paper, we extend our DL approach for determining the EOS of dense matter and associated nuclear properties using mass-radius M (R) measurements of neutron stars.In particular, we pay special attention to deducing the slope, curvature, and skewness of the nuclear symmetry energy, in addition to the EOS.Our results demonstrate that DL algorithms can accurately and reliably extract the NS EOS and nuclear matter properties from observational data.Moreover, we find that our DL approach can successfully reconstruct realistic EOSs and nuclear matter properties, which brings us one step closer to revealing the true nature of the EOS of dense, neutron-rich matter.
In the present work, we have structured our discussion as follows.In the first section, we have provided a brief introduction.In Section II, we have presented the main aspects of our formalism.This encompasses a comprehensive overview of the essential characteristics of the EOS employed in our analysis, along with the details of our DL algorithms, such as data generation, neural network architectures, and training methodologies.Subsequently, in Section III, we have put forth our results and their implications.Finally, in Section IV, we have summarized our findings and provided future research directions.

II. FORMALISM
In this section, we present the methodologies utilized in our study.First, we provide a comprehensive overview of the key characteristics and specifications of the EOS used.Subsequently, we briefly outline the procedure for solving the static NS structure equations.In Section II C, we discuss the DL approach adopted in mapping the NS mass-radius M (R) observations to the EOS, as well as the procedure for mapping the reconstructed EOS to selected nuclear matter properties.

A. Equation of State
The equation of state plays a critical role in determining the properties of neutron stars, such as mass M and radius R. To determine the nuclear matter EOS, two main theoretical approaches are commonly used: phenomenological and microscopic methods.Phenomenological approaches rely on effective interactions to describe the ground state of finite nuclei.These methods, including those based on Skyrme interactions [122,123] and relativistic mean-field (RMF) models [124], have been widely used in the study of low-density nuclear systems.However, they are not well-suited for high isospin asymmetry systems and at large densities, where experimental data are unavailable to constrain such interactions, predictions based on these methods can be far from realistic behavior [125].On the other hand, microscopic approaches use realistic two-body and three-body nucleon forces to describe the behavior of nucleons.These interactions can be based on meson-exchange theory [126,127] or more recent χEFT [42,[128][129][130].Microscopic many-body methods, such as the Brueckner-Hartree-Fock (BHF) approach [131], the Dirac-Brueckner-Hartree-Fock (DBHF) theory [132,133], the variational approach [134], the Quantum Monte Carlo technique and its derivatives [135,136], the self-consistent Green's function technique [137], χEFT [45], and the V low;k approach [138] are based on these interactions.The major challenge for these methods is the treatment of the short-range repulsive core of the nucleon-nucleon interaction, which distinguishes the different techniques from each other.
The nucleonic component of the EOS can be described by two quantities: the binding energy of symmetric nuclear matter, E SN M (ρ), and the symmetry energy, E sym (ρ) (see Eq. ( 1)).These two quantities can be expanded as Taylor series around ρ 0 as given by: where x ≡ (ρ − ρ 0 )/3ρ 0 .The coefficients of these expansions can be related to various physical properties of nuclear matter and can be experimentally constrained.They have the following meanings [139]:  [140].Several of these parameters have rather moderate uncertainty.For example, the binding energy E 0 is estimated to be −15.9± 0.4 MeV, while the magnitude of the symmetry energy S 0 is 31.7 ± 3.2 MeV.However, many of them still have significant uncertainty, such as the curvature of the symmetry energy K sym , which could range from −400 MeV to 100 MeV, and the higher order coefficients, J 0 and J sym , with even wider uncertainty ranges.
Although the Taylor expansions given by Equations ( 2) and (3) are known to diverge at higher densities [141], these expressions can also be viewed as parameterizations with free parameters [140].This duality means that, for systems with low isospin asymmetries, the Taylor expansions are valid near saturation density, while for highly neutron-rich systems at supra-saturation densities, Equations ( 2) and (3) should be treated as parameterizations [140].For further information on the relationship between the Taylor expansions and the parameterizations, we refer the reader to Ref. [140].These expressions are frequently utilized in modeling the NS EOS, and they have been applied, for instance, in solving the inverse structure problem of NSs and constraining high-density symmetry energy through NS observational data [140,142].In addition, the NS EOS metamodel has been used in Bayesian analyses to determine the most likely values of high-density EOS parameters through inference from NS data [143].Compared to the widely used piecewise polytropes, these parameterizations have the advantage of including isospin dependence and composition information throughout the density range, while still allowing for modeling a wide range of EOSs from various many-body approaches.This feature is particularly important for deducing the high-density E sym (ρ), as the parameterizations separate clearly the contribution of the symmetry energy to the EOS.For example, these parameterizations were instrumental in our previous work [121] for extracting the nuclear symmetry energy directly from NS observational data via deep neural networks.
In this study, we utilize an EOS metamodel to facilitate the extraction of the EOS and selected nuclear matter properties from NS observational data.By varying the parameters of the EOS, we can generate numerous EOSs and the corresponding sequences of mass and radius (M − R) by solving the NS structure equations.The matter in the core of the neutron star is modeled as a mixture of protons, neutrons, electrons, and muons in beta-equilibrium (referred to as the npeµ-model).We use the expressions for E SN M (ρ) and E sym (ρ) from Equations ( 2) and (3) to calculate E(ρ, δ) through Equation ( 1).The pressure of the neutron star matter in β-equilibrium can then be determined from the energy density, , where M N is the average nucleon mass and ε l (ρ, δ) is the lepton energy density.Further details on calculating ε l (ρ, δ) can be found in, e.g., Ref. [144].
When the density of the neutron star matter falls below approximately 0.07 f m −3 , the core EOS is complemented by a crustal EOS, which is more appropriate for lower density regions.For the inner crust, we use the EOS provided by Pethick et al. [145] and for the outer crust, the EOS by Haensel and Pichon [146].
In our analysis, we use Equations ( 2) and (3) as parameterizations together with the parabolic approximation of the nucleonic EOS given by Equation ( 1).The values of E 0 and S 0 are fixed at their most probable current values, which have been obtained through a combination of nuclear laboratory experiments and theoretical calculations.We subsequently vary the rest of the parameters, K 0 , J 0 , L, K sym , and J sym to generate many samples of the EOS.With the expanded parameter space, compared to the one considered in our previous work [121], we are able to model a wider class of EOSs as predicted by various many-body approaches, and models of the nuclear interaction.The effect of varying the individual parameters is shown in Figure 1.While, in principle, these parameters are absolutely free, the asymptotic boundary conditions of the EOS near ρ 0 and δ = 0 provide some prior knowledge of the ranges of these parameters.Their ranges are further restricted by imposing the requirement that the EOSs must satisfy causality and the microscopic stability condition, and the resultant NS models can support a maximal mass of at least 2.14 M of the heaviest pulsar observed so far [84].The ranges of E SN M and E sym (ρ) satisfying all constraints are shown in Figure 2.

B. Structure Equations of Static Neutron Stars
In this section, we briefly revisit the procedure for calculating the mass M and radius R of static neutron stars.For a spherically symmetric relativistic star, Einstein's field equations can be simplified to the Tolman-Oppenheimer-Volkoff (TOV) equation [147], as follows: where the mass within a sphere of radius r is determined by To solve the above equations, one needs to supplement them with the EOS in the form P (ε).Starting with the initial conditions m(r = 0) = 0 and ε c = ε(r = 0) at the NS center (r = 0), integration of Equations ( 5) and ( 6) is carried out until the pressure P reaches zero, marking the edge of the star.Some care should be taken at r = 0 since the above equations are singular at the center.The point r = R where P vanishes determines the NS radius and For a given EOS, there is a unique relationship between the stellar mass and the central density ε c .Thus, for a particular EOS, there is a unique sequence of NSs parameterized by the central density (or equivalently the central pressure P c = P (0)).In Figure 3 we show the range of possible EOSs (P (ρ)) satisfying all constraints (left window), and the resultant M − R NS sequences (right window).

C. Artificial Neural Networks
In this section, we briefly discuss the basic setup, structure, and workflow associated with implementing DNNs for our specific application.For more extensive discussions, the reader is referred to a number of machine learning articles [96,148] and textbooks [97,149].
We apply a combination of two DNNs with similar architectures to first extract the EOS of dense neutron-rich matter from a set of mass-radius NS measurements, and then deduce selected nuclear matter properties from the β-equilibrium NS EOS.We refer to the first neural network as EOS DNN (equation of state deep neural network), and the second one NuPRO DNN (nuclear matter properties deep neural network).The procedure of using these DNNs to extract the EOS and selected nuclear matter properties is illustrated schematically in Figure 4.

EOS Network (EOS DNN)
To extract the EOS, in this analysis, we apply a supervised DL approach and formulate a regression problem, where the input to the DNN consist of M (R) sequences (sets of points representing pairs of NS mass-radius measurements), while the output consist of EOS (P (ε)) estimates.Accordingly, the datasets used for the training, validation, and testing of EOS DNN consist of M (R) sequences and P (ε) samples.We use the EOS metamodel discussed in Section II A, and vary the parameters in Equations ( 2) and (3) to generate many samples of the EOS, and subsequently by solving the NS structure equations, corresponding M − R sequences.Specifically, we set E 0 = 15.9MeV and S 0 = 31.7 MeV, and vary the rest of the parameters by randomly sampling their values from their respective ranges: K 0 = 240 ± 20 MeV, −300 ≤ J 0 ≤ 400 MeV, L = 58.7 ± 28.1 MeV, −400 ≤ K sym ≤ 100 MeV, −200 ≤ J sym ≤ 800 MeV.Recently, the latest results of the PREX collaboration suggested a rather high value of L with an upper limit at 143 MeV [150].Examining the effect of higher L values is left to following works.The resultant EOSs P (ε) are checked regarding whether they satisfy (i) the microscopic stability condition, i.e., dP dε ≥ 0, and (ii) the causality condition, i.e., the speed of sound c s ≡ dP dε ≥ c.In addition, the resultant NS models must be able to sustain a maximal mass of at least 2.14 M [84].These constraints restrict the values of K 0 , J 0 , L, K sym and J sym , and the final EOS samples.To simulate NS observational data, from a given genuine M − R sequence, we randomly choose 50 points in the range of 1M to M max supported by the given EOS [167].Then, each input sample is an array of dimension 2 × 50 consisting  of 50 pairs of (M , R) values.The values of M and R are scaled by dividing them by 3 and 20 respectively to ensure that the input data are in the (0, 1) range.Similarly, each output sample is an array of dimension 2 × 50 consisting of 50 pairs of estimated (P , ε) values, representing the EOS in the density range from ∼ 0.4ρ 0 to 5ρ 0 .In this respect, the DNN maps an input M (R) sequence to an output EOS, P (ε).
In supervised learning, the data are divided into training, validation, and testing data sets.The training data set is used by the DNN to learn from, the validation data are used to verify whether the network is learning correctly, and the testing data are used to assess the performance of the trained model.Here, the training dataset consist of 120,000 independent M (R) sequences, representing the DNN inputs, and 120,00 matching EOS samples, P (ε), representing the DNN outputs.From each M (R) sequence we further draw 50 ensembles, each containing 50 randomly selected (M , R) pairs.In this way, each EOS sample in the training data set is represented by 50 different random ensembles drawn from the same genuine M (R) curve.The final training data set therefore consist of 6 × 10 6 samples.Similarly, the final validation data set consist 250,000 samples, where each of 5,000 independent output EOS samples is represented by 50 different random ensembles drawn from the same M (R) sequence.Finally, the testing data set consist of 5,000 unique input and output samples, not used in the training and validation of the DNN.
EOS DNN is a feedforward neural network with 5 hidden, dense, fully connected layers of dimension 500, and ReLU activation functions, followed by a dense linear layer of dimension 100.The first layer corresponds to the input to the neural network, which, in this case, is a 2 × 50 array containing the NS M and R values drawn from a given M (R) curve.At the end, there is a linear output layer of dimension 2 × 50 returning the estimated EOS, P (ε) (pairs of (P , ε) points).The network design was optimized by fine-tuning multiple hyper-parameters, which include here the number and type of network layers, the number of neurons in each layer, and the type of activation function.The optimal network architecture was determined through multiple experiments and tuning of the hyper-parameters.The architecture of the EOS DNN is summarized in Table I.
To build and train the neural network, we used the Python toolkit Keras (https://keras.io),which provides a high-level application programming interface (API) to the TensorFlow [151] (https://www.tensorflow.org)deep learning library.We applied the technique of stochastic gradient descent with an adaptive learning rate with the ADAM method [152] with the AMSgrad modification [153].To train the DNN, we used an initial learning rate of 0.001 and chose a batch size of 1000.During each training session, the number of epochs was limited to 5000, or until the validation error stopped decreasing.By employing checkpoint-callback, we selected the model achieving the lowest loss value on the validation dataset.The training of the DNN was performed on an NVIDIA Tesla V100 GPU and the size of the mini-batches was chosen automatically depending on the specifics of the GPU and data sets.We used the mean absolute error (MAE) as a cost (or loss) function: where n is the number of samples, y i is the "true" value and ŷi is the DNN model prediction.The NuPRO DNN architecture comprises an input layer with 50 dimensions, corresponding to the 50 equally spaced points of the EOS P (ρ).This is followed by five dense fully connected layers of varying dimensions, culminating in an output layer that returns the estimated nuclear matter parameters K0, J0, L, Ksym, and Jsym.The total number of trainable parameters in this DNN model is 118,555.Further information about this architecture can be found in the text.

Nuclear Matter Properties Network (NuPRO DNN)
Once the β-stable matter EOS becomes available, we can proceed with the extraction of chosen nuclear matter properties.In order to achieve this goal, we have trained another DNN, which we refer to as the NuPRO DNN.The input to the DNN consists of EOS samples represented as P (ρ), which are sets of 50 equally spaced points within the interval ρ = [0.08 − 0.8]f m −3 , where the input data was converted to decimal logarithm values.On the other hand, the output corresponds to estimations of selected nuclear matter properties, with respect to each input EOS sample.Our aim is to learn the mapping y(x) with x i = [P β (ρ 1 ), P β (ρ 2 ), ..., P β (ρ 50 )] and y i = [K 0 , J 0 , L, K sym , J sym ] being the corresponding set of parameters.
In our work, we utilized a training dataset comprising 120,000 samples of the equation of state P (ρ), each with an accompanying set of parameters (K 0 , J 0 , L, K sym , and J sym ) matching the specific EOS realization.Additionally, we constructed separate validation and testing datasets, each consisting of 5,000 data samples.The architecture of the NuPRO DNN model that we have implemented involves a feedforward structure with five hidden layers, each being dense with dimensions 200, 200, 200, 100, and 50, respectively.We have utilized the ReLU activation functions for each of these layers.The neural network architecture of the EOS DNN was optimized through an iterative process involving multiple experiments and tuning of the hyper-parameters.The neural network's input layer has a dimension of 50 and corresponds to the 50 uniformly distributed data points representing the equation of state, P (ρ), while the output layer has a dimension of 5, which returns the estimated nuclear matter parameters.A summary of the network architecture is provided in Table II.
We used Keras and TensorFlow to develop and train our neural network.Similarly to before, we employed stochastic gradient descent with an adaptive learning rate through the ADAM method [152], which was further modified with the AMSgrad technique [153].For training the DNN, we selected a batch size of 1000 and initialized the learning rate at 0.001.Additionally, we set a limit of 5000 epochs for each training session or until the validation error reached a minimum.Using a checkpoint-callback, we selected the model with the lowest loss value on the validation dataset.The cost function used for this task is the mean-squared error (MSE), which is defined as the sum of the squared differences between the predicted values of the DNN model, ŷi , and the actual or "true" values, y i , divided by the number of samples, n: III. RESULTS

A. Extracting the EOS
We first examine the ability of the DNN to reconstruct the EOS, P (ε), from a set of mass and radius M − R measurements that may result from electromagnetic observations of neutron stars, such as those from the NICER mission, for instance.In particular, we apply the trained EOS DNN, described in the previous section, to a test dataset containing ∼5000 simulated M (R) sequences, and compare the corresponding estimated output EOS with the exact EOS for each sample.In Figure 5, we show results for five representative examples from the test dataset.It is seen that the EOS (broken colored lines) for each input M − R sequence matches almost exactly the "true" EOS (solid black lines) over the entire density range considered here.For the purpose of presentation, the EOS is shown in the form P (ρ).The results are very similar for the rest of the test data samples.Quantitatively, the mean absolute error over the whole test dataset is 0.5 MeV fm −3 with standard deviation of 1.3 MeV fm −3 .Choosing different ensembles of randomly selected points from the genuine M (R) curves does not alter appreciably the accuracy with which the EOS is estimated.
Realistic NS observations inevitably carry uncertainties, which result in corresponding uncertainties in the estimated EOS.To investigate the effect of the observational uncertainties on extracting the EOS via our DL approach, we prepared a test dataset assuming that the NS mass and radius measurements are subject to a measurement error.After randomly selecting N points (M i , R i ) from a genuine M (R) curve, with {i = 1, 2, ..., N ; N = 50}, we draw the actual M and R values from normal distributions, N (µ M , σ M ) and N (µ R , σ R ) respectively.Specifically, we examined the effect of smaller and larger observational errors by choosing σ M = 0.02 M and σ R = 0.1 km to simulate smaller uncertainties, and σ M = 0.1 M and σ R = 1 km to model larger uncertainties, where µ M = M i and µ R = R i .To quantify the effect of observational uncertainties for each sample in the test dataset, for each case we draw 100 ensembles from the respective normal distribution and calculate the mean absolute errors.In Figure 6, we illustrate the effect of measurement errors for a representative example from the test dataset.It is seen that for smaller observational uncertainties (σ M = 0.02 M , σ R = 0.1 km) the estimated EOS P (ρ) (broken green line) matches almost exactly the "true" EOS (solid black line).The greenish shaded band represents the corresponding mean absolute errors.For larger observational uncertainties (σ M = 0.1 M , σ R = 1 km), we see that at higher densities, of ∼ ρ/ρ 0 ≥ 4, the reconstructed EOS (broken red line) starts to moderately diverge from the ground-truth EOS, but it is still within the reconstruction errors represented by the reddish shaded band (corresponding to the mean absolute errors).The results follow a very similar trend for the rest of the test data samples, and suggest that the accuracy of the EOS estimation is mainly affected by the magnitude of the assumed observational errors.
We emphasize that the observational uncertainties are introduced in the analysis by directly subjecting only the test dataset to an assumed level of "noise", and the trained DNN model does not have prior knowledge of similar uncertainties in the training data.Nevertheless, these results show that the neural network is able to accurately reconstruct the EOS from moderately noisy NS M (R) observational data.The performance of the DNN can be further improved by introducing the measurement errors also in the training dataset.Since the systematic study of the measurement uncertainty effects is not the major focus of this work, detailed investigations, and also studying the effect of introducing "noise" in the training data, are left for following articles.

B. Application to realistic EOSs
To test further the performance of the trained EOS DNN model, we apply it to several realistic EOSs from the Com-pOSE repository [154] (https://compose.obspm.fr).CompOSE is an online tool that provides data tables containing state-of-the-art EOSs that can be readily used for various applications in astrophysics and nuclear physics.For the purpose of our analysis we chose several EOSs within the range of the parameter space of the DNN training dataset: APR [134], BL [155], QMC-RMF2 [156], QMC-RMF3, [156], SK255 [157], and SK272 [157].We generated the required input data following the procedure outlined in Section II C 1.For each EOS (Figure 7, left window), we calculated the M − R relation (Figure 7, right window), and then drew 10 5 random ensembles of 50 (M i , R i ) points to determine the mean and MAE of the reconstructed EOSs.The results of this test are shown in Figure 8, and demonstrate the ability of the EOS DNN to reconstruct realistic EOSs from M (R) data.In all frames, the solid blue lines denote the ground-truth EOS and the red dot characters represent the mean of the DNN predictions respectively.The error bars represent the MAE with which the EOS is reconstructed in each case due to the random drawing of the M (R) data, and in order to clearly separate the error contribution of this effect alone, they do not include the effect of assumed observational uncertainties.Among the realistic EOSs we have considered, it is seen that the predicted P (ρ) relation matches almost exactly the ground-truth values for the APR, BL, QMC-RMF2 and QMC-RMF3 models, while the the SK255 and SK272 EOSs are reconstructed less accurately, but still within the reconstruction errors.
Here we briefly recall the main features of the realistic EOSs used in our analysis.The APR EOS is calculated using variational approaches with the A18 + delta v + UIX* interaction [134].The BL EOS is obtained using realistic two-body and three-body nuclear interactions derived in the framework of χEFT and including the ∆(1232) isobar intermediate state [155].This EOS has been derived using the Brueckner-Bethe-Goldstone quantum manybody theory in the Brueckner-Hartree-Fock (BHF) approximation with the continuous choice for the auxiliary single particle potential.The QMC-RMF2 and QMC-RMF3 EOSs are computed using a relativistic mean-field (RMF) theory constrained by χEFT calculations of pure neutron matter (from 0.08 fm −3 to 0.32 fm −3 ) and by properties of isospin-symmetric nuclear matter around ρ 0 [156].The SK255 and SK272 EOSs are unified models by Gulminelli and Raduta [157] computed with the SK255 and SK272 effective interactions [158].The APR and BL EOSs are microscopic while the rest of the EOSs are based on phenomenological models.
These results clearly demonstrate that a DNN, trained on a relatively simple dataset generated with the EOS metamodel discussed in Section II A, is able to generalize the task of reconstructing the EOS and predict accurately realistic EOSs.

Performance on the Test Dataset
In the following analysis, we examine the effectiveness of the trained NuPRO DNN model in extracting particular nuclear matter properties from the EOS of β-equilibrium NS matter.After the model is trained and the optimal architecture is determined (as shown in Table II), we assess its final performance by evaluating it on a test dataset composed of 5,000 samples of P (ρ), and corresponding sets of selected nuclear matter properties matching each EOS sample.To evaluate the model's performance, we compute the standard deviation, σ i , of the residuals, i = Q DN N i − Q i , for each of the nuclear matter parameters: K 0 , J 0 , L, K sym , and J sym .Here Q i is one of the selected 5 nuclear matter properties.By examining the standard deviation, we can determine the degree of accuracy and precision of the model's predictions for each of these parameters.
In Figure 9, we present the results of our evaluation of the performance of the trained NuPRO DNN model in extracting selected nuclear matter properties from the EOS of β-equilibrium NS matter.We provide scatter plots of the distribution of the residuals for each of the nuclear matter parameters, along with the numerical values for the mean FIG.9: Residuals of the model for each of the selected nuclear matter parameters along with the numerical values for the mean, µ, and standard deviation, σ.As observed, for the lower-order parameters (K0 and L), the mean values of the residuals are less than 0.5 MeV (with |µ| ≈ 0.2 MeV for both cases).Additionally, it can be seen that the standard deviation is comparatively smaller for the lower-order parameters.This can be attributed to the fact that the range of possible values for the lower order parameters is smaller, resulting in better interpolation precision.On the other hand, the higher-order parameters exhibit larger values of σ, owing to the larger range of possible values.For further information, refer to the text.
µ and the standard deviation σ .These results clearly demonstrate that the trained NuPRO DNN model achieved a high degree of accuracy in extracting the nuclear matter parameters.In particular, we observed that the lower-order terms were extracted with higher accuracy, which can be attributed to the smaller range of possible values compared with the higher-order terms.The better interpolation precision of the model associated with the smaller range of possible values allowed for more accurate extraction of the lower-order terms.These results are also summarized in Table III.

Reconstructing Esym(ρ)
Next, we focus on evaluating the ability of the trained NuPRO DNN model to accurately extract nuclear symmetry energy parameters, namely L, K sym , and J sym , which are used to reconstruct E sym (ρ).The nuclear symmetry energy is a crucial yet uncertain component of the high-density equation of state, and it is imperative to explore whether DL  techniques could offer a viable means of deducing it from astrophysical observations of neutron stars.
In our previous work [121], we pioneered the use of DL methods for extracting the nuclear symmetry energy from a set of neutron star observations in the M − R or M − Λ planes.In our "proof-of-concept" study [121], our main focus was on the extraction of E sym (ρ), and thus we generated our datasets by holding all parameters in Equation ( 2), representing the energy of symmetric nuclear matter, constant and varying only the nuclear symmetry energy parameters L and K sym in Equation (3).We demonstrated that, under the given model assumptions, DNNs could extract E sym (ρ) effectively and accurately directly from astronomical observations of neutron stars.In our present investigation, we have advanced our deep learning methodology for extracting the nuclear symmetry energy E sym (ρ) by significantly enlarging the parameter space of our neural network training dataset.To achieve this, we have kept only the parameters E 0 and S 0 fixed at their most probable values while varying the other parameters, namely K 0 , J 0 , L, K sym , and J sym , in Equations ( 2) and (3), to generate numerous samples of P (ρ).As depicted in Figure 7, the augmented parameter space of the neural network training datasets also permits the modeling of predictions of modern realistic equations of state, which satisfy the constraints from recent mass-radius observations of neutron stars.
In Figure 10 we show results for five selected instances from the test dataset.It can be seen that the reconstructed nuclear symmetry energy (shown as broken colored lines) for each input P (ρ) sample agrees nearly perfectly with the true E sym (ρ) (depicted by solid black lines).Similar outcomes are obtained for the remaining test data samples.We emphasize that the reconstructed nuclear symmetry energy is deduced by substituting the estimated values of the nuclear symmetry energy parameters, namely L, K sym , and J sym , predicted by the NuPRO DNN, into Equation (3).Moreover, we assume that the β-equilibrium NS EOS, P (ρ), is already known with a certain level of precision, such as being extracted from mass-radius observations of neutron stars by using the EOS DNN trained model.

Model Uncertainty
Let us not forget that realistic observations of neutron stars unavoidably harbor uncertainties, which in turn give rise to uncertainties in the inferred EOS of β-stable NS matter, and consequently in the extracted nuclear matter parameters and the symmetry energy E sym (ρ).To investigate the impact of errors in the reconstruction of the EOS on the inferred nuclear matter parameters and symmetry energy, we have incorporated "noise" into the P (ρ) data samples that portray the EOS, and evaluated the nuclear matter parameters and E sym (ρ).We have conducted experiments by varying the level of noise and investigated the resulting effect on the accuracy of the extracted nuclear matter parameters and symmetry energy.As shown in Figure 11, we elucidate the effect of introducing a 20% uncertainty to the input P (ρ) data samples on the reconstructed E sym (ρ).In the left window, we show the exact EOS (represented by the solid blue line) and an EOS data sample containing 20% uncertainty (indicated by the red broken line).The reddish colored band denotes the uncertainty of P (ρ).In order to assess the uncertainty in determining the nuclear matter parameters and the symmetry energy, we generate 10 5 random sets of 50 equally spaced points in ρ, P (ρ i ), where i = 1, 2, ..., 50, lying within the uncertainty band.We subsequently compute the mean and standard   deviation for each of the nuclear matter parameters for each set.Thereafter, with every estimated set of nuclear matter parameters, we determine E sym (ρ) through Equation (3).The reconstructed symmetry energy is illustrated in the right window of Figure 11.The reddish colored band represents the mean absolute error (MAE) in deducing E sym (ρ), the solid line depicts the exact symmetry energy, and the red broken line indicates the mean symmetry energy.As anticipated, since the inferred symmetry energy is reconstructed via Equation (3), it closely follows the exact E sym (ρ) in a qualitative manner.Quantitatively, the estimated values begin to deviate moderately from the exact ones at approximately ρ ≥ 2ρ 0 , however, they remain within the range specified by the mean absolute errors of the model for the assumed uncertainty of the input EOS.The mean, standard deviation, and MAE for each of the nuclear matter parameters for the specific example shown in Figure 11 are presented in Table IV.The results are highly analogous for the remainder of the data samples from our test dataset.
It is important to note that the uncertainties presented in our analysis are solely introduced to the test dataset, and the trained DNN model does not possess any prior knowledge of uncertainties in the training data.Despite this, our findings demonstrate that the neural network is capable of accurately extracting the nuclear matter parameters and reconstructing E sym (ρ), even when faced with moderately noisy input EOS data.In order to further improve the performance of the DNN, it may be beneficial to introduce uncertainties to the training dataset as well.Although the impact of measurement uncertainties is not the primary focus of the current study, the potential effects of such uncertainties can be studied in detail in future works.This includes investigating the effects of introducing "noise" to the training data, as well as conducting systematic studies on the impacts of measurement uncertainties.These investigations may provide further insight into the behavior of the DNN and help to enhance its performance in future applications.Having established this, we now proceed to applying the model to a set of realistic EOSs, which were previously discussed in Section III B. However, before we discuss the results, it is important to highlight the complexity of the inference task and the limitations of our model assumptions.Firstly, the DNN model was trained on a dataset that assumes the EOS of nuclear matter depends on the matter isospin asymmetry via a quadratic dependence only, as specified in Equation (1).Secondly, for the hadronic component of the EOS, we used parameterizations given by Equations ( 2) for symmetric nuclear matter and (3) for the nuclear symmetry energy E sym (ρ) in the density range from approximately 0.04 f m −3 to 0.8 f m −3 .It is important to note that beyond the saturation density ρ 0 , these expressions should be regarded solely as parameterizations, and not as Taylor expansions.Thirdly, we made the assumption that neutron star matter is composed of nucleons, electrons, and muons in β-equilibrium.This assumption is made to simplify the problem, and it may not accurately represent the composition of matter in neutron stars, which may include other exotic particles, such as hyperons, or quark matter.
For the purpose of our analysis, in Figure 12, we present the residuals of the nuclear matter parameters K 0 , L, and K sym of the trained NuPRO DNN model for each realistic EOS considered in this study.These residuals correspond to the differences between the predicted values of the parameters and their true values obtained from the CompOSE repository.We note that K sym values were not available for the APR and BL EOSs, and hence, we do not show the residuals for these models.The standard deviations, σ i , of the residuals were also calculated to assess the uncertainty associated with the estimation of the nuclear matter parameters from a real β-equilibrium NS EOS.The standard deviations for K 0 , L, and K sym are 30.18MeV, 11.22 MeV, and 19.09 MeV, respectively.These values are smaller than the reported uncertainties in the literature [159].
It is important to emphasize that our model assumptions and limitations should be taken into account when interpreting these results.Therefore, the applicability of our results to other types of matter, such as hyperonic matter or quark matter, remains an open question.Nonetheless, our results demonstrate the potential of using DNN models to extract nuclear matter parameters from astrophysical observations of neutron stars.
Precise measurements of the masses and radii of a sufficient number of neutron stars would ultimately allow for the accurate determination of the EOS of β-stable matter through converting the M (R) curve, via various methods, to the underlying EOS [119].However, extracting the nuclear matter properties from the β-equilibrium EOS poses another challenge itself since the interior composition of a neutron star is unknown, and even the determination of the proton fraction is highly challenging [119].For instance, in a Bayesian approach presented in Ref. [160], the authors were unable to deduce the nuclear matter properties from the β-stable matter EOS.Similarly, the authors of Ref. [161] demonstrated the existence of multiple solutions for the determination of the NS interior composition from the β-stable matter EOS, owing to the high level of degeneracy.Furthermore, the determination of the nuclear symmetry energy from the β-equilibrium EOS requires an accurate knowledge of the EOS of symmetric nuclear matter [162], which is necessary for determining the proton fraction in the NS interior.These considerations underscore the importance and potential of the DL methods presented in this work as they provide a model-independent avenue to deducing the EOS of β-stable matter, and in turn, the nuclear matter parameters and E sym (ρ).

IV. SUMMARY AND OUTLOOK
In this study, we have demonstrated the feasibility of using a DL approach to directly extract the EOS of dense neutron-rich matter from observational data of neutron stars.Through analysis of simulated mass and radius measurements of neutron stars, we have shown that deep neural networks can accurately extract the EOS of β-stable NS matter.Furthermore, we have illustrated the ability of a trained DNN model to deduce selected nuclear matter properties, including the E sym (ρ).Most importantly, we have demonstrated that our DL approach can accurately extract realistic EOSs and nuclear matter properties from NS observational data.These results represent an important step towards the ultimate goal of determining the EOS of dense nuclear matter, and highlight the potential of DL-based techniques in the era of multi-messenger astrophysics, where a growing volume of NS observational data is rapidly becoming available.
In the near future, we plan to systematically examine the uncertainties associated with the NS observational data and the DNN model, and their impact on the model's performance.By understanding the effect of these uncertainties, we aim to explore potential approaches to further enhance the model's reliability and performance.In particular, in order to apply our approach practically, it is essential to consider the empirical errors and uncertainties and incorporate them consistently into the formalism.A possible strategy to achieve this is to recast the regression problem of extracting the EOS and nuclear matter properties into a probabilistic framework.Specifically, in subsequent research, we intend to use Bayesian neural networks to perform the inference task.In this paradigm, instead of obtaining deterministic values, the weights of the network are characterized by probability distributions by placing a prior over the network weights [163].In future studies, we also plan to apply our DL approach to real observational data of neutron stars, which would enable the extraction of a model-independent EOS, nuclear matter properties, and symmetry energy.Finally, we also plan to investigate likelihood-free inference methods using normalizing flows [164].These techniques are able to model complex posteriors by applying nonlinear transformations to a simple posterior shape, such as a multivariate Gaussian, without evaluating the likelihood directly.This approach has already generated considerable interest in the scientific community, and it has been successfully applied in multiple research domains.For instance, a recent study [165] applied a likelihood-free inference method using normalizing flows to rapidly estimate the parameters of eight GW binary black hole (BBH) events in the first LIGO Gravitational Wave Transient Catalog, GWTC-1 [166].With the next-generation of space telescopes and GW detectors, which will be sensitive enough to detect and observe compact binary collisions and neutron stars throughout the history of the universe, identifying over a million events per year, including thousands of BNS and NSBH detections per year, it will be crucial to process the incoming observational data quickly and accurately.In this context, it is important to emphasize that traditional Bayesian inference methods are not scalable to the study of thousands of BNS and NSBH events per year, and modern normalizing flow models, and similar approaches, could play a critical role in accurately and promptly extracting important NS parameters.
In the end, with the increasing number of observed events involving neutron stars, these contemporary data-driven techniques will enable us to rapidly process the growing volume of neutron star observational data and accurately determine the equation of state of dense nuclear matter and the nuclear symmetry energy.
Data Availability: Codes and data from this analysis are available upon request from the author.

FIG. 1 :
FIG. 1: (Upper left) Energy per particle of SNM as a function of the reduced density ρ/ρ0 for various values of K0, with E0 = 15.9MeV and J0 = 0 MeV.(Upper middle) Same as the upper left window but for various values of J0, with E0 = 15.9MeV and K0 = 240 MeV.(Upper right) Symmetry energy Esym as a function of ρ/ρ0 for various values of L, with S0 = 31.7 MeV, Ksym = 0 MeV, and Jsym = 0 MeV.(Lower left) Same as the upper right window but for various values of Ksym, with S0 = 31.7 MeV, L = 58.7 MeV and Jsym = 0 MeV.(Lower right) Same as the previous two windows but for various values of Jsym, with S0 = 31.7 MeV, L = 58.7 MeV and Ksym = 0 MeV.See text for details.

FIG. 2 :
FIG. 2: Range of the energy of symmetric nuclear matter ESNM (left window) and the nuclear symmetry energy Esym (right window).The ESNM and Esym are plotted as functions of the reduced density ρ/ρ0.

FIG. 3 :
FIG. 3: (Left window) Range of the EOS incorporating all constraints: Total pressure P as a function of the reduced density ρ/ρ0.(Right window) Range of mass-radius relation: Corresponding M − R sequences of the NS models computed with the EOSs considered in this study.The mass ranges of the three heaviest pulsars known at present [82-84] are indicated in the right window.

FIG. 5 :
FIG. 5: Example input M (R) sequences (left window) and corresponding estimated P (ρ) (right window).The input samples consist of 50 randomly selected points, denoted by the "o" characters, from the genuine M (R) curves, denoted by the solid lines, in the range of 1-Mmax M .The output data samples consist of 50 P (ρ) points in the range of ∼0.4-5 ρ0.Broken colored lines in the right window denote the estimated EOS and the solid lines represent the "true" EOS.Same curve colors in both windows denote pairs of input M (R) sequences and corresponding output EOS samples.

FIG. 6 :
FIG. 6: Example EOS (P (ρ)) predictions of the trained DNN model illustrating the effect of observational uncertainties included in the test dataset (simulated M (R) measurements with errors).The broken green and red lines represent the mean of the extracted EOS for smaller and larger measurement errors respectively, while the greenish and reddish shaded bands denote the corresponding mean absolute errors.The magnitude of the errors introduced in the input M (R) measurements is controlled by the values of σM and σR defining the normal distributions from which M and R are drawn.Specifically, to model smaller uncertainties, we choose σM = 0.02 M and σR = 0.1 km.Similarly, to model larger uncertainties, we choose σM = 0.1 M and σR = 1 km.The solid black line denotes the ground-truth EOS.See text for details.

FIG. 10 :
FIG. 10: Reconstructed Esym(ρ) from the β-equilibrium NS EOS, P (ρ), for several representative samples from our training dataset.The black solid curves represent the ground-truth symmetry energy and the broken colored lines denote the DNN predictions respectively.The predicted symmetry energy is obtained through Equation (3) with the parameters L, Ksym and Jsym estimated by the trained NuPRO EOS model.

FIG. 11 :
FIG. 11: (Left window) EOS data sample, P (ρ), with added 20% uncertainty.The ground-truth equation of state, P (ρ), is represented by the blue solid line, while the red dashed line shows an instance of data with random noise added within the range of uncertainty.The uncertainty of the input equation of state is illustrated by the reddish band.(Right window) Estimated nuclear symmetry energy, Esym(ρ).The precise value of nuclear symmetry energy is indicated by the blue solid line, while the red dashed line represents the average value of the derived Esym(ρ), and the reddish colored band indicates the mean absolute error (MAE).The calculation of Esym(ρ) is done using Equation (3) with the nuclear matter parameters L, Ksym, and Jsym extracted through NuPRO DNN.Further information can be found in the text.

FIG. 12 :
FIG. 12: NuPRO DNN model residuals for K0 (left window), L (middle window), and Ksym (right window) for the EOSs considered in our analysis.Note that Ksym values are not available for the APR and BL EOSs and therefore residuals for these models are not shown in the figure.See text for details.

TABLE I :
The EOS DNN architecture comprises an input layer followed by 6 dense fully connected layers.The output layer returns the estimated equation of state, P (ε).The model has 1,102,600 trainable parameters.Further details can be found in the text.

TABLE IV :
Values for the exact, predicted, mean µi, standard deviation σi, and mean absolute error (MAE) of the nuclear matter parameters Qi=[K0, J0, L, Ksym, Jsym] for the example illustrated in Figure11.All values are given in MeV.Please see the text for further details.