Decoding rotating neutron stars: Role of the symmetry energy slope

For non-rotating neutron stars, the equations of the hydrostatic equilibrium are the so-called Tolman-Oppenheimer-Volkof (TOV) equations (Oppenheimer & Volkoff, 1939), which reads:

where $\epsilon$ is the energy density and $p$ is the pressure, which must be obtained from an equation of state (EOS), $M(r)$ is the mass enclosed in radius $r$. Moreover, here I use G = c = 1.

In Fig. 1, we display the mass-radius relation coming from the TOV solution for different values of the symmetry energy slope. For all parameterizations, we use the BPS EOS (Baym et al., 1971b) for the neutron star’s outer crust and the BBP EOS (Baym et al., 1971a) for the inner crust. We use the BPS+BBP EoS up to the density of 0.0089 fm^-3 for all values of $L$, and from this point on, we use the QHD EoS, as suggested in ref. Glendenning (2000). Imposing $p_{crust}=p_{core}$ at the crust-core transition, I found that the core EoS begins around 0.03 fm^-3 for $L=44$ MeV, up to around 0.05 fm^-3 for $L=116$ MeV. In ref. Fortin et al. (2016), the authors compare the BPS+BBP crust EoS with a unified EoS. They show that for the canonical star, there is a variation in the radius of $60\mathrm{m}<~{}R_{1.4}~{}<150\mathrm{m}$. For a radius of 13 km, this implies an uncertainty around 1$\%$. It is also worth pointing out that unified EOS predict that the crust-core phase transition happens at higher densities, around 0.08 fm^-3 (see ref. Carreau et al. (2019); Grams et al. (2022) to additional discussion). The neutron stars’ core EOS derived from the QHD are chemically stable and electrically neutral. We also discuss a couple of observational constraints coming from the NICER X-ray telescope. The first, and maybe the most important one, is the PSR J0740+6620, whose mass and radius lie in the range of $M=2.08\pm 0.07~{}M_{\odot}$, and 11.41 km $<~{}R~{}<$ 13.70 km, respectively (Riley et al., 2021). The other constraint is related to the radius of the canonical $M=1.4~{}M_{\odot}$ star. Two NICER results constrain the radius of the canonical star between 11.52 km $<~{}R_{1.4}~{}<$ 13.85 km (Riley et al., 2019) and between 11.96 and 14.26 km (Miller et al., 2019). Here we use the union set of both constraints as a more conservative approach. Explicitly, we use 11.52 km $<~{}R_{1.4}~{}<$ 14.26 km as a constraint. The correlation between $L$ and the radii of neutron stars (i.e., that models with larger $L$ have larger radii) has been well-established in the literature (Cavagnoli et al., 2011; Dexheimer et al., 2019; Lopes et al., 2023). Moreover, as can be seen, the increase of the slope causes a small increase in the maximum mass.

In Newtonian formalism, the mass shedding limit, which is the maximum angular velocity a gravitational bounded object can rotate is given by (Hartle, 1967; Hartle & Thorne, 1968):

As the gravity in general relativity is stronger than his Newtonian counterpart, one can expect that a relativistic star can rotate even faster than the Newtonian limit. However, it was shown in ref. Weber & Glendenning (1992), that it is the opposite. The Kepler frequency is actually lower due to the drag of the inertial frame. A more precise expression for the Kepler frequency is given by the empirical formulae:

As pointed out in ref. Glendenning (2000); Weber & Glendenning (1992), the great advantage of eq. 5 is the fact that the relativistic Kepler frequency of a neutron star rotating at its mass shedding limit can be estimated from the mass and radius of the corresponding non-rotating limiting-mass star. In Fig 2, we display the Kepler limit in function of the neutron stars’ mass for different values of $L$. Moreover, in Tab 3, we display the inferior mass limit for four different values of $\Omega$: 1500 $s^{-}1$, 3000 $s^{-}1$, 4500 $s^{-}1$, and 5000 $s^{-}1$. We can notice that as we increase the angular velocity, we also increase the inferior limit of the neutron star mass ($M_{lim}$), where, in the Hartle-Throne approximation, only the results for neutron stars with $M~{}>M_{limt}$ are valid. Moreover, larger values of $L$ present higher values of $M_{lim}$ for a fixed $\Omega$. This was expected, due to the relation between the radii and the slope. For $L$ = 116 MeV, we see that the inferior limit is close to the canonical mass for $\Omega=5000~{}s^{-1}$. An interesting feature is the fact that for $\Omega=1500~{}s^{-1}$, the lower values of $L$ present the higher value of $M_{lim}$. The reason is for such low angular velocity, we have a very low neutron star mass with a low central density. And, as shown in Fig. 2 in ref. Lopes et al. (2023), for low densities, the EOS is stiffer for low values of $L$. Then it is reversed for high densities.

The detailed formalism about slowly rotating neutron stars can be found in the Hartle and Thorne original papers (Hartle, 1967; Hartle & Thorne, 1968; Hartle, 1973). Additional discussion can also be found in ref. Glendenning (2000); Weber & Glendenning (1992); Chubarian et al. (2000); Pattersons & Sulaksono (2021). Here we present only a small discussion and the necessary equations to build a numerical code. This is especially useful to help readers who are not familiar with the Hartle-Thorne (HT) formalism. The metric of a non-rotating, spherical symmetric star has the Schwarzschild form:

where:

Therefore, we can write the metric for a slowly rotating star as a perturbed Schwarzschield metric:

where the perturbative terms can be expanded up to the second-order:

where $P_{2}(\cos\theta)$ is the Legendre polynomial of order 2. Due to the symmetry $k_{0}(r)$ = 0. Moreover, as will become clear at the end of the calculations, the terms $h_{0}(r)$, and $m_{2}(r)$ do not contribute to the mass correction nor to the star deformation in our approximation, and will not be taken into account any further (they expression are nevertheless presented in ref. Hartle & Thorne (1968)). We also define $\nu_{2}~{}\equiv~{}h_{2}+k_{2}$ as suggested in ref. Hartle (1967). The parameter $\omega=\omega(r)=d\phi/dt$ is the angular velocity acquired by a free-falling observer due to the drag of the inertial frame. The angular velocity of the star relative to the local inertial frame, $\bar{\omega}=\Omega-\omega$ is obtained by solving:

with

The angular momentum ($J$) of a spherical symmetrical star with radius $R$ is given by:

Now, if $p$, $\epsilon$, and $n$ are the pressure, the energy density, and the number density at some given ($r,\theta)$ for a static star, the correspondent values of the pressure, energy density, and number density at the same given ($r,\theta)$ that is momentarily moving with the fluid in a rotating star are $p+\Delta p$, $\epsilon+\Delta\epsilon$, and $n+\Delta n$. The perturbation terms can also be expanded up to second-order:

where, $p_{0}^{*}$ and $p_{2}^{*}$ are dimensionless functions of $r$. Now, the perturbative terms, $m_{0}$, $p_{0}^{*}$, $\nu_{2}$, $h_{2}$, and $p_{2}^{*}$ are obtained by applying Einstein’s field equations. We obtain:

with the boundary conditions: $m_{0}(0)=p_{0}^{*}(0)=h_{2}(0)=\nu_{2}(0)=0$. Finally, we have:

which also implies $p_{2}^{*}(0)$ = 0. The mass correction is given by:

and the deformation of the star is represented by:

where

Therefore, if a non-rotating spherical symmetric neutron star with a central density $n_{c}$ presents a mass $M$ and a radius $R$, for the same central density, a slowly rotating neutron star will be a spheroid with mass $M+\delta M$ and a radius $R(\theta)=R+\delta R(R,\theta)$. Nevertheless, they are not the same star, as they present different baryon numbers, as pointed out in ref. Glendenning (2000). Finally, we can express the polar ($R_{p}$) and equator $(R_{e})$ radii as (Pattersons & Sulaksono, 2021):

this allows us to define the eccentricity:

As the measurement of the angular velocity is an easy task, once it is related to the period: $\Omega=2\pi/T$, where $T$ is the period, and the eccentricity does not depend on the absolute values of the radii, this quantity can play a very special role to fix the symmetry energy slope.

Once I established the basic equations of the Hartle-Thorne formalism, I now solve these coupled equations for four different values of $L$: 41, 76, 100, and 116 MeV. And for each value of $L$ we study four different values of angular velocities: $\Omega$: 1500 $s^{-}1$, 3000 $s^{-}1$, 4500 $s^{-}1$, and 5000 $s^{-}1$. The focus here is to investigate how different values of the slope change the main properties of neutron stars, I especially focus on the maxim mass, the correspondent central density, and the variation of the radius of some fixed mass. The first one is the canonical 1.4$M_{\odot}$ and the second is the 2.01 $M_{\odot}$ star, which is not only the most probable mass value of PSR J0348+0432 (Antoniadis et al., 2013) but also the lower limit of PSR J0740+6620, whose gravitational mass is 2.08 $\pm$ 0.07 $M_{\odot}$ (Miller et al., 2021; Riley et al., 2021).

In Fig. 3 I display the equator radius in function of the stellar mass for rotating neutron stars, for different values of $L$ and $\Omega$, while in fig 4 I show the eccentricity of the star in function of the mass for the same $L$ and $\Omega$. From a qualitative point of view, the effect of rotation is the same for all values of $L$. As we increase the angular velocity, we also increase the equator radius for all values of a fixed $M$. We also increase the maximum mass. The same is true of the eccentricity.

From the quantitative point of view, we see that the equator radius of the canonical star increases from 12.58 km for a non-rotating neutron star up to 13.72 km for a rotating neutron star with $\Omega=5000~{}(s^{-1})$ in the case of $L=44$ MeV, which represents an increase of 1.14 km or $9.1\%$. For the case of $L=116$ MeV, we see that the equator radius of the canonical star increases from 14.30 km up to 16.51 km. An increase of 2.21 km or 15$\%$. Therefore, not only the absolute values of $R_{eq}$ are dependent on the slope but also the relative ones. This means that for the same angular velocity, the star will appear flatter if the slope is higher. This can easily be seem from the eccentricity. While for $L$ = 44 MeV we have $e=0.51$ for $\Omega=5000~{}(s^{-1})$, for $L$ = 116 MeV we have $e=0.61$ for the same angular velocity. The same is true for $M=2.01~{}M_{\odot}$ but with more modest results. The equator radius for a 2.01 solar mass star increases from 12.40 km up to 13.07 km for $L=44$ MeV and from 13.82 km up to 15.03 km for $L$ = 166 MeV. Both at $\Omega=5000~{}s^{-1}$. This implies an increase of 0.67 km or $5.4\%$ for $L$ = 44 MeV and 1.48 km or $11\%$ for $L$ = 116 MeV. The eccentricities are 0.39 and 0.46 for $L$ = 44 and 116 MeV respectively.

I now analyze the variation of the maximum mass and the correspondent central density for $\Omega=5000~{}s^{-1}$. For $L$ = 44 the maximum mass increase from 2.31 $M_{\odot}$ up to 2.37 $M_{\odot}$. A change of 0.06 $M_{\odot}$ or $2.6\%$. The central density drops from 0.939 fm^-3 to 0.899 fm^-3, or a decrease of $4.3\%$. On the other hand, the maximum mass increases from 2.39 $M_{\odot}$ up to 2.47 $M_{\odot}$, an increment of $3.3\%$, while the central density decreases from 0.872 fm^-3 to 0.823 fm^-3, which represent a decrease of $5.6\%$. All relevant quantities are summarized in Tab. 4. The percentage is expressed with two significant digits. It is also important to emphasize that these results will be slightly different within a unified EOS.

For a fixed mass value and $\Omega$, the higher the symmetry energy, the higher will be the eccentricity. This fact can be very important in a future analysis of the 33 recently discovered millisecond pulsars (Smith et al., 2023). As the measurement of the angular velocity is easy, and the eccentricity does not depend on the absolute value of the neutron stars’ radii, its measure is a powerful tool, which combined with other techniques can improve the constraints of the symmetry energy slope.