PSR J0030+0451, GW170817, and the Nuclear Data: Joint Constraints on Equation of State and Bulk Properties of Neutron Stars

Parameterized representations of the EoS can reasonably describe the main properties of the dense matters in NSs (see Baiotti 2019 for a review). Currently there are two widely adopted methods to parameterize the EoS, namely, the piecewise polytropic expansion method (Vuille & Ipser 1999; Read et al. 2009; Lackey & Wade 2015) and the spectral decomposition method (Lindblom 2010; Carney et al. 2018).

In this work we take the same piecewise polytropic expansion method and parameter settings used in Jiang et al. (2019), along with the same pressure-based spectral decomposition method used in Annala et al. (2018) to parameterize the EoS. The piecewise polytropic expansion method uses four pressures $\{{P}_{1{{\rm{n}}}_{{\rm{s}}}},{P}_{1.85{{\rm{n}}}_{{\rm{s}}}},{P}_{3.7{{\rm{n}}}_{{\rm{s}}}},{P}_{7.4{{\rm{n}}}_{{\rm{s}}}}\}$ located, respectively, at four different rest mass densities $\{{\rho }_{1{{\rm{n}}}_{{\rm{s}}}},{\rho }_{1.85{{\rm{n}}}_{{\rm{s}}}},{\rho }_{3.7{{\rm{n}}}_{{\rm{s}}}},{\rho }_{7.4{{\rm{n}}}_{{\rm{s}}}}\}$ to parameterize the EoS, where ${n}_{{\rm{s}}}$ is the nuclear saturation density. Between each pair of adjoining bounds of densities, we approximately take a polytropic form. While the spectral decomposition method uses four expansion parameters {γ₀, γ₁, γ₂, γ₃} to describe the relation between the pressure and the adiabatic index, and thus uniquely determine an EoS.

Nuclear experiments and gravitational-wave data have been proved to be essential in constraining the EoS (e.g., Lattimer & Steiner 2014; Tews et al. 2017; Abbott et al. 2018; Annala et al. 2018; De et al. 2018; Most et al. 2018). As summarized in Section 2.4 of Jiang et al. (2019), the current nuclear data have set interesting constraints on both ${P}_{1{{\rm{n}}}_{{\rm{s}}}}$ and ${P}_{1.85{{\rm{n}}}_{{\rm{s}}}}$, which are $4.70\geqslant {P}_{1{{\rm{n}}}_{{\rm{s}}}}/({10}^{33}\mathrm{dyn}\,{\mathrm{cm}}^{-2})\geqslant 3.12$ and ${P}_{1.85{{\rm{n}}}_{{\rm{s}}}}$/$({10}^{34}\mathrm{dyn}\,{\mathrm{cm}}^{-2})\geqslant 1.21$, respectively. We also follow that paper (Section 2.5 therein) to take into account the constraints set by the gravitational-wave data of GW170817. However, in this work we do not consider the low-mass X-ray binary data further because the latest NICER measurement of PSR J0030+0451 is much more direct and suffers from significantly fewer systematic uncertainties.

The NICER teams have measured the mass and the radius of PSR J0030+0451 through pulse profile modeling methods. Riley et al. (2019a) reported a mass (radius) ${1.34}_{-0.16}^{+0.15}{M}_{\odot }$ $({12.71}_{-1.19}^{+1.14}\mathrm{km})$, while Miller et al. (2019a) reported a mass (radius) ${1.44}_{-0.14}^{+0.15}{M}_{\odot }$ $({13.02}_{-1.06}^{+1.24}\mathrm{km})$, which are consistent with each other. We mimic the mass–radius posterior distribution of this source by a two-dimensional Gaussian kernel density estimation (KDE) to constrain the EoS:

$$\begin{eqnarray}&&P(M,R)=\mathrm{KDE}(M,R| {\boldsymbol{S}}),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{S}}$ is posterior samples taken from the best fitting three ovals case of Miller et al. (2019b) and the ST+PST case of Riley et al. (2019b). Unless specified, we only show the results with the data of Riley et al. (2019b), because the data of Miller et al. (2019b) yield rather similar results, as shown in Table 2. Besides, we do not consider the impact of spin on the radius because Raaijmakers et al. (2019) have shown that the spin effect of this source is small in comparison to the rather large systematic uncertainties.

In addition to satisfying the above bounds, the EoS should meet the following general requests: (i) the causality condition, i.e., the speed of sound of the dense matter can never exceed the speed of light c; (ii) the microscopical stability condition, i.e., the pressure cannot be smaller in denser matters; (iii) the maximal gravitational mass of nonrotating NSs (M_TOV) should be above all those accurately measured. Note that for the rapidly rotating NSs, the gravitational mass has been enhanced. For PSR J0740+6620, the enhancement of the gravitational mass is about 0.01M_⊙ (Breu & Rezzolla 2016; Ma et al. 2018). Thus we have a robust lower limit of ${M}_{\mathrm{TOV}}\geqslant 2.04{M}_{\odot }$ based on the mass measurement of PSR J0740+6620 (Cromartie et al. 2020). Additionally, the adiabatic index Γ(p) for the spectral decomposition method are limited in the range [0.6, 4.5] (Abbott et al. 2018).

As mentioned in Section 2.1, two parameterization methods are adopted in this work. To investigate the impact of M_TOV on constraining the EoS, we consider the possible regions of (2.04, 2.40) M_⊙ (i.e., the low M_TOV case) and (2.40, 2.90) M_⊙ (i.e., the high M_TOV case), respectively.³ With these considerations we carry out four tests:

• (i)  
  Using the piecewise method to parameterize the EoS, assuming ${M}_{\mathrm{TOV}}\in (2.04,2.40){M}_{\odot }$.
• (ii)  
  The same as (i), except assuming ${M}_{\mathrm{TOV}}\in (2.40,2.90){M}_{\odot }$.
• (iii)  
  Using the spectral method to parameterize the EoS, assuming the low M_TOV case.
• (iv)  
  The same as (iii), except assuming the high M_TOV case.

To investigate the role of prior in each test, we also construct the prior for piecewise method and spectral method, where we just take into account the causality condition, the microscopical stability condition, and the lower limit M_TOV > 2.04M_⊙, while in the spectral case the adiabatic index is assumed to be within the range [0.6, 4.5].

Since we combine the data of GW170817 and PSR J0030+0451 to do the analyis, the likelihood in each test will take the form

$$\begin{eqnarray}&&L\propto \exp [-2{\int }_{0}^{\infty }\displaystyle \frac{| \tilde{d}(f)-\tilde{h}(f;{{\boldsymbol{\theta }}}_{\mathrm{gw}}){| }^{2}}{{S}_{n}(f)}\,{df}]\times P(M,R),\end{eqnarray} \tag{ 2 }$$

where ${{\boldsymbol{\theta }}}_{\mathrm{gw}}=\{{{ \mathcal M }}_{{\rm{c}}},q,{\chi }_{1},{\chi }_{2},{\theta }_{\mathrm{jn}},{t}_{{\rm{c}}},{\rm{\Psi }},{{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2}\}$ are the gravitational-wave parameters. The ${{ \mathcal M }}_{{\rm{c}}}$, q, χ_i, θ_jn, t_c, Ψ, and Λ_i are chirp mass, mass ratio, spin of the ith neutron star, inclination angle, coalescence time, polarization, and tidal deformability of the ith neutron star, respectively. The $\tilde{d}(f)$, $\tilde{h}(f)$, and S_n(f) are frequency domain gravitational-wave data of GW170817, frequency domain waveform, and power spectral density of the GW170817 data, respectively. The two tidal deformabilities of source NSs of GW170817 and the mass/radius of PSR J0030+0451 are determined by

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Lambda }}}_{1} & = & {{\rm{\Lambda }}}_{1}({{\boldsymbol{\theta }}}_{\mathrm{eos}},{M}_{1}),\,\,\,{{\rm{\Lambda }}}_{2}={{\rm{\Lambda }}}_{2}({{\boldsymbol{\theta }}}_{\mathrm{eos}},{M}_{2}),\\ M & = & M({{\boldsymbol{\theta }}}_{\mathrm{eos}},{p}_{{\rm{c}}}),\,\,\,R=R({{\boldsymbol{\theta }}}_{\mathrm{eos}},{p}_{{\rm{c}}}),\end{array}\end{eqnarray} \tag{ 3 }$$

where ${{\boldsymbol{\theta }}}_{\mathrm{eos}}$ are parameters needed to describe an EoS, in the case of the piecewise method, ${{\boldsymbol{\theta }}}_{\mathrm{eos}}=\{{P}_{1{{\rm{n}}}_{{\rm{s}}}},{P}_{1.85{{\rm{n}}}_{{\rm{s}}}},{P}_{3.7{{\rm{n}}}_{{\rm{s}}}},{P}_{7.4{{\rm{n}}}_{{\rm{s}}}}\}$, whereas in the case of spectral method, ${{\boldsymbol{\theta }}}_{\mathrm{eos}}=\{{\gamma }_{0},{\gamma }_{1},{\gamma }_{2},{\gamma }_{3}\}$. The p_c is the central pressure of PSR J0030+0451. The M₁ and M₂ are evaluated from ${{ \mathcal M }}_{{\rm{c}}}$ and q. Since we marginalize the distance and phase in analyzing the gravitational wave, there are only 12 parameters in total, namely, ${\boldsymbol{\theta }}={{\boldsymbol{\theta }}}_{\mathrm{gw}}\cup {{\boldsymbol{\theta }}}_{\mathrm{eos}}\cup {p}_{c}$, and we sample these parameters using the PyMultiNest code (Buchner et al. 2014; Buchner 2016) implemented in Bilby (Ashton et al. 2019a, 2019b), while the likelihood contribution of the gravitational wave is calculated using PyCBC (The PyCBC Team 2018; Biwer et al. 2019).

We construct the P − ρ relation for each posterior sample of ${{\boldsymbol{\theta }}}_{\mathrm{eos}}$ (Figure 1(a)). The narrow pressure uncertainty at ρ_sat is mainly governed by the nuclear constraint, and the higher density part between ρ_sat and 3ρ_sat is mainly determined by the gravitational-wave data and the M − R data of PSR J0030+0451, while the highest density region is mainly constrained by the limits on M_TOV. Meanwhile, for the two investigated M_TOV regions, the EoSs at high densities are slightly different, because the more massive the compact star, the stiffer the EoS for the dense matter.

The character of dense matter can also be described by the sound speed ${c}_{{\rm{s}}}=\sqrt{{dp}/d\epsilon }$, for which there are two interesting limits ${c}_{{\rm{s}}}\leqslant c$ and ${c}_{{\rm{s}}}\leqslant c/\sqrt{3}$ predicted by causality and asymptotically free theories like QCD, respectively. Our results are presented in Figure 1(b), where the megenta region shows the central density of PSR J0030+0451, and the dashed purple line represents the asymptotic limit of ${c}_{{\rm{s}}}=c/\sqrt{3}$. It is unclear whether such a limit has been reached by PSR J0030+0451 because of the relatively large uncertainties of both ρ_c and c_s (see also Greif et al. 2019; Raaijmakers et al. 2019). Future observation of NSs by NICER and Advanced LIGO/Virgo are necessary to robustly solve this issue.

The M − R and M−z_g relations are also constructed (Figures 1(c) and (d)), and the credible regions are consistent with that of Raaijmakers et al. (2019). Interestingly, the measured gravitational redshift of the isolated NS RX J0720.4-3125 (${z}_{g}={0.205}_{-0.003}^{+0.006};$ Hambaryan et al. 2017) is nearly the same as that of PSR J0030+0451 (${z}_{{\rm{g}}}\simeq {0.206}_{-0.017}^{+0.014};$ Riley et al. 2019a). If the NSs share the same EoS and the spin effect is neglected, then the masses of these two isolated objects are expected to be nearly the same, because the radii of NSs are almost unchanged in a relative narrow mass range. For z_g ≃ 0.206, using results shown in Figure 1(c) we find that the mass of PSR J0030+0451 is consistent with that of binary neutron star (BNS) systems (Kiziltan et al. 2013), suggesting no evidence for experiencing significant accretion of these isolated objects (see also Tang et al. 2020).

For a given EoS parameter ${{\boldsymbol{\theta }}}_{\mathrm{eos}}$, the canonical properties of the NSs (or the bulk properties at M = 1.4M_⊙) can be obtained through optimizing the central pressure p_c to reach a gravitational mass M = 1.4M_⊙. With a group of posterior samples of ${{\boldsymbol{\theta }}}_{\mathrm{eos}}$, we can deduce the probability distributions of bulk properties, such as Λ_1.4, R_1.4, and ${z}_{{\rm{g}},1.4}$. As shown in Figure 2, large ${M}_{\mathrm{TOV}}$ will boost the ${{\rm{\Lambda }}}_{1.4}$ and R_1.4 to higher values in both parameterization methods. This is understandable, because larger M_TOV can lead to stiffening of EoS. We also notice that the spectral method is more easily affected by the M_TOV condition than the piecewise method, because smooth parameterization models (e.g., spectral method) usually couple the high-density EoS, which sets M_TOV, from the low-density EoS, that determines R_1.4 (Capano et al. 2019; Tews et al. 2019). The canonical gravitational redshift ${z}_{{\rm{g}},1.4}$, which is connected to the R_1.4 by one-to-one mapping, mainly reflects the variation contrary to the R_1.4. While the combined dimensionless tidal deformability $\tilde{{\rm{\Lambda }}}$ of GW170817 simply keeps the same trend as Λ_1.4. This is mainly because both of the masses of NSs in GW170817 are near the region of 1.4M_⊙ (Abbott et al. 2019). Besides, in comparison to using the data of Riley et al. (2019b), the incorporation of the data of Miller et al. (2019b) yields slightly harder EoS in the low M_TOV case, while in high M_TOV case the results are remarkably consistent.

With the posterior samples of three pairs of M − Λ, namely, $\{{M}_{1},{{\rm{\Lambda }}}_{1},{M}_{2},{{\rm{\Lambda }}}_{2},{M}_{3},{{\rm{\Lambda }}}_{3}\}$, it is possible to deduce some bulk properties of NSs, as done in Landry & Kumar (2018) and Abbott et al. (2018). Here we adopt three universal relations to transform our posterior samples into constraints of tidal deformability, moment of inertia, and binding energy of NSs in the mass range (1.1–1.7)M_⊙. These relations (see Yagi & Yunes 2017, and references therein) read

$$\begin{eqnarray}\begin{array}{rcl}\lambda (M) & \simeq & {\lambda }_{\mathrm{ref}}+{\lambda }^{1}(M-{M}_{\mathrm{ref}})/{M}_{\odot },\\ {\mathrm{log}}_{10}\bar{I} & = & \sum _{n=0}^{4}{a}_{{\rm{n}}}{\left({\mathrm{log}}_{10}{\rm{\Lambda }}\right)}^{n},\\ \mathrm{BE}/M & = & \sum _{n=0}^{4}{b}_{{\rm{n}}}{\bar{I}}^{-n},\end{array}\end{eqnarray} \tag{ 4 }$$

where $\lambda {(M)\equiv {\rm{\Lambda }}(M)({GM}/{c}^{2})}^{5}$ is the tidal deformability (its dimensionless form is Λ), $\bar{I}\equiv {c}^{4}I/{G}^{2}{M}^{3}$ is the dimensionless moment of inertia (its dimensional form is I), BE is the binding energy, and G is the Newton's gravitational constant. We take the coefficients a_n and b_n from Landry & Kumar (2018) and Steiner et al. (2016), respectively. For each reference mass M_ref, we can obtain a best fit of λ_ref and λ¹ with a single posterior sample $\{{M}_{1},{{\rm{\Lambda }}}_{1},{M}_{2},{{\rm{\Lambda }}}_{2},{M}_{3},{{\rm{\Lambda }}}_{3}\}$. Repeating this process for a population of posterior samples, we get a sample of λ_ref and λ¹. The λ_ref is inserted into the $\bar{I}-{\rm{\Lambda }}$ relation to get $\bar{I}$. Then, the resulting samples of $\bar{I}$ are inserted into the ${BE}-\bar{I}$ relation to obtain samples of BE. We finally get constraints on Λ, I, and BE for each reference mass in the range of (1.1–1.7) M_⊙. We find that all four tests give results that are consistent with each other and shrink the prior of EoS to a softer region, with high M_TOV cases showing slightly larger Λ, larger I, and lower BE due to the stiffening of EoS in these cases (as shown in Figure 1(a) and Table 2). For PSR J0737-3039A, which has an accurately measured mass of 1.338M_⊙ and is expected to have a precise determination of moment of inertia in the next few years via radio observations, our Test (i) predicts its moment of inertia $I={1.35}_{-0.14}^{+0.26}\times {10}^{38}\mathrm{kg}\cdot {{\rm{m}}}^{2}$. This value is higher than that of Landry & Kumar (2018), which is mainly caused by the fact that they adopted a group of Λ_1.4 smaller than ours.