Decoding Neutron Star Observations: Revealing Composition through
  Bayesian Neural Networks
We exploit the great potential offered by Bayesian Neural Networks (BNNs) to
directly decipher the internal composition of neutron stars (NSs) based on
their macroscopic properties. By analyzing a set of simulated observations,
namely NS radius and tidal deformability, we leverage BNNs as effective tools
for inferring the proton fraction and sound speed within NS interiors. To
achieve this, several BNNs models were developed upon a dataset of $\sim$ 25K
nuclear EoS within a relativistic mean-field framework, obtained through
Bayesian inference that adheres to minimal low-density constraints. Unlike
conventional neural networks, BNNs possess an exceptional quality: they provide
a prediction uncertainty measure. To simulate the inherent imperfections
present in real-world observations, we have generated four distinct training
and testing datasets that replicate specific observational uncertainties. Our
initial results demonstrate that BNNs successfully recover the composition with
reasonable levels of uncertainty. Furthermore, using mock data prepared with
the DD2, a different class of relativistic mean-field model utilized during
training, the BNN model effectively retrieves the proton fraction and speed of
sound for neutron star matter.
Valéria Carvalho, Márcio Ferreira, Tuhin Malik, Constança Providência
# Decoding Neutron Star Observations: Revealing Composition through Bayesian Neural Networks

###### Abstract

We exploit the great potential offered by Bayesian Neural Networks (BNNs) to directly decipher the internal composition of neutron stars (NSs) based on their macroscopic properties. By analyzing a set of simulated observations, namely NS radius and tidal deformability, we leverage BNNs as effective tools for inferring the proton fraction and sound speed within NS interiors. To achieve this, several BNNs models were developed upon a dataset of \(\sim\) 25K nuclear EoS within a relativistic mean-field framework, obtained through Bayesian inference that adheres to minimal low-density constraints. Unlike conventional neural networks, BNNs possess an exceptional quality: they provide a prediction uncertainty measure. To simulate the inherent imperfections present in real-world observations, we have generated four distinct training and testing datasets that replicate specific observational uncertainties. Our initial results demonstrate that BNNs successfully recover the composition with reasonable levels of uncertainty. Furthermore, using mock data prepared with the DD2, a different class of relativistic mean-field model utilized during training, the BNN model effectively retrieves the proton fraction and speed of sound for neutron star matter.

## I Introduction

The extreme matter conditions inside neutron stars (NSs) are impossible to recreate in terrestrial laboratories, making the equation of state (EoS) of dense and asymmetric nuclear matter (realized inside NSs) an interesting and still unknown quantity. Modeling NS matter is restricted by constraints coming from the observations of massive NSs: PSR J1614-2230 [1; 2; 3] with \(M=1.908\pm~{}0.016M_{\odot}\), PSR J0348 - 0432 with \(M=2.01\pm~{}0.04~{}M_{\odot}\)[4], PSR J0740+6620 with \(M=2.08\pm~{}0.07~{}M_{\odot}\)[5], and J1810+1714 with \(M=2.13\pm~{}0.04~{}M_{\odot}\)[6]. Additionally, theoretical calculations, such as chiral effective field theory (cEFT), are applicable only at very low densities, while perturbative quantum chromodynamics (pQCD) is reliable for extremely high densities. Recently, multi-messenger astrophysics has become an exciting field, allowing to understand NSs physics by connecting information carried by different sources, such as gravitational waves, photons, and neutrinos. The detection by LIGO/Virgo collaboration of compact binary coalescence events, such as GW170817 [7] and GW190425 [8], allowed us to further constrain the EoS of NS matter. Recent results from NICER (Neutron star Interior Composition ExploreR) on the PSR J0030+045 pulsar mass [9; 10] and PSR J0740+6620 radius [11; 12; 13] were relevant in restricting the possible neutron stars physics. Future expected observations from experiments such as the enhanced X-ray Timing and Polarimetry mission (eXTP) [14; 15], the (STROBE-X) [16], and Square Kilometer Array [17] telescope will allow for the determination of NSs radii and masses with a few % uncertainty.

Numerous statistical methods have been extensively explored to determine the most probable EoS based on observational data of NSs. These methods include Bayesian inference [18; 19] and Gaussian processes [20]. However, even if the EoS is known with high precision, the challenge remains in constraining the composition of neutron star matter. Previous studies have highlighted the impossibility of recovering nuclear matter properties solely from the \(\beta\)-equilibrium EoS without knowledge of the compositions (or symmetry energy at high densities) [21; 22; 23] or without information about the EoS of symmetric nuclear matter in conjunction with compositions [24]. However, these studies were either limited to meta-models or based on restricted models. Motivated by this, we aim to construct an artificial neural network that directly maps NS observational properties to EoS composition using a large set of EoS derived from the Relativistic Mean Field (RMF) approach.

Deep learning is another field that has become a buzzword in trying to solve the dense matter EoS problems, and all kinds of physics problems [25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43; 44]. The inference problem of determining the EoS from observational data can be roughly divided into two main categories: the reconstructing of the EoS, pressure or speed of sound, from either mass radius or tidal deformability, [25; 26; 27; 28; 29; 30; 31; 32; 34; 35; 37; 40; 41; 44], or focus directly on the nuclear matter saturation properties [33; 36; 38; 39; 42; 43]. As an example, Fujimoto et. al. [25; 26; 27] explored a framework based on neural networks (NNs) where an observational set of NS mass, radius and respective variances was used as input and the speed of sound squared as output. A similar perspective was followed in [30]. The works [33; 36] fall into the second category: determining specific saturation properties of nuclear matter, in this case, the density dependence of the nuclear symmetry energy, directly from observational NS data.

However, the majority of these NNs based models face a considerable drawback, namely, the lack of uncertainty quantification. Questions such as _how confident is a model about its predictions?_ is the main focus of the present work, in which the uncertainty modeling is explored by implementing
a very appealing approach called Bayesian Neural Networks (BNNs), [45]. BNNs have already started being used in different fields of physics [46; 47], and have been shown useful in uncertainty quantification. Our goal is to implement an inference framework that gives a prediction uncertainty to any model prediction. We analyze the density dependence of the proton fraction and speed of sound inside NS matter. For that, several synthetic datasets made of observational data will be constructed and the impact of adding information on the tidal deformability into the model predictions will be analyzed. Let us make clear an important distinction between our work and the majority of studies: instead of applying a widely used approach of parameterizing the EoS, e.g., with polytropes, we used a specific family of nuclear models to construct a set of possible EoS. Despite their flexibility in exploring the entire region of possible EoS, the generic agnostic parametrizations of the EoS are unable to model and track the different degrees of freedom inside NS. The use of a microscopic model has the crucial advantage of accessing the density dependence of all degrees of freedom and thus unable us to analyze the proton fraction.

The paper is organized as follows. A basic introduction to BNNs is presented in Sec II. The family of nuclear models chosen is presented in Sec. III and also the Bayesian inference framework employed to construct the EoS dataset. The generation of the synthetic observation datasets is explained in Sec. IV. The model results for the proton fraction and speed of sound are discussed in Sec. V, and lastly, the conclusions are drawn in Sec. VII.

## II Bayesian neural networks

Despite the great success of (feedforward) neural networks (NNs) in different fields of science, they come with some drawbacks that require special attention. NNs are susceptible to over-fitting and are unable to access the uncertainty of its predictions, which may lead to overconfident predictions. Bayesian Neural Networks (BNNs) is a Bayesian approach framework that introduces stochastic weights to NNs making them uncertainty-aware models [48].

NNs are capable of representing arbitrary functions and are composed, in their simplest architecture, of a sequence of blocks where a linear transformation is followed by a nonlinear operation (activation functions). To simplify the notation, let us denote a NN by \(\mathbf{y}=f_{\mathbf{\theta}}(\mathbf{x})\), where \(\mathbf{\theta}=(\mathbf{W},\mathbf{b})\) represent all NN weights. The vectors \(\mathbf{W}\) and \(\mathbf{b}\) denote, respectively, the connections (weights) and bias of all linear transformations of the network, which define completely the NN model. Training the NN consists in determining the numerical procedure (back-propagation algorithm) of finding the \(\mathbf{\theta}^{*}\) that minimizes a chosen cost function on the training data. This traditional approach of estimating a single model defined by \(\mathbf{\theta}^{*}\) ignores all other possible parametrizations \(\mathbf{\theta}\).

BNNs simulate multiple possible NNs models by introducing stochastic weights. These networks operate by first choosing a functional model, i.e., a network architecture, and then the stochastic model, i.e., the probability distributions for the weights. Bayesian inference is then required to train the network by defining the likelihood function of the observed data, \(P(D|\mathbf{\theta})\), and the prior probability distribution over the model parameters, \(P(\mathbf{\theta})\). It is then possible to employ Bayes theorem and obtain the posterior probability distribution, i.e. the probability of the model parameters given the data:

\[P(\mathbf{\theta}|D)=\frac{P(D|\mathbf{\theta})P(\mathbf{\theta})}{P(D)} \tag{1}\]

where \(P(D)=\int_{\mathbf{\theta}^{\prime}}P(D|\mathbf{\theta}^{\prime})P(\mathbf{\theta}^{\prime })d\mathbf{\theta}^{\prime}\) is the evidence. Having a distribution on the weights, the BNNs predictions become a Bayesian model average: the probability distribution of some unknown \(\mathbf{y}^{*}\) given an input \(\mathbf{x}^{*}\) is

\[P(\mathbf{y}^{*}|\mathbf{x}^{*},D)=\int_{\mathbf{\theta}}P(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{ \theta})P(\mathbf{\theta}|D)d\mathbf{\theta}. \tag{2}\]

\(P(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{\theta})\) is considered to be the likelihood of our data, the distribution that comes out of the network and captures the noise present in our data, and \(P(\mathbf{\theta}|D)\) is the posterior distribution of our weights, that brings up the uncertainty on the model. Another advantage of using these networks is that they capture two types of uncertainty, aleatoric uncertainty, uncertainty on tainty on the data, and epistemic uncertainty, uncertainty on the model estimation defined as \(P(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{\theta})\) and \(P(\mathbf{\theta}|D)\) respectively. However, solving Eq. 2 is a very complex task because the posterior \(P(\mathbf{\theta}|D)\) depends on the evidence \(P(D)=\int_{\mathbf{\theta}^{\prime}}P(D|\mathbf{\theta}^{\prime})P(\mathbf{\theta}^{\prime })d\mathbf{\theta}^{\prime}\), which is a non-analytic expression that requires marginalizing over all model parameters. Fortunately, there are multiple ways of tracking it by using either Markov Chain Monte Carlo or variational inference (more information can be found in [48]). We are going to implement the variational inference method that is presented in the following.

### Varitional inference formalism

The variational inference method aims to approximate a variational posterior \(q_{\mathbf{\phi}}(\mathbf{\theta})\) to the real posterior \(P(\mathbf{\theta}|D)\) by using the Kullback-Leibler (KL) divergence. The KL divergence is a measure of dissimilarity between two probability distributions. It approaches zero when the variational posterior and the true posterior are identical and is positive otherwise. Fundamentally, we want to find the variational posterior that corresponds to the minimum value of the KL divergence between the variational posterior and the true posterior:

\[q_{\phi^{*}}=\operatorname*{arg\,min}_{q_{\phi}}\ \text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{ \theta}|D)), \tag{3}\]

where the KL divergence is defined by

\[\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}|D)) =\mathbb{E}_{q_{\theta}(\mathbf{\theta})}\left[\log\left(\frac{q_{ \phi}(\mathbf{\theta})}{P(\mathbf{\theta}|D)}\right)\right] \tag{4}\] \[=\int_{\mathbf{\theta}}q_{\phi}(\mathbf{\theta})\log\left(\frac{q_{\phi}( \mathbf{\theta})}{P(\mathbf{\theta}|D)}\right)d\mathbf{\theta} \tag{5}\]
In order for the dependence on the true posterior to disappear, the last equation can be rewritten, with the help of Bayes rule Eq. 1, as

\[\text{KL} (q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}|D))=\int_{\mathbf{\theta}}q_{\phi}( \mathbf{\theta})\log\left(\frac{q_{\phi}(\mathbf{\theta})P(D)}{P(D|\mathbf{\theta})P(\mathbf{ \theta})}\right)d\mathbf{\theta}\] \[=\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}))-\mathbb{E}_{q_{ \phi}(\mathbf{\theta})}(\log P(D|\mathbf{\theta}))+\log P(D)\] \[=F(D,\theta)+\log P(D), \tag{6}\]

where \(F(D,\theta)=\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}))-\mathbb{E}_{q_{ \phi}(\mathbf{\theta})}(\log P(D|\mathbf{\theta}))\) is called the variational free energy. We end up with

\[\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}|D))=F(D,\theta)+\log P(D). \tag{7}\]

As the last term, \(\log P(D)\), does not depend on the variational posterior and his parameters, which we are optimizing, minimizing \(\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}|D))\) with respect to \(\phi\) is the same as minimizing \(F(D,\theta)\). ELBO is another important quantity, which stands for evidence lower bound, and it is defined as the negative free energy, i.e., \(\text{ELBO}=-F(D,\theta)\). Equation 7 can then be rewritten as

\[\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}|D))=-\text{ELBO}+\log P(D), \tag{8}\]

or

\[\text{ELBO}=-\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}|D))+\log P(D). \tag{9}\]

ELBO is called the lower bound of the evidence because \(\text{ELBO}\leq\log P(D)\). In other words, \(\text{KL}(q_{\phi}(\mathbf{\theta})||P(\mathbf{\theta}|D))\) is minimized by maximizing the evidence lower bound.

In the end, our optimization objective resumes to

\[q_{\phi^{\star}} =\operatorname*{arg\,min}_{q_{\phi}}\ \text{KL}(q_{\phi}(\mathbf{ \theta})||P(\mathbf{\theta}|D))\] \[=\operatorname*{arg\,max}_{q_{\phi}}\ \text{ELBO}=\operatorname*{arg \,min}_{q_{\phi}}F(D,\theta)\] \[=\operatorname*{arg\,min}_{q_{\phi}}\left[\text{KL}(q_{\phi}( \mathbf{\theta})||P(\mathbf{\theta}))-\mathbb{E}_{q_{\phi}(\mathbf{\theta})}(\log P(D| \mathbf{\theta}))\right].\]

The above general formalism is applied to our specific case, where we have chosen a multivariate Gaussian for the variational posterior, \(q_{\phi}(\mathbf{\theta})=\mathcal{N}(\mathbf{\mu}_{q},\mathbf{\Sigma}_{q})\), and a multivariate Gaussian with diagonal covariance matrix for the prior, \(P(\mathbf{\theta})=\mathcal{N}(\mathbf{0},\mathbf{I})\). The final loss function we are trying to minimize uses Monte Carlo sampling to obtain the expected values, where \(\mathbf{\theta}^{(n)}\) is being sampled from the variational posterior, \(q_{\phi}(\mathbf{\theta})\), and for our specific model, also the exact value of the KL divergence between two Gaussian's with covariance matrix is used,

\[F(D,\phi) =\frac{1}{2D}\left[(-\log\det(\mathbf{\Sigma}_{q}))-k+\text{tr}\left( \mathbf{\Sigma}_{q}\right)+(\mathbf{\mu}_{q})^{T}(\mathbf{\mu}_{q})\right]\] \[-\frac{1}{B}\sum_{i=1}^{B}\frac{1}{N}\sum_{n=1}^{N}\log P(\mathbf{y}_ {i}|\mathbf{x}_{i},\mathbf{\theta}^{(n)}), \tag{10}\]

where \(B\) is the number of points of the dataset, \(k\) is the dimension of the identity matrix of the prior and \(N\) is the number of samples we take of the variational posterior (we used \(N=10^{4}\)). One aspect that stands out in these networks is how back-propagation works. Without going into much detail, the back-propagation updates \(\mathbf{\mu}_{q}\) and \(\mathbf{\Sigma}_{q}\) for each of the networks's parameters (more information can be found in [49]). Once the network is trained and the best mean \(\mathbf{\mu}_{q}\) and covariance matrix \(\mathbf{\Sigma}_{q}\) are obtained for the distribution of the parameters, Eq. 2 becomes solvable and predictions are obtained using Monte Carlo estimations

\[P(\mathbf{y}^{\star}|\mathbf{x}^{\star},D)= \int_{\mathbf{\theta}}P(\mathbf{y}^{\star}|\mathbf{x}^{\star},\mathbf{\theta})q_{ \phi}(\mathbf{\theta})d\mathbf{\theta} \tag{11}\] \[= \frac{1}{N}\sum_{n=1}^{N}P(\mathbf{y}^{\star}|\mathbf{x}^{\star},\mathbf{ \theta}^{(n)}),\quad\mathbf{\theta}^{(n)}\sim q_{\phi}(\mathbf{\theta}). \tag{12}\]

The mean \(\hat{\mathbf{\mu}}\) and variance \(\hat{\mathbf{\sigma}}^{2}\) vectors of the predicting distribution \(P(\mathbf{y}^{\star}|\mathbf{x}^{\star},D)\) can be calculated, for a fixed \(\mathbf{x}^{\star}\), by applying the law of total expectation and total variance. From

\[\mathbb{E}[\mathbf{y}^{\star}|\mathbf{x}^{\star},D]=\mathbb{E}_{q_{\phi}(\mathbf{\theta})} \left[\mathbb{E}[\mathbf{y}^{\star}|\mathbf{x}^{\star},\mathbf{\theta}]\right] \tag{13}\]

and

\[\operatorname{Var}\left[\mathbf{y}^{\star}|\mathbf{x}^{\star},D\right]= \,\mathbb{E}_{q_{\phi}(\mathbf{\theta})}\left[\operatorname{Var}\left[ \mathbf{y}^{\star}|\mathbf{x}^{\star},\mathbf{\theta}\right]\right]\] \[-\operatorname{Var}_{q_{\phi}(\mathbf{\theta})}\left[\mathbb{E}[\mathbf{y}^ {\star}|\mathbf{x}^{\star},\mathbf{\theta}]\right], \tag{14}\]

we obtain

\[\hat{\mathbf{\mu}}=\frac{1}{N}\sum_{n=1}^{N}\hat{\mathbf{\mu}}_{\mathbf{\theta}_{n}}, \tag{15}\]

and

\[\hat{\mathbf{\sigma}}^{2}=\underbrace{\frac{1}{N}\sum_{n=1}^{N}\hat{\mathbf{\sigma}}_{ \mathbf{\theta}_{n}}^{2}}_{\text{Aleator uncertainty}}+\underbrace{\frac{1}{N}\sum_{n=1}^{N}( \hat{\mathbf{\mu}}_{\mathbf{\theta}_{n}}-\hat{\mathbf{\mu}})\odot(\hat{\mathbf{\mu}}_{\mathbf{ \theta}_{n}}-\hat{\mathbf{\mu}})}_{\text{Epistemic uncertainty}}, \tag{16}\]

where \(\odot\) denotes element-wise multiplication. The predicted variance captures both epistemic and aleatoric uncertainties [50].

## III Nuclear models

A field theoretical approach is adopted to calculate a dataset of nuclear equations of state (EoSs). The approach incorporates self-interactions and mixed meson terms within a relativistic mean field (RMF) description. A wide and reasonable region of the parameter space is considered, providing an accurate representation of presently known nuclear properties. Including non-linear terms is crucial for determining the density dependence of the EoS. In this treatment, the nucleons interact through the exchange of scalar-isoscalar mesons (\(\sigma\)), vector-isoscalar mesons (\(\omega\)), and vector-isovector mesons (\(\rho\)). The Lagrangian governing the baryonic degrees of freedom
can be expressed as follows: \(\mathcal{L}=\mathcal{L}_{N}+\mathcal{L}_{M}+\mathcal{L}_{NL}\) with

\[\mathcal{L}_{N}= \bar{\Psi}\Big{[}\gamma^{\mathbf{\mu}}\left(i\partial_{\mathbf{\mu}}-g_{ \omega}\omega_{\mathbf{\mu}}-g_{\varrho}\mathbf{t}\cdot\mathbf{\varrho}_{\mathbf{\mu}}\right)\] \[-\left(m-g_{\sigma}\sigma\right)\Big{]}\Psi\] \[\mathcal{L}_{M}= \frac{1}{2}\left[\partial_{\mathbf{\mu}}\sigma\partial^{\mathbf{\mu}} \sigma-m_{\sigma}^{2}\sigma^{2}\right]\] \[-\frac{1}{4}F_{\mathbf{\mu}\nu}^{(\omega)}F^{(\omega)\mathbf{\mu}\nu}+ \frac{1}{2}m_{\omega}^{2}\omega_{\mathbf{\mu}}\omega^{\mathbf{\mu}}\] \[-\frac{1}{4}\mathbf{F}_{\mathbf{\mu}\nu}^{(\varrho)}\cdot\mathbf{F}^{(\varrho )\mathbf{\mu}\nu}+\frac{1}{2}m_{\varrho}^{2}\mathbf{\varrho}_{\mathbf{\mu}}\cdot\mathbf{ \varrho}^{\mathbf{\mu}}\] \[\mathcal{L}_{NL}= -\frac{1}{3}b\,m\;g_{\sigma}^{3}(\sigma)^{3}-\frac{1}{4}cg_{ \sigma}^{4}(\sigma)^{4}+\frac{\xi}{4!}g_{\omega}^{4}(\omega_{\mathbf{\mu}}\omega^{ \mathbf{\mu}})^{2}\] \[+\Lambda_{\omega}g_{\varrho}^{2}\mathbf{\varrho}_{\mathbf{\mu}}\cdot\mathbf{ \varrho}^{\mathbf{\mu}}g_{\omega}^{2}\omega_{\mathbf{\mu}}\omega^{\mathbf{\mu}},\]

where the field \(\Psi\) denotes the Dirac spinor that describes the nucleon doublet (neutron and proton) with a bare mass \(m\), \(\gamma^{\mathbf{\mu}}\) are the Dirac matrices, and \(\mathbf{t}\) is the isospin operator. The vector meson tensors are defined as \(F^{(\omega,\cdot\cdot)\mathbf{\mu}\nu}=\partial^{\mathbf{\mu}}A^{(\omega,\cdot\cdot)\bm {\nu}}-\partial^{\nu}A^{(\omega,\cdot)\mathbf{\mu}}\). The couplings of the nucleons to the meson fields \(\sigma\), \(\omega\), and \(\varrho\) are denoted by \(g_{\sigma}\), \(g_{\omega}\), and \(g_{\varrho}\), respectively. The meson masses are given by \(m_{\sigma}\), \(m_{\omega}\), and \(m_{\varrho}\). More information on the specifics of the model can be found in [19] and references therein. The parameters \(g_{\sigma}\), \(g_{\omega}\), \(g_{\rho}\), \(b\), \(c\), \(\xi\), and \(\Lambda_{\omega}\) are systematically sampled within a Bayesian framework, adhering to minimal constraints imposed by several nuclear saturation properties. Furthermore, these parameters are subject to the conditions of the neutron star maximum mass exceeding \(2M_{\odot}\), as well as the EoS for low-density pure neutron matter, which is meticulously generated through a precise N3LO calculation in chiral effective field theory. A detailed discussion on these aspects will be presented in the subsequent subsection.

### The Bayesian setup

Based on observed or fitted data, a prior belief (expressed as a prior distribution) is updated using Bayesian inference. The posterior distribution is derived according to Bayes' theorem [51]. In order to establish a Bayesian parameter optimization system, four key components must be defined: the prior, the likelihood function, the fit data, and the sampler.

_The Prior:-_ A broad range of nuclear matter saturation properties is carefully considered in the prior domain of the adopted RMF model. As a result of Latin hypercube sampling, we determine the prior range in our Bayesian setup. Uniform priors are chosen for each parameter, as described in Table 1.

_The Fit Data:-_ The fit data, presented in Table 2, include the nuclear saturation density \(\rho_{0}\), the binding energy per nucleon \(\epsilon_{0}\), the incompressibility coefficient \(K_{0}\), and the symmetry energy \(J_{\text{sym},0}\), all evaluated at \(\rho_{0}\). Additionally, we incorporate the pressure of pure neutron matter (PNM) at densities of 0.08, 0.12, and 0.16 fm\({}^{-3}\) from N\({}^{3}\)LO calculations in chiral effective field theory (\(\chi\)EFT) [52], accounting for 2\(\times\) N\({}^{3}\)LO data uncertainty. Furthermore, the likelihood also includes the requirement of the neutron star maximum mass exceeding 2.0 \(M_{\odot}\) with uniform probability.

_The Log-Likelihood:-_ We optimize a log-likelihood function as a cost function for the given fit data in Table 2. Equation 17 represents the log-likelihood function, taking into account the uncertainties \(\sigma_{j}\) associated with each data point \(j\). The maximum mass of neutron stars is treated differently, using step function probability,

\[Log(\mathcal{L})\propto-\sum_{j}\left\{\left(\frac{d_{j}-m_{j}(\mathbf{\theta})}{ \sigma_{j}}\right)^{2}+Log(2\pi\sigma_{j}^{2})\right\}. \tag{17}\]

To populate the six-dimensional posterior, we employ the nested sampling algorithm [53], specifically the PyMultinest sampler [54; 55], which is well-suited for low-dimensional problems. The EoS dataset for subsequent analyses will be generated using the full posterior, which contains 25287 EoS.

\begin{table}
\begin{tabular}{c c c c} \hline \hline  & \multicolumn{2}{c}{Constraints} & \\ Quantity & & Value/Band & Ref \\ \hline  & \(\rho_{0}\) & \(0.153\pm 0.005\) & [56] \\  & \(\epsilon_{0}\) & \(-16.1\pm 0.2\) & [57] \\  & \(K_{0}\) & \(230\pm 40\) & [58; 59] \\  & \(J_{\text{sym},0}\) & \(32.5\pm 1.8\) & [60] \\ \multicolumn{2}{c}{PNM} & & \\
[MeV fm\({}^{-3}\)] & \(P(\rho)\) & \(2\times\) N\({}^{3}\)LO & [52] \\  & \(dP/d\rho\) & \(>0\) & \\  & [\(M_{\odot}\)] & \(M_{\text{max}}\) & \(>2.0\) & [5] \\ \hline \end{tabular}
\end{table}
Table 2: The Bayesian inference imposes constraints on various quantities to generate sets of models. These constraints include the binding energy per nucleon \(\epsilon_{0}\), incompressibility \(K_{0}\), and symmetry energy \(J_{\text{sym},0}\) at the nuclear saturation density \(\rho_{0}\), each with a 1\(\sigma\) uncertainty. Additionally, the pressure of pure neutron matter (PNM) is considered at densities of 0.08, 0.12, and 0.16 fm\({}^{-3}\), obtained from a \(\chi\)EFT calculation [52]. The likelihood incorporates a 2\(\times\) N\({}^{3}\)LO uncertainty for the PNM pressure, noting that it increases with density. Furthermore, the maximum mass of neutron stars is constrained to be above \(2M_{\odot}\).

\begin{table}
\begin{tabular}{c c c c} \hline \hline  & & \multicolumn{2}{c}{_Set 0_} \\ \cline{3-4} No & Parameters & min & max \\ \hline
1 & \(g_{\sigma}\) & 6.5 & 15.5 \\
2 & \(g_{\omega}\) & 6.5 & 15.5 \\
3 & \(g_{\varrho}\) & 6.5 & 16.5 \\
4 & \(B\) & 0.5 & 9.0 \\
5 & \(C\) & -5.0 & 5.0 \\
6 & \(\xi\) & 0.0 & 0.04 \\
7 & \(\Lambda_{\omega}\) & 0 & 0.12 \\ \hline \end{tabular}
\end{table}
Table 1: We use a uniform prior range for the parameters of the RMF models. Specifically, B and C are \(b\times 10^{3}\) and \(c\times 10^{3}\), respectively. Distribution minimums and maximums are indicated by ’min’ and ’max’ respectively.
Dataset

Our goal is to train BNNs models to predict the speed of sound and proton fraction of NS matter from a given set of NS mock observations. To understand the effect of different NS properties on the prediction uncertainty, we generate different datasets using the 25287 EoS that were obtained through the Bayesian analysis formalism (see Sec. III.1).

### Structure

Our BNN model, \(P(\mathbf{Y}|\mathbf{X},\mathbf{\theta})\), attributes a probability distribution to \(\mathbf{Y}\) (output space) given a set of NS mock observations \(\mathbf{X}\) (input space), where \(\mathbf{Y}\) denotes the different NS matter properties under study, i.e., the speed of sound \(\mathbf{v_{s}^{2}}(\mathbf{n})\) and proton fraction \(\mathbf{y_{p}}(\mathbf{n})\). We have chosen to characterize each element of \(\mathbf{Y}\) at \(15\) fixed baryonic densities \(n_{k}\), e.g., \(\mathbf{y_{p}}(\mathbf{n})=[y_{p}(n_{1}),y_{p}(n_{2}),...,y_{p}(n_{15})]\). The density points are equally spaced between \(n_{1}=0.15\) fm\({}^{-3}\) and \(n_{15}=1.0\) fm\({}^{-3}\), \(n_{k}=\{0.15,0.21,0.27,...,1.0\}\) fm\({}^{-3}\). The number of points (\(N_{Y}=15\)) was selected as a trade-off between computational training time and the interpolation accuracy, i.e. low residuals between the interpolation and the real values. We illustrate this discretization process of the output space elements for the proton fraction in Fig. 1.

Regarding the structure of the input space \(\mathbf{X}\), two different element structures are studied: i) \(\mathbf{X}=[M_{1},...,M_{5},R_{1},...,R_{5}]\) corresponding to five \(M_{i}(R_{i})\) simulated observations and ii) \(\mathbf{X}=[M_{1},...,M_{5},R_{1},...,R_{5},M_{1}^{{}^{\prime}},...,M_{5}^{{}^{ \prime}},\Lambda_{1},....,\Lambda_{5}]\) corresponding to five \(M_{i}(R_{i})\) and five \(\Lambda_{j}(M_{j}^{{}^{\prime}})\) simulated observations. In summary, the output elements \(\mathbf{Y}_{i}\) of our datasets are specified by 15-dimensional vectors and the input space elements \(\mathbf{X}_{i}\) by 10-dimensional or 20-dimensional vectors, depending on the dataset type under study. The statistical procedure for generating the different synthetic observational datasets is presented in the following.

### Generation

The first step of the generation of the datasets consists in randomly splitting the total number of EoS into train and test sets in a proportion of 80%/20%, i.e., the train set contains 22758 EoS while the test has 2529 EoS. Secondly, we generate two types of datasets that share the \(\mathbf{Y}_{i}\) structure but with different \(\mathbf{X}_{i}\) structures:

\[\mathbf{X_{i}} =[M_{1},...,M_{5},R_{1},...,R_{5}]\] \[\mathbf{X_{i}^{{}^{\prime}}} =[M_{1},...,M_{5},R_{1},...,R_{5},M_{1}^{{}^{\prime}},...,M_{5}^{ {}^{\prime}},\Lambda_{1},....,\Lambda_{5}].\]

These different structures allows us to compare and access how informative is the tidal deformability on the model predictions. The statistical generating procedure is composed of the following steps. For each EoS, we randomly select 5 NS mass values, \(M_{i}^{(0)}\), from a uniform distribution between \(1.0M_{\odot}\) and \(M_{\text{max}}\). Then, the radius \(R_{i}\) is sampled from a Gaussian distribution centred at the TOV solution, denoted as \(R(M_{i}^{(0)})\), and with a standard deviation of \(\sigma_{R}\). Finally, we sample the final NS mass from a Gaussian distribution centered at \(M_{i}^{(0)}\) and a standard deviation of \(\sigma_{M}\). The above process can be summarized in the following equations:

\[M_{i}^{(0)} \sim\mathcal{U}[1,M_{\text{max}}]\quad(\text{in units of M}_{\odot}) \tag{18}\] \[R_{i} \sim\mathcal{N}\left(R\left(M_{i}^{(0)}\right),\sigma_{R}^{2}\right)\] (19) \[M_{i} \sim\mathcal{N}\left(M_{i}^{(0)},\sigma_{M}^{2}\right),\quad i=1,...,5 \tag{20}\]

The final generated elements consists of \(\mathbf{X_{i}}=[M_{1},...,M_{5},R_{1},...,R_{5}]\), and each one is a possible realization (_observation_) that characterizes the \(M(R)\) diagram of the specific EoS. This procedure is similar to the one used in [25], where a Gaussian noise was applied to 15 values from \(M(R)\) curve, and then shifted from the original mass radius curve: \(M_{i}=M_{i}^{(0)}+\mathcal{N}\left(0,\sigma_{M}^{2}\right)\) and \(R_{i}=R\left(M_{i}^{(0)}\right)+\mathcal{N}\left(0,\sigma_{R}^{2}\right)\), for \(i=1,...,15\). The second kind of datasets includes the tidal deformability and has the additional steps

\[M_{j}^{{}^{\prime}} \sim\mathcal{U}[1,M_{\text{max}}]\quad(\text{in units of M}_{\odot}) \tag{21}\] \[\Lambda_{j} \sim\mathcal{N}\left(\Lambda(M_{j}^{{}^{\prime}}),\sigma_{ \Lambda}^{2}(M_{j}^{{}^{\prime}})\right)\quad j=1,...,5, \tag{22}\]

where \(\Lambda(M_{j}^{{}^{\prime}})\) is given by the \(\Lambda(M)\) relation of the specific EoS and \(\sigma_{\Lambda}(M_{j}^{{}^{\prime}})\) describes an overall dispersion around the mean value. Samples values with \(\Lambda<0\) were discarded. The functional form \(\sigma_{\Lambda}(M_{j}^{{}^{\prime}})\) should reflect our expectation on the mock observational uncertainty as a function of the NS mass. In the present work, we considered \(\sigma_{\Lambda}(M_{j}^{{}^{\prime}})=\text{constant}\times\hat{\sigma}(M_{j}^ {{}^{\prime}})\), where \(\hat{\sigma}(M_{j}^{{}^{\prime}})\) is the standard deviation of \(\Lambda(M)\) determined from dataset of EoSs. The generated point is \(\mathbf{X_{i}}=[M_{1},...,M_{5},R_{1},...,R_{5},M_{1}^{{}^{\prime}},...,M_{5}^{{ }^{\prime}},\Lambda_{1},....,\Lambda_{5}]\) (a similar approach can be found in [37]). In the above procedures, there is an additional parameter, which we denote by \(n_{\text{s}}\), that specifies the number of mock observations for
each EoS, i.e., the number of times the above procedures are applied to each EoS. For instance, choosing \(n_{s}=20\) would mean running the above procedures 20 times for each EoS (20 _observations_), and thus obtaining \(\{\mathbf{X_{1}},\mathbf{X_{2}},...,\mathbf{X_{20}}\}\).

Applying the above formalism to both the train and test sets, we have generated a total of 4 datasets whose properties are displayed in Table 3. Sets 1 and 2 only contain information about the NS radii (input space \(\mathbf{X}\) is 10-dimensional) while sets 3 and 4 also include the tidal deformability (input space \(\mathbf{X}\) is 20-dimensional). The analysis of sets 1 and 2 allows to understand how a decrease in the spread of the mock observations around the TOV solution affects the predictions and uncertainties. In the same manner, sets 3 and 4 aim to understand possible effects on the model predictions arising from an increase of simulated observations scattering around the mean value on the tidal deformability. We use 60 mock observations, \(n_{s}=60\), on the training sets for each EoS while \(n_{s}=1\) was employed for the test sets. This key difference tries to simulate a real case scenario in which we only have access to a single mock observation of the _true_ EoS. Here, by single mock observation, we mean \(n_{s}=1\) that corresponds to five \(M_{i}(R_{i})\) mock observations (sets 1 and 2) or five \(M_{i}(R_{i})\) mock observations and five \(\Lambda_{j}(M_{j})\) mock observations (sets 3 and 4). To illustrate the dataset generation, Fig. 2 displays the 60 mock observations for two distinct EoSs that belong to the generated datasets (dataset 1 and 2 on the right and left figures, respectively). Note that, for each EoS, there are 300 points in the \(M(R)\) diagram, i.e., consisting of 60 EoS simulated observations for 5 NS mock observations each. It is clear the increase of \(\sigma_{R}\) from dataset 1 to 2, highlighting the differences between the two datasets. Furthermore, we provide Figure 3 to depict the tidal deformability values for datasets 3 and 4.

### Training procedure

To assess the response of Bayesian Neural Networks (BNNs) to varying input noises and output targets, we have conducted experiments involving the training of diverse functional and stochastic models (as explained in Section II). These BNN models were trained using distinct datasets generated as outlined in Section IV.2. During the training stage, a subset of the training data was randomly selected as validation step, i.e., the training data was split into 80% actually for training and 20% for validation. Moreover, the input data \(\mathbf{X}\) was standardized.

Defining the functional models involves adjusting the number of neurons, layers, and activation functions. Table 4 shows the best functional models for each dataset mentioned in Table 3. For the hidden layers, we explore hyperbolic tangent, softplus, and sigmoid activation functions, while utilizing a linear activation function for the output layer. As the input vector sizes differ (10 for sets 1 and 2, and 20 for sets 3 and 4), we employ more neurons per layer for larger input spaces. This is due to the increased complexity demanded by a greater number of parameters. Specifically, we use 15 and 10 neurons for sets 1 and 2, and 20 and 25 neurons for sets 3 and 4, respectively. The output layer consistently contains 30 neurons, with 15 representing the mean and 15 representing the standard deviation of the output probability distribution function. It is worth noting that we deliberately excluded the use of correlation in the output layer due to the inferior performance observed when attempting to incorporate it. As a result, the output layer is solely focused on capturing the mean and standard deviation information of the output distribution. The architecture employed in this study involved

Figure 3: The \(n_{s}=60\) mock observations generated in the \(\Lambda-M\) diagram for two EoSs in dataset 3 (left) and dataset 4 (right). The grey area represents the extremes of our EoS dataset. The two EoS coincide with the ones used in Fig. 1.

\begin{table}
\begin{tabular}{c c c c} \hline \hline Dataset & \(\sigma_{M}\)\([M_{\odot}]\) & \(\sigma_{R}\) [km] & \(\sigma_{\Lambda}(M_{j})\) \\ \hline
1 & 0.05 & 0.15 & — \\ \hline
2 & 0.1 & 0.3 & — \\ \hline
3 & 0.1 & 0.3 & 0.5\(\hat{\sigma}(M_{j})\) \\ \hline
4 & 0.1 & 0.3 & 2\(\hat{\sigma}(M_{j})\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Generation parameters for each dataset. \(\hat{\sigma}(M_{j})\) denotes the standard deviation of \(\Lambda(M)\) calculated on the train set.

Figure 2: The \(n_{s}=60\) mock observations generated for two EoSs in dataset 1 (left) and dataset 2 (right). The grey area represents the extremes of our EoS dataset. The two EoS coincide with the ones used in Fig. 1.
utilizing two to three hidden layers. The number of neurons within each hidden layer remained consistent, but it varied depending on the size of the input vector, as explained earlier. During the grid search process, we systematically explored four different architectures for each output. However, we narrowed our focus to datasets 1 and 3, as we specifically aimed to identify the most suitable architecture for the two different input sizes. The best outcomes are obtained by employing sigmoid as the activation function in the hidden layers, ensuring minimal loss and preventing divergence. Across all eight dataset configurations, i.e., two outputs (\(v_{s}^{2}\) and \(y_{p}\)) and four datasets, see Table 3, it was found that the optimal number of hidden layers is two. For sets 1 and 2, the best performance was achieved with 15 neurons in each hidden layer, while for sets 3 and 4, 25 neurons were utilized in each hidden layer. Detailed information on these configurations can be found in Table 4 for the two output variables. During training, we employ a learning rate of 0.001 and utilize the ADAM optimizer [61] with the AMSgrad improvement [62]. The models are trained for 4000 epochs, with a mini-batch size of 768.

Regarding the stochastic model, we adopt a Gaussian prior with mean zero and standard deviation of one as mentioned in Section II.1. While this prior choice lacks a specific theoretical justification, it serves as a reasonable default prior, as discussed in [48]. Future research could delve further into investigating the impact of prior parameters, similar to the approach taken in reference [46]. Additionally, we select a multivariate normal distribution as the variational posterior as explained in Section II.1, initialized with mean 0 and a diagonal covariance matrix, where the standard deviation is equal to \(\log(1+\exp{0})=0.693\). Furthermore, we opt for a deterministic output layer instead of a probabilistic one, as it has demonstrated improved results in our experiments. This decision is motivated by the fact that the deterministic output layer aligns better with the specific requirements of our model architecture and the nature of the problem we are addressing. All BNNs models were coded using Tensorflow library [63], more specifically we use Keras [64], an high-level API of the TensorFlow.

## V Results: neutron star matter properties

In the following, we discuss the results for the speed of sound squared \(v_{s}^{2}(n)\) and proton fraction \(y_{p}(n)\). To analyze how the observational uncertainty in \(R\) and \(\Lambda\) affects the model confidence, we are going to compare the results from the different models of Tab. 4, which were trained on the different datasets presented in Tab. 3. For the sake of simplicity, whenever we want to distinguish the different 8 BNNs models, we will, hereafter, refer to their training (data)sets, 1 to 4, and predicting quantity, \(v_{s}^{2}(n)\) or \(y_{p}(n)\). For instance, when we say the results of dataset 2 for \(y_{p}(n)\), we mean the BNN model trained on dataset 2 that has the proton fraction as the target quantity.

### Speed of sound

Let us start by illustrating the different BNNs predictions for the speed of sound squared on a randomly selected EoS from the test set. The results are displayed in Figure 4, where the top panel shows sets 1 (blue) and 2 (orange) models while the below figure shows sets 3 (green) and 4 (purple). Model predictions are computed using Eq. 2, from which we show the mean values (solid lines) and \(2\sigma\) regions (color regions). The upper plot clearly shows that the prediction uncertainty, characterized by the distribution's standard deviation \(\sigma\), is smaller on set 1 than on set 2, and most importantly, the predicted mean values are close to the real values. The same pattern is seen in the lower figure: set 3 model (purple) has a lower prediction variance than set 4. While a deeper understanding of the overall behaviour requires that we analyze the whole test set, Fig. 4 already indicates that the BNNs models are able to capture the characteristics of the different datasets (see Table 3): the increased dispersion of the NS mock observations around the true values translate into a larger uncertainty when inferring the corresponding EoS properties.

To investigate how the models predictions behave over the entire test set, we define the normalized residuals' predictions as \(\Gamma(n_{k})=\left(v_{s}^{2}(n_{k})-v_{s}^{2}(n_{k})^{\text{true}}\right)/ \sigma(n_{k})\) and the dispersion by \(\Sigma(n_{k})=\sigma(n_{k})\) at each of the prediction densities, i.e., \(k=1,...,15\). A summary statistics of both quantities over all EoS of the test set is shown in Fig. 5. The distribution of \(\Gamma(n)\) (top panel), for all 4 datasets, has 50 % of the values near zero, indicating that the median of the prediction values is unbiased. Furthermore, at the 2.3 % and 97.7% of the cumulative percentage, we can see that the distribution lies between -2\(\sigma\) and 2\(\sigma\) respectively, indicating that the prediction mean deviates from the true value less than 2\(\sigma\) 95% of the times. The fact that the distribution properties of \(\Gamma(n)\) are similar across all datasets and independent of the density reveals that the BNNs models are correctly modelling the dispersion of the predictions considering the corresponding mean residuals at each density value. The value of \(\Sigma(n)\) (bottom panel) grows with increasing density reflecting the training set statistics: the EoS dataset was

\begin{table}
\begin{tabular}{c|c|c|c} \multirow{2}{*}{**Layers**} & \multirow{2}{*}{**Activation**} & \multicolumn{2}{c}{**Neurons**} \\ \cline{3-4}  & & **Dataset 1 \& 2** & **Dataset 3 \& 4** \\ \hline Input & N/A & 10 & 20 \\ \hline Hidden Layer 1 & Sigmoid & 15 & 25 \\ \hline Hidden Layer 2 & Sigmoid & 15 & 25 \\ \hline Output & Linear & 30 & 30 \\ \hline \end{tabular}
\end{table}
Table 4: Structures of the final BNN models. The \(v_{s}^{2}(n)\) and \(y_{p}(n)\) models have the same structure.
generated by Bayesian inference where saturation properties were imposed, leading to a wider uncertainty at higher densities while low density regions are strongly constrained. From the \(\Sigma(n)\) plot, we see that the whole distribution of the BNN model trained on set 1 (blue line) shifts to lower values, specifically on \(2.3\%\) and \(50\%\) of the values, showing that there is a considerable decrease in uncertainty when the dispersion of NS mock observations is reduced by a factor of two, from (\(\sigma_{M}=0.1M_{\odot},\sigma_{R}=0.3\) km) to (\(\sigma_{M}=0.05M_{\odot},\sigma_{R}=0.15\) km) (see Table 3).

To estimate the overall performance of the different BNNs models, we show the coverage probability of each model in Fig. 6. The coverage probability quantifies how the model is perceiving the distribution of the data, by determining if the percentage of values contained in \(1\sigma\) of the output distribution, i.e., the number of values in a specific interval divided by the total number of values, corresponds to 68 % of the number of values we are using for the test set. The same analogy is repeated for \(2\sigma\) (95%) and \(3\sigma\) (99.7%), this was implemented for each of the 15 values of the output independently, obtaining the respective coverage probabilities relative to the output densities, and then we estimated the mean of this 15 coverage probabilities represented in Fig. 6 for the 4 sets. The overall results show that the model is correctly estimating the distribution of the data, since the bars are very close to the three percentages, a little fluctuation can be seen in the 68%, where the 4 sets are overestimating the uncertainty of the results, meaning, we have a percentage

Figure 4: The BNNs predictions for \(v_{s}^{2}(n)\) using one EoS of the test. The models trained on datasets 1 (blue) and 2 (orange) are in the upper figure while datasets 3 (purple) and 4 (green) models are in the lower figure. The prediction mean values (solid lines) and \(2\sigma\) confidence intervals are shown. The true values are shown in black dots and the range of \(v_{s}^{2}(n)\) from the train set is indicated by the grey region.
of data bigger than 68% in \(1\sigma\) of the model, so the \(\sigma\) of the model should be smaller, this can be seen even more for set 3.

To quantify how the increase in the mock observational scattering of \(\left(R,\Lambda\right)\) affects the model prediction uncertainties, let us define the following quantity

\[\eta\left[a,b\right]\left(n_{k}\right)=\frac{1}{T}\left(\sum_{i=1}^{T}\frac{ \sigma_{i}^{b}(n_{k})-\sigma_{i}^{a}(n_{k})}{\sigma_{i}^{a}(n_{k})}\right) \times 100, \tag{23}\]

where \(T\) is the total number of EoS in the test set, and \(k=1,...,15\). This quantity defines the percentage uncertainty deviation between models \(b\) and \(a\) at density \(n_{k}\). Figure 7 shows the results, where we plotted 5 different comparisons. The first conclusion, looking at the result for \(\eta\left[2,1\right]\) (cyan), is that the prediction uncertainty increased when we consider the BNN model of dataset 1 compared with the one from dataset 2. Furthermore, the \(\eta\left[2,1\right]\) reached the maximum value of 20% at \(n_{4}=0.33\) fm\({}^{-3}\). In other words, when the (synthetic) observational data scattering doubles, from \(\left(\sigma_{M}=0.05M_{\odot},\sigma_{R}=0.15\text{ km}\right)\) to \(\left(\sigma_{M}=0.1M_{\odot},\sigma_{R}=0.3\text{ km}\right)\), the uncertainty increases of the order 4.9-20%, having the largest value at the densities where the NS radius, \(R(M)\), and \(v_{s}^{2}(n)\) have the highest correlation (see Anne A for details and ref. [65]). The region of larger sensitivity to the uncertainty of the mock observational data coincides with the density interval where the speed of sound increases steadily, and, in many agnostic approaches attains a maximum followed by a decrease or flattening at larger densities [66; 67; 68]. The nucleonic EoS that has been used to train the model also shows similar behavior for the speed of sound [19]. This behavior of the speed of sound for densities below three times saturation density is dictated by the two solar mass constraints.

The second important conclusion is the impact of adding the \(\Lambda\) information on the inference properties. This point is clear when analysing the dependence and values of \(\eta\left[3,2\right]\) (the difference between datasets 2 and 3 is that the latter contains information on the tidal deformability, see Tab. 3). The negative values reflect the fact that the prediction uncertainty decreases when the tidal information is added to the training procedure - the tidal deformability is informative of the \(v_{s}^{2}(n)\) of neutron star matter. The maximum value of the uncertainty decrease is 7%, occurring at \(n_{6}=0.45\) fm\({}^{-3}\). Similarly, when considering \(\eta\left[4,2\right]\) (blue), it also exhibits negative values throughout, indicating a decrease in uncertainty compared to dataset 2. However, beyond the seventh density, the two values become very close to each other, suggesting that no additional information is being perceived by including more dispersion on the tidal deformability. Looking at \(\eta\left[4,3\right]\) and comparing it with \(\eta\left[2,1\right]\), we see that for the first density, they differ very little, but from there on \(\eta\left[2,1\right]\) is approximately always more than 5 times bigger. Noting an important consideration here, the proportion of input values being altered: dataset 3 and 4 involves changing only the uncertainty of 5 quantities out of the 20 input values, i.e. changing a quarter of our input vector, while in dataset 1 and 2 we have changed all the input values. So we could anticipate an at least four-fold increase in \(\eta\) for datasets 1 and 2, however, this percentage is the majority of times even bigger. This implies that the dispersion of the mass-radius pairs has a more significant impact on the model compared to the mass-tidal deformability pairs. By acknowledging this difference in the proportion of modified input values, we gain insight into how the model's understanding and response to dispersion changes are influenced.

### Proton fraction

Let us now analyze the model's predictions for the proton fraction \(y_{p}(n)\). Using a specific EoS from the test dataset (for illustration purposes), we show in Fig. 8 the models' prediction for each dataset: 1 (blue) and 2 (orange) on the upper

Figure 6: Coverage probability calculated on the test set of the \(v_{s}^{2}(n)\) BNNs models.

Figure 7: Prediction uncertainty deviation \(\eta\left[a,b\right]\) between the \(v_{s}^{2}(n)\) BNN models \(a\) and \(b\) (see text for details).
panel and 3 (purple) and 4 (green) in the lower panel. The range of \(y_{p}(n)\) from the train set is indicated by the grey region and the dashed grey line displays the \(99.9\%\) percentage of data/probability line. The train statistics is important to point out that the upper region between the \(99.9\%\) probability and maximum boundary lines is caused by the presence of just one _extreme_ EoS.

The conclusion drawn from Fig. 8 is similar to the \(v_{s}^{2}(n)\) results (see Fig. 4): it is evident that dataset 1 exhibits the narrowest prediction uncertainty, whereas dataset 2 demonstrates a considerable increase in uncertainty; the prediction uncertainties are similar for datasets 3 and 4. Definite conclusions require some statistics over a set of EoS, and, for that purpose, we are going to analyze the whole test set once again.

Figure 9 shows the following quantities over the four datasets: the model residuals \(\delta(n_{k})=y_{p}(n_{k})-y_{p}(n_{k})^{\text{true}}\) (left), the standard deviation \(\Sigma(n_{k})=\sigma(n_{k})\) (center), and the normalized model residuals \(\Gamma(n_{k})=\left(y_{p}(n_{k})-y_{p}(n_{k})^{\text{true}}\right)/\sigma(n_{k})\) (right). We calculate these quantities, at each density \(n_{k}\), for each of the 2529 EoS from the test set. We see that the model has a broader residual spread around 0.4-0.5 fm\({}^{-3}\) that is counterbalanced by larger standard deviation values in this region - this is a direct conclusion when looking at the normalized model residuals steadiness. In other words, the model correctly captures and models the data statistics: the BNN models capture larger prediction uncertainties in regions where \(y_{p}(n)\) has a larger dispersion, as expected. The overall quality of the models prediction is seen in the density independence of the \(\Gamma(n)\) (right panel) statistics, the models residuals are \(95\%\) of the times between \(2\sigma\). The insight of the train set statistics, in Fig. 10, entails interpretability for \(\delta(n)\) and particularly \(\Sigma(n)\) behaviour, since it becomes clearer that \(\sigma(n)\) of the train set has a non-monotonic behaviour, reaching a maximum value around 0.5 fm\({}^{-3}\) decreasing for lower/larger densities. The coverage probability for the \(y_{p}(n)\) models are similar to the ones obtained for \(v_{s}^{2}(n)\) (see Fig. 6) and thus the models are correctly estimating the distribution of the test set data.

The dataset comparison for \(\eta\left[a,b\right](n)\) (see Eq. 23) is displayed in Fig. 11. Firstly, we observe that \(\eta\left[2,1\right]\) exhibits a behavior similar to \(v_{s}^{2}\). However, it reaches its maximum value earlier at \(n_{3}\). Interestingly, this is precisely where \(R(M)\) and \(y_{p}(n)\) demonstrate the highest correlation, as explained in Annex A. When comparing \(\eta[2,1]\) with \(\eta[4,3]\), the behaviour is consistent with the one observed in the speed of sound, albeit this time is even more outstanding, which comes back again to the correlation, that is almost zero for \(\Lambda(M)\) and \(y_{p}(n)\). Based on this observation, we can speculate that \(\Lambda\) does not contribute with significant information to the model, especially since this time the ratio between \(\eta[4,3]\) and \(\eta[2,1]\) deviates even further from the expected one-quarter proportion, that comes from the change in the input vector, than it did on the \(v_{s}^{2}\). Increasing the model's complexity without adding substantial information can lead to increased confusion and greater uncertainty, as evident in \(\eta[3,2]\) and \(\eta[4,2]\). These two cases exhibit positive values of \(\eta\), indicating higher uncertainty for sets 3 and 4 when compared with set 2.

### BNNs epistemic and aleatoric uncertainties

Let us briefly analyze how the prediction uncertainty of BNNs is modelled and its different components. We are going to focus on the \(y_{p}(n)\) BNN models. However, all the features discussed below are also seen for the \(v_{s}^{2}(n)\) models. We have seen in Eq. 16 of Sec. II that the prediction variance \(\hat{\mathbf{\sigma}}^{2}\) is a combination of two terms,

\[\hat{\mathbf{\sigma}}^{2}=\hat{\mathbf{\sigma}}^{2}_{\text{alea}}+\hat{\mathbf{\sigma}}^ {2}_{\text{epist}},\]

Figure 8: The BNNs predictions for \(y_{p}(n)\) using one EoS of the test. The models trained on datasets 1 (blue) and 2 (orange) are in the upper figure while datasets 3 (purple) and 4 (green) models are in the lower figure. The predicted mean values (solid lines) and \(2\sigma\) confidence intervals are shown. The true values are shown in black dots.
where the aleatoric uncertainty \(\hat{\mathbf{\sigma}}_{\text{dea}}^{2}\) measures the mean variance of the models' ensemble, while the epistemic uncertainty \(\hat{\mathbf{\sigma}}_{\text{epist}}^{2}\) measures the spread of the models around the ensemble mean \(\hat{\mathbf{\mu}}\). The epistemic uncertainty arises from limited information or data and is encoded on the posterior probability \(P(\mathbf{\theta}|D)\), i.e., the model distribution. On the other hand, aleatoric uncertainty is due to the inherent randomness of the dataset and is encoded into the data likelihood \(P(\mathbf{y}^{\star}|\mathbf{x}^{\star},\mathbf{\theta})\). While the epistemic uncertainty decreases when more data is available, the aleatoric uncertainty value does not depend on the amount of data has it is a property of the generating data process.

To analyze the proportions of both uncertainty types on the total prediction variance, we calculate the epistemic percentage as \(f_{\text{epist}}=(\hat{\mathbf{\sigma}}_{\text{epist}}^{2}/\hat{\mathbf{\sigma}}^{2}) \times 100\%\). Using the BNNs models for \(y_{p}(n)\), we show in the left panel of Fig. 12 the mean (dashed lines) and 68% confidence interval region (colored regions) of \(f_{\text{epist}}\) across the entire test (2 and 3) sets. Additionally, the right panel shows \(f_{\text{epist}}\) using BNN models trained on set 1 but with different numbers of mock observations: \(n_{s}=20\) and \(n_{s}=60\) (the value used in this work). The following conclusions can be drawn: i) the prediction variance \(\hat{\mathbf{\sigma}}^{2}\) is composed mainly of aleatoric uncertainty (left panel), around 95%, and it is due to the already high number of mock observations \(n_{s}=60\); ii) the right panel shows that decreasing the number of mock observations \(n_{s}\), and thus the total number of training points, increases the epistemic uncertainty. The epistemic uncertainty converges to zero when the data points go to infinity; iii) left panel also shows that \(f_{\text{epist}}\) is smaller for set 2 (orange) than set 3 (purple) because the input dimensions increase from 10 to 20, which is reflected on the posterior \(P(\mathbf{\theta}|D)\). Lastly, let us argue why the epistemic uncertainty is larger at densities 0.2 - 0.4 fm \({}^{-3}\). When constructing the predicting ensemble, \(P(\mathbf{y}^{\star}|\mathbf{x}^{\star},D)=\frac{1}{N}\sum_{n=1}^{N}P(\mathbf{y}^{\star}| \mathbf{x}^{\star},\mathbf{\theta}^{(n)})\), by sampling from

Figure 11: Prediction uncertainty deviation \(\eta\left[a,b\right]\) between the \(y_{p}(n)\) BNN models \(a\) and \(b\) (see text for details).

Figure 10: Some statistics of \(y_{p}(n)\) calculated from the train dataset.
the variational posterior, \(\mathbf{\theta}^{(n)}\sim q_{\phi}(\mathbf{\theta})\), the density points \(n_{k}\) with larger correlation with \(y_{p}(n)\) are much more sensitive to model sampling than other density points where correlations are much weaker.

## VI Prediction for the DD2 nuclear model

As a final test, we applied the BNN model (trained on set 1, see Tab. 3) to a nuclear model with different properties from the ones used to train, in particular, obtained within a different microscopic description of nuclear matter. We select the DD2 model which is a generalized relativistic mean-field model (RMF) model with density-dependent couplings [69], which has been calibrated to describe properties of finite nuclei. One key difference between the DD2 and the RMF family we used to generate the set of EoS consists of the high density behavior of the symmetry energy. In DD2 model, the coupling to the \(\rho\)-meson that defines the isovector channel of the EoS goes to zero at sufficiently high densities, favoring very asymmetric matter. One of the main consequences is that nucleonic direct Urea processes inside NS are not predicted by DD2 [70; 71]. Another noticeable difference between the DD2 class of models and the class of models used to train the BNN is the behavior of the speed of sound with density: for DD2 like models the speed of sound increases monotonically, although it remains always well below \(c\), while for the class of models used to train BNN the speed of sound flattens or even decreases above \(\sim 3\rho_{0}\). These two differences will be reflected in the performance of the BNN model.

After selecting the DD2 EoS, following the statistical procedure described in Sec. IV, we generated one mock observation (\(n_{s}=1\)) using the dataset 1 properties, which is the one with lower \(\sigma_{R}\) and does not contain the information about \(\Lambda\). The BNN model predictions for the speed of sound (top panel) and proton fraction (lower panel) are shown in Fig. 13. Despite the DD2 lying outside the training values (grey region) for the speed of sound, the model prediction uncertainty extends beyond the training maximum values and almost contains completely the DD2 results. The \(y_{p}(n)\) prediction is quite good, being the prediction of the mean value close to the real one. The DD2 proton fraction reflects the property described above concerning its favoring large neutron-proton asymmetries due to the behavior of the isovector channel. However, our BNN model was able to capture this behavior. Despite the results being quite compelling, there are some crucial points that we would like to point out. During the above test stage, we generated just one mock observation (\(n_{s}=1\)) from the DD2 \(M(R)\) curve to simulate a real case scenario, where a very limited number of NS observations are accessible. Being the sampling procedure of generating one mock observation (\(n_{s}=1\)), i.e., five \(M_{i}(R_{i})\) values, a random process, implies that different samples will originate different predictions (\(v_{s}^{2}(n)\) is much more sensitive than \(y_{p}(n)\), since the DD2 target is completely inside the training values region). This is a somehow expected behaviour since we are trying to characterize the whole \(M(R)\) with only 5 random \((M,R)\) values. While a given sample may be sufficient to inform the general dependence of the \(M(R)\) (like the one we generated), others might be as well almost uninformative of the actual \(M(R)\) curve, i.e., a sample where all five mock observations \(M_{i}(R_{i})\) cluster around the same \(M\) value. This is a general problem that shows up regardless of the inference model or framework: inferring the EoS from a very limited number of NS observations. The BNN performance assessment for the DD2 EoS would get more reliable as the number of points \(M_{i}(R_{i})\), which compose each mock observation (5 in the present work), increase since a random sample, in that case, would be much more informative of the true \(M(R)\) curve.

## VII Conclusions

We have explored Bayesian Neural Networks (BNNs), which is a probabilistic machine learning model, to predict the proton fraction and speed of sound of neutron star matter from a set of NS mock observations. This method is based upon the usual neural networks but with the crucial advantage of attributing an uncertainty measurement to its predictions. Our EoS dataset was generated from a relativistic mean field approach through a Bayesian framework, where constraints from nuclear matter properties and NS observations were applied. The choice of a specific microscopic nuclear model, instead of a more flexible EoS parameterization, as the ones discussed for instance in [72], is justified because we want to analyze the possibility of inferring the neutron star composition, specifically, the proton fraction, from NS observations. From the set of 25287 EoS, four different mock observational sets, simulating four different scenarios of mock observational uncertainties, were generated. Two of them are only composed of \(M(R)\) simulated observations and the other two have also information regarding \(\Lambda(M)\). In the end, 8 different BNNs were trained to predict the \(v_{s}^{2}(n)\) and \(y_{p}(n)\) in each of the 4 datasets.

With this study, we have shown that using BNNs, the measurements of the mass and radius of five neutron stars allow

Figure 12: The pdf of \(f_{\text{pixel}}=(\hat{\mathbf{\sigma}}_{\text{pixel}}^{2}/\hat{\mathbf{\sigma}}^{2}) \times 100\%\) for \(y_{p}(n)\) BNN models in sets 2 and 3 (left) and for set 1 model but trained in different training datasets with different number of mock observations \(n_{s}\) (right).
us to recover information from the equation of state of nuclear matter with associated uncertainty, not only for a quantity that is more connected with the isoscalar behavior of the EoS, the speed of sound, but also for the proton fraction, a property that is determined by the isovector behavior of the EOS. In several recent works, the attempt to determine the proton fraction from the mass radius measurements was unsuccessful [21; 22; 23]. In all these descriptions a polynomial expansion of the EoS until the third or fourth order has been considered. In [23], this was attributed to the existence of multiple solutions. The authors of [22] identify the correlations among higher-order parameters as a difficulty. The BNN approach allows the model to learn the full density dependence of the EoS avoiding the short comes of the density expansion with a finite number of terms. It was shown that the uncertainty associated with the predicted quantities is particularly sensitive to the precision of the observational data if some kind of correlation exists between the data and the property that is being calculated. For the speed of sound this was reflected in a larger sensitivity for densities below three times saturation density, where the NS radius is strongly correlated with the speed of sound as discussed in [65]. It was also shown that adding extra observational mock data, in particular, the tidal deformability, could decrease the uncertainty associated with the prediction, but not always. There was a clear improvement for the speed of sound but not for the proton fraction. Too scattered data does not bring an improvement on the uncertainty determination due to the increase of the complexity of the model compared with the quality of the data. It is important to point out that the improvement attained with the tidal deformability data with a smaller uncertainty worked for the speed of sound because it shows a correlation with the tidal deformability for densities of the order of twice saturation density similar to the one with the radius. This correlation does not exist between the proton fraction and the tidal deformability so no improvement in the proton fraction prediction was attained when introducing the tidal deformability observation. The proton fraction has shown some sensitivity at twice saturation density to the radius uncertainty and this can be traced back to the existing correlation of low mass star radius with the symmetry energy slope [73], a quantity that strongly determines the proton fraction. This correlation weakens quickly with the increase of the NS mass, and is much weaker with the tidal deformability. We have also tested the BNN model with a mock measurement obtained from the DD2 EoS generated with a microscopic framework different from the one used to generate the EoS used to train the BNN model. The results have confirmed the validity of the model and its predicting power.

We have been very conservative concerning the uncertainties attached to the observations. In the future, observatories such as STROBE-X [16] and eXTP [14] may give us radius measurements with uncertainties as small as 2%-5% and this will improve the predictions as demonstrated in the present study.

There are several potential paths for further improvement and exploration in this work. One possibility is to extend the analysis to include other properties of neutron stars and investigate their relationship with observable quantities. Another possibility for improvement as discussed in the results obtained for the DD2 model, is increasing the number of observable pairs used as input can enhance the model's performance. However, it is worth noting that in the case of Bayesian neural networks (BNNs), expanding the number of pairs introduces a greater increase in the model parameters compared to traditional architectures, which is why we used a lesser amount of pairs compared with previous articles using conventional neural networks, as demonstrated in studies like [25] and related articles. Furthermore, for the stochastic model, it would be interesting to improve the prior as mentioned in section II.1.

###### Acknowledgements.

 This work was partially supported by national funds from FCT (Fundacao para a Ciencia e a Tecnologia, I.P, Portugal) under Projects No. UIDP/04564/2020, No. UIDB/04564/

Figure 13: The BNN model predictions, \(v_{s}^{2}\) (upper) and \(y_{p}\) (lower), for one mock observation (\(n_{s}=1\)) of the DD2 EoS, the blue area represents the 95% confidence interval, and the solid line the mean.
2020 and 2022.06460.PTDC.

## Appendix A Correlation between NS properties and EoS

Figure 14 shows the Pearson correlation coefficient between \(v_{s}^{2}(n)\) and \(R(M)\) (left panel) and \(\Lambda(M)\) (right panel) for specific NS masses (colors) and the average value, by considering \(M/M_{\odot}\in[1,2.2]\). The same is performed for the proton fraction in Fig. 15. The Pearson correlation is calculated as \(\text{Corr}(a,b)=\text{Cov}(a,b)/(\sigma_{a}\sigma_{b})\), where \(a\) consists of \(v_{s}^{2}\) and \(y_{p}\) and \(b\) of \(R\) and \(\Lambda\). Note, however, that this correlation measure is only sensitive to linear dependencies, and higher order ones can be missed. These correlations have been discussed in [65].
