aainstitutetext: Institute for Nuclear Theory, University of Washington, Box
351550, Seattle, WA 98195, USAbbinstitutetext: Department of Physics, The
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japanccinstitutetext: Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto 606-8502 Japanddinstitutetext: Department of Physics, Tokyo
Metropolitan University, Hachioji 192-0397, Japan

# Uncertainty quantification in the machine-learning inference from neutron
star probability distribution to the equation of state

Yuki Fujimoto b Kenji Fukushima b Syo Kamata c,d Koichi Murase<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We discuss the machine-learning inference and uncertainty quantification for
the equation of state (EoS) of the neutron star (NS) matter directly using the
NS probability distribution from the observations. We previously proposed a
prescription for uncertainty quantification based on ensemble learning by
evaluating output variance from independently trained models. We adopt a
different principle for uncertainty quantification to confirm the reliability
of our previous results. To this end, we carry out the MC sampling of data to
infer an EoS and take the convolution with the probability distribution of the
observational data. In this newly proposed method, we can deal with arbitrary
probability distribution not relying on the Gaussian approximation. We
incorporate observational data from the recent multimessenger sources
including precise mass measurements and radius measurements. We also quantify
the importance of data augmentation and the effects of prior dependence.

††preprint: INT-PUB-24-001, YITP-24-05

##  1 Introduction

Unveiling properties of QCD (Quantum Chromodynamics) matter in the regime of
finite baryon density is a vital challenge in modern nuclear physics. Despite
various efforts over several decades, many aspects of finite-density QCD have
remained elusive, such as the nature of QCD transitions (see refs.
Fukushima:2010bq ; Baym:2017whm ; Aarts:2023vsf  for comprehensive reviews).
The equation of state (EoS) of dense matter is of central importance in
describing the properties of QCD matter. Ample knowledge of the EoS, that is,
the behavior of thermodynamic quantities in extreme environments, would help
us reveal the properties of states emerging from QCD. This approach has indeed
worked successfully in the case of QCD at small density and high temperature,
where the EoS from numerical simulations of first-principles lattice QCD is
the bedrock for the crossover from the hadronic state to a quark-gluon plasma
Borsanyi:2013bia ; HotQCD:2014kol ; HotQCD:2018pds . By contrast, the finite-
density state is inaccessible by the Monte-Carlo (MC) algorithm due to the
sign problem, though there are a few useful exceptions and one can perform the
finite-density simulations in lattice QCD Alford:1998sd ; Son:2000xc ;
Iida:2022hyy ; Brandt:2022hwy ; Abbott:2023coj ; Fujimoto:2023unl .

Neutron stars (NSs) by far provide the most reliable and robust constraints on
the dense-matter EoS in the region, $n_{B}\gtrsim 2n_{0}$, where $n_{B}$ is
the baryon density and $n_{0}\simeq 0.16~{}\text{fm}^{-3}$ is the nuclear
saturation density. The NS observables include relationships between the mass
$M$, the radius $R$, the tidal deformability $\Lambda$, etc., while the EoS is
expressed as a relation between the pressure, $p$, and the energy density,
$\varepsilon$. The one-to-one correspondence between the NS observables and
the EoS via the general relativistic equation called the TOV (Tolman-
Oppenheimer-Volkoff) equation Tolman:1939jz ; Oppenheimer:1939ne ;
1992ApJ…398..569L  can constrain possible EoS candidates. Several studies have
argued that it is feasible to anchor the EoS in the regions from $n_{B}\simeq
n_{0}$ up to $\sim 2n_{0}$ by microscopic calculations based on the chiral
effective field theory Hebeler:2010jx ; Hebeler:2013nza ; Gandolfi:2019zpj ;
Drischler:2020hwi ; Drischler:2020fvz ; Drischler:2021bup . One could also fix
the very high-density part of the EoS around $n_{B}\gtrsim 40n_{0}$ by the
perturbative QCD (pQCD) calculations Freedman:1976ub ; Baluni:1977ms ;
Vuorinen:2003fs ; Kurkela:2009gj ; Gorda:2018gpy ; Gorda:2021znl ;
Gorda:2021kme . With these setups, from the current lower/upper bounds for the
heaviest known $M$’s Demorest:2010bx ; Antoniadis:2013pzd ; Fonseca:2016tux ;
NANOGrav:2017wvv ; NANOGrav:2019jur ; Fonseca:2021wxt , $R$’s Riley:2019yda ;
Miller:2019cac ; Miller:2021qha ; Riley:2021pdl , and $\Lambda$ in the merging
binary NSs LIGOScientific:2017vwq ; LIGOScientific:2018cki ;
LIGOScientific:2018hze , a large portion of the EoS has been strictly excluded
Kurkela:2014vha ; Annala:2017llu ; Annala:2019puf ; Drischler:2020fvz ;
Annala:2021gom ; Drischler:2021bup .

To take one step further to obtain a likely candidate, one can perform
statistical inference of the EoS from the observed NS data incorporating the
uncertainty of astronomical observations by using various approaches. The
standard framework used in the community is Bayesian inference with a
parametrized EoS developed in refs. Steiner:2010fz ; Steiner:2012xt ;
Ozel:2015fia . In a similar but more direct method discussed in earlier papers
Ozel:2009da ; Ozel:2010fw , the probability distribution of the observables is
directly transformed into the probability of the EoS by the MC integration
with a Jacobian marginalizing over the observables. These Bayesian approaches
with a parametrized EoS have been extended to a non-parametric formulation
employing the Gaussian process, as introduced in ref. Landry:2018prl . Machine
learning (ML) using a deep neural network (NN) is another independent stream
of research initiated in the series of works Fujimoto:2017cdo ;
Fujimoto:2019hxv ; Fujimoto:2021zas  by three of the present authors. The
problem of the EoS inference is attracting a growing interest in the nuclear
astrophysics community, and different methods are beneficial for improving
systematics. Indeed there are many papers published even within the past few
years along the lines of these approaches; an incomplete list of such studies
includes refs. Alvarez-Castillo:2016oln ; Raaijmakers:2019dks ; Greif:2020pju
; Al-Mamun:2020vzu ; Huth:2021bsp ; Altiparmak:2022bke ; Lim:2022fap ;
Chimanski:2022wzi ; Marczenko:2022jhl ; Brandes:2022nxa ; Ecker:2022dlg ;
Jiang:2022tps ; Takatsy:2023xzf ; Annala:2023cwx ; Brandes:2023hma ;
Musolino:2023edi ; Richter:2023zec ; Providencia:2023rxc ; Cai:2023pkt ;
Pang:2023dqj ; Fan:2023spm ; Tsang:2023vhh ; Tang:2023owf ; Komoltsev:2023zor
; Imam:2024gfh  for Bayesian inference, refs. Essick:2019ldf ; Landry:2020vaw
; Essick:2020flb ; Essick:2020ghc ; Essick:2021kjb ; Legred:2021hdx ;
Essick:2021ezp ; Legred:2022pyp ; Gorda:2022jvk ; Annala:2023cwx ;
Essick:2023fso ; Mroczek:2023zxo ; Komoltsev:2023zor  for Gaussian process,
and refs. Ferreira:2019bny ; Morawski:2020izm ; Traversi:2020dho ;
Krastev:2021reh ; Soma:2022qnv ; Thete:2022eif ; Farrell:2022lfd ;
Soma:2022vbb ; Ferreira:2022nwh ; Goncalves:2022smd ; Chatterjee:2023ecc ;
Krastev:2023fnh ; Farrell:2023ojk ; Zhou:2023cfs ; Carvalho:2023ele ;
Carvalho:2024kgf  for machine learning.

The results in most literature suggest that the EoS exhibits a rapid
stiffening at a few times $n_{0}$, i.e. a sharp increase in $p$ with
increasing $\varepsilon$. The stiffness of the EoS is characterized by the
squared sound speed defined by the slope of the EoS, i.e., $c_{s}^{2}:={\rm
d}p/{\rm d}\varepsilon$, and it approaches $c_{s}^{2}\to 1/3$ in the conformal
limit at asymptotically high density. At a few times $n_{0}$, $c_{s}^{2}$ may
exceed $1/3$ such that the EoS is stiff enough to account for the existence of
the heavy pulsars Bedaque:2014sqa ; Alford:2015dpa ; Tews:2018kmu . It should
be noted, however, that $c_{s}^{2}>1/3$ can be consistent with conformality
that should be better characterized by the vanishing conformal anomaly
$\varepsilon-3p$ as $c_{s}^{2}$ is a derivative quantity Fujimoto:2022ohj .
Some of the literature claims that this rapid stiffening is followed by
softening, resulting in a peak structure in $c_{s}^{2}$ as a function of
$\varepsilon$ Fujimoto:2019hxv ; Fujimoto:2021zas ; Altiparmak:2022bke ;
Gorda:2022jvk . Anchoring the pQCD results at $n_{B}\simeq 40n_{0}$ will help
reduce $c_{s}^{2}$ at far higher densities out of the range of NS matter;
whether pQCD significantly impacts the EoS in the NS density range is still
under intense debates Komoltsev:2021jzg ; Somasundaram:2022ztm ; Gorda:2022jvk
; Brandes:2023hma ; Mroczek:2023zxo ; Komoltsev:2023zor . Such a peak
structure in $c_{s}^{2}$ is generally expected if the EoS is bridged between
the hadronic and the quark regions without a sharp phase transition
Masuda:2012ed ; Masuda:2012kf ; Kojo:2014rca ; McLerran:2007qj ;
McLerran:2018hbz ; Jeong:2019lhv ; Kojo:2021ugu ; Fujimoto:2022xhv ;
Fujimoto:2023ioi , while purely nucleonic descriptions with the repulsive
three-body forces Akmal:1998cf ; Hebeler:2009iv ; Gandolfi:2011xu  or the
vector mean field Serot:1997xg  tend to increase stiffness monotonically with
the density. Thus, a peak structure in $c_{s}^{2}$ may signal the onset of
quark degrees of freedom in the hadron-quark crossover picture. If the peak
structure is confirmed from the NS data without any theoretical bias, it
strongly signifies the nontrivial change in underlying physics (see also ref.
Tan:2021nat  for possible detection through the binary Love relation).

Previously, we reported a rapidly stiffening EoS using the deep ML inference
Fujimoto:2019hxv ; Fujimoto:2021zas . It should be noted that in what follows
we use “machine learning” or ML to refer to our method specifically in terms
of the NNs, as explained later, though Bayesian analysis can also be regarded
as an ML technology in a general sense. The sound speed from our previous ML
method rises to $c_{s}^{2}\simeq 0.61\pm 0.16$ (with 68% credence) around
$\varepsilon/\varepsilon_{0}\simeq 4$–$5$, where $\varepsilon_{0}$ is the
energy density at $n_{B}=n_{0}$. The high-density part around
$\varepsilon/\varepsilon_{0}\simeq 6$–$8$ reads $c_{s}^{2}\simeq 0.55\pm
0.10$, i.e., $c_{s}^{2}$ was slightly decreased as in Fig. 4 of ref.
Fujimoto:2019hxv . One may thus identify a peak structure in $c_{s}^{2}$ as a
function of $\varepsilon$, although the statistical significance was not
enough. Going into more detail, we first developed the method in ref.
Fujimoto:2017cdo  and extended it in ref. Fujimoto:2019hxv ; we trained an NN
with a sufficiently large dataset of mock observations increased by a few
orders of magnitude with data augmentation. For a given input of $(M,R)$ pairs
of the NSs and their standard deviations, the NN outputs a set of parameters
of the EoS, for which we employ a piecewise $c_{s}^{2}$. Since there was no
established method of uncertainty quantification for deep ML regression, the
treatment of the input and the output uncertainties in our previous analysis
assumed ad hoc simplification. Regarding the input uncertainties in the NS
data, we extracted the standard deviations of the probability distributions of
$(M,R)$ along the $M$ and $R$ directions by marginalizing the distributions.
In effect, we made the Gaussian approximation and used the standard deviations
as the uncertainty quantifying parameters in the NN input. For the output
uncertainties, we attempted two different methods. One is based on the
comparison between the validation data for which we know the ground truth, and
the other is based on the ensemble learning method, i.e., the method called
bagging specifically. We also recognized an issue with the high-density part
of the inferred EoS being not well-constrained by the data at the time of
refs. Fujimoto:2019hxv ; Fujimoto:2021zas .

The purpose of this paper is to verify the nontrivial behavior of $c_{s}^{2}$
in our previous work by using an alternative method based on the convoluted
analysis of the ML reconstruction of the EoS and the MC integration over the
NN input space. The present analysis can be regarded as complementary to our
previous ones, and on top of that, the novel facets are found as follows: (a)
the MC integration handles the full shape of the probability distribution
including the input uncertainty, and accordingly, it quantifies the output
uncertainty, (b) we take into account the new NS data including the heavier
ones than those used in our previous studies, and (c) we numerically check the
effects of our prior assumptions on the inferred results. In the following
paragraphs, we will give further background to the above three points in order
and motivate the present analysis. We shall elaborate the technical details in
the following sections.

First, we make full use of the entire information of the NS data, instead of
using only the standard deviations of the marginalized distributions. For
this, we incorporate the original shape of the joint probability distribution
of $M$ and $R$ into the analysis by the MC integration. Each $(M,R)$ is
sampled for the NN input and thus weighted by the probability; through the MC
integration over $(M,R)$-space, the uncertainties of the input data flow to
the output data in the NN, that is, the EoS parameters in our setup. In this
way, our method proposed in this work is an alternative formulation of the
uncertainty quantification for the ML inference. In particular, we shall
numerically confirm whether the input uncertainty propagates to the output
uncertainty as it should in the proper uncertainty quantification. We note
that this way of uncertainty quantification is somehow similar to that in
refs. Ozel:2009da ; Ozel:2010fw , which additionally requires the explicit
evaluation of the Jacobian of the variable transformation from $(M,R)$ to the
EoS. Treating the original shape of the probability distribution properly is
demanded also by the fact that the joint probability distribution of $M$ and
$R$ is significantly skewed. In some astronomical measurements, moreover, the
correlation between $M$ and $R$ is essential, e.g. the data from NICER
(Neutron star Interior Composition Explorer) has a positive $R$-$M$
correlation because of the measurement of the compactness $M/R$. In our
previous prescription, we discarded the information on the $R$-$M$
correlation, as we dropped the finite skewness of the probability
distributions. These points comprise the main results of this work.

Secondly, as a straightforward extension, we bring our NS data up to date. A
major update since the last work is the radius measurement of heavy mass
neutron star J0740+6620 from NICER Miller:2021qha ; Riley:2021pdl .
Integrating the heavy NS data significantly enriches the information content
in the high-density part of the EoS, which was lacking in our previous works.
We also take the heavy mass pulsars explicitly into account to inform the NN
of the existence of such objects and incorporate the $M$-$R$ probability
distribution of binary NSs in the GW170817 event which is inferred from the
EoS-insensitive relation LIGOScientific:2018hze . The data we used in ref.
Fujimoto:2019hxv  is the $M$-$R$ relation of NSs in photospheric radius
expansion X-ray bursts (which we abbreviate as the PREs) and in the low-mass
X-ray binaries in quiescent period (which we abbreviate as the qLMXBs)
Ozel:2015fia ; Bogdanov:2016nle . These sets of data are obtained from
_spectral_ measurements of X-rays and are known to contain a systematic
uncertainty associated with the atmospheric composition of the star (see e.g.
refs. Miller:2013tca ; Miller:2016pom  for a review), while the NICER performs
_timing_ measurements of X-rays, so it is thought to be free from such
systematic errors. Thus, including a different source of data is in general
beneficial in the viewpoint of controlling the systematic errors in the
predicted EoS.

Finally, as our extended analysis, we shall examine and intentionally adjust
several prior assumptions in the training dataset generation that may
potentially affect the inferred results. In our previous paper
Fujimoto:2021zas , we have found that the inferred results in the density
range unconstrained by the input data depend on the prior values; it is partly
responsible for decreasing behavior in $c_{s}^{2}$ at high densities, so it is
crucial to check the dependence of the results on the prior values in our new
analysis here. In particular, we will consider three major factors: (a) the
EoS parametrization, (b) the $M$-$R$ data sampling in the input layer in the
training dataset, and (c) the data augmentation. We will elaborate procedures
in secs. 3.2 and 3.3, and we will discuss related issues in section 5.

The paper is organized as follows. In section 2, we review our previous method
and its subtleties, and then conceptually introduce our new approach. In
section 3, we elaborate on the technicalities of the new method of mapping the
$M$-$R$ data to the EoS using the NN and the MC integration. Section 4
presents our main results. In section 5, we perform numerical assessments for
the validity of our results including the prior dependence checks in our NN
model. Throughout the paper, we use the natural unit: $c=G=1$.

##  2 Preview of the Present Work

We shall overview our idea in advance to technical details about concrete
analyses and results. For the scientific application, it is crucial how to
treat the uncertainties correctly for the input and output data. In physics,
any experimentally observed data carries an uncertainty in principle, which is
usually estimated by the standard error $\sigma$, namely the standard
deviation of the error. For the inference problem from the input of the
observational data to the output of the theoretical parameters, as we address
in this work, we need to assess the standard errors in the output in response
to the given uncertainty in the input. In this section, we will revisit the
previously adopted method and propose a new approach for the error estimate.

###  2.1 Previous method

In our previous works Fujimoto:2017cdo ; Fujimoto:2019hxv ; Fujimoto:2021zas ,
we proposed a method for the error estimate based on the bagging technique. In
this technique randomly chosen subsets of the training data are used to build
parallel models for the purpose of independent training processes. It should
be noted that we use the term, bagging, in a loose sense as remarked in ref.
Fujimoto:2021zas . Then, the final output is given as a function (e.g., mean)
of all the outputs from the parallel models. This bagging technique, often
referred to as bootstrap aggregating (see figure 1), is useful particularly
for coping with noisy data. At the same time, we can characterize fluctuations
associated with the parallel models by measuring the variance or the standard
deviation of the outputs. Figure 1 schematically summarizes this idea. If the
physical results are reasonably constrained, we should anticipate that the
outputs from these parallel models should be consistent with each other, which
leads to small fluctuations and thus a small value of the standard deviation.
So far, this argument to identify the standard deviation in the bagging as the
error bar is plausible; however, we must be careful of the interpretation as
we discuss in what follows here.

Figure 1: Previously adopted method: Bagging to estimate the uncertainty using
the fluctuation from the parallel models.

In addition, it is more nontrivial how to take account of the error bar in the
input side of observational data in a way that is desirable for the physics
problem. The natural tendency of the estimated error should be as follows: the
more uncertain the input observational data is, the larger error bars the
output theoretical parameters should have. Below, we argue that it is not
obvious whether this expected tendency of the error is satisfied or not.

###  2.2 Subtleties in the previous method and a sketch of the new method

Figure 2: Generation of the training data in the previous works. The set of
$(M_{i},R_{i})$ is sampled on the $M$-$R$ curve obtained from the EoS with
random errors $(\sigma_{M,i},\sigma_{R,i})$ corresponding to the observation
errors. Then, randomly shifted data is generated to train the model so that
the model can learn that the input data may have deviation within
$(\sigma_{M,i},\sigma_{R,i})$.

To clarify our motivation in this work, we shall explain the method we adopted
in the previous works Fujimoto:2017cdo ; Fujimoto:2019hxv ; Fujimoto:2021zas
more closely. Figure 2 shows how to take account of the observational
uncertainty in preparing the training data for supervised learning. First, we
generate many EoSs randomly. We will explain the parametrization of the EoS in
the next subsection. Then, we solve the TOV equation to obtain a family of the
NS mass $M$ and the radius $R$ as a function of the central pressure
$p_{\mathrm{c}}$, which is called the $M$-$R$ curve. Next, we randomly sample
$(M_{i},R_{i})$ on the $M$-$R$ curve as well as the errors
$(\sigma_{M,\,i},\sigma_{R,\,i})$ as sketched in the left panel of figure 2.
In reality, the data for each NS may have a far more non-trivial distribution
than the Gaussian distribution, and we denote the real data distribution by
$\mathcal{P}_{\mathrm{obs}}$. The ML model should “know” that the given data
is subject to deviation from the genuine $M$-$R$ relation according to
$\mathcal{P}_{\mathrm{obs}}$. Thus, as depicted in the right panel of figure
2, for each sampled point we augment the data by shifting $(M_{i},R_{i})$ to
$(M_{i}+\Delta M_{i},R_{i}+\Delta R_{i})$ in accord with the Gaussian
distribution with the standard deviations, $(\sigma_{M,\,i},\sigma_{R,\,i})$.
In this way, we make noised copies of the data so that the ML model can be
tolerant of uncertainty in the input and correctly guess the answer avoiding
the over-fitting problem. The error bar in the output is estimated by the
bagging as schematically illustrated in figure 1.

Let us reiterate the above strategy in a symbolic way. Here, we assume that
there are $n_{\mathrm{NS}}(\in\mathbb{N})$ observed NSs, and we can generalize
the mass and the radius as “properties” of each observation specified by
$P(\in\mathbb{N})$ parameters such as the mass, the radius, the tidal
deformability, the moment of inertia, etc. Then, the observation data has
$n_{X}:=n_{\mathrm{NS}}P$ parameters in total and is collectively expressed as
$\bm{X}=\\{(X_{1}^{(i)},\dots,X_{P}^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}\in\mathbb{R}^{n_{X}}$,
where $X_{p}^{(i)}$ is the $p$-th property of the $i$-th NS.111Generally
speaking, we can also take $\bm{X}$ as spectra from the actual telescope
observation as discussed in refs. Farrell:2022lfd ; Farrell:2023ojk . In this
case, however, the relation between the telescope spectra and the EoS is not
as clear-cut as the relation between the $M$-$R$ curve and the EoS that is
firmly based on general relativity. Therefore, one needs to introduce the
nuisance parameters in the former case and marginalization over these
parameters. In the present work, specifically, $(X_{1}^{(i)},X_{2}^{(i)})$
stands for the mass and the radius of the $i$-th NS, i.e.,
$(M^{(i)},\,R^{(i)})$. As already discussed at the beginning of this section,
astronomical observations do not directly give a point $\bm{X}$ because of
observational uncertainties. Thus, the input from the observation is
inherently represented by a joint probability density
$\mathcal{P}_{\mathrm{obs}}(\bm{X})$ over $\mathbb{R}^{n_{X}}$. In our
previous works Fujimoto:2019hxv ; Fujimoto:2021zas , we tried to characterize
$\mathcal{P}_{\mathrm{obs}}$ by the mean
$\\{(\mu_{M}^{(i)},\mu_{R}^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}$ and the width
$\\{(\sigma_{M}^{(i)},\sigma_{R}^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}$. We then
proposed to build an NN model, $\Psi_{\mathrm{prev}}$, that directly predicts
a likely EoS:

$\displaystyle\begin{array}[]{cccc}\Psi_{\mathrm{prev}}:&\mathbb{R}^{n_{X}}\times\mathbb{R}^{n_{X}}&\to&\mathcal{F}(\mathbb{R})\\\
&\rotatebox[origin={c}]{90.0}{$\in$}&&\rotatebox[origin={c}]{90.0}{$\in$}\\\
&(\bm{\mu}_{\bm{X}},\bm{\sigma}_{\bm{X}})&\quad\mapsto&\bm{Y}\end{array}$ (4)

Here, the output EoS of our central interest is a relation, $p(\varepsilon)$,
between the pressure $p$ and the energy density $\varepsilon$, which is
collectively denoted by $\bm{Y}=p(\cdot)\in\mathcal{F}(\mathbb{R})$, where
$\mathcal{F}(\Omega)$ denotes the function space $\\{\Omega\to\mathbb{R}\\}$.

There are, however, several subtle points in this approach. The observational
data is not necessarily a simple product of the Gaussian distributions
characterized by $(\sigma_{M,\,i},\sigma_{R,\,i})$ only. The data for each NS
could have a rather complicated distribution as shown by the shaded blobs in
the top-left panel of figure 3, which schematically represents
$\mathcal{P}_{\mathrm{obs}}$. The problem is how to capture arbitrary shapes
of the distribution $\mathcal{P}_{\mathrm{obs}}$. As sketched in figure 3 (and
figure 4 for the real data), the distribution of the observational data may be
significantly distorted and could not be well approximated by
$(\sigma_{M,\,i},\sigma_{R,\,i})$.

There is another subtle problem, i.e., the uncertainty quantification in the
final output. To understand the problem better, let us assume a hypothetical
parameter that is not constrained by the observational data at all. If such a
parameter exists, the uncertainty must be large. In the actual estimate,
however, such an unconstrained parameter could be dominantly chosen by the
prior distribution, which is essentially the distribution in the training
dataset in the present context. For a fixed prior distribution, consequently,
the standard deviation of the unconstrained parameter may happen to be small.
Roughly speaking, the undetermined output parameter may take close values that
are far from the true answer.

The present paper focuses on the problem of the $\mathcal{P}_{\mathrm{obs}}$
treatment. In contrast, it is generally difficult to overcome the latter
problem, and we should judiciously scrutinize the prior dependence. In this
paper, we will take a close look at different properties of the estimated
uncertainties and defer the systematic investigation of the prior dependence
to the future study. Nonetheless, we will present exploratory analyses on the
sensitivity of our results to selected prior choices.

Figure 3: Our proposed procedures to deal with the error bars. The
observational data is supposed to have a distribution
$\mathcal{P}_{\mathrm{obs}}$ which spans on the $M$-$R$ plane in the present
case (top-left). Sets of data points are randomly sampled according to
$\mathcal{P}_{\mathrm{obs}}$ (top-right). The ML model is trained to infer
sets of the theoretical parameters, i.e., the sound velocities at various
densities in the present case (bottom-right). Finally, the ensemble average is
taken with the weight $\mathcal{P}_{\mathrm{obs}}$ on the input side, leading
to the estimate for the probability distribution $\mathcal{P}_{\mathrm{EoS}}$
on the output side (bottom-left).

Figure 3 is a cartoon to explain the new procedures we propose in the present
work. It is not our intention to claim that the present method supersedes the
previous one. Rather, we aim to establish an independent estimate for treating
the error bars, which would be useful for the reliability check. Here, for
convenience, we use the terminology in a specific field of physics, but the
problem of uncertainty quantification is ubiquitous, and it can be applied to
a broad spectrum of problems. In the top-left panel of figure 3,
$\mathcal{P}_{\mathrm{obs}}$ represents the probability distribution of the
observational data. Given an ML model from the observational $M$-$R$ points to
the theoretical EoS points, let us carry out the MC sampling of the $M$-$R$
data points following $\mathcal{P}_{\mathrm{obs}}$. For the points sampled on
the $M$-$R$ plane as in the top-right panel of figure 3, we can deduce the EoS
parameters, which are depicted in the bottom-right panel of figure 3. Since
the $M$-$R$ points are randomly sampled, we should take the ensemble average
in the end, which is nothing but making the convolution with
$\mathcal{P}_{\mathrm{obs}}$. We can finally reach the probability
distribution of the EoS parameters denoted by $\mathcal{P}_{\mathrm{EoS}}$,
from which we can immediately read out the mean and the variance.

###  2.3 Strategy of the present method

We shall go into the details of the present method using the symbolic notation
as used in the previous subsection. Our goal is to make constraints on the NS
EoS using the observations of static NS properties, for which we adopt only
the mass $M$ and the radius $R$ in this work. The problem is that the exact
$\bm{Y}$ cannot be identified by the observational input because the
observation data takes a distribution $\mathcal{P}_{\mathrm{obs}}(\bm{X})$,
and thus the corresponding EoS should naturally be given by a distribution as
well. Even if the observational data comes along the genuine $M$-$R$ line
without uncertainties, in principle, the observation of a finite number of the
NSs is generally insufficient to pin down the exact EoS in an infinite-
dimensional function space. Therefore, we should consider the constraint to
the NS EoS given in the form of a probability distribution,
$\mathcal{P}_{\mathrm{EoS}}(\bm{Y})$. In the general term, we may express our
ultimate goal as a regression $\mathscr{R}_{\mathrm{inv}}$ of “a probability
distribution of EoS” $\mathcal{P}_{\mathrm{EoS}}$ on “a probability
distribution of the NS properties” $\mathcal{P}_{\mathrm{obs}}$:

$\displaystyle\begin{array}[]{cccc}\mathscr{R}_{\mathrm{inv}}:&\mathcal{F}(\mathbb{R}^{n_{X}})&\to&\mathcal{F}(\mathcal{F}(\mathbb{R}))\\\
&\rotatebox[origin={c}]{90.0}{$\in$}&&\rotatebox[origin={c}]{90.0}{$\in$}\\\
&\mathcal{P}_{\mathrm{obs}}(\cdot)&\quad\mapsto&\mathcal{P}_{\mathrm{EoS}}(\cdot)\end{array}$
(8)

One general and well-defined framework to perform this is the Bayesian
analysis based on the MC integration, but it would be numerically demanding if
we test the method in various observational scenarios and try to understand
the characteristics of the method itself. As an alternative approach, it would
be interesting to think about a simpler ML-based setup.

In this work, we examine a deep ML approach. Our central idea is to render the
problem (8) into the pointwise regression with a supervised ML:

$\displaystyle\begin{array}[]{cccc}\Psi_{\mathrm{inv}}:&\bm{X}&\quad\mapsto&\bm{Y}\end{array}$
(10)

In practice, because of the limitation of computational power, we shall
restrict $\mathcal{F}(\mathbb{R})$ to an $N$-dimensional EoS subspace that
spans over the parametrized form of $p(\varepsilon)$ with the EoS parameters
$\bm{Y}\in\mathbb{R}^{N}$. We represent $\Psi_{\mathrm{inv}}$ by an NN model
and train the model using a dataset of (input, output) pairs,
$(\bm{X},\bm{Y})$, for supervised learning. Once $\Psi_{\mathrm{inv}}$ is
obtained, we introduce the following Ansatz for the EoS distribution to model
$\mathcal{P}_{\mathrm{EoS}}(\bm{Y})\simeq\mathcal{P}_{\mathrm{EoS}}^{\prime}(\bm{Y})$:

$\displaystyle\mathscr{R}^{\prime}_{\mathrm{inv}}:\quad\mathcal{P}_{\mathrm{EoS}}^{\prime}(\bm{Y})$
$\displaystyle=\int\delta(\bm{Y}-\Psi_{\mathrm{inv}}(\bm{X}))\mathcal{P}_{\mathrm{obs}}(\bm{X}){\rm
d}\bm{X}\,.$ (11)

Using the distribution (11), we will discuss the most likely EoS and
corresponding uncertainties. Equation (11) is applicable for arbitrary shapes
of the observational distribution $\mathcal{P}_{\mathrm{obs}}$ without
updating the trained model, $\Psi_{\mathrm{inv}}$. Also, we can discuss the
uncertainties of the resulting EoS using
$\mathcal{P}_{\mathrm{EoS}}^{\prime}$.

Figure 4: Distribution of the actual observational data. The $M$-$R$ region
from the inferred EoS distribution (shown by the red band) is overlaid with
the uncertainty estimate based on the method we propose in this work.

For $\mathcal{P}_{\mathrm{obs}}$, we use the collection of the data summarized
in figure 4. These probability distributions are assembled from various
sources Ozel:2015fia ; Bogdanov:2016nle ; Demorest:2010bx ; Antoniadis:2013pzd
; Fonseca:2016tux ; Miller:2019cac ; Miller:2021qha ; LIGOScientific:2018cki .
In figure 4, we plot the contour of the 68% credible region of the posterior
probability density function in Bayesian statistics except for the GW170817,
J1614-2230, and J0348+0432 data. The uncertainties of two NSs in the GW170817
data are intentionally shown by error bars because the distributions of the
two NSs observed in a single binary NS merger event are highly correlated with
each other through the tight constraint on the total mass of the system. The
rectangular bands for the J1614-2230 and J0348+0432 data show that the
$R$-direction is uniformly distributed while the $M$-direction is limited
within the 68% confidence interval. The NSs marked with thick lines in figure
4 are the data not incorporated in the previous analysis Fujimoto:2019hxv  and
newly considered in this work. The red band is the $M$-$R$ curve with
uncertainties corresponding to the EoS predicted by our NN method as shown in
figure 7), which constitutes the main results of this paper. The details will
be elaborated in section 4.

##  3 Mapping $M$-$R$ distributions to EoS using neural network

We summarize the setup of our inference to predict the EoS parameters,
$\bm{Y}$, from the NS properties, $\bm{X}$.

Figure 5: Schematic figure of our strategy to infer the distribution,
$\mathcal{P}_{\mathrm{EoS}}\simeq\mathcal{P}_{\mathrm{EoS}}^{\prime}$ of the
EoS parameters, $\bm{Y}={\bm{c}_{s}^{2}}$ (and $p_{c}$), from the experimental
distribution, $\mathcal{P}_{\mathrm{obs}}(\bm{X})$.

###  3.1 Synopsis of our method

Our strategy is illustrated in figure 5, and we describe our method along this
schematic figure in this subsection.

We consider
$\bm{X}=\\{(M^{(i)},R^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}\in\mathbb{R}^{2n_{\mathrm{NS}}}$
for the NS properties, where $M^{(i)}$ and $R^{(i)}$ are the mass and the
radius, respectively, of the $i$th NS. For the EoS parameters $\bm{Y}$, we
consider a piecewise polytropic form of the EoS and represent the EoS with a
set of the speed of sound,
$\bm{Y}:=\bm{c_{s}^{2}}=\\{c_{s,n}^{2}\\}_{n=1}^{N}$ (see section 3.2 for
detail). The joint probability distribution
$\mathcal{P}_{\mathrm{obs}}(\bm{X})$ is specified by the observational data
itemized in section 3.5.

We first generate the dataset,
$\mathfrak{D}=\\{(\bm{X}_{i},\bm{Y}_{i})\\}_{i=1}^{I}$, for supervised
learning (see the yellow box on the top left corner in figure 5). We randomly
sample EoS parameters $\bm{c_{s}^{2}}=\bm{Y}$. The corresponding NS properties
$\bm{X}$ are obtained by adding noises to the solutions of the TOV equation
$\bm{X}_{\ast}$ (the red arrow on the top left corner). The noises are sampled
from the normal distribution
$\mathcal{N}(0,(\sigma_{M}^{(i)})^{2})\times\mathcal{N}(0,(\sigma_{R}^{(i)})^{2})$,
and the standard deviations $\\{(\sigma_{M}^{(i)},\sigma_{R}^{(i)})\\}$ are
determined so as to match those of $\mathcal{P}_{\mathrm{obs}}(\bm{X})$ (the
bottom left corner of figure 5). Further details on the generation of the
dataset $\mathfrak{D}$ will be elaborated in section 3.3. We then design the
NN architecture of $\Psi_{\mathrm{inv}}$ (explained in section 3.4) and train
it with the dataset $\mathfrak{D}$ (the green arrow in figure 5). We finally
analyze the “posterior” distribution $\mathcal{P}_{\mathrm{EoS}}^{\prime}$
obtained by the MC integration over the $\bm{X}$-space convoluted with
$\mathcal{P}_{\mathrm{obs}}(\bm{X})$ as given in eq. (11) (the bottom right
corner of figure 5).

To infer a “likely” EoS using $\mathcal{P}_{\mathrm{EoS}}^{\prime}$, we
consider the Bayes estimator of the EoS parameters:
$\bm{Y}^{*}:=\langle\bm{Y}\rangle_{\mathcal{P}_{\mathrm{EoS}}^{\prime}}$,
where $\langle\cdot\rangle_{\mathcal{P}_{\mathrm{EoS}}^{\prime}}$ denotes the
average by the distribution $\mathcal{P}_{\mathrm{EoS}}^{\prime}$. We then
estimate its uncertainties by the covariances
$\langle(\bm{Y}-\bm{Y}^{*})(\bm{Y}-\bm{Y}^{*})^{\mathrm{T}}\rangle_{\mathcal{P}_{\mathrm{EoS}}^{\prime}}$.
Instead of explicitly constructing $\mathcal{P}_{\mathrm{EoS}}^{\prime}$, we
perform the MC integration to directly evaluate the averages
$\langle\cdot\rangle_{\mathrm{\mathcal{P}_{\mathrm{EoS}}^{\prime}}}$:

$\displaystyle\langle\mathcal{O}[\bm{Y}]\rangle_{\mathcal{P}_{\mathrm{EoS}}^{\prime}}$
$\displaystyle=\int\mathcal{O}[\Psi_{\mathrm{inv}}(\bm{X})]\mathcal{P}_{\mathrm{obs}}(\bm{X}){\rm
d}\bm{X}\simeq\frac{1}{n_{\mathrm{MC}}}\sum_{i=1}^{n_{\mathrm{MC}}}\mathcal{O}[\Psi_{\mathrm{inv}}(\bm{X}_{i})]\,,$
(12)

where $\mathcal{O}[\bm{Y}]$ is a quantity written in terms of the EoS
parameters $\bm{Y}$. Also, $\bm{X}_{i}$ is the $i$th MC sampling point
following the distribution $\mathcal{P}_{\mathrm{obs}}(\bm{X})$, and
$n_{\mathrm{MC}}$ is the number of the sample points.

###  3.2 EoS parametrization

In principle, an EoS can be defined through the speed of sound
$c_{s}^{2}(\varepsilon)=\partial p/\partial\varepsilon$ as
$p(\varepsilon)=\int_{0}^{\varepsilon}c_{s}^{2}(\varepsilon^{\prime})d\varepsilon^{\prime}$.
In this study, we employ the piecewise polytropic expression for an EoS, and
it is parametrized by a set the speed of sound,
$\bm{c_{s}^{2}}=\\{c_{s,n}^{2}\\}_{n=1}^{N}$. We split the range of considered
density $[\varepsilon_{0},\varepsilon_{N}]$ into $N$ ($\in\mathbb{N}$)
segments, i.e., $[\varepsilon_{n-1},\varepsilon_{n}]$ ($n=1,\dots,N$), and
specify a value of the average speed of sound:
$c_{s,n}^{2}=\int_{\varepsilon_{n-1}}^{\varepsilon_{n}}c_{s}^{2}(\varepsilon^{\prime})d\varepsilon^{\prime}/(\varepsilon_{n}-\varepsilon_{n-1})$
in each segment. The density range is equally partitioned into segments in the
logarithmic scale as

$\displaystyle\ln\varepsilon_{n}:=\ln\varepsilon_{0}+\frac{n}{N}(\ln\varepsilon_{N}-\ln\varepsilon_{0})\qquad(n=1,\dots,N)\,,$
(13)

where we chose $\varepsilon_{0}=150~{}\text{MeV}/\text{fm}^{3}$ to approximate
the energy density at the nuclear saturation density $n_{B}=n_{0}$, and we set
$\varepsilon_{N}:=8\varepsilon_{0}$. For instance, a choice of $N=5$ leads to
$(\varepsilon_{0},\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4},\varepsilon_{5})\simeq(1.0,1.5,2.3,3.5,5.2,8.0)\varepsilon_{0}$.
The low-density EoS, $p(\varepsilon<\varepsilon_{0})$, is constrained by
existing models of nuclear matter, for which we employ SLy4 Douchin:2001sv .
The pressure $p_{n}$ at the segment boundary $\varepsilon_{n}$ is obtained as

$\displaystyle p_{n}$
$\displaystyle:=p_{n-1}+\int_{\varepsilon_{n-1}}^{\varepsilon_{n}}c_{s}^{2}(\varepsilon^{\prime})d\varepsilon^{\prime}$
$\displaystyle=p_{n-1}+c^{2}_{s,n}(\varepsilon_{n}-\varepsilon_{n-1})\qquad(n=1,\dots,N)\,,$
(14)

where $p_{0}=p(\varepsilon_{0})$ is fixed by SLy4. The EoS inside each segment
is interpolated by a polytrope function, i.e.,
$p=p_{\mathrm{pol}}(\varepsilon)=K_{n}\varepsilon^{\Gamma_{n}}$ for
$\varepsilon_{n-1}<\varepsilon<\varepsilon_{n}$, where $(\Gamma_{n},K_{n})$ is
determined by the boundary conditions,
$p_{n-1}=p_{\mathrm{pol}}(\varepsilon_{n-1})$ and
$p_{n}=p_{\mathrm{pol}}(\varepsilon_{n})$, which yields:

$\displaystyle\Gamma_{n}=\frac{\ln p_{n}-\ln
p_{n-1}}{\ln\varepsilon_{n}-\ln\varepsilon_{n-1}}\,,\qquad K_{n}$
$\displaystyle=p_{n}^{1/(1-\ln\varepsilon_{n}/\ln\varepsilon_{n-1})}p_{n-1}^{1/(1-\ln\varepsilon_{n-1}/\ln\varepsilon_{n})}\,.$
(15)

In this paper, we present results for $N=5$ unless otherwise specified.

###  3.3 Dataset generation

Here, we describe the procedure to generate the dataset used for training and
validation. Specifically, we take the following steps to prepare our dataset,
which is mostly in line with ref. Fujimoto:2021zas  except for step 4 in the
following:

  1. 1.

We first generate an EoS, $\bm{Y}=\\{c_{s,n}^{2}\\}_{n=1}^{N}$, by randomly
sampling each $c_{s,n}^{2}$ by the uniform distribution whose support is
$[\delta,1-\delta]$ with $\delta=0.01$.

  2. 2.

We then obtain the corresponding $M$-$R$ line by solving the TOV equations
described below. If the maximum mass $M_{\mathrm{max}}$ does not reach
$2.01M_{\odot}$, we reject the EoS. If the EoS is accepted, we proceed to the
next step.

  3. 3.

We randomly sample $n_{\mathrm{NS}}$ NSs,
$\\{(M_{*}^{(i)},R_{*}^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}$, on the $M$-$R$ line
by the uniform distribution with respect to $M_{*}^{(i)}$ in the range
$[M_{\mathrm{min}},M_{\mathrm{max}}]$, where we take
$M_{\mathrm{min}}=1.0M_{\odot}$.

  4. 4.

We sample $n_{\mathrm{NS}}$ observational data points,
$\bm{X}=\\{(M^{(i)},R^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}$, from the Gaussian
distributions:
$(M^{(i)},R^{(i)})\sim\mathcal{N}(M_{*}^{(i)},(\sigma_{M}^{(i)})^{2})\times\mathcal{N}(R_{*}^{(i)},(\sigma_{R}^{(i)})^{2})$.
Here, the standard deviations,
$\\{(\sigma_{M}^{(i)},\sigma_{R}^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}$, are
matched with those of $\mathcal{P}_{\mathrm{obs}}(\bm{X})$.

  5. 5.

A pair of $(\bm{X},\bm{Y})$ forms one data point. We repeat steps 1–4 to
generate many pairs to create a dataset. To increase the dataset size, we may
sample multiple sets of $n_{\mathrm{NS}}$ NSs for a single step in step 3, and
multiple sets of observational points for a single set of $n_{\mathrm{NS}}$
NSs in step 4.

We note that in the previous method Fujimoto:2019hxv ; Fujimoto:2021zas , at
step 4, we had to sample the standard deviations,
$\\{(\sigma_{M}^{(i)},\sigma_{R}^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}$, randomly
by the uniform distributions denoted by $U(0,M_{\odot})$ and
$U(0,5\,\mathrm{km})$. For step 2, $M$ and $R$ of NSs can be obtained from the
EoS through the TOV equation Tolman:1939jz ; Oppenheimer:1939ne . The TOV
equation describes the hydrostatic equilibrium condition of spherically
symmetric configuration in general relativity, which takes the following form:

$\displaystyle\frac{{\rm d}{p}(r)}{{\rm d}{r}}$
$\displaystyle=-\frac{m(r)\varepsilon(p(r))}{r^{2}}\left[1+\frac{p(r)}{\varepsilon(p(r))}\right]\left[1+\frac{4\pi
r^{3}p(r)}{m(r)}\right]\left[1-\frac{2m(r)}{r}\right]^{-1}\,,$ (16)
$\displaystyle\frac{{\rm d}{m}(r)}{{\rm d}{r}}$ $\displaystyle=4\pi
r^{2}\varepsilon(p(r))\,,$ (17)

where $r$ is the radial distance from the NS center, $p(r)$ is the pressure
inside the NS at $r$, and $m(r)$ is the mass enclosed within the sphere of
radius $r$. We can numerically solve eqs. (16) and (17) with the boundary
conditions,

$\displaystyle p(r=0)=p_{c}\,,\qquad m(r=0)=0,$ (18)

using a standard method for ordinary differential equations outward from $r=0$
until the pressure vanishes. The NS radius, $R$, is defined by the condition
of vanishing pressure, and the NS mass, $M$, is defined accordingly as

$\displaystyle R:\ p(R)=0\,,\qquad M:=m(R)\,.$ (19)

Therefore, the resulting pair of $M$ and $R$ depends on the parameter $p_{c}$
in eq. (18). We can draw the $M$-$R$ curve by changing the parameter $p_{c}$.

###  3.4 Neural-network models

To represent the fitting function $\Psi_{\mathrm{inv}}(\bm{X})$, we employ a
fully connected deep neural network. The input
$\bm{X}\in\mathbb{R}^{d_{\mathrm{in}}}$ is normalized for pre-processing:

$\displaystyle\bar{M}^{(i)}:=\frac{M^{(i)}}{M_{\rm
norm}}\,,\qquad\bar{R}^{(i)}:=\frac{R^{(i)}}{R_{\rm
norm}}\qquad(i=1,\dots,n_{\mathrm{NS}})\,,$ (20)

with $M_{\mathrm{norm}}=3M_{\odot}$ and $R_{\mathrm{norm}}=30~{}\text{km}$ in
our choice. The normalized NS properties are plugged into the input layer:
$\bm{x}^{(0)}=\bm{X}$. The NN model iteratively calculates the values of the
subsequent layers as

$\displaystyle
x_{i}^{(l)}=a^{(l)}\biggl{(}\sum_{j=1}^{n_{l-1}}W^{(l)}_{ij}x_{j}^{(l-1)}+b_{i}^{(l)}\biggr{)}\qquad(l=1,\dots,n_{\mathrm{layer}})\,,$
(21)

where $\bm{x}^{(l)}\in\mathbb{R}^{d^{(l)}}$ is the state of the $l$-th layer
with $d^{(l)}$ being the number of neurons in the $l$-th layer. The activation
function for the $l$-th layer is $a^{(l)}(x)$, and $W_{ij}^{(l)}$ and
$b_{j}^{(l)}$ are the trainable parameters. The output in the last layer gives
$\bm{Y}=\bm{x}^{(n_{\mathrm{layer}})}$, which is the EoS parameters,
$\bm{Y}=\bm{c_{s}^{2}}$, in our model. It should be noted that our extended
analysis in section 5.4 will instead adopt
$\bm{Y}=(\bm{c_{s}^{2}},\\{p_{c}^{(i)}\\}_{i=1}^{n_{\mathrm{NS}}})$ to
additionally predict the central pressure $p_{c}^{(i)}$ of the $i$-th NS. For
the motivation of analysis with and without $p_{c}$, we will closely discuss
it later in section 5.4. The NN structure used in this work is summarized in
table 1. The dimensions of the input and output and the number of trainable
parameters of the NN models in our analysis are summarized in table 2, which
are estimated from $\dim[W^{(l)}]=d^{(l-1)}d^{(l)}$ and
$\dim[b^{(l)}]=d^{(l)}$.

Layer |  $l=0$ (Input) | $l=1$ | $l=2$ | $l=3$ | $l=4$ |  $l=n_{\mathrm{layer}}$ (Output)
---|---|---|---|---|---|---
$d^{(l)}$ | $d_{\rm in}$ | $4d_{\rm in}$ | $4d_{\rm in}$ | $2d_{\rm in}$ | $2d_{\rm in}$ | $d_{\rm out}$
$a^{(l)}$ | $-$ | ReLU | ReLU | ReLU | ReLU | Sigmoid
Table 1: Our NN architecture. The constants $d_{\rm in/out}$ represent the dimensions (the number of neurons) of the input/output layer.  |  | $n_{\mathrm{NS}}=14$ | $n_{\mathrm{NS}}=6$ | $n_{\mathrm{NS}}=20$
---|---|---|---|---
Model w/o $p_{c}$ | $(d_{\rm in},\,d_{\rm out})$ | (28, 5) | (12, 5) | (40, 5)
| # params | 25,709 | 4,877 | 52,085
Model with $p_{c}$ | $(d_{\rm in},\,d_{\rm out})$ | (28, 19) | (12, 11) | (40, 25)
| # params | 26,507 | 5,027 | 53,705
Table 2:  The numbers of trainable parameters. In the input layer, $d_{\rm
in}=2n_{\mathrm{NS}}$ corresponds to our choice of $M$ and $R$ for the NS
properties. In the output layer, $d_{\rm out}=N$ for the model without $p_{c}$
and $d=N+n_{\mathrm{NS}}$ for the model with $p_{c}$. We take $N=5$ and
$n_{\mathrm{NS}}=14$, $6$, and $20$ (see the text for details).

For our numerics, we use the Keras package software:Keras  for the NN
construction and employ Adam DBLP:journals/corr/KingmaB14  and msle for the
optimizer and the loss function, respectively. We generate $200,000$
independent $M$-$R$ lines, repeating 20 times of sampling $n_{\mathrm{NS}}$
observational data points (corresponding to step 3 in section 3.3) per one
generated $M$-$R$ line, and then repeat 10 times of augmenting the data points
with Gaussian noises per one set of the data points (corresponding to step 4
in section 3.3). Hence, we prepare the dataset consisting of $200,000\times
20\times 10=40,000,000$ observation data. We split it into $40,000,000\times
0.9=36,000,000$ and $40,000,000\times 0.1=4,000,000$ observation data as the
datasets used for training and validation, respectively. We take $200$ epochs
with the batch size of $100,000$ for the optimization222 We checked that
reducing the batch size to 10,000, for example, does not change the results
within the error bars. .

###  3.5 Observed properties of neutron stars

We explain how to prepare $\mathcal{P}_{\mathrm{obs}}(\bm{X})$ in the present
analysis in what follows in this subsection. All the sources of data are given
in the form of the posterior probability density function in Bayesian
statistics; see figure 4 for the 68% contour of their credible regions.

We use the following 14 sources of spectral measurement in X-ray, which will
be collectively referred to as “PREs/qLMXBs” hereafter. They were also used in
our previous analyses Fujimoto:2019hxv ; Fujimoto:2021zas .

  * •

Photospheric radius expansion X-ray bursts (PREs): 4U 1608-52, 4U 1724-207, KS
1731-260, EXO 1745-248, SAX J1748.9-2021, and 4U 1820-30 Ozel:2015fia
($\mbox{\\# distributions}=6$)

  * •

Quiescent low-mass X-ray binaries (qLMXBs): M13, M28, M30, NGC 6304, NGC 6397,
$\omega$ Cen Ozel:2015fia , X5, and X7 Bogdanov:2016nle  ($\mbox{\\#
distributions}=8$)

On top of them, we further incorporate the following 5 new data into the
analyses in this work333Precisely speaking, the NICER data from J0030+0451 was
already taken into account in our previous analysis Fujimoto:2021zas ., which
will be collectively referred to as “HP/NICER/GW” later.

  * •

Heavy pulsars (HP) from precise radio astronomy: PSR J1614-2230
Demorest:2010bx ; Fonseca:2016tux , and PSR J0348+0451 Antoniadis:2013pzd
($\mbox{\\# distributions}=2$)

  * •

Radii from the timing measurements in the X-ray band in NICER: PSR J0030+0451
Miller:2019cac , and PSR J0740+6620 Miller:2021qha  ($\mbox{\\#
distributions}=2$)

  * •

Gravitational wave event from LIGO-Virgo collaboration (GW): GW170817
LIGOScientific:2018cki
($\mbox{\\# distributions}=1$)

The heavy pulsars are the precise radio measurements of masses of the heaviest
NS known to date except for the black widow pulsar Romani:2022jhd . We did not
adopt the black widow pulsar data because of a larger uncertainty compared to
the two pulsars listed above. We limited the radii of the heavy pulsars within
$8.5$–$16.5~{}\text{km}$. The edge values correspond to the smallest and
largest possible radii Drischler:2020fvz , respectively, extrapolated from the
crust EoS assumed up to half the saturation density.

We define $\mathcal{P}_{\mathrm{obs}}(\bm{X})$ by a product of all those
distributions. It should be noted that two pairs of $(M,R)$ from GW170817 in
the binary system of two NSs are correlated to each other. Hence, the numbers
of NSs for the two categories are $n_{\mbox{\tiny PREs/qLMXBs}}=14$ and
$n_{\mbox{\tiny HP/NICER/GW}}=6$, respectively, and the total number amounts
to $n_{\mathrm{NS}}=20$.

##  4 Results for the Equation of State

Our main results of the EoS, in particular the behavior of $c_{s}^{2}$, are
discussed in this section.

Figure 6: Speed of sound squared, $c_{s,n}^{2}$, as a function of the energy
density. The results inferred from the full, PREs/qLMXBs, and HP/NICER/GW
distributions are shown by the red, blue, and green bars, respectively, which
represent the standard deviation of $\bm{Y}$ over
$\mathcal{P}_{\mathrm{EoS}}^{\prime}$.

Figure 6 shows $c_{s,n}^{2}$ estimated by the full distribution of
$\mathcal{P}_{\mathrm{obs}}$. We also compare it with the results from partial
distributions limited to only the PREs/qLMXBs and the HP/NICER/GW data. The
mean values $\mu[\bm{Y}]$ and the standard deviation $\sigma[\bm{Y}]$ are
obtained by

$\displaystyle\mu[Y_{i}]$ $\displaystyle=\langle
Y_{i}\rangle_{\mathcal{P}_{\mathrm{EoS}}^{\prime}}\,,\qquad\sigma[{Y}_{i}]=\sqrt{\langle(Y_{i}-\mu[Y_{i}])^{2}\rangle_{\mathcal{P}_{\mathrm{EoS}}^{\prime}}}\,.$
(22)

The behavior of PREs/qLMXBs significantly differs from that of HP/NICER/GW. As
seen in figure 6, the PREs/qLMXBs data generates a soft EoS with
$c_{s}^{2}\sim 0.1$ in the low $\varepsilon$ regime and the EoS starts to be
stiff only around the middle segments at $\varepsilon/\varepsilon_{0}\gtrsim
2$, whereas the EoS inferred from the HP/NICER/GW data starts to be stiff
already in the second segment; the latter EoS has a larger sound speed around
$\varepsilon/\varepsilon_{0}\simeq 2$ compared to the typical value, e.g., the
ab initio calculation from the chiral effective theory Drischler:2020fvz ;
Drischler:2020hwi  although they are still consistent within the uncertainty.

The high-density part of the EoS inferred from the NP/NICER/GW data,
particularly the fifth segment, may be unconstrained by the data as it is
apparently set by the average of the prior distribution $c_{s,5}^{2}\simeq
0.5$. Further, the small uncertainty signifies that this may be the biased
output from the NN and does not necessarily reflect the actual statistical
errors. We will come back to this point later in section 5.4 when we include
the central pressure of NSs in the analysis.

Figure 7: EoS with the uncertainty estimate based on the method we propose in
this work. The colored lines and bands are the mean values and the errors,
respectively. The upper and lower lines of the bands are obtained by
substituting $\bm{c_{s}^{2}}=\mu[\bm{Y}]\pm\sigma[\bm{Y}]$ in eqs. (14) and
(15).

In figure 7, we show the pressure $p$ obtained from the constrained
$c_{s,n}^{2}$. The pressure is calculated as a piecewise polytropic form of a
function of $\varepsilon$ by substituting $\bm{c_{s}^{2}}=\mu[\bm{Y}]$ in eqs.
(14) and (15). The EoS marked as “This work” in figure 7 corresponds to
$c_{s,n}^{2}$ marked as “Full” in figure 6. We compared our EoS results with
our previous estimates from ref. Fujimoto:2019hxv . The inferred EoS in this
work exhibits a slightly steeper transition from a soft EoS at low density to
a stiff EoS at middle density compared to that in the previous work.

(a) PREs/qLMXBs (b) HP/NICER/GW
---
(c) Full
Figure 8: Correlation coefficients of $c_{s,n}^{2}$ ($n=1,\dots,5$) for
$n_{\mathrm{NS}}=14$ (PREs/qLMXBs), $6$ (HP/NICER/GW), and $20$ (Full
distributions). The red and blue colored cells correspond to positive and
negative values, respectively. The diagonal components are unity by definition
of eq. (23).

The EoS presented by the thick solid line in figure 7 is based on the point
estimation, but it is just one possibility among various EoSs. As already
discussed in sections 1 and 2, it is crucial to estimate the uncertainty in
constraining physics corresponding to the uncertainty in observational data.
The standard deviation over $\mathcal{P}_{\mathrm{EoS}}^{\prime}$ is presented
by the band in figure 7 as one measure of the uncertainty. Moreover, it should
be noted that those estimated values of $c_{s,n}^{2}$ at different $n$ are not
independent of each other. They can be correlated with each other through the
nontrivial shape of $\mathcal{P}_{\mathrm{EoS}}^{\prime}$ generated by
$\mathcal{P}_{\mathrm{obs}}$. Such correlation affects the shape of the
possible EoSs. As the simplest measure for the correlation among
$c_{s,n}^{2}$’s, we consider the Pearson correlation coefficients:

$\displaystyle r_{ij}$
$\displaystyle:=\frac{\langle(Y_{i}-\mu[Y_{i}])(Y_{j}-\mu[Y_{j}])\rangle_{\mathcal{P}_{\mathrm{EoS}}^{\prime}}}{\sigma[Y_{i}]\sigma[Y_{j}]}\,.$
(23)

Figure 8 shows the correlation coefficients $r_{ij}$ corresponding to figure
6. Panels (a) and (b) in figure 8 show a sharp contrast particularly in
$r_{23}$ and $r_{34}$, while the qualitative features of (c) are similar to
those of (a). One interpretation for this difference is as follows. We see a
general tendency for the inferred EoS to undergo a rapid stiffening at some
threshold point of $\varepsilon$ possibly to support heavy NSs. If the EoS is
stiffer before the threshold, the EoS after the threshold tends to be softer,
which is also favored to set the maximum mass and radius reasonably small. In
figure 8 (a) for the PREs/qLMXBs distribution, $r_{34}$ is negative reflecting
that the stiffening threshold is found around $c_{s,n=3}^{2}$ in view of the
blue data in figure 6. In the same way, negative $r_{23}$ in figure 8 (b) for
the HP/NICER/GW distribution is understood from the stiffening threshold at
$c_{s,n=2}^{2}$ as shown by the green data in figure 6.

##  5 Discussions

We discuss consistency checks of our results $\mathcal{P}_{\mathrm{obs}}$ by
manipulating our NN procedures. To discuss uncertainty propagation, in section
5.1, we perform test analysis with artificially rescaled
$\mathcal{P}_{\mathrm{obs}}$. We then argue parametrization dependence for the
discretization of $\varepsilon$ in section 5.2 and distribution for sampling
dataset and its prior dependence in section 5.3. In section 5.4, we extend the
NN architecture with additional physical information of the central pressure
$p_{c}$, which is useful to manifest the correlation between the NS properties
and our results. We also use the NN with $p_{c}$ to discuss noise dependence
in the data augmentation in section 5.5. Finally, we comment on the direct
evaluation of Jacobian as a possible alternative framework in section 5.6.

###  5.1 Test analysis with rescaled observational distributions

Figure 9: Speed of sound squared, $c_{s,n}^{2}$, resulting from rescaled
$\mathcal{P}_{\mathrm{obs}}$ in the MC integration using the full distribution
with $n_{\mathrm{NS}}=20$. The error bars show the standard deviation of
$\bm{Y}$ over the modified distribution by replacing
$\mathcal{P}_{\mathrm{obs}}$ in eq. (11) with $\mathcal{P}_{s}$ defined in eq.
(24).

We naturally expect that the uncertainty on the output side should be enhanced
if the input has a larger uncertainty. Here, we verify whether this expected
behavior is properly realized in our NN method by modifying
$\mathcal{P}_{\mathrm{obs}}$ in the MC integration. Specifically, we rescale
$\mathcal{P}_{\mathrm{obs}}(\bm{X})$ around its mode,
$\hat{\bm{X}}=\arg\max_{\bm{X}}\mathcal{P}_{\mathrm{obs}}(\bm{X})$ and
introduce a scaling parameter $s\in\mathbb{R}_{>0}$. We define a modified
distribution as

$\displaystyle\mathcal{P}_{s}(\bm{X}){\rm
d}\bm{X}:=s\mathcal{P}_{\mathrm{obs}}(\hat{\bm{X}}+s^{-1}\Delta\hat{\bm{X}}){\rm
d}\bm{X}\,,\qquad\Delta\hat{\bm{X}}:=\bm{X}-\hat{\bm{X}}\,.$ (24)

Figure 9 shows the results for $s=0.5$, $1.0$, and $2.0$. We obtained the
expected $s$-dependence, namely, the error bars become larger (smaller) with
increasing (decreasing, respectively) $s$. The uncertainty on the input side
characterized by $\mathcal{P}_{s}$ reasonably propagates to the uncertainty on
the output side quantified by the error bars in figure 9. We also observe a
trend of increasing mean value with increasing $s$. This can be understood
from the difference between the average and mode in
$\mathcal{P}_{\mathrm{obs}}$.

###  5.2 Dependence on the number of segments in the EoS parametrization

Figure 10: Comparison of $c_{s,n}^{2}$ estimated by $N$-segment polytropes
with $N=5$ (blue) and $N=8$ (red). The error bars show the standard deviation
of $\bm{Y}$ over $\mathcal{P}_{\mathrm{EoS}}^{\prime}$.

As explained in section 3.2, we parametrize our EoS by dividing the range of
density into $N$ segments in the logarithmic scale, i.e., by a $N$-segment
polytrope. Throughout this work, we fix $N=5$ unless otherwise noted. To check
the bias on the specific choice of the EoS parametrization, we can check the
robustness of the results on different $N$. We show a comparison between the
cases with $N=5$ and $N=8$ in figure 10. Overall, the $N=8$ results are fairly
consistent with the $N=5$ results, interpolating the $N=5$ points. This
implies that our NN architecture could be safely extended for a larger $N$ to
capture more detailed features of EoS when the available data is increased in
both quantity and quality in the future. The EoS rapidly rises around
$\varepsilon/\varepsilon_{0}\simeq 2$–$4$ for both $N=5$ and $N=8$, and we can
safely conclude that the rapid increase in $c_{s}^{2}$ is physical and not
caused by the EoS parametrization artifact. The discussion on the optimal EoS
parametrization will be more closely reported elsewhere.

###  5.3 Training dataset sampling and prior dependence check

Figure 11: Dependence of $c_{s,n}^{2}$ on the distribution assumed for
sampling $n_{\mathrm{NS}}$ points from the $M$-$R$ curve. For descriptions of
the assumed distributions, see the text. The error bars show the standard
deviation of $\bm{Y}$ over $\mathcal{P}_{\mathrm{EoS}}^{\prime}$.

In the training data generation, as explained in step 3 in section 3.3, we
sample $n_{\mathrm{NS}}$ neutron stars from the uniform distribution in the
range $[M_{\min},M_{\max}]$, which is denoted by $U(M_{\min},M_{\max})$. We do
not have a priori knowledge about the distribution of NSs on the $M$-$R$
curve. Our choice of the uniform distribution in mass is just a possible
scenario. This is a factor to comprise our prior assumption in the analysis.
We should avoid any biases by this distribution choice, i.e., the results
should not depend on our prior assumption if the EoS is well-constrained by
the NS data. To check this, we discuss the results with different choices for
the distributions of sampled NS masses.

Here, we compare three different distributions. The first choice is the
uniform distribution, $U(M_{\min},M_{\max})$, as adopted in the main analysis
of this work. Instead of mass, one can also specify a point on the $M$-$R$
curve through the NS central pressure $p_{c}$ as explained in eq. (18). Here,
it would be reasonable to take $U(p_{c,{\min}},p_{c,{\max}})$ as our second
choice, where $p_{c,\mathrm{min/max}}$ is fixed by the condition,
$M(p_{c,\mathrm{min/max}})=M_{\mathrm{min/max}}$. We also make a comparison
with the log uniform distribution, namely, $U(\log p_{c,{\min}},\log
p_{c,{\max}})$ as the third choice.

In figure 11, we compare the results from different sampling distributions. We
observe that the results are consistent with each other within the error bars
up to the fourth segment. In the fifth (highest density) segment, $c_{s}^{2}$
depends on the sampling distributions relatively more strongly than the lower-
density regions. This implies that the fifth segment may not be sufficiently
constrained by the current analysis. One may think that the available data has
no sensitivity to such high-density regions, but this is not necessarily the
case as we demonstrate in the next subsection.

###  5.4 Extension with the NS central pressure

Figure 12: Central pressure of the NSs from the NN extended with $p_{c}$. The
results inferred from the full, PREs/qLMXBs, and HP/NICER/GW distributions are
shown by the red, blue, and green bars, respectively, which represent the
standard deviation of $\bm{Y}$ over $\mathcal{P}_{\mathrm{EoS}}^{\prime}$.

As discussed in section 2, the uncertainty in the EoS chiefly comes from two
factors, i.e., the finiteness of the observed NS number and the limitation of
observational constraints on the observed NS properties. Some randomness is
inevitable in determining both the true NS properties and the observation
process. Although our primary interest is the true EoS inferred in principle
from $\\{(M_{*}^{(i)},R_{*}^{(i)})\\}_{i=1}^{n_{\mathrm{NS}}}$, the other NS
properties, if any, should also be physical information, which would be
precisely determined by the data if the observation had unlimited resolution.
To have an insight into the randomness involved in the EoS determination, it
is useful to check how the other NS properties are estimated and how they are
correlated to the estimated EoS within our method. Here, for the demonstration
purpose, we choose a single parameter, i.e., the central pressure $p_{c}$ to
represent other NS properties on top of the EoS parameters.

Here, we extend our NN model in section 3.4 to include the NS central
pressures, $\\{p_{c}^{(i)}\\}_{i=1}^{n_{\mathrm{NS}}}$, in the output layer.
Specifically, the output data is

$\displaystyle\bm{Y}$
$\displaystyle:=\\{c_{s,n=1}^{2},\dots,c_{s,n=5}^{2},\,p_{c}^{(1)},\dots,p_{c}^{(n_{\mathrm{NS}})}\\}\,,$
(25)

where $p^{(i)}_{c}$ is given in unit of $\text{MeV}/\text{fm}^{3}$ and
determined by $(M_{*}^{(i)},R_{*}^{(i)})$ in step 3 in section 3.3. We
optimize this model in the same way as previously but with the dataset
extended with $\\{p_{c}^{(i)}\\}_{i=1}^{n_{\mathrm{NS}}}$.

Figure 12 shows the results for the full distribution ($n_{\mathrm{NS}}=20$ by
the red bars), the PREs/qLMXBs distribution ($n_{\mathrm{NS}}=14$ by the blue
bars), and the HP/NICER/GW distribution ($n_{\mathrm{NS}}=6$ by the green
bars). The mean values and errors are evaluated by eq. (22). In the full
results, $p_{c}$’s for the HP/NICER/GW NSs indicated by the rightmost six red
bars in figure 12 are larger than $p_{c}$’s from the HP/NICER/GW distribution
only indicated by the green bars in figure 12. This enhancement of $p_{c}$ for
the HP/NICER/GW NSs in the full case is caused by the inclusion of the
PREs/qLMXBs NSs. This also explains the difference in the $c_{s}^{2}$ behavior
inferred from the full and the partial HP/NICER/GW distributions shown in
figure 6.

Figure 12 also tells us the highest density range that each NS can reach. From
figure 7, one can read out the central energy density $\varepsilon_{c}$
corresponding to $p_{c}$ of the red bars in figure 12. The largest $p_{c}$ is
realized for J0740+6620, and the corresponding energy density turns out to be
$\varepsilon_{c}\simeq 5.6\varepsilon_{0}$, which is larger than the boundary
of the fifth segment at $\varepsilon\simeq 5.2\varepsilon_{0}$. Therefore, on
the one hand, in the current analysis with the full distribution, we can say
that the observational data constrains the EoS up to the fifth segment. On the
other hand, in the previous analysis, we only used the PREs/qLMXBs NS data.
The central density $\varepsilon_{c}$ corresponding to the largest $p_{c}$
read from the maximum value of the blue bars in figure 12 is smaller than the
boundary of the fifth segment. Therefore, we confirm that the fifth segment of
the EoS in our previous works Fujimoto:2019hxv ; Fujimoto:2021zas  was not
well constrained by the NS data, which is consistent with our previous
discussion in ref. Fujimoto:2021zas .

A similar conclusion can be drawn for the EoS inferred from the HP/NICER/GW NS
data only. The central density $\varepsilon_{c}$ corresponding to the maximum
$p_{c}$ of the green bars in figure 12 is smaller than the boundary of the
fifth segment. Thus, the fifth segment of the EoS may not be constrained. This
explains why the green error bar in the fifth segment in figure 6 is small
around $c_{s,n}^{2}\simeq 0.5$.

###  5.5 Noise dependence in data augmentation

(a)

(b)

Figure 13: Noise dependence of (a) $c_{s,n}^{2}$ and (b) $p_{c}$ for the full
distribution with $n_{\mathrm{NS}}=20$. The noises are generated by the
Gaussian distribution with the standard deviation $s\sigma$ modified by a
parameter $s$. The plotted results are for $s=0.1$ (black), $0.5$ (orange),
$1.0$ (blue), $2.0$ (green), and $4.0$ (red). The error bars show the standard
deviation of $\bm{Y}$ over $\mathcal{P}_{\mathrm{EoS}}^{\prime}$.

We discuss noise dependence in data augmentation using the NN model extended
with $p_{c}$. The inclusion of $p_{c}$ is useful to clarify the effect of the
data augmentation, which is important in the training data generation. As
explained in step 4 in section 3.3, we randomly sample $n_{\mathrm{NS}}$
observational points from the Gaussian distribution. A concern with this
method is that there is no inherent way to fix the noise magnitude. Below, we
quantitatively study the effect of changing the noise magnitude on the
inferred results.

We vary the noise size by changing the training dataset with
$s\in\mathbb{R}_{>0}$ as

$\displaystyle({M}^{(i)},R^{(i)})\;\sim\;\mathcal{N}(M_{*}^{(i)},(s\sigma_{M}^{(i)})^{2})\times\mathcal{N}(R_{*}^{(i)},(s\sigma_{R}^{(i)})^{2})\,.$
(26)

We show the results for $s=0.1$, $0.5$, $1.0$, $2.0$, and $4.0$ in figure 13.
In panel (a) of figure 13, we see that the error bars are larger for smaller
$s$. The behavior of the sound velocity also does not appear to be reasonable
when $s$ is too small $\sim 0.1\sigma$; the sound velocity is mostly vanishing
in the third segment, while it is mostly the light velocity in the fourth and
fifth segments. This is because the NN trained with small $s$ does not
properly learn the existence of the observational errors and expects the
points $(M^{(i)},R^{(i)})$ in the input to align on a single curve. The NN
with small $s$ cannot properly handle the unaligned input $\bm{X}$ generated
by $\mathcal{P}_{\mathrm{obs}}$ and starts to produce biased results.

The results do not seem sensible either in the opposite case of large $s\sim
4$; the error bars are unnaturally small. This is because the NN learns that
the input is not trustable when it is trained by large $s$ so tends to predict
a constant output without using the input. In other words, the NN trained by
large $s$ outputs a constant prediction mainly determined by the prior
distribution, to which the structures of the input distribution
$\mathcal{P}_{\mathrm{obs}}$ are not reflected. In particular, this loss of
the input distribution can be observed in the results of the central
pressures, $p_{c}$, as shown in the panel (b) of figure 13. The results of
$p_{c}$ for $s=4.0$ are located around $p_{c}\sim 100\;{\rm MeV/fm^{3}}$ only.
This implies that data augmentation with too large noise unreasonably
undervalues the meaningful information of NS features in the observation data.
This noise dependence could be relevant to other prior dependence. We will
further investigate this issue in our future work.

###  5.6 Alternative framework: direct evaluation of Jacobian

Finally, we mention a possible alternative framework. In section 2.3, we
constructed the inverse mapping of the TOV equation $\Psi_{\mathrm{inv}}$ in
eq. (10) to obtain the most likely EoS. We already know the forward mapping
$\Psi_{\mathrm{TOV}}:\bm{Y}\mapsto\bm{X}$ by solving the TOV equation, so that
we can construct the _pullback_ of $\mathcal{P}_{\mathrm{obs}}$ by
$\Psi_{\mathrm{TOV}}$ to obtain $\mathcal{P}_{\mathrm{EoS}}$, i.e.,

$\displaystyle\mathcal{P}_{\mathrm{EoS}}(\bm{Y}){\rm
d}\bm{Y}=(\Psi_{\mathrm{TOV}}^{*}{\rm
d}\mathcal{P}_{\mathrm{obs}})(\bm{Y})=\mathcal{P}_{\mathrm{obs}}(\Psi_{\mathrm{TOV}}(\bm{Y}))\biggl{|}\frac{{\rm
d}\Psi_{\mathrm{TOV}}(\bm{Y})}{{\rm d}\bm{Y}}\biggr{|}{\rm d}\bm{Y}\,.$ (27)

In reality, the above formula should be generalized to take account of
randomness in the observation uncertainty:

$\displaystyle\mathcal{P}_{\mathrm{EoS}}(\bm{Y}){\rm d}\bm{Y}\sim\int{\rm
d}\mathcal{P}(\omega)\mathcal{P}_{\mathrm{obs}}(\Psi_{\mathrm{TOV}}(\bm{Y},\omega))\biggl{|}\frac{{\rm
d}\Psi_{\mathrm{TOV}}(\bm{Y},\omega)}{{\rm d}\bm{Y}}\biggr{|}{\rm d}\bm{Y}\,,$
(28)

where $\omega$ symbolically denotes the random variables involved in sampling
$M$ with observational errors. In this framework, the inversion is
incorporated in the pullback, so an explicit construction of the inverse
mapping is unnecessary. Historically speaking, this pullback method, in which
one must explicitly compute the Jacobian, was the first approach taken in the
realm of Bayesian inference Ozel:2009da ; Ozel:2010fw .

Indeed, if the inversion were exact (i.e.,
$\Psi_{\mathrm{inv}}\circ\Psi_{\mathrm{TOV}}=\mathrm{id}_{\mathrm{EoS}}$ and
$\Psi_{\mathrm{TOV}}\circ\Psi_{\mathrm{inv}}=\mathrm{id}_{MR}$), the above
definition by the pullback exactly reproduces eq. (11). In practice, however,
the evaluation of the Jacobian is a nontrivial numerical task, and our method
presented in this work is computationally more straightforward.

##  6 Conclusions

We developed a new method to properly incorporate the shape of the joint
probability distribution of NS masses and radii, which we denoted
$\mathcal{P}_{\mathrm{obs}}$, by performing the MC integration. The new method
naturally incorporates the uncertainty quantification in the final output from
the NN, which had been more like a black box in our previous method. The
results in the present work imply that the previous evaluation of the output
error by bootstrap aggregating gives an estimate consistent with the current
method.

The chief result of this work is presented in figure 6. The inclusion of new
NS data upon the older ones, particularly the $M$-$R$ data from the NICER and
GW170817 event, led to a sharp change in the stiffness of the EoS as displayed
in figure 7 (see also the error bars marked as “Full ($n_{\mathrm{NS}}=20$)”
in figure 6). This is quantified by the weaker correlation between the second
and third segments of the piecewise polytrope shown in panel (c) of figure 8
compared to panel (a). Our previous analysis of observational data suggested a
peak structure in the speed of sound squared, $c_{s}^{2}$, as a function of
energy density, and the current analysis including new data strengthens this
claim, though the statistical significance is not yet achieved with the
current accuracy of the NS data.

It is an important question up to which density the observational data can in
principle constrain. In the present study, the decrease of the sound speed at
the largest energy density range, in particular, is certainly a feature under
the constraint. This contrasts with our previous work. This decrease in the
previous work may have been attributed to the prior assumptions due to the
lack of constraining data. In the new analysis, we substantiated that the
$c_{s}^{2}$ behavior in this high-density region is not a mere reflection of
the prior assumptions but is also constrained by the new data: First, we
changed the uncertainty in the input data and verified that the input
uncertainty is properly propagated to the output uncertainty as one can see in
figure 9. This fact implies that the last segment is also constrained by the
new data. Indeed, as shown in figure 12, the pulsar J0740+6620 is predicted to
have the central pressure $p_{c}\approx 300~{}\text{MeV/fm}^{3}$ in our
analysis, which roughly corresponds to the energy density in the last segment
in our parametrization.

Furthermore, we examined three artificially adjusted prior assumptions to see
the effects on the EoS. We demonstrated that the prior assumptions do not
affect the results within the uncertainties. First, we changed the number of
segments in the piecewise polytrope. We verified in figure 10 that the results
are insensitive to the number of segments within the uncertainty. Second, for
the input to the NN, we changed the distributions for sampling
$n_{\mathrm{NS}}$ points along the $M$-$R$ curve. This change may lead to
substantial prior dependence, but as seen in figure 11, the results are still
consistent with each other within the uncertainties. Finally, we changed the
magnitude of the noise used in the data augmentation. In figure 13, we found
that too small noise leads to physically unacceptable results, but too large
noise also leads to unconstrained results with underestimated output
uncertainty. The current choice of noise comparable to the observational error
($1.0\sigma$) turns out to be the most reasonable choice. Related to this, we
also note that our new method requires a smaller size of the training dataset
because this MC integration can lessen the magnitude of the data augmentation
necessary in the uncertainty analysis. A mathematically founded method to
optimize the data augmentation procedure to properly incorporate the effect of
uncertainty in the NN is yet to be explored. More detailed investigations of
these prior effects will be presented elsewhere.

###### Acknowledgements.

We thank Len Brandes, Tetsuo Hatsuda, and Shuhei Minato for useful
discussions. The work of Y.F. is supported by the Japan Society for the
Promotion of Science (JSPS) through the Overseas Research Fellowship and by
the INT’s U.S. DOE Grant No. DE-FG02-00ER41132. This work was supported by
JSPS KAKENHI Grant Nos. 22H01216 (K. F), 22H05118 (K. F and S. K), and
23K13102 (K.M.).

## References

  * (1) K. Fukushima and T. Hatsuda, _The phase diagram of dense QCD_ , _Rept. Prog. Phys._ 74 (2011) 014001 [1005.4814].
  * (2) G. Baym, T. Hatsuda, T. Kojo, P.D. Powell, Y. Song and T. Takatsuka, _From hadrons to quarks in neutron stars: a review_ , _Rept. Prog. Phys._ 81 (2018) 056902 [1707.04966].
  * (3) G. Aarts et al., _Phase Transitions in Particle Physics: Results and Perspectives from Lattice Quantum Chromo-Dynamics_ , _Prog. Part. Nucl. Phys._ 133 (2023) 104070 [2301.04382].
  * (4) S. Borsanyi, Z. Fodor, C. Hoelbling, S.D. Katz, S. Krieg and K.K. Szabo, _Full result for the QCD equation of state with 2+1 flavors_ , _Phys. Lett. B_ 730 (2014) 99 [1309.5258].
  * (5) HotQCD collaboration, _Equation of state in ( 2+1 )-flavor QCD_ , _Phys. Rev. D_ 90 (2014) 094503 [1407.6387].
  * (6) HotQCD collaboration, _Chiral crossover in QCD at zero and non-zero chemical potentials_ , _Phys. Lett. B_ 795 (2019) 15 [1812.08235].
  * (7) M.G. Alford, A. Kapustin and F. Wilczek, _Imaginary chemical potential and finite fermion density on the lattice_ , _Phys. Rev. D_ 59 (1999) 054502 [hep-lat/9807039].
  * (8) D.T. Son and M.A. Stephanov, _QCD at finite isospin density_ , _Phys. Rev. Lett._ 86 (2001) 592 [hep-ph/0005225].
  * (9) K. Iida and E. Itou, _Velocity of sound beyond the high-density relativistic limit from lattice simulation of dense two-color QCD_ , _PTEP_ 2022 (2022) 111B01 [2207.01253].
  * (10) B.B. Brandt, F. Cuteri and G. Endrodi, _Equation of state and speed of sound of isospin-asymmetric QCD on the lattice_ , _JHEP_ 07 (2023) 055 [2212.14016].
  * (11) R. Abbott, W. Detmold, F. Romero-López, Z. Davoudi, M. Illa, A. Parreño et al., _Lattice quantum chromodynamics at large isospin density: 6144 pions in a box_ , 2307.15014.
  * (12) Y. Fujimoto and S. Reddy, _Bounds on the Equation of State from QCD Inequalities and Lattice QCD_ , 2310.09427.
  * (13) R.C. Tolman, _Static solutions of Einstein’s field equations for spheres of fluid_ , _Phys. Rev._ 55 (1939) 364.
  * (14) J.R. Oppenheimer and G.M. Volkoff, _On massive neutron cores_ , _Phys. Rev._ 55 (1939) 374.
  * (15) L. Lindblom, _Determining the Nuclear Equation of State from Neutron-Star Masses and Radii_ , _Astrophys. J._ 398 (1992) 569.
  * (16) K. Hebeler, J.M. Lattimer, C.J. Pethick and A. Schwenk, _Constraints on neutron star radii based on chiral effective field theory interactions_ , _Phys. Rev. Lett._ 105 (2010) 161102 [1007.1746].
  * (17) K. Hebeler, J.M. Lattimer, C.J. Pethick and A. Schwenk, _Equation of state and neutron star properties constrained by nuclear physics and observation_ , _Astrophys. J._ 773 (2013) 11 [1303.4662].
  * (18) S. Gandolfi, J. Lippuner, A.W. Steiner, I. Tews, X. Du and M. Al-Mamun, _From the microscopic to the macroscopic world: from nucleons to neutron stars_ , _J. Phys. G_ 46 (2019) 103001 [1903.06730].
  * (19) C. Drischler, R.J. Furnstahl, J.A. Melendez and D.R. Phillips, _How Well Do We Know the Neutron-Matter Equation of State at the Densities Inside Neutron Stars? A Bayesian Approach with Correlated Uncertainties_ , _Phys. Rev. Lett._ 125 (2020) 202702 [2004.07232].
  * (20) C. Drischler, S. Han, J.M. Lattimer, M. Prakash, S. Reddy and T. Zhao, _Limiting masses and radii of neutron stars and their implications_ , _Phys. Rev. C_ 103 (2021) 045808 [2009.06441].
  * (21) C. Drischler, S. Han and S. Reddy, _Large and massive neutron stars: Implications for the sound speed within QCD of dense matter_ , _Phys. Rev. C_ 105 (2022) 035808 [2110.14896].
  * (22) B.A. Freedman and L.D. McLerran, _Fermions and Gauge Vector Mesons at Finite Temperature and Density. 3. The Ground State Energy of a Relativistic Quark Gas_ , _Phys. Rev. D_ 16 (1977) 1169.
  * (23) V. Baluni, _Nonabelian Gauge Theories of Fermi Systems: Chromotheory of Highly Condensed Matter_ , _Phys. Rev. D_ 17 (1978) 2092.
  * (24) A. Vuorinen, _The Pressure of QCD at finite temperatures and chemical potentials_ , _Phys. Rev. D_ 68 (2003) 054017 [hep-ph/0305183].
  * (25) A. Kurkela, P. Romatschke and A. Vuorinen, _Cold Quark Matter_ , _Phys. Rev. D_ 81 (2010) 105021 [0912.1856].
  * (26) T. Gorda, A. Kurkela, P. Romatschke, S. Säppi and A. Vuorinen, _Next-to-Next-to-Next-to-Leading Order Pressure of Cold Quark Matter: Leading Logarithm_ , _Phys. Rev. Lett._ 121 (2018) 202701 [1807.04120].
  * (27) T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi and A. Vuorinen, _Soft Interactions in Cold Quark Matter_ , _Phys. Rev. Lett._ 127 (2021) 162003 [2103.05658].
  * (28) T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi and A. Vuorinen, _Cold quark matter at N3LO: Soft contributions_ , _Phys. Rev. D_ 104 (2021) 074015 [2103.07427].
  * (29) P. Demorest, T. Pennucci, S. Ransom, M. Roberts and J. Hessels, _Shapiro Delay Measurement of A Two Solar Mass Neutron Star_ , _Nature_ 467 (2010) 1081 [1010.5788].
  * (30) J. Antoniadis et al., _A Massive Pulsar in a Compact Relativistic Binary_ , _Science_ 340 (2013) 6131 [1304.6875].
  * (31) E. Fonseca et al., _The NANOGrav Nine-year Data Set: Mass and Geometric Measurements of Binary Millisecond Pulsars_ , _Astrophys. J._ 832 (2016) 167 [1603.00545].
  * (32) NANOGrav collaboration, _The NANOGrav 11-year Data Set: High-precision timing of 45 Millisecond Pulsars_ , _Astrophys. J. Suppl._ 235 (2018) 37 [1801.01837].
  * (33) NANOGrav collaboration, _Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar_ , _Nature Astron._ 4 (2019) 72 [1904.06759].
  * (34) E. Fonseca et al., _Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620_ , _Astrophys. J. Lett._ 915 (2021) L12 [2104.00880].
  * (35) T.E. Riley et al., _A $NICER$ View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation_, _Astrophys. J. Lett._ 887 (2019) L21 [1912.05702].
  * (36) M.C. Miller et al., _PSR J0030+0451 Mass and Radius from $NICER$ Data and Implications for the Properties of Neutron Star Matter_, _Astrophys. J. Lett._ 887 (2019) L24 [1912.05705].
  * (37) M.C. Miller et al., _The Radius of PSR J0740+6620 from NICER and XMM-Newton Data_ , _Astrophys. J. Lett._ 918 (2021) L28 [2105.06979].
  * (38) T.E. Riley et al., _A NICER View of the Massive Pulsar PSR J0740+6620 Informed by Radio Timing and XMM-Newton Spectroscopy_ , _Astrophys. J. Lett._ 918 (2021) L27 [2105.06980].
  * (39) LIGO Scientific, Virgo collaboration, _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101 [1710.05832].
  * (40) LIGO Scientific, Virgo collaboration, _GW170817: Measurements of neutron star radii and equation of state_ , _Phys. Rev. Lett._ 121 (2018) 161101 [1805.11581].
  * (41) LIGO Scientific, Virgo collaboration, _Properties of the binary neutron star merger GW170817_ , _Phys. Rev. X_ 9 (2019) 011001 [1805.11579].
  * (42) A. Kurkela, E.S. Fraga, J. Schaffner-Bielich and A. Vuorinen, _Constraining neutron star matter with Quantum Chromodynamics_ , _Astrophys. J._ 789 (2014) 127 [1402.6618].
  * (43) E. Annala, T. Gorda, A. Kurkela and A. Vuorinen, _Gravitational-wave constraints on the neutron-star-matter Equation of State_ , _Phys. Rev. Lett._ 120 (2018) 172703 [1711.02644].
  * (44) E. Annala, T. Gorda, A. Kurkela, J. Nättilä and A. Vuorinen, _Evidence for quark-matter cores in massive neutron stars_ , _Nature Phys._ 16 (2020) 907 [1903.09121].
  * (45) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis et al., _Multimessenger Constraints for Ultradense Matter_ , _Phys. Rev. X_ 12 (2022) 011058 [2105.05132].
  * (46) A.W. Steiner, J.M. Lattimer and E.F. Brown, _The Equation of State from Observed Masses and Radii of Neutron Stars_ , _Astrophys. J._ 722 (2010) 33 [1005.0811].
  * (47) A.W. Steiner, J.M. Lattimer and E.F. Brown, _The Neutron Star Mass-Radius Relation and the Equation of State of Dense Matter_ , _Astrophys. J. Lett._ 765 (2013) L5 [1205.6871].
  * (48) F. Ozel, D. Psaltis, T. Guver, G. Baym, C. Heinke and S. Guillot, _The Dense Matter Equation of State from Neutron Star Radius and Mass Measurements_ , _Astrophys. J._ 820 (2016) 28 [1505.05155].
  * (49) F. Ozel and D. Psaltis, _Reconstructing the Neutron-Star Equation of State from Astrophysical Measurements_ , _Phys. Rev. D_ 80 (2009) 103003 [0905.1959].
  * (50) F. Ozel, G. Baym and T. Guver, _Astrophysical Measurement of the Equation of State of Neutron Star Matter_ , _Phys. Rev. D_ 82 (2010) 101301 [1002.3153].
  * (51) P. Landry and R. Essick, _Nonparametric inference of the neutron star equation of state from gravitational wave observations_ , _Phys. Rev. D_ 99 (2019) 084049 [1811.12529].
  * (52) Y. Fujimoto, K. Fukushima and K. Murase, _Methodology study of machine learning for the neutron star equation of state_ , _Phys. Rev. D_ 98 (2018) 023019 [1711.06748].
  * (53) Y. Fujimoto, K. Fukushima and K. Murase, _Mapping neutron star data to the equation of state using the deep neural network_ , _Phys. Rev. D_ 101 (2020) 054016 [1903.03400].
  * (54) Y. Fujimoto, K. Fukushima and K. Murase, _Extensive Studies of the Neutron Star Equation of State from the Deep Learning Inference with the Observational Data Augmentation_ , _JHEP_ 03 (2021) 273 [2101.08156].
  * (55) D. Alvarez-Castillo, A. Ayriyan, S. Benic, D. Blaschke, H. Grigorian and S. Typel, _New class of hybrid EoS and Bayesian M-R data analysis_ , _Eur. Phys. J. A_ 52 (2016) 69 [1603.03457].
  * (56) G. Raaijmakers et al., _Constraining the dense matter equation of state with joint analysis of NICER and LIGO/Virgo measurements_ , _Astrophys. J. Lett._ 893 (2020) L21 [1912.11031].
  * (57) S.K. Greif, K. Hebeler, J.M. Lattimer, C.J. Pethick and A. Schwenk, _Equation of state constraints from nuclear physics, neutron star masses, and future moment of inertia measurements_ , _Astrophys. J._ 901 (2020) 155 [2005.14164].
  * (58) M. Al-Mamun, A.W. Steiner, J. Nättilä, J. Lange, R. O’Shaughnessy, I. Tews et al., _Combining Electromagnetic and Gravitational-Wave Constraints on Neutron-Star Masses and Radii_ , _Phys. Rev. Lett._ 126 (2021) 061101 [2008.12817].
  * (59) S. Huth et al., _Constraining Neutron-Star Matter with Microscopic and Macroscopic Collisions_ , _Nature_ 606 (2022) 276 [2107.06229].
  * (60) S. Altiparmak, C. Ecker and L. Rezzolla, _On the Sound Speed in Neutron Stars_ , _Astrophys. J. Lett._ 939 (2022) L34 [2203.14974].
  * (61) Y. Lim and J.W. Holt, _Neutron Star Radii, Deformabilities, and Moments of Inertia from Experimental and Ab Initio Theory Constraints of the 208Pb Neutron Skin Thickness_, _Galaxies_ 10 (2022) 99 [2204.09000].
  * (62) E.V. Chimanski, R.V. Lobato, A.R. Goncalves and C.A. Bertulani, _Bayesian Exploration of Phenomenological EoS of Neutron/Hybrid Stars with Recent Observations_ , _Particles_ 6 (2023) 198 [2205.01174].
  * (63) M. Marczenko, L. McLerran, K. Redlich and C. Sasaki, _Reaching percolation and conformal limits in neutron stars_ , _Phys. Rev. C_ 107 (2023) 025802 [2207.13059].
  * (64) L. Brandes, W. Weise and N. Kaiser, _Inference of the sound speed and related properties of neutron stars_ , _Phys. Rev. D_ 107 (2023) 014011 [2208.03026].
  * (65) C. Ecker and L. Rezzolla, _Impact of large-mass constraints on the properties of neutron stars_ , _Mon. Not. Roy. Astron. Soc._ 519 (2022) 2615 [2209.08101].
  * (66) J.-L. Jiang, C. Ecker and L. Rezzolla, _Bayesian Analysis of Neutron-star Properties with Parameterized Equations of State: The Role of the Likelihood Functions_ , _Astrophys. J._ 949 (2023) 11 [2211.00018].
  * (67) J. Takatsy, P. Kovacs, G. Wolf and J. Schaffner-Bielich, _What neutron stars tell about the hadron-quark phase transition: A Bayesian study_ , _Phys. Rev. D_ 108 (2023) 043002 [2303.00013].
  * (68) E. Annala, T. Gorda, J. Hirvonen, O. Komoltsev, A. Kurkela, J. Nättilä et al., _Strongly interacting matter exhibits deconfined behavior in massive neutron stars_ , 2303.11356.
  * (69) L. Brandes, W. Weise and N. Kaiser, _Evidence against a strong first-order phase transition in neutron star cores: impact of new data_ , 2306.06218.
  * (70) C. Musolino, C. Ecker and L. Rezzolla, _On the maximum mass and oblateness of rotating neutron stars with generic equations of state_ , 2307.03225.
  * (71) J. Richter and B.-A. Li, _Empirical radius formulas for canonical neutron stars from bidirectionally selecting EOS features in extended Bayesian analyses of observational data_ , 2307.05848.
  * (72) C. Providência, T. Malik, M.B. Albino and M. Ferreira, _Neutron star equation of state: identifying hadronic matter characteristics_ , 2307.05086.
  * (73) B.-J. Cai, B.-A. Li and Z. Zhang, _Central Speed of Sound, Trace Anomaly and Observables of Neutron Stars from Perturbative Analyses of Scaled TOV Equations_ , 2307.15223.
  * (74) P.T.H. Pang, L. Sivertsen, R. Somasundaram, T. Dietrich, S. Sen, I. Tews et al., _Probing Quarkyonic Matter in Neutron Stars with the Bayesian Nuclear-Physics Multi-Messenger Astrophysics Framework_ , 2308.15067.
  * (75) Y.-Z. Fan, M.-Z. Han, J.-L. Jiang, D.-S. Shao and S.-P. Tang, _Maximum gravitational mass $M_{\rm TOV}=2.25^{+0.08}_{-0.07}M_{\odot}$ inferred at about $3\%$ precision with multimessenger data of neutron stars_, 2309.12644.
  * (76) C.Y. Tsang, M.B. Tsang, W.G. Lynch, R. Kumar and C.J. Horowitz, _Determination of the Equation of State from Nuclear Experiments and Neutron Star Observations_ , 2310.11588.
  * (77) S.-P. Tang, M.-Z. Han, Y.-J. Huang, Y.-Z. Fan and D.-M. Wei, _Mass and Radius of the Most Massive Neutron Star: the Probe of the Equation of State and Perturbative QCD_ , 2311.13805.
  * (78) O. Komoltsev, R. Somasundaram, T. Gorda, A. Kurkela, J. Margueron and I. Tews, _Equation of state at neutron-star densities and beyond from perturbative QCD_ , 2312.14127.
  * (79) S.M.A. Imam, T. Malik, C. Providência and B.K. Agrawal, _Implications of comprehensive nuclear and astrophysics data on the equations of state of neutron star matter_ , 2401.06018.
  * (80) R. Essick, P. Landry and D.E. Holz, _Nonparametric Inference of Neutron Star Composition, Equation of State, and Maximum Mass with GW170817_ , _Phys. Rev. D_ 101 (2020) 063007 [1910.09740].
  * (81) P. Landry, R. Essick and K. Chatziioannou, _Nonparametric constraints on neutron star matter with existing and upcoming gravitational wave and pulsar observations_ , _Phys. Rev. D_ 101 (2020) 123007 [2003.04880].
  * (82) R. Essick, I. Tews, P. Landry, S. Reddy and D.E. Holz, _Direct Astrophysical Tests of Chiral Effective Field Theory at Supranuclear Densities_ , _Phys. Rev. C_ 102 (2020) 055803 [2004.07744].
  * (83) R. Essick and P. Landry, _Discriminating between Neutron Stars and Black Holes with Imperfect Knowledge of the Maximum Neutron Star Mass_ , _Astrophys. J._ 904 (2020) 80 [2007.01372].
  * (84) R. Essick, I. Tews, P. Landry and A. Schwenk, _Astrophysical Constraints on the Symmetry Energy and the Neutron Skin of Pb208 with Minimal Modeling Assumptions_ , _Phys. Rev. Lett._ 127 (2021) 192701 [2102.10074].
  * (85) I. Legred, K. Chatziioannou, R. Essick, S. Han and P. Landry, _Impact of the PSR J0740+6620 radius constraint on the properties of high-density matter_ , _Phys. Rev. D_ 104 (2021) 063003 [2106.05313].
  * (86) R. Essick, P. Landry, A. Schwenk and I. Tews, _Detailed examination of astrophysical constraints on the symmetry energy and the neutron skin of Pb208 with minimal modeling assumptions_ , _Phys. Rev. C_ 104 (2021) 065804 [2107.05528].
  * (87) I. Legred, K. Chatziioannou, R. Essick and P. Landry, _Implicit correlations within phenomenological parametric models of the neutron star equation of state_ , _Phys. Rev. D_ 105 (2022) 043016 [2201.06791].
  * (88) T. Gorda, O. Komoltsev and A. Kurkela, _Ab-initio QCD Calculations Impact the Inference of the Neutron-star-matter Equation of State_ , _Astrophys. J._ 950 (2023) 107 [2204.11877].
  * (89) R. Essick, I. Legred, K. Chatziioannou, S. Han and P. Landry, _Phase transition phenomenology with nonparametric representations of the neutron star equation of state_ , _Phys. Rev. D_ 108 (2023) 043013 [2305.07411].
  * (90) D. Mroczek, M.C. Miller, J. Noronha-Hostler and N. Yunes, _Nontrivial features in the speed of sound inside neutron stars_ , 2309.02345.
  * (91) M.a. Ferreira and C. Providência, _Unveiling the nuclear matter EoS from neutron star properties: a supervised machine learning approach_ , _JCAP_ 07 (2021) 011 [1910.05554].
  * (92) F. Morawski and M. Bejger, _Neural network reconstruction of the dense matter equation of state derived from the parameters of neutron stars_ , _Astron. Astrophys._ 642 (2020) A78 [2006.07194].
  * (93) S. Traversi and P. Char, _Structure of Quark Star: A Comparative Analysis of Bayesian Inference and Neural Network Based Modeling_ , _Astrophys. J._ 905 (2020) 9 [2007.10239].
  * (94) P.G. Krastev, _Translating Neutron Star Observations to Nuclear Symmetry Energy via Deep Neural Networks_ , _Galaxies_ 10 (2022) 16 [2112.04089].
  * (95) S. Soma, L. Wang, S. Shi, H. Stöcker and K. Zhou, _Neural network reconstruction of the dense matter equation of state from neutron star observables_ , _JCAP_ 08 (2022) 071 [2201.01756].
  * (96) A. Thete, K. Banerjee and T. Malik, _Inferring the dense matter equation of state from neutron star observations via artificial neural networks_ , 2208.13163.
  * (97) D. Farrell, P. Baldi, J. Ott, A. Ghosh, A.W. Steiner, A. Kavitkar et al., _Deducing neutron star equation of state parameters directly from telescope spectra with uncertainty-aware machine learning_ , _JCAP_ 02 (2023) 016 [2209.02817].
  * (98) S. Soma, L. Wang, S. Shi, H. Stöcker and K. Zhou, _Reconstructing the neutron star equation of state from observational data via automatic differentiation_ , _Phys. Rev. D_ 107 (2023) 083028 [2209.08883].
  * (99) M. Ferreira, V. Carvalho and C. Providência, _Extracting nuclear matter properties from the neutron star matter equation of state using deep neural networks_ , _Phys. Rev. D_ 106 (2022) 103023 [2209.09085].
  * (100) G. Gonçalves, M. Ferreira, J.a. Aveiro, A. Onofre, F.F. Freitas, C. Providência et al., _Machine-learning Love: classifying the equation of state of neutron stars with transformers_ , _JCAP_ 12 (2023) 001 [2210.08382].
  * (101) S. Chatterjee, H. Sudhakaran and R. Mallick, _Analyzing the speed of sound in neutron star with machine learning_ , 2302.13648.
  * (102) P.G. Krastev, _A Deep Learning Approach to Extracting Nuclear Matter Properties from Neutron Star Observations_ , _Symmetry_ 15 (2023) 1123 [2303.17146].
  * (103) D. Farrell, P. Baldi, J. Ott, A. Ghosh, A.W. Steiner, A. Kavitkar et al., _Deducing neutron star equation of state from telescope spectra with machine-learning-derived likelihoods_ , _JCAP_ 12 (2023) 022 [2305.07442].
  * (104) W. Zhou, J. Hu, Y. Zhang and H. Shen, _Nonparametric Model for the Equations of State of a Neutron Star from Deep Neural Network_ , _Astrophys. J._ 950 (2023) 186 [2305.03323].
  * (105) V. Carvalho, M. Ferreira, T. Malik and C. Providência, _Decoding neutron star observations: Revealing composition through Bayesian neural networks_ , _Phys. Rev. D_ 108 (2023) 043031 [2306.06929].
  * (106) V. Carvalho, M. Ferreira and C. Providência, _From NS observations to nuclear matter properties: a machine learning approach_ , 2401.05770.
  * (107) P. Bedaque and A.W. Steiner, _Sound velocity bound and neutron stars_ , _Phys. Rev. Lett._ 114 (2015) 031103 [1408.5116].
  * (108) M.G. Alford, G.F. Burgio, S. Han, G. Taranto and D. Zappalà, _Constraining and applying a generic high-density equation of state_ , _Phys. Rev. D_ 92 (2015) 083002 [1501.07902].
  * (109) I. Tews, J. Carlson, S. Gandolfi and S. Reddy, _Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations_ , _Astrophys. J._ 860 (2018) 149 [1801.01923].
  * (110) Y. Fujimoto, K. Fukushima, L.D. McLerran and M. Praszalowicz, _Trace Anomaly as Signature of Conformality in Neutron Stars_ , _Phys. Rev. Lett._ 129 (2022) 252702 [2207.06753].
  * (111) O. Komoltsev and A. Kurkela, _How Perturbative QCD Constrains the Equation of State at Neutron-Star Densities_ , _Phys. Rev. Lett._ 128 (2022) 202701 [2111.05350].
  * (112) R. Somasundaram, I. Tews and J. Margueron, _Perturbative QCD and the neutron star equation of state_ , _Phys. Rev. C_ 107 (2023) L052801 [2204.14039].
  * (113) K. Masuda, T. Hatsuda and T. Takatsuka, _Hadron–quark crossover and massive hybrid stars_ , _PTEP_ 2013 (2013) 073D01 [1212.6803].
  * (114) K. Masuda, T. Hatsuda and T. Takatsuka, _Hadron-Quark Crossover and Massive Hybrid Stars with Strangeness_ , _Astrophys. J._ 764 (2013) 12 [1205.3621].
  * (115) T. Kojo, P.D. Powell, Y. Song and G. Baym, _Phenomenological QCD equation of state for massive neutron stars_ , _Phys. Rev. D_ 91 (2015) 045003 [1412.1108].
  * (116) L. McLerran and R.D. Pisarski, _Phases of cold, dense quarks at large N(c)_ , _Nucl. Phys. A_ 796 (2007) 83 [0706.2191].
  * (117) L. McLerran and S. Reddy, _Quarkyonic Matter and Neutron Stars_ , _Phys. Rev. Lett._ 122 (2019) 122701 [1811.12503].
  * (118) K.S. Jeong, L. McLerran and S. Sen, _Dynamically generated momentum space shell structure of quarkyonic matter via an excluded volume model_ , _Phys. Rev. C_ 101 (2020) 035201 [1908.04799].
  * (119) T. Kojo, _Stiffening of matter in quark-hadron continuity_ , _Phys. Rev. D_ 104 (2021) 074005 [2106.06687].
  * (120) Y. Fujimoto, K. Fukushima, K. Hotokezaka and K. Kyutoku, _Gravitational Wave Signal for Quark Matter with Realistic Phase Transition_ , _Phys. Rev. Lett._ 130 (2023) 091404 [2205.03882].
  * (121) Y. Fujimoto, T. Kojo and L.D. McLerran, _Momentum Shell in Quarkyonic Matter from Explicit Duality: A Solvable Model Analysis_ , 2306.04304.
  * (122) A. Akmal, V.R. Pandharipande and D.G. Ravenhall, _The Equation of state of nucleon matter and neutron star structure_ , _Phys. Rev. C_ 58 (1998) 1804 [nucl-th/9804027].
  * (123) K. Hebeler and A. Schwenk, _Chiral three-nucleon forces and neutron matter_ , _Phys. Rev. C_ 82 (2010) 014314 [0911.0483].
  * (124) S. Gandolfi, J. Carlson and S. Reddy, _The maximum mass and radius of neutron stars and the nuclear symmetry energy_ , _Phys. Rev. C_ 85 (2012) 032801 [1101.1921].
  * (125) B.D. Serot and J.D. Walecka, _Recent progress in quantum hadrodynamics_ , _Int. J. Mod. Phys. E_ 6 (1997) 515 [nucl-th/9701058].
  * (126) H. Tan, V. Dexheimer, J. Noronha-Hostler and N. Yunes, _Finding Structure in the Speed of Sound of Supranuclear Matter from Binary Love Relations_ , _Phys. Rev. Lett._ 128 (2022) 161101 [2111.10260].
  * (127) S. Bogdanov, C.O. Heinke, F. Özel and T. Güver, _Neutron Star Mass-Radius Constraints of the Quiescent Low-mass X-ray Binaries X7 and X5 in the Globular Cluster 47 Tuc_ , _Astrophys. J._ 831 (2016) 184 [1603.01630].
  * (128) M.C. Miller, _Astrophysical Constraints on Dense Matter in Neutron Stars_ , _Astrophys. Space Sci. Libr._ 461 (2020) 1 [1312.0029].
  * (129) M.C. Miller and F.K. Lamb, _Observational Constraints on Neutron Star Masses and Radii_ , _Eur. Phys. J. A_ 52 (2016) 63 [1604.03894].
  * (130) F. Douchin and P. Haensel, _A unified equation of state of dense matter and neutron star structure_ , _Astron. Astrophys._ 380 (2001) 151 [astro-ph/0111092].
  * (131) F. Chollet, “Keras: Deep learning library for theano and tensorflow.” https://github.com/fchollet/keras, 2015.
  * (132) D.P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , _CoRR_ abs/1412.6980 (2014) [1412.6980].
  * (133) R.W. Romani, D. Kandel, A.V. Filippenko, T.G. Brink and W. Zheng, _PSR J0952 $-$0607: The Fastest and Heaviest Known Galactic Neutron Star_, _Astrophys. J. Lett._ 934 (2022) L17 [2207.05124].