Decoding Rotating Neutron Stars: Role of the Symmetry Energy Slope

In 2023 December, the Fermi Large Area Telescope Catalog announced the discovery of 33 new millisecond pulsars. Motivated by that, in this work, I study how different values of the symmetry energy slope (L) affect the properties of static and slowly rotating neutron stars. For fixed values of angular velocity, I study how the slope influences the increase of the maximum mass, the radii of the canonical 1.4 M ⊙, its eccentricity, as well the same quantities for the 2.01 M ⊙ stars. I show that different slope values cause different variations not only in the absolute quantities but also in relative ones. Indeed, different slope values predict different values for the eccentricity, which does not depend on the absolute value of the neutron stars’ radii. Therefore, this quantity can be a powerful tool to constrain the symmetry energy slope.


INTRODUCTION
Pulsars are among the fastest-spinning objects in nature, and their period can be as low as a few milliseconds.The first millisecond pulsar was discovered over forty years ago, the PSR B1937+21, with a frequency of 642 Hz (Ω = 4034 s −1 ), (Backer et al. 1982).Currently, PSR J1748-2446ad is the fastest-spinning pulsar known, at 716 Hz (Ω = 4499 s −1 ) (Hessels et al. 2006).Another remarkable example is the black widow pulsar PSR J0952-0607, which is not only the second fastest-spinning pulsar, rotating at 707 Hz (Ω = 4442 s −1 ), but also it is potentially the heaviest neutron star ever known, whose mass lies between 2.35 ± 0.17 M ⊙ (Romani et al. 2022).A better understanding of the effects of rotation in neutron stars become even more important in the last weeks.In December 2023, the third Fermi Large Area Telescope Catalog of gamma-ray pulsars announced the discovery of 33 new millisecond pulsars, among these, at least a dozen had an inferred angular velocity above 3500 s −1 (Smith et al. 2023).
The modern theory of pulsars considers them as rotating neutron stars.In the neutron star's crust, where the densities are usually below half of the saturation density, the matter is inhomogeneous and presents a high degree of complexity.In this regime clusters coexisting with neutron (super)fluid.Moreover, the competition between attractive nuclear force and repulsive Coulomb interaction can turn the nuclear matter into a frustrated system, i.e., the system presents more than one low-energy configuration, which can cause the onset of unusual nuclear shapes with different geometries.This is called nuclear pasta phase (Ravenhall et al. 1983;Lorenz et al. 1993).Although these properties are relevant for rotating neutron stars since they are believed to affect pulsar glitches, their studies are beyond the scope of the present work.I refer the interested reader to Chamel & Haensel (2008) and the references therein.On the other hand, in the neutron stars' core, where the densities can reach several times the nuclear saturation density, pressure is so strong that nucleons can no longer form clusterized matter, and neutrons, protons, and leptons form a homogeneous liquid in beta-equilibrium.The tool to describe a neutron star's interior is therefore many-body nuclear physics.A good model must be able to reproduce at least five well-constrained parameters at the saturation density: the nucleon effective mass (M * N /M N ), the binding energy per baryon (E/A), the compressibility (K), the symmetry energy (S 0 ), and the saturation density itself (n 0 ).
A sixth important parameter is the symmetry energy slope, or simply slope (L), However, unlike the other five parameters, the slope is not yet well constrained.In the earlier 2010s, most studies pointed to a relatively low value a llopes@cefetmg.comfor L. For instance, in refs.Paar et al. (2014); Lattimer & Steiner (2014); Lattimer & Lim (2013) upper limits of 54.6, 61.9, and 66 MeV respectively were suggested.However, in the last couple of years, the situation has changed and new experiments have pointed to a significantly higher upper limit.For instance, in a study about the spectra of pions in intermediate energy collisions, an upper limit of 117 MeV was obtained (Estee et al. 2020), while in one of the PREX II analyses (Reed et al. 2020) an upper limit of 143 MeV was suggested.All these conflicting results have been well summarized in a recent paper (Tagami et al. 2022): the CREX group points to a slope in the range 0 < L < 51 MeV, while PREX II results point to 76 MeV < L < 165 MeV.The CREX and PREX II results do not overlap.That is a huge problem that must be solved.
It is also well-known that the symmetry energy slope affects the neutron star properties (Lattimer & Steiner 2014;Lattimer & Lim 2013;Cavagnoli et al. 2011;Lopes & Menezes 2014).Therefore, neutron star observations can indeed be used to constrain the symmetry energy slope.However, most works on this subject use the static approximation.As neutron stars can spin very fast, the effects of their rotation can influence their macroscopic properties.Quantifying the effects of the slope in rotating neutron stars is therefore our main goal.
In this work, I study how different values of the symmetry energy slope affect neutron stars that rotate at different angular velocities.To accomplish this task, I use quantum hadrodynamics (QHD) with the traditional σ − ω − ρ mesons (Glendenning 2000;Serot 1992).Moreover, to keep the symmetry energy fixed while varying the slope, we use two extensions of QHD: to reduce the slope, we add the non-linear ω − ρ coupling as presented in the IUFSU model (Fattoyev et al. 2010;Cavagnoli et al. 2011;Dexheimer et al. 2019), while to increase the slope, we add the scalar-isovector δ meson (Kubis & Kutschera 1997;Liu et al. 2002).The rotation is introduced via Hartle and Thorne's slowly-rotating formalism (Hartle 1967;Hartle & Thorne 1968;Hartle 1973).In general, most of the studies related to rotating neutron stars compare the static results with those in which the star rotates at the mass shedding limit, also called Keplerian limit, Ω K (Weber & Glendenning 1992;Chubarian et al. 2000;Dhiman et al. 2007;Jha et al. 2008;Haensel et al. 2009;Most et al. 2020;Rather et al. 2021).Here, I try a more realistic approach, and study rotating neutron stars with angular velocity close to those observed in nature.As pointed out earlier, the fastest-spinning pulsar has an angular velocity of Ω = 4499 s −1 .As even fast-spinning neutron stars may exist, I use Ω = 5000 s −1 , as the maximum angular velocity, which corresponds to around 10% above the angular velocity of the PSR J1748-2446ad.This implies that most of the neutron stars studied here are significantly below the mass shedding limit, which guarantees Hartle and Thorne's slowly rotating approach as a good approximation.

NUCLEAR MODEL
To describe the nuclear interaction, I use here an extended version of the quantum hadrodynamics (QHD) (Glendenning 2000;Serot 1992), which includes the ω − ρ non-linear coupling (Fattoyev et al. 2010;Cavagnoli et al. 2011;Dexheimer et al. 2019), as well the scalar-isovector δ meson (Kubis & Kutschera 1997;Liu et al. 2002).Considering a pure nucleonic matter, the Lagrangian density reads: Here, ψ N is the Dirac field of the nucleons, and σ, ω µ , δ, and ρ µ are the mesonic fields.The g's are the Yukawa coupling constants that simulate the strong interaction, M N is the nucleon mass, and m s , m v , m δ , and m ρ are the masses of the σ, ω, δ, and ρ mesons, respectively.The self-interaction term U (σ) was introduced in ref.Boguta & Bodmer (1977) to fix the incompressibility, is given by: The Pauli matrices are denoted by τ , the antisymmetric mesonic field strength tensors are given by their usual expressions: . The γ µ are the Dirac matrices and Λ ωρ is a non-linear isoscalar-isovector mixing coupling that provides a simple and efficient method of softening the symmetry energy without compromising the success of the model in reproducing well-determined ground-state observables (Fattoyev et al. 2010).The detailed calculation of the EOS in the mean field approximation can be found in refs.Glendenning (2000); Serot (1992)  In this work, I follow ref.Lopes et al. (2023) and use a slightly modified version of the L3ωρ parametrization proposed in ref.Lopes (2022).The model parameters and their predictions for symmetric nuclear matter are presented in Table 1.The five well-known nuclear constraints at saturation density discussed above are taken from two extensive review articles: Dutra et al. (2014); Oertel et al. (2017), and are also included in Table 1.The parameters (g ρ /m ρ ) 2 , (g δ /m δ ) 2 , and Λ ωρ are chosen in order to fix the symmetry energy at the saturation point S 0 = 31.7 MeV while varying the slope L. The values of these parameters are reported in Table 2.As the parameterizations presented in this work are the same as presented in ref.Lopes et al. (2023), additional discussion about the symmetry energy and its slope, as well complementary results for non-rotating neutron stars can be found in ref.Lopes et al. (2023).Furthermore, in Table 2 we also calculate the second derivative of the symmetry energy (Zhang et al. 2018): (2) This quantity was not calculated in previous works for the present parameterizations.As can be seen, K sym is negative for lower values of L, and positive for higher ones.From the field theory point of view, K sym is negative for models that present non-linear ω − ρ coupling and positive for models that incorpore the scalar-isovector δ meson.

Static Limit
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 ǫ 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 60m < R 1.4 < 150m.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 ± 0.07 M ⊙ , 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 ⊙ 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 Ω: 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 Ω.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 Ω = 5000 s −1 .An interesting feature is the fact that for Ω = 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.

Rotating neutron stars
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 L. L. Lopes (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 θ) 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 ν 2 ≡ h 2 + k 2 as suggested in ref.Hartle (1967).The parameter ω = ω(r) = dφ/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, ω = Ω − ω is obtained by solving: The angular momentum (J) of a spherical symmetrical star with radius R is given by: Now, if p, ǫ, and n are the pressure, the energy density, and the number density at some given (r, θ) for a static star, the correspondent values of the pressure, energy density, and number density at the same given (r, θ) that is momentarily moving with the fluid in a rotating star are p + ∆p, ǫ + ∆ǫ, and n + ∆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 , ν 2 , h 2 , and p * 2 are obtained by applying Einstein's field equations.We obtain: with the boundary conditions: 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 +δM and a radius R(θ) = R+δR(R, θ).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: Ω = 2π/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: Ω: 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.4M ⊙ and the second is the 2.01 M ⊙ 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 ± 0.07 M ⊙ (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 Ω, while in fig 4 I show the eccentricity of the star in function of the mass for the same L and Ω.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 Ω = 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 Ω = 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 ⊙ 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 Ω = 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 Ω = 5000 s −1 .For L = 44 the maximum mass increase from 2.31 M ⊙ up to 2.37 M ⊙ .A change of 0.06 M ⊙ 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 ⊙ up to 2.47 M ⊙ , 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 Ω, 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.

CONCLUSION
In this work, I study the effects of the symmetry energy slope on rotating neutron stars.For four different values of L, I choose four values of Ω, and see how some macroscopic properties are affected: the maximum mass, the central number density of the maximum mass, the equator radius, and the eccentricity of the canonical 1.4M ⊙ star and the 2.01 M ⊙ star.We can see that the rotation affects all the neutron stars' properties.Moreover, the higher the slope, the higher the perturbation induced by the rotation, in both: absolute and relative terms.For a fixed value of Ω, we see that the effects are small in the maximum mass, are more significant at the 2.01 M ⊙ star, but strongly affect the canonical star.Furthermore, as the eccentricity does not depend on the absolute value of the star's radius, and the measurement of angular velocity is not a hard task, this quantity can potentially help to constraint the symmetry

Figure 1 .
Figure 1.Mass-radius relation for different values of L for non-rotating neutron stars.The hatched areas are the constraints discussed in the text.

Figure 3 .
Figure 3. (Color online) Equator radius as a function of the stellar mass for different values of L and different values of Ω.

Figure 4 .
Figure 4. Eccentricity of rotating neutron stars as a function of the stellar mass for different values of L and different values of Ω.

Table 1 .
Oertel et al. (2017)herein.Model parameters used in this study and their predictions for symmetric nuclear matter at saturation density.The parametrization was taken from ref.Lopes et al. (2023), and the phenomenological constraints were taken from refs.Dutra et al. (2014);Oertel et al. (2017).

Table 3 .
Kepler frequency in function of the neutron stars mass for different values of L. The vertical line corresponds to the canonical M = 1.4M⊙ star.Inferior limit of the neutron stars' masses for a fixed Ω.