The equation of state of neutron star matter based on the G -matrix and observations

G -matrix theory, which has a sound basis both theoretically and phenomenologically in this region. In addition, a phenomenological third-order term of baryon density is introduced to control the stiffness of the EOS in this density region. In the higher-density region, we introduce a parameterized EOS. The adjustable parameters are ﬁxed utilizing the statistical method developed by Steiner et al. [Astro-phys. J. 722 , 33 (2010)] to be consistent with observational data on neutron stars. As a result, we ﬁnd that an EOS softened by the additional third-order term of baryon density in the density region lower than 3.5–4.0 times the normal density and a stiff EOS in the higher-density region are preferable. The resultant EOS is similar to an EOS with the assumption of the hadron–quark crossover proposed by Masuda et


Introduction
The average densities of typical neutron stars are about two times the nuclear density (≈ 3.0 × 10 14 g/cm 3 ). Therefore, an explanation of the equation of state (EOS) of neutron stars is one of the fundamental problems in nuclear physics. In fact, many nuclear physicists have attempted to explain the masses of neutron stars with an EOS based on realistic baryon-baryon interactions [1][2][3]. In this context, the hyperon mixing in neutron stars causes too soft an EOS and fails to explain the masses of typical neutron stars with 1.4 solar mass (M ). Further, the recent discovery of the heavy neutron stars PSR J1614+2230 [4] with mass M = 1.97 ± 0.04 M and PSR J0348+0432 [5] with mass M = 2.01 ± 0.04 M has caused great difficulties with this problem. In addition, the appearance of new degrees of freedom softens the EOS and the maximum mass of neutron stars is reduced considerably. Some authors call this difficulty the hyperon puzzle [6]. The hyperon puzzle can be solved in principle by constructing the EOS to give the maximum mass of neutron stars around 2 M . The solution of the hyperon puzzle can be achieved in frameworks based on the model with quarkmeson coupling [7][8][9] and the vector baryon-meson coupling model [10]. The density-dependent relativistic mean-field model involving meson-hadron coupling constants and meson masses also leads to the stiffening of the EOS and a large neutron star mass [11]. In high-density matter, physical effects such as many-body forces, boson condensations, or effects of quark degrees of freedom are expected to be important.

Equation of state
The EOS of nuclear matter is an important ingredient to determine the M-R relation for neutron stars. Using realistic baryon-baryon interactions, we determine the EOS of high-density β-stable baryonic matter and calculate the M-R relation for neutron stars by solving the TOV equation. We consider three regions in neutron stars as illustrated in Fig. 1.
The first region is the crust of neutron stars. Although the mass of the neutron star crust constitutes only ≈ 1% of the neutron star mass, and its thickness is typically less than one-tenth of the star radius, the neutron star crust plays an important role in determining the M-R relation. In this region, we use the Bogomol'nyi-Prasad-Sommerfield (BPS) EOS [16] and its extrapolation up to the transition baryon density ρ crust = ρ 0 /2, where ρ 0 = 0.17 fm −3 is the nuclear saturation baryon density.
The second region is the theoretical EOS region defined by the baryon densities from ρ 0 /2 to nρ 0 , where n is a variable. This region is assumed to be dominated by baryon-baryon interactions. In this Fig. 1. Pressure p as a function of energy density . In the high-density region, we use a parameterized EOS, which is the general piecewise linear function. We take ρ crust = ρ 0 /2, 2.5 ≤ n ≤ 5, and 0 < ν 1,...,4 ≤ 1.

2/18
Downloaded from https://academic.oup.com/ptep/article-abstract/2016/7/073D02/1752908 by guest on 27 July 2018 region, we assume an energy density given by = theor (ρ) + 3 ρ 3 for ρ crust < ρ < nρ 0 (1) where and ρ are the energy density and the baryon density, respectively. theor is determined by the G-matrix theory with baryon-baryon interactions and the 3 term is introduced as a phenomenological third-order term of baryon density. The 3 term does not directly mean the three-body force effect. The aim of the 3 term is not to construct a consistent model describing both symmetric nuclear matter and neutron star matter, but to control the stiffness of the neutron star matter EOS in the second region. In this work, we assume that the saturation properties of symmetric nuclear matter are not affected by the 3 term. It is well known that, for a given neutron star mass, the neutron star radius depends on how stiff or soft the EOS is. Steiner et al. estimated that the radius of a neutron star with mass 1.4 M is between 10.4 and 12.9 km [17]. Because neutron star radii are mainly decided by the theoretical EOS in the second region, the 3 term plays a role in controlling neutron star radii (see Sect. 3).
To determine theor , we perform the G-matrix calculation [18] for β-stable baryonic matter with two models of baryon-baryon interactions, NSC97e [19] and fg2014 [20]. Both models have similar properties for NN and N interactions but give very different predictions for N interactions. The former provides an attractive − -nuclear matter interaction, which causes − mixing in neutron star matter at relatively low densities (around 2 ρ 0 ). The latter provides a repulsive interaction and does not cause − mixing in the second region (n < 5). Recent experimental and theoretical knowledge in hypernuclear physics support the repulsive − -nuclear matter interaction. As a result, NSC97e gives a softer EOS than fg2014 at densities higher than two times the nuclear density. Because NSC97e and fg2014 give relatively soft EOSs, both of them produce neutron stars with small radii. To compare with these two models, we also consider the EOS model K=240 (GM3) [21], which gives larger radii for neutron stars than the NSC97e and fg2014 models. If we employ other models to determine the theoretical EOS in the second region, we may obtain quantitatively different results. This means that it may be possible for us to decide the EOS in this region by using a more accurately observed M-R relation in the future.
To determine an adequate theoretical EOS of neutron star matter, we must use the baryon-baryon interaction, which reproduces the symmetry energy at around the saturation density (ρ 0 ). Therefore, we calculate the saturation baryon density (ρ s ), the binding energy (B), the symmetry energy (S v ), and its derivative (L), which are defined by where x = ρ p /ρ denotes the proton fraction (x = 1/2 corresponds to symmetric nuclear matter) and m p (m n ) is the proton (neutron) rest mass. Using the NSC97e and fg2014 models, we obtain (ρ s [fm − where Using NSC97e and fg2014, we obtain (ρ, 0) as shown in Fig. 2 (30.8, 34.6) for NSC97e and fg2014, respectively. S 0 and L 0 are not the same as S v and L defined at x = 1/2. Since neutron star matter is very similar to pure neutron matter, we employ S 0 and L 0 as criteria of the properties at the saturation density. When the contribution from the 3 term is taken into account, we obtain the result shown in Fig. 3. From this figure, to ensure that the symmetry energy does not become too small, we adopt the condition 3 > −0.2 and 3 > −0.1 for NSC97e and fg2014, respectively. From Figs. 2 and 3, we find that NSC97e gives slightly a stiffer EOS than fg2014 around the saturation density. As mentioned above, at higher densities (a few times the saturation density), NSC97e gives a softer EOS than fg2014 because of the effect of hyperon mixing. If we extend these EOSs to high densities, they cannot support a neutron star of 2 M . Various effective approaches have been used to reproduce 2 M neutron stars, such as assuming the hadron-quark crossover [12] and the universal repulsive three-body force effect [13,14]. In this paper, we extend the EOS to higher densities (the third region) using piecewise linear functions.
The third region is that with densities higher than n times the saturation density. This region is described by an EOS satisfying the causality condition (dp/d < 1). However, the inner structure of neutron stars is still unknown due to the theoretical and observational uncertainties. Theoretically, the inner cores of neutron stars with very high density are expected to be composed of various exotic particles, such as pions, kaons, hyperons [21], and strange quark states [24,25]. However, the EOSs  of dense matter beyond the nuclear density are still quite uncertain in particle and nuclear physics. Therefore, for energy densities above 0 = (nρ 0 ), we assume the pressure p( ) as a parameterized piecewise linear function of the energy density given by where ν i is a constant slope. We use four linear functions (N = 4) with slopes ν 1,...,4 , which make it possible to vary the stiffness of the EOS. In numerical calculations, we have confirmed that N = 4 is robust and the results are unchanged even if additional linear functions are introduced. Note that we parameterize the high-density EOS as a function of the energy density . We choose the transition baryon density nρ 0 = (2.5-5.0) ρ 0 and p 0 is determined by the continuity of p( ) at the transition density 0 = (nρ 0 ). We vary the constant slope parameters over the ranges ν i−1 ≤ ν i ≤ 1. Thus we have a total of six EOS parameters n, 3 , ν 1,...,4 . The parameter is large enough to ensure that the EOS in the high-density region could be parameterized reasonably well. It is convenient to remember that the densities 3-5 ρ 0 correspond to the energy densities 450-800 MeV fm −3 . For ≤ 0.3 fm −4 , it is impossible to parameterize the EOS in the high-density region. For ≥ 0.5 fm −4 , we obtained the result with similar values of ν 3 and ν 4 near the causality limit. Therefore, ν 4 becomes redundant because this is the same as the three-line-segment case. Hence, we choose = 0.4 fm −4 . By extending the EOS to higher densities in this way, we show that it is possible to derive constraints on the EOS and on the M-R relation.

Mass-radius relation
For a given EOS, masses and radii of neutron stars can be determined as functions of central pressure (or central energy density) by solving the TOV equation [26,27] using the EOS. The TOV equation is given by This is a system of simultaneous first-order ordinary differential equations for the pressure p(r), including the mass m(r) with the EOS p( ). We can solve it by integrating from r = 0 with the initial values where c is the central energy density and p c is the central pressure. The neutron star radius R is given by the condition of p(R) = 0, and the neutron star mass M is defined by M = m(R). The piecewise linear function described in the previous section and 3 are used to generate the EOS. Using each EOS, we determine M and R as functions of p c . Plotting M and R for various p c , we obtain an M-R relation for each EOS. If we make the theoretical EOS stiffer in the second region, we obtain larger radii but the maximum mass does not change much, as shown in Fig. 4. This means that accurate determination of neutron star radii is very important to clarify the properties of the EOS around ρ 0 < ρ < nρ 0 (the second region). In Fig. 5, we show the dependence of the M-R relation on the parameterized EOS. By assuming that the parameterized EOS is stiffer, we can obtain a larger maximum mass for neutron stars, but the radii do not change much.

Application of statistical methods to constraints on the EOS
The mass-radius (M-R) relation is important to determine the EOS of neutron star matter. However, at present, constraints on masses and radii of neutron stars are uncertain. Moreover, the number of neutron stars for which the mass and radius have been measured is very limited. Therefore, we apply a statistical method to obtain the constraints on the EOS of neutron star matter. In this work, we employ the Bayesian statistical method introduced by Steiner et al. [15] for using neutron star observational data to probe suitable EOSs. Based on these analyses, it will be possible to conclude which EOS is the most probable.

Bayesian analysis
The Bayes theorem shows us how to obtain the quantity of the conditional probability P(M|D) of the model M given the data D. The Bayes theorem [28] is written as follows: In Eq. (11) the masses and radii of neutron stars. Therefore, we need an additional model parameter, e.g., the central pressure p c to uniquely determine the mass and radius for each neutron star. In this treatment, p c for each l is a function of M l and the six EOS parameters. Therefore, M(six EOS parameters, p c for l = 1, . . . , 9) is equivalent to M(six EOS parameters, M 1 , . . . , M 9 ). On the other hand, p c depends strongly on the EOS used and its lower and upper bounds are unclear. Therefore, the p c for l = 1, . . . , 9 are unsuitable as model parameters. This is the reason why we choose to treat the {M l } as model parameters. Substituting in Eq. (12), we have where N = N p + N M = 6 + 9 = 15 is the dimension of our model space, where N p is the total number of EOS parameters and N M is the total number of neutron stars in our data set. Table 1 shows the masses and radii of the 9 observed neutron stars that we use in our paper. In order to apply Eq. (13) to our problem, we assume that P(D|M) is proportional to the product over the probability distributions D l evaluated at the masses that are determined in the model and at the radii that are determined from the model M, i.e., In our calculation, Eq. (13) can be rewritten as where we assume that the prior probability P(M) is uniform under several conditions for model parameters: 2.5 < n < 5, ν i ≤ ν i+1 < 1, 3 > −0.2 for NSC97e (−0.1 for fg2014), and supporting 2 M . For the data D l , we use the probability distributions D l (l = 1, . . . , 9) derived from the 9 neutron star observations listed in Table 1. In this work, we assume that none of the probability distributions D l (M , R) of the 9 neutron stars has probability outside of the ranges M low < M < M high

Observed masses and radii of neutron stars
Observations of mass and radius information were obtained from astrophysical observations of X-ray bursts and thermal emissions from quiescent low-mass X-ray binaries (LMXBs). When neutron stars pull material away from companion stars, they can become much brighter. Using observation of Xrays at different wavelengths, combined with theoretical models of neutron star atmospheres, we can estimate the relationship between the radius and mass of neutron stars. This work has been performed by Heinke [32], Webb and Barret, as explained in Ref. [33], and Guillot in Ref. [34]. All of these observations were done for neutron star binaries in globular clusters. Because of thermonuclear explosions on their surfaces, the atmospheres of neutron stars expand. If observers catch one of these bursts, they can calculate its surface area based on the cooling process. After that, when this calculation is combined with an independent estimate of the distance to the neutron star, the mass and radius of this star can be estimated. Ozel and Guver have applied this technique in Refs. [29][30][31]35].
The papers referred to above provide information about the neutron star's M-R relation and we use this information to construct probability distributions D l . In detail, the probability distributions D l (M , R) for 4U 1608-52 (l = 1) and 4U 1820-30 (l = 3) are described as the Gaussian distribution with the values given in Table 2.
For EXO 1745-248 (l = 2), we use where the values of the parameters are given in Table 3. The values of A l in Eqs. (16) and (17)    we estimated the probability distributions shown in Figs. 6, 7, and 8. Our probability distributions are similar to those given by Steiner et al. [15] but not the same.

Results from the statistical analysis
In Table 4, we summarize the results of the conditional probability − log P[M|D] and the EOS parameters. A smaller − log P[M|D] corresponds to a better fit, i.e., the EOS used is more probable.
We can see that, for all models of the EOS, good fits are obtained with small 3 and large slopes ν 3 , ν 4 . This means a soft EOS in the theoretical region, stiff in the high-density region. For each model of the theoretical EOS, corresponding to the best cases of − log P[M|D], we draw the M-R relation to compare with probability distributions of observed neutron stars in Figs. 9, 10, and 11. Our results support the suggestion that neutron stars with mass 2.0 M have small radii 9-10 km. Figures 12 and 13 show the EOS bands consisting of EOSs that fulfill the following three conditions: (1) the EOS supports neutron stars with masses larger than 2.0 M . In the region = 150-600 MeV/fm 3 , the upper limit of the EOS band is determined by condition (2) and the lower limit by condition (3). Condition (1) provides a stiff EOS in the region ≥ 650 MeV/fm 3 . Because of bad results, we discard and do not draw the EOS band for the K=240 model. 10     The reason why the K=240 model has bad results is because most of the neutron star observations used in this work have small radii (see Figs. 9, 10, and 11). Because the K=240 model gives a stiff EOS in the second region, this model supports 1.4 M neutron stars with large radii (∼ 13 km). On the other hand, most of the neutron star observations that we use have small radii (∼ 10 km). In the future, if the radii of neutron stars precisely determined in observations are larger than those used here, this result will surely change. We also discard all the cases with n ≥ 4.5 because they cannot support neutron stars with M ≥ 2.0 M unless large values of 3 (≥ 1), which give bad − log P[M|D], are used. When looking at the contours of the distribution probabilities D l (M , R), most of them support neutron stars with small radii, except X7 in 47 Tuc. Because the differences between the radii of X7 and the others are nearly 5 km, it is impossible to obtain an M-R relation that fits all 9 distribution probabilities. Also for this reason, the − log P [M|D] values are somewhat bad (the best value that we found is ∼ 20). Moreover, the widths of the contours are large and two separate peaks for KS 1731-260, U24, EXO 1745-248 are found. These uncertainties allow many cases of EOSs with similar − log P[M|D].

12/18
Downloaded from https://academic.oup.com/ptep/article-abstract/2016/7/073D02/1752908 by guest on 27 July 2018     Finally, we consider the case of n = 4. We can see that, with negative values of 3 , it is impossible to support neutron stars with mass 2 M though a very stiff parameterized EOS is used. To compare with observational data, the M-R relations are shown in Figs. 9, 10, and 11. These M-R curves suggest that a relatively small radius (9-10 km) is consistent with our above remarks on the soft theoretical EOS. The EOS, which is stiff in both the theoretical and the parameterized EOS regions, may not be denied. However, this type of EOS cannot give high probability (− log P[M|D] > 25) because of large radii and is disfavored by the present 9 observational pieces of data used in this work. More neutron star observations with mass and radius constraints would enable us to improve our results.
The most probable EOS bands are shown in Figs. 12 and 13. In these figures, it is remarkable that, for both NSC97e and fg2014, the EOS bands rapidly become stiff at the energy density ∼ 600-650 MeV/fm 3 (around 3.5-4.0 times saturation density). Various effective approaches have been used to reproduce 2 M neutron stars, such as assuming the hadron-quark crossover [12] and the universal repulsive three-body force effect [13,14]. If the G-matrix EOS is adopted around the normal nuclear density, a stiff EOS is preferable in the high-density region with 2 ρ 0 so as to account for the existence of 2 M neutron stars. As well as making the EOS stiffer at high 15/18 Downloaded from https://academic.oup.com/ptep/article-abstract/2016/7/073D02/1752908 by guest on 27 July 2018 densities, many phase transitions have been considered, such as superfluid transitions [36], kaon or pion condensates [37,38], hyperon matter [39], etc. Eventually, phase transitions from nuclear matter to quark matter [40,41] are expected at very high densities. In general, the EOS becomes soft after the phase transition. The behavior of our EOS is in contrast to the general phase transitions but is similar to that assuming the hadron-quark crossover [12], which leads to a stiffening of the EOS. However, we note that our results are obtained based on Bayesian analysis and do not depend on specific physical assumptions about neutron star matter. Therefore, this behavior of our EOS must be confirmed in the future. Below the transition density, our most probable EOS bands are softer than the EOS bands suggested by Steiner et al. [15] but, above the transition density, our EOS bands are stiffer. Our EOS bands support neutron stars with radii 9-10 km whereas Steiner et al.'s support neutron stars with radii 11-12 km [15]. This difference comes from the differences between our input probability distributions and those of Steiner et al. While Steiner et al. performed their calculations with six observational pieces of data, we add data on NGC2808, U24, and KS 1731-260 to our calculation. Our probability distributions of 4U 1608-52 and 4U 1820-30 are constructed based on central masses, central radii, and their uncertainties, which were determined by Guver et al. [29,30]. On the other hand, using the Monte Carlo method, Steiner et al. [15] constructed them based on their own calculations of Eddington fluxes, angular areas assuming different parameters from Guver et al. For this reason, in the work by Steiner et al., the 4U 1608-52 and 4U 1820-30 probability distributions consist of 1 or 2 peaks (depending on the photospheric radius assumption) while, in our work, they consist of only one peak. It is important to note that, if the probability distributions as input change, the results would change accordingly. In Fig. 14, we show calculated profiles of neutron stars with masses M = 1.48 M and 2.0 M . For heavy neutron stars with mass 2.0 M , we find that a small radius (R < 10 km) implies very high energy densities in the central region, strongly depending on the EOS in the second region (ρ 0 < ρ < nρ 0 ). In the case of the canonical

Summary
By combining nuclear matter theory and astrophysical observations, we have constructed the EOS of neutron star matter. By using realistic baryon-baryon interactions, we determined the EOS of highdensity β-stable baryonic matter and calculated the M-R relation of neutron stars by solving the TOV equation. We have introduced three regions in neutron stars. The first region is the crust of neutron stars. The second region is dominated by baryon-baryon interactions with densities lower than n times the saturation density, where n is a variable. In this region, we introduce a phenomenological thirdorder term of baryon density. The third is the region with densities higher than n times the saturation density. This region is described by a parameterized EOS satisfying the causality condition. The adjustable parameters are fixed utilizing the statistical method developed by Steiner et al. [15] to be consistent with observational data on neutron stars.
The recent discovery of neutron stars with 2 M poses many challenges to theoretical physics. In general, the EOS based on the G-matrix theory is thought to be too soft to account for even 1.4 M neutron stars. However, we have indicated that neutron stars with mass 2 M can be reproduced by stiffening of the EOS at high densities. By applying the Bayesian statistical method, we have found that an EOS softened by the additional third-order term of baryon density in the second region and a stiff EOS in the third region are preferable. These EOSs lead to a constraint on the symmetry energy of S 0 ≤ 32.08 MeV for the NSC97e model (31.47 MeV for fg2014) and its density derivative L 0 ≤ 44.08 MeV for the NSC97e model (39.39 MeV for fg2014). In addition, based on the most probable EOS bands, we predict a rapid change of stiffness around 3.5-4.0 times the saturation density. The behavior of our EOS is in contrast to the general phase transitions but is similar to that assuming the hadron-quark crossover [12], which leads to a stiffening of the EOS.