Exploring universal characteristics of neutron star matter with relativistic ab initio equations of state

Neutron stars are one of the most compact objects in the universe: their central densities can reach as high as 5 to 10 times the saturation density of nuclear matter, $\rho_{0}\approx 0.16\ \text{fm}^{-3}$ Lattimer and Prakash (2004), which is far beyond what can be achieved in terrestrial laboratories. Therefore, neutron stars are ideal laboratories for studying ultradense matter, and have established close connections among nuclear physics, particle physics, and astrophysics.

The astrophysical observations of the global properties of neutron stars provide important constraints for the equation of state (EOS) of dense matter Lattimer and Prakash (2007); Özel and Freire (2016); Lattimer and Prakash (2016); Burgio et al. (2021), which is the only ingredient needed to unveil the structure of neutron stars theoretically. The high-precision mass measurements of massive neutron stars constitute nowadays one of the most stringent astrophysical constraints on the nuclear EOS; such measurements include PSR J1614-2230 ($1.928\pm 0.017M_{\odot}$) Demorest et al. (2010); Fonseca et al. (2016), PSR J0348+0432 ($2.01\pm 0.04M_{\odot}$) Antoniadis et al. (2013), and PSR J0740+6620 ($2.08\pm 0.07M_{\odot}$) Cromartie et al. (2020); Fonseca et al. (2021). Recently, the Neutron star Interior Composition Explorer (NICER) mission reported two independent Bayesian parameter estimations of the mass and equatorial radius of the millisecond pulsar PSR J0030+0451: $1.34_{-0.16}^{+0.15}\ M_{\odot}$ and $12.71_{-1.19}^{+1.14}\ \text{km}$ Riley et al. (2019) as well as $1.44_{-0.14}^{+0.15}\ M_{\odot}$ and $13.02_{-1.06}^{+1.24}\ \text{km}$ Miller et al. (2019). In combination with constraints from radio timing, gravitational wave (GW) observations, and nuclear physics experiments, these posterior distributions have been used to infer the properties of the dense matter EOS (see Ref. Riley et al. (2021) and references therein). Moreover, two independent Bayesian estimations of the radius for the massive millisecond pulsar PSR J0740+6620 have also been reported Riley et al. (2021); Miller et al. (2021).

Another unique probe for studying the properties of dense matter was extracted from the recent observation of GW signals emitted from a binary neutron star merger, i.e., GW170817 Abbott and et al. (2017). The tidal deformability, which denotes the mass quadrupole moment response of a neutron star to the strong external gravitational field induced by its companion Damour et al. (1992); Hinderer (2008); Flanagan and Hinderer (2008); Damour and Nagar (2009); Postnikov et al. (2010), can be inferred from the GW signals. The limits on the tidal deformability have been widely used to constrain the neutron star radius Fattoyev et al. (2018); Annala et al. (2018); Most et al. (2018); Tews et al. (2018), the asymmetric nuclear matter EOS Malik et al. (2018); Zhang et al. (2018); Tong et al. (2020), and hence the neutron skin thickness of ²⁰⁸Pb Fattoyev et al. (2018).

Besides, as rotating objects, the internal structures of neutron stars are strongly constrained by the moment of inertia, which can be determined from the measurements of spin-orbit coupling in double pulsar systems Lyne et al. (2004). Such a measurement of the moment of inertia for neutron star would have crucial implications for delimiting the EOS significantly Lattimer and Schutz (2005) and can be used to distinguish neutron stars from quark stars Yagi and Yunes (2013a). Special attention has been attracted by the system PSR J0737-3039 Burgay et al. (2003); Lyne et al. (2004); Kramer and Wex (2009), which is the only currently known double pulsar system. It is hoped that the moment of inertia of the 1.338$M_{\odot}$ primary component in this system, i.e., PSR J0737-3039A, will be measured eventually to within 10% Kramer et al. (2021), and could be used to impose new constraints on the EOS Lattimer (2021).

The global properties of neutron stars, like masses, radii, tidal deformabilities, and moments of inertia are highly sensitive to the EOS for neutron star matter Lim et al. (2019); Li and Sedrakian (2019); Tong et al. (2022). Nevertheless, it has been shown Yagi and Yunes (2013a, b) that, for slowly rotating neutron stars, there exist universal relations between the moment of inertia $I$, the tidal deformability $\Lambda$ (or Love number), and the quadrupole moment $Q$ of neutron stars, i.e., the so-called $I$-Love-$Q$ relations, which are approximately independent of the internal composition and the EOS for neutron star matter. Wide attention has been attracted by these universal relations (see Ref. Yagi and Yunes (2017) for a review). Although so far the reasons for these universal behaviors are not well understood Yagi et al. (2014); Jiang and Yagi (2020), attempts have been made to combine the universal relations with GW detections to infer neutron star properties Landry and Kumar (2018); Kumar and Landry (2019).

The robustness of the $I$-Love-$Q$ relations has been extensively studied with EOSs constructed from a variety of nuclear models (see Refs. Yagi and Yunes (2017); Wei et al. (2019) and references therein), including the relativistic Brueckner-Hartree-Fock (RBHF) theory Müther et al. (1987); Brockmann and Machleidt (1990); Li et al. (1992); Engvik et al. (1996); Gross-Boelting et al. (1999); Katayama and Saito (2015). Since the 1980s the RBHF theory has played an important role in understanding the properties of dense nuclear matter from realistic nucleon-nucleon ($NN$) interactions Anastasio et al. (1980, 1981). In the RBHF theory, the single-particle motion of the nucleon in nuclear matter is described with the Dirac equation, where the medium effects are absorbed into the single-particle potentials. In principle, the scalar and the vector components of the single-particle potentials should be determined in the full Dirac space Nuppenau et al. (1989), i.e., by considering the positive-energy states (PESs) and negative-energy states (NESs) simultaneously. However, to avoid the difficulties induced by NESs, the RBHF calculations are primarily performed in the Dirac space without NESs Müther et al. (1987); Brockmann and Machleidt (1990); Li et al. (1992); Engvik et al. (1996); Horowitz and Serot (1984); Gross-Boelting et al. (1999).

Recently, a self-consistent RBHF calculation in the full Dirac space was achieved for symmetric nuclear matter (SNM) Wang et al. (2021, 2022a) and asymmetric nuclear matter (ANM) Wang et al. (2022b). By decomposing the matrix elements of single-particle potential operator in the full Dirac space, the momentum-dependent scalar and vector components of the single-particle potentials are determined uniquely Wang et al. (2021). The long-standing controversy about the isospin dependence of the effective Dirac mass in relativistic ab initio calculations of ANM is also clarified Wang et al. (2022b).

The RBHF theory in the full Dirac space has been applied to neutron stars Wang et al. (2022b); Tong et al. (2022), where the mass, radius, and tidal deformability are calculated with realistic $NN$ interactions Bonn A, B, and C Machleidt (1989). The maximum mass of a neutron star is found less than $2.4M_{\odot}$ and the neutron star radius for $1.4M_{\odot}$ is predicted about 12 km, which are consistent with the astrophysical observations of massive neutron stars and simultaneous mass-radius estimations by NICER Miller et al. (2019). The tidal deformabilities for a $1.4M_{\odot}$ neutron star are predicted as 376, 473, and 459 for the three parametrizations of $NN$ interactions respectively, and all lie in the region $\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ inferred from the revised analysis by LIGO and Virgo Collaborations Abbott and et al. (2018).

In this work, we employ the RBHF theory in the full Dirac space to study other global properties of neutron stars, including the gravitational redshift, moment of inertia, and quadrupole moment under the slow-rotation and small-tidal-deformation approximation. The main focus will be the relation between the moment of inertia and the compactness parameter, as well as the universal $I$-Love-$Q$ relations. This paper is organized as follows. In Sec. II, the theoretical framework of the RBHF theory and structure equations for neutron star properties are briefly described. The obtained results and discussions are presented in Sec. III. Finally, a summary is given in Sec. IV.

In the left panel of Fig. 1, the speed of sound squared for neutron star matter from the RBHF theory with potential Bonn A is depicted as a function of density. The results calculated with the projection method and momentum-independence approximation show slowly increasing tendencies with the increase of density. In the full Dirac space, the speed of sound squared increases quickly with the density, and reaches 0.77 at $0.57\ \text{fm}^{-3}$. Above this density, the amplitude of scalar potential $U_{S}$ exceeds the nucleon rest mass, and the RBHF iteration in the full Dirac space is very difficult to achieve. This fact might be related to the extension of the Bonn potential to the full Dirac space, which was determined with PESs only. When the central density of a neutron star is fixed at $\rho_{c}=0.57\ \text{fm}^{-3}$, our calculations could support a neutron star with mass equal to $1.9M_{\odot}$. To further explore the maximum mass of a neutron star from the RBHF theory in the full Dirac space, we continue with an EOS where the speed of sound is equal to the speed of light Rhoades and Ruffini (1974); Gandolfi et al. (2012); Wang et al. (2022b) (red dotted line). This would provide an upper bound on the maximum mass of a neutron star.

Based on the EOS from the RBHF theory, the mass-radius relations of a neutron star can be calculated from the TOV equations, which are shown in the right panel of Fig. 1. The 68% and 95% contours of the joint probability density distribution of the mass and radius of PSR J0030+0451 Miller et al. (2019) and PSR J0740+6620 Miller et al. (2021) from the NICER analysis are also shown. It can be found that the results obtained in the full Dirac space, with the projection method, and with the momentum-independence approximation are consistent with the recent constraints by NICER. The maximum masses $M_{\text{max}}$ predicted by the three methods are $2.43M_{\odot}$, $2.31M_{\odot}$, and $2.18M_{\odot}$ respectively, which are consistent with the available astrophysical constraints from massive neutron star observations, such as PSR J1614-2230 ($1.928\pm 0.017M_{\odot}$) Demorest et al. (2010); Fonseca et al. (2016), PSR J0348+0432 ($2.01\pm 0.04M_{\odot}$) Antoniadis et al. (2013), and PSR J0740+6620 ($2.08\pm 0.07M_{\odot}$) Cromartie et al. (2020); Fonseca et al. (2021). The radii $R_{1.4M_{\odot}}$ of a $1.4M_{\odot}$ neutron star from the three methods are $11.97\ \text{km}$, $12.38\ \text{km}$, and $12.35\ \text{km}$, respectively.

In Table 1 we summarize the maximum mass of a neutron star and the radius and the central density for a $1.4M_{\odot}$ neutron star obtained by the RBHF theory in the full Dirac space, together with those obtained with the projection method and momentum-independence approximation. The smallest value for $R_{1.4M_{\odot}}$ found in the full Dirac space corresponds to the softest EOS for neutron star matter below a density of about 0.4 fm^-3. This fact can be further related to nuclear matter properties, which are also summarized in Tab. 1, including the binding energy per nucleon $E/A$, the symmetry energy $E_{\text{sym}}$, and its slope $L$. It is found that the symmetry energy and its slope at the saturation density obtained in the full Dirac space are the smallest among the three methods. This explains the reason for smallest $R_{1.4M_{\odot}}$ and shows how important it is to take both the PESs and NESs into account.

Figure 2 shows the relations between the gravitational redshift $z$ of a neutron star with its mass $M$ and radius $R$ obtained by the RBHF theory in the full Dirac space with the potentials Bonn A, B, and C. It can be seen that the $z$-$M$ relation is not strongly affected by the $NN$ interactions. With the increase of mass, the gravitational redshift shows a monotonically increasing behavior. A backbending phenomenon is found for the gravitational redshift with the decrease of radius, similar to the case for $M$-$R$ relations as shown in Fig. 1.

The clear one-to-one correspondence relation for gravitational redshift and mass established in the left panel in Fig. 2 can be used to infer the mass of a isolated neutron star, when the observation of the gravitational redshift is provided. In Fig. 2, astrophysical observations of gravitational redshift for isolated neutron stars RX J0720.4-3125 Hambaryan, V. et al. (2017), RBS 1223 Hambaryan et al. (2014), and RX J1856.5-3754 Hambaryan et al. (2014) are shown as shadow bands. The predicted masses by combining these observations and the theoretical calculations from the RBHF theory in the full Dirac space with Bonn A potential are listed in the third column of Table 2. The uncertainties are from the gravitational redshift measurements. In Table 2 the corresponding radii are also shown in the last column.

In Fig. 3 we display the ratio of the moment of inertia $I$ to $MR^{2}$ as a function of the compactness parameter $M/R$ obtained by the RBHF theory in the full Dirac space with $NN$ interactions Bonn A, B, and C. Lattimer et al. Lattimer and Schutz (2005) showed that, in the absence of phase transition and other effects that strongly soften the EOS at supranuclear densities, there is a relatively unique relation between the quantity $I/MR^{2}$ and $M/R$:

This relation is shown as the purple band in Fig. 3. It is found that our results are consistent with the universal relations obtained in Ref. Lattimer and Schutz (2005) for the range where $M/R>0.08\ M_{\odot}/\text{km}$. The derivation for smaller compactness is unimportant, since the observational evidence for neutron star masses and radii lie in the ranges of $1.2\ M_{\odot}\leq M\leq 2.2\ M_{\odot}$ and $9\ \text{km}\leq R\leq 15\ \text{km}$ respectively, which leads to the range of compactness as $0.08\ M_{\odot}/\text{km}<M/R<0.24\ M_{\odot}/\text{km}$. Lim and collaborators Lim et al. (2019) have investigated neutron star moments of inertia from Bayesian posterior probability distributions of the nuclear EOSs that incorporate information from microscopic many-body theory and empirical data of finite nuclei. The probability distribution for $I/MR^{2}$ is shown as in Fig. 3. They found that over the entire range of neutron star compactness $0<C\leq 0.25\ M_{\odot}/\text{km}$, their results can be well fitted with the following formula:

This result is also shown as the orange dashed line in Fig. 3. It can be seen that our results are very close to that obtained by Lim et al., especially for the neutron stars with small compactness.

Although the ratio of the neutron star moment of inertia $I$ to $MR^{2}$ has a universal function of the compactness parameter $M/R$, the moment of inertia itself depends sensitively on the neutron star’s internal structure. It has been suggested Lattimer and Schutz (2005) that a measurement accuracy of 10% for $I$ is sufficient to place strong constraints on the EOS.

In Table 3, we show the momenta of inertia $I_{1.338M_{\odot}}$ and radius $R_{1.338M_{\odot}}$ for PSR J0737-3039A predicted by the RBHF theory in the full Dirac space with three parametrizations for $NN$ interactions. The results obtained by the projection method Gross-Boelting et al. (1999); van Dalen et al. (2004) and momentum-independence approximation Brockmann and Machleidt (1990) are also shown. The RBHF theory in the full Dirac space leads to minimum values compared to the approximations in the Dirac space without NESs. This is understandable since the RBHF theory in the full Dirac space gives the minimum radius of a neutron star for the fixed canonical mass, as shown in Table 1.

The moment of inertia for PSR J0737-3039A predicted by the RBHF theory in the full Dirac space with Bonn A is $1.356\times 10^{45}\ \mathrm{g~{}cm^{2}}$, which is very close to the most probable value $1.36\times 10^{45}\ \mathrm{g~{}cm^{2}}$ obtained from Bayesian analysis (95% credibility) Lim et al. (2019). The result from the nonrelativistic ab initio variational calculations Akmal et al. (1998) is also shown in Table 3, where the Argonne $v18$ interaction (AV18) Wiringa et al. (1995) is used, together with the relativistic boost corrections to the two-nucleon interaction as well as three-nucleon interactions modeled with the Urbana force Pudliner et al. (1995). The nonrelativistic ab initio calculation leads to a moment of inertia smaller that what we obtain, similarly to the case for radius. In Ref. Landry and Kumar (2018), by using well-known universal relations among neutron star observables, the reported 90% credible bound on the tidal deformability $\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ from GW170817 Abbott and et al. (2018) has been translated into a direct constraint on the moment of inertia of PSR J0737-3039A, giving $I_{1.338M_{\odot}}=1.15^{+0.38}_{-0.24}\times 10^{45}\ \mathrm{g~{}cm^{2}}$. It can be seen that the results with the three methods for RBHF theory are consistent with this constraint.

Let us now compare the EOSs obtained by the relativistic ab initio calculations, i.e., the RBHF theory in the full Dirac space, projection method, and momentum-independence approximation with Bonn potentials, to the universal $I$-Love-$Q$ relations. The $I$-Love as well as $Q$-Love relations and $I$-$Q$ relations are shown in the top panels of Figs. 4 and 5, respectively. The dimensionless moment of inertia $\bar{I}$ and dimensionless quadrupole moment $\bar{Q}$ are defined in Eq. (15). A single parameter along the curve is the mass or compactness, which increases to the left of the plots. Similarly to Ref. Yagi and Yunes (2017), we only show data with the mass of an isolated, nonrotating configuration in the range $1M_{\odot}<M<M_{\text{max}}$ with $M_{\text{max}}$ representing the maximum mass for such a configuration. One observes that the universal relations hold very well. Since the relations are insensitive to EOS, one can construct a single fit (black solid curves) given by Yagi and Yunes (2013b, 2017)

where coefficients are listed in Table 4. These coefficients are very close to that in Ref. Yagi and Yunes (2017), where a large number of EOSs are considered. The bottom panels of Figs. 4 and 5 show the absolute fractional difference between all the data and the fit, which is less than $1\%$ in the whole range.

The universal relations between $\bar{I}$ and $\Lambda$ allows one to extract the momenta of inertia of a neutron star with $1.4$ solar mass, $\bar{I}_{1.4M_{\odot}}$, from the tidal deformability $\Lambda_{1.4M_{\odot}}$ from GW170817. The revised analysis from LIGO and Virgo Collaborations, $\Lambda_{1.4M_{\odot}}=190^{+390}_{-120}$ Abbott and et al. (2018), corresponds to $\bar{I}_{1.4M_{\odot}}=10.30^{+3.39}_{-2.10}$ as shown in the left panel of Fig. 4. From $\bar{I}_{1.4M_{\odot}}$ and the relation $\bar{I}=I/M^{3}$ we obtain $I_{1.4M_{\odot}}=1.22^{+0.40}_{-0.25}\times 10^{45}\mathrm{g\ cm^{2}}$. These values are consistent with the results $\bar{I}_{1.4M_{\odot}}=11.10^{+3.64}_{-2.28}$ and $I_{1.4M_{\odot}}=1.15^{+0.38}_{-0.24}\times 10^{45}\mathrm{g\ cm^{2}}$ in Ref. Landry and Kumar (2018), where the $I$-Love relation is obtained by using a large set of candidate neutron star EOSs based on relativistic mean-field and Skyrme-Hartree-Fock theory.

In summary, the RBHF theory in the full Dirac space has been employed to study the gravitational redshift, moment of inertia, and quadrupole moment of neutron star under the slow-rotation and small-tidal-deformation approximation. The one-to-one correspondence relation for gravitational redshift and mass is established and used to infer the masses of isolated neutron stars by combining gravitational redshift measurements. The ratio of the moment of inertia $I$ to $MR^{2}$ as a function of the compactness $M/R$ is obtained, and is consistent with the universal relations shown by Lattimer et al. Lattimer and Schutz (2005) and that from Bayesian posterior probability distributions by Lim et al. Lim et al. (2019). Using $NN$ interactions Bonn A, B, and C, the moment of inertia for $1.338M_{\odot}$ pulsar PSR J0737-3039A is predicted to be $1.356\times 10^{45}$, $1.381\times 10^{45}$, and $1.407\times 10^{45}\ \mathrm{g~{}cm^{2}}$, which are consistent with the constraint translated from the tidal deformability deduced from GW170817 with universal relations among neutron star observables. The EOSs constructed by the RBHF theory in the full Dirac space, together with those by the projection method and momentum-independence approximation, are compared successfully to universal $I$-Love-$Q$ relations. By combing the tidal deformability $\Lambda_{1.4M_{\odot}}$ from GW170817 and the numerical fitting for these universal relations from relativistic ab initio EOSs, the moment of inertia of neutron star with 1.4 solar mass is deduced as $I_{1.4M_{\odot}}=1.22^{+0.40}_{-0.25}\times 10^{45}\mathrm{g\ cm^{2}}$.