A Bayesian Analysis on Neutron Stars within Relativistic Mean Field Models

Abstract

We present a Bayesian analysis on the equation of state of neutron stars based on a class of relativistic mean field models. The priors on the equation of state are related to the properties of nuclear matter at saturation and the posteriors are obtained through the Bayesian procedure by exploiting recent astrophysical constraints on the mass–radius relations of neutron stars. We find indications of a tension (within the adopted model) between the prior on the nuclear incompressibility and its posterior which in turn seems to suggest a possible phase transition at about twice saturation density to a phase where the nucleon effective mass is strongly reduced. A possible relation with the chiral phase transition in dense matter is also discussed.

1. Introduction

Theoretical calculations of the Equation of State (EoS) of a fluid of interacting particles is in general a very difficult task and it is even more difficult when the gas is degenerate (quantum correlations are important) and particles are strongly interacting. In this respect, theory needs experimental inputs as a guide for choosing the approximations which are necessary to provide a relation e.g., between pressure and energy density or particle density. The case of nucleonic matter (i.e., matter composed by neutrons and protons) is particularly challenging due to the strong interaction and at the same time phenomenologically very relevant since such mixture of particles is present in finite nuclei but also in the core of neutron stars. In both cases, the gas is highly degenerate (the temperature is negligible with respect to the chemical potential) but a first basic difference is the isospin asymmetry, with nuclei being almost symmetric and with the matter of neutron stars being beta-stable (proton fraction at maximum of the order of $0.1$). Second and more importantly, in neutron stars, gravity compresses matter up to densities which can reach ten times the nuclear saturation density. Observational data are available both from experimental nuclear physics and from astrophysical data on the structure of neutron stars. In particular, the very recent detections of the gravitational waves signals from the merger of two neutron stars offer a new EoS sensitive quantity, namely the average tidal deformability of the binary which adds to the more “traditional” measurements of masses and radii of neutron stars [1,2].

A very powerful approach for using data to improve our theoretical understanding of the EoS of nucleonic matter relies on the Bayesian approach. Since the seminal paper [3], a huge effort has been put in trying to infer the properties of the EoS of neutron stars by using astrophysical data. There are several different approaches, see [4] for a complete list, which make use of several different parametrizations of the EoS: e.g., piecewise-polytropic, constant speed of sound EoSs, Taylor expansion around the nuclear matter point, etc. In [4], a relativistic mean field model parametrization has been adopted which is based on the non linear Walecka model proposed in [5,6]. The main result that we review and further analyze in this contribution is the possibility of obtaining mass–radius relations which are qualitatively very similar to the ones obtained when considering hybrid stars i.e., neutron stars in which, at some density in the inner core, a phase transition to quark matter occurs. In the relativistic mean field model we use, there are no quark degrees of freedom and therefore we cannot speak about deconfinement but we do see, for EoSs optimized by using the Bayesian analysis, a strong reduction of the nucleon effective mass which is reminiscent of a chiral phase transition at finite density.

2. Theoretical Model and Data Set

The mean field Lagrangian that we adopt here contains five parameters which correspond to properties of nuclear matter at saturation. These parameters have priors set by nuclear physics data that will be turned into posteriors by using the astrophysical constraints, within the framework of Bayesian analyses. The five parameters are: the binding energy per nucleon ($E_0$), the saturation density ($n_0$), the symmetry energy (S), the incompressibility (K), and the effective mass of nucleon ($m^{*}$), see [7] and there are five corresponding parameters in the Lagrangian which are the coupling of the $\sigma$, $\omega$, and $\rho$ mesons with nucleons and two additional self-interaction parameters for the sigma field: $g_{\sigma }/m_{\sigma }$, $g_{\omega }/m_{\omega }$, $g_{\rho }/m_{\rho }$, b, c. Interestingly, in the relativistic mean field model that we are using, the relations between couplings and nuclear properties at saturation are analytic (see [7]) and therefore one does not need to generate EoSs of symmetric matter but just the ones needed for neutron stars i.e., beta stable EoSs.

Relativistic mean field models have, as all other models, advantages and drawbacks. The main drawback is related to the limited and somehow arbitrary number and structure of interactions terms, five in the model adopted in this paper (three minimal couplings and two self interactions terms). More terms can be added with couplings fixed by other properties of saturation such as, for instance, the symmetry energy slope L. An improvement would be to adopt models with interactions “inspired” by properties of QCD such as chiral symmetry within relativistic chiral models such as the one studied in Ref. [8] or the parity doublet chiral model studied in Refs. [9,10]. Another point concerns the stability of the mean field potentials: in the model adopted here the parameter c associated with the quartic self interaction can become negative in some parametrizations thus limiting the range of applicability of the model. A standard viewpoint is to consider these models as phenomenological models that need validation from experiments and observations. The predictive power of this class of models is clearly limited, but it is anyway more reliable w.r.t. the predictions and extrapolations of models for the EoS based for instance on piecewise polytropic expansions. Indeed, the main advantage of relativistic mean field models consists of treating on equal footing both symmetric and beta-stable matter, on the possibility to extend the calculation to finite temperatures (which in piecewise polytropic models would instead be based on an adiabatic thermal index arbitrarily chosen) and finally to take into account the possibility of forming hyperons, deltas, and meson condensates. Additionally, one can derive not only the relation between pressure and energy density (which sets the structure of compact stars) but also response functions (heat capacity) and transport properties (shear and bulk viscosity for instance). It is clear then that our approach aims at improving our knowledge on the zero temperature EoS from structural properties of compact stars and then use this knowledge to study, for instance, the merger of compact stars or the supernovae explosions whose simulations require tabulated EoS with a wide range of density, temperature, and proton fraction.

Let us introduce now the priors for our Bayesian analysis. We refer to the recent analyses presented in [11,12] where a complete overview on experimental and theoretical results on nuclear matter properties has been presented. We extract from those two studies a uniform prior and a Gaussian prior with ranges specified in

Table 1. The ranges for the nuclear matter properties for the two priors adopted in this work and extracted from Refs. [11,12].

Priors	Range	${\mathit{m}}^{*}$	K (MeV)	${\mathit{n}}_0$ (${\mathbf{fm}}^{-3}$)	S (MeV)	${\mathit{E}}_0$ (MeV)
Marg_unif	Min	$0.64$	219	$0.145$	$31.19$	$-16.35$
	Max	$0.71$	355	$0.153$	$38.71$	$-16.12$
Marg_Gauss	Mean	$0.67$	268	$0.1494$	$35.11$	$-16.24$
	$\sigma$	$0.02$	34	$0.0025$	$2.63$	$0.06$

.

The aim of the Bayesian analysis is to update the priors (in our case the five properties of nuclear matter at saturation) with new data. The astrophysical data we will exploit are the following. The mass–radius constraints on the thermonuclear bursters and on the quiescent low-mass-X-ray binaries of Ref. [13] (11 sources in total, see Figure 1 and Figure 2, upper panels). We further use the mass–radius constraint obtained from the X-ray burst cooling spectra of 4U 1702–429 Nättilä et al. [14] for which $M=1.9\pm 0.3M_{\odot }$ and $R=12.4\pm 0.4$ km, and the mass radius constraint derived by the NICER collaboration from the pulse profile modelling of PSR J0030+0451: $M={1.34}_{-0.16}^{+0.15}{M}_{\odot }$, $R={12.71}_{-1.19}^{+1.14}$ km [15] (One could have used also the constraints of Ref. [16] as in Ref. [17]. We plan to do it in a future work.). For both sources we use, following [18], a bivariate Gaussian distribution (see Figure 1 and Figure 2 lower left panels). Finally, we use the mass–radius constraint as extracted from the gravitational waves signal GW170817 [1] (see Figure 1 and Figure 2 lower right panels). Let us briefly review how does the Bayesian analysis work, see Refs. [19,20,21,22,23,24] for previous applications of this kind of tool to gravitational waves (GW170817) and to electromagnetic signals from neutron stars as detected by NICER. We indicate with $P\left(\theta \right|M)$ the prior probability of the parameters set (i.e., the uniform or Gaussian distributions of the five parameters listed in Table 1). For each parameter set, we can generate first the couplings of the mean field model, then the beta stable EoSs and finally the mass–radius relations. We construct then the posterior probability $P\left(\theta \right|D,M)$ (once a specific data set has been used) as

$$P\left(\theta \right|D,M)\propto P(D|\theta ,M)P\left(\theta \right|M),$$

(1)

where $P\left(D\right|\theta ,M)$ is the likelihood. From a practical point of view, one has to sample this multidimensional distribution. For this aim, we use a Markov-Chain Monte Carlo (MCMC) tool within the python emcee package with a stretch-move algorithm [25].

3. Results

In Figure 1 and Figure 2, we show the results of the Bayesian analysis for the mass–radius relation of neutron stars. We display, for the two priors adopted in this work, the bunch of M-R relations corresponding to the $68\%$ Confidence Interval (CI). The black dashed lines refer to the most probable EoS, namely the one with the highest joint posterior probability. There are two remarkable facts that one can notice when comparing the case of uniform prior and Gaussian prior. In the former, it is clearly difficult to obtain compact configurations: the radius of the $1.4M_{\odot }$ stars, $R_{1.4}$, is larger than about $13.5$ km and this makes the agreement (even at $90\%$) with the constraints obtained by the gravitational waves emitted by GW170817 difficult (see lower right panel of Figure 1). Similarly, it is also difficult to reconcile this bunch of EoSs with the data on bursters (upper right panel of Figure 1). It is instead possible to fulfill the constraints of qLMXBs (upper left panel, which have however large error bars) and the large stellar configurations inferred for the two sources indicated in the lower left panel. In

Table 2. Most probable empirical parameters of the joint posterior distribution for the two priors. We also report the values of the slope of the symmetry energy at saturation L.

Models	${\mathit{m}}^{*}$	K (MeV)	${\mathit{n}}_0$ (fm${}^{-3}$)	S (MeV)	${\mathit{E}}_0$ (MeV)	L (MeV)
Marg_unif	$0.710$	$219.1$	$0.152$	$31.3$	$-16.26$	$90.9$
Marg_Gauss	$0.713$	$163.8$	$0.150$	$34.5$	$-16.23$	$100.6$

, we report the most probable values of the five nuclear physics parameters. In the case of uniform prior, the need to obtain massive enough compact stars (EoSs with maximum mass smaller than $2M_{\odot }$ are excluded from the sample) forces the effective mass $m*$ to be of the order of $0.7$ (in units of nucleonic mass), whereas the need to obtain also small configurations (such as the ones associated by GW170817) forces the incompressibility K to sit at the lower edge of the uniform prior distribution, $K\sim 220$ MeV. The importance of $m*$ for increasing the maximum mass and of K for increasing the compactness of small mass stars have been reported in [26,27], respectively. When allowing for a wider exploration of the parameter space within the Gaussian prior, we notice a much satisfying agreement with the observational constraints. Indeed, now $R_{1.4}$ can be as small as ∼12.5 km while the maximum mass is still above $2M_{\odot }$. By comparing with the previous case, now K has dropped to ∼164 MeV while $m*$ is still ∼0.7. While a small value of K is inferred also from the analysis of the data of the KaoS collaboration [28], the value obtained within our analysis is $3\sigma$ below the central value as obtained in the Gaussian distribution of [11,12]. It is very possible that this result is an indication of a tension between our theoretical model and astrophysical data that could be eventually solved by using different terms in the Lagrangian. However, there is another point that is actually suggesting that this tension could have a physical explanation. One can notice in the mass–radius relations of Figure 2 (see the dashed line) that a smooth kink appears in the curve which is very similar to the one obtained when allowing for the formation of a mixed phase between nucleons and a new phase, e.g., the quark phase. Qualitatively, that curve is very similar to the curve of hybrid stars with an early onset of the phase transition and with a mixed phase in which pasta structures appear, see, for instance, [29]. We remark that this is just a similarity, here there are not quark degrees of freedom. This interpretation is supported by the analysis on the behavior of the nucleon effective mass as a function of the density, displayed in the upper panel of Figure 3. By comparing the red and green lines corresponding to the uniform and Gaussian priors respectively, one notices that in the latter case there is a fast drop of $m^{*}$ at a density of about twice the saturation density. This fast drop affects both the mass radius relations of Figure 2 and the relation between the mass and the central density shown in the lower panel of Figure 3. The drop of $m^{*}$ flattens this curve, which then starts to increase with a larger slope above a density of about three times the saturation density. A smooth phase transition to a chirally restored phase would display a similar behavior. Notice that such a nonlinear behavior of the EoS cannot be obtained within analyses based on piecewise polytropes or on Taylor expansions around saturation. In this respect, using a nonlinear relativistic mean field model allows for capturing possible phase transitions in dense matter associated e.g., with the appearance of new degrees of freedom or to a change of the order parameter associated with chiral symmetry. Notice that this behavior depends on the choice of the parameters and not on the choice of the Lagrangian which is the same for the two cases discussed here.

4. Conclusions

In this contribution, we have presented a first attempt to perform a Bayesian analysis on the EoS of neutron stars based on a class of relativistic mean field models. The most striking feature of this model is its nonlinearity in the $\sigma$ field which allows for a phenomenological analysis of possible phase transitions in the same spirit of the Landau–Ginzburg approach. Indeed, we find that the astrophysical data, which seem to suggest in many cases the existence of quite compact stellar objects, “force” our Monte Carlo to prefer solutions with a rapidly decreasing nucleon effective mass at increasing density. It is therefore tempting to infer a possible phase transition to a chirally restored phase. Indeed, the chiral condensate can be used as an order parameter in dense matter and as the density increases, see the predictions of the chiral parity doublet model [9,10], it reduces and leads to a first order or to a smoother (second order or crossover) chiral phase transition which could in principle have an effect on the structure of compact stars.

This preliminary study needs to be extended by including new and/or different interaction terms which allow for better constraining the density dependence of symmetry energy [30] and the value of the incompressibility [31]. In this respect, a further analysis using density dependent couplings models would be valuable [32,33]. Another important extension of this work, partially studied in [4], is the inclusion of heavier baryons such as hyperons and deltas since their possible formation in the core of neutron stars has been confirmed by a number of theoretical studies; see, e.g., [34,35,36,37,38,39]. Finally, the possibility of a strong first order phase transition to quark matter, in the twin star scenario [40,41,42,43,44] or the two-families scenario [45,46,47,48,49], should also be studied within a Bayesian approach with the aim of comparing the “evidence” (or Bayes factor) of the standard one-family scenario and the scenarios of disconnected mass–radius relations [50,51].