© The Authors 2026

1. Introduction

A millisecond pulsar (MSP) is a neutron star (NS) with a rotational period less than about 10 milliseconds (ms), exhibiting significantly faster spin rates than ordinary pulsars. Millisecond pulsars are generally older than standard pulsars. They have formed through a prolonged period of accretion from a companion star, which spins them up to periods less than about 10 ms (Nagase 1989). This mass transfer process causes the NS to spin faster and leads to a reduction in its magnetic field. A significant proportion of MSPs are observed in binary star systems (Bhattacharyya & Roy 2022). Their masses typically range from 1 to 2 M_⊙ (M_⊙ ≃ 1.989 × 10³³ g is the mass of the Sun), while their radii generally fall within the range of 10 to 14 km. Recently, a large population of MSPs has been discovered in globular clusters (GCs) (Pan et al. 2021; Ridolfi et al. 2021; DeCesar et al. 2015; Wang et al. 2020), which have been shown to undergo a process of accretion acceleration. The first MSP, PSR B1937+21, was discovered in 1982 (Backer et al. 1982), with a period of 1.558 ms (642 Hz), once known as the second fastest-spinning pulsar. Currently, PSR J1748-2446ad is recognized as the fastest-spinning pulsar known, with a frequency of 716 Hz (1.397 ms) (Hessels et al. 2006). M14A is a black widow pulsar with a rotational period of approximately 1.98 ms (Pan et al. 2021). It is the second fastest-spinning black widow pulsar in GCs, after the one found in GC Terzan 5 (1.67 ms) (Ransom et al. 2005). Another remarkable MSP is the black widow pulsar PSR J0952-0607, which is not only the second fastest-spinning pulsar now, with a period of 1.41 ms (707 Hz), but it may also be the heaviest NS ever known, with a mass of 2.35 ± 0.17 M_⊙ (Bassa et al. 2017; Romani et al. 2022).

FAST (Five-hundred-meter Aperture Spherical radio Telescope) is the most sensitive single-antenna radio telescope in the world today (Li et al. 2018). It has identified a lot of new pulsars since its commissioning, including MSPs (Han 2025). Specific work includes the Galactic Plane Pulsar Snapshot (GPPS) survey, which discovered 201 pulsars, of which 40 are MSPs (Han et al. 2021). More than 30 GC pulsars were discovered between 2017 and 2021 (Pan et al. 2020; Wang et al. 2020). PSR J1717+4307A was the first pulsar discovered in the GC Messier 92 (M92), with a spin period of 3.16 ms (316.5 Hz) (Pan et al. 2020). PSR J1853-0842A, found through high-sensitivity searches using FAST, has a period of 2.15 ms and is part of a black widow system (Yan et al. 2021). Motivated by these findings, in this work we study the effects of nuclear symmetry energy and its slope on the bulk properties of MSPs.

The reasonable range of symmetry energy and its slope have been widely discussed (Tews et al. 2017; Holt & Kaiser 2017; Reed et al. 2021; Essick et al. 2021; Tagami et al. 2022), incorporating both nuclear experimental constraints and astrophysical observations. However, this remains an open problem, particularly since the PREX II results (Reed et al. 2021) do not overlap with those from other constraints. It is well known that the nuclear symmetry energy and its slope have significant effects on the global properties of NSs (Cavagnoli et al. 2011; Baldo & Burgio 2016; Margueron et al. 2018; Wu & Shen 2019; Hu et al. 2020; Wu et al. 2021; Wang et al. 2024), but most studies on this topic focus on static NSs. Cavagnoli et al. (2011) demonstrated that for NSs with masses exceeding 1 M_⊙, the radius exhibits a linear relationship with the symmetry energy slope. Utilizing a Bayesian statistical approach, Margueron et al. (2018) found that symmetry energy slope and curvature play an important role in determining the radius of NSs. Tsaloukidis et al. (2019) examined the effects of symmetry energy on the transition density and r-mode instability of a rotating NS. Li et al. (2020) investigated how the isoscalar skewness coefficient and the symmetry energy slope influence the masses of rapidly rotating NSs. Lopes (2024) studied how the symmetry energy slope of the core EOS affect the properties of static and slowly rotating NSs using non-unified EOSs with a fixed crust segment. In most calculations of NS properties, a non-unified EOS is employed, i.e., the core EOS is matched to a crust EOS obtained from different models (Read et al. 2009; Koliogiannis & Moustakidis 2021; Lopes 2024). However, it should be more advantageous to calculate the structures of static NSs and MSPs using a unified EOS, whereby both the core and crust EOS are based on the same theoretical framework. The interior of a NS has three qualitatively different domains: the outer crust, the inner crust, and the core, separated by phase transition boundaries, namely the neutron drip density (interface between the outer crust and the inner crust) and the crust-core interface. The outer crust is mainly determined by masses of neutron-rich nuclei (Baym et al. 1971; Haensel et al. 2007; Fortin et al. 2016), and its properties are described in this work by the EOS of Baym et al. 1971 (Baym et al. 1971), referred to as BPS. The inner crust has to be studied using phenomenological models due to the presence of inhomogeneous nuclear matter. The unified EOS is obtained from the same Lagrangian to describe the neutron drip density, the inner crust, and the core. The use of a unified EOS allows for a self-consistent treatment of the core-crust transition. It has been reported that the EOS of the crust significantly impacts NS properties (Fortin et al. 2016; Ji et al. 2019; Suleiman et al. 2021; Huang et al. 2024). In this work, we compare the results using unified EOSs with those using non-unified EOSs, and explore the effects of the symmetry energy and its slope using the parameter family TM1 introduced in Bao et al. (2014), in which the self-consistent Thomas-Fermi approximation has been employed to calculate the inner crust structure. A similar calculation using the family IUFSU (Bao et al. 2014) was also performed to check the symmetry energy dependence of our results. The original TM1 and IUFSU parameter sets describe the properties of finite nuclei and the saturation properties of nuclear matter very well (Sugahara & Toki 1994; Fattoyev et al. 2010). The TM1 and IUFSU families predict similar ground-state properties of finite nuclei, whereas they affect the interface densities between the crust and core and the crust composition due to their different symmetry energy slopes. The unified EOS utilizes a consistent model for both the crust and core, while for the non-unified EOS, we adopt the same TM1e or original IUFSU EOS for the core (Shen et al. 2020; Fattoyev et al. 2010), but the crust features different symmetry energy slopes within the same parameter family.

The NS/MSP observations are crucial for understanding symmetry energy and its slope, especially the NS radius. The mass measurement of MSP J1614+2230 (P = 3.15 ms) indicates a mass of M = 1.97 ± 0.04 M_⊙ (Demorest et al. 2010), with an updated result of M = 1.908 ± 0.016 M_⊙ (Arzoumanian et al. 2018). This finding raised the observational limit of NS mass to approximately 2 M_⊙. In 2013, PSR J0348+0432 (P = 39.12 ms) was reported to have a mass of M = 2.01 ± 0.04 M_⊙ (Antoniadis et al. 2013). Later, in 2019, PSR J0740+6620 with a spin period of about 2.88 ms was found to have a mass of M = 2.08 ± 0.07 M_⊙ (Cromartie et al. 2020; Fonseca et al. 2021), which is recognized as the most massive MSP measured precisely through pulse timing observations and detailed assessments of general relativistic effects. The Neutron Star Interior Composition Explorer (NICER) has enabled simultaneous measurements of mass and radius. For PSR J0740+6620, the measurements yield a mass of M = 2.08 ± 0.07 M_⊙, and a radius of $R=13.7_{-1.5}^{+2.6}\mathrm{km}$ (Miller et al. 2021). An alternatively result is a mass of $2.{072}_{-0.066}^{+0.067}{M}_{\odot }$ with a radius of $12.{39}_{-0.98}^{+1.30}$ km (Riley et al. 2021). For PSR J0030+0451, its spin period is known to be 4.87 ms, and its mass was determined to be $M=1.{44}_{-0.14}^{+0.15}{M}_{\odot }$ with a radius of $R=13.{02}_{-1.06}^{+1.24}\mathrm{km}$ (Miller et al. 2019), or alternatively, $1.{34}_{-0.16}^{+0.15}{M}_{\odot }$ with a radius of $12.{71}_{-1.06}^{+1.24}$ km (Miller et al. 2019; Riley et al. 2019).

To infer the radii of NSs, NICER observes the X-ray pulse profiles emitted by hot spots on the surfaces of rotating NSs (Raaijmakers et al. 2019). A critical ingredient in the pulse profile model is the oblate shape of the star. The shape eccentricity introduced by rotation could impact the X-ray light curve from rotating pulsars (Morsink et al. 2007). To reduce the number of parameters used to describe the shape and multipoles of rotating NSs, Bauböck et al. (2013) established a universal relation between the eccentricity and the compactness. Gao et al. (2023) proposed a universal relation that connects the eccentricity to the radius and moment of inertia of rotating NSs. In this work, we also discuss the dependence of the eccentricity and moment of inertia of MSPs on the symmetry energy slope, L, and the spin period, P.

This paper is organized as follows. Section 2 provides a brief description of the relativistic mean-field (RMF) model and the unified EOSs, while the non-unified connection of the crust and core EOSs is also introduced. Additionally, the Hartle-Thorne slow-rotation approximation is reviewed briefly. In Section 3, we present our main results regarding the effects of symmetry energy on the properties of rotating MSPs. Finally, Section 4 is devoted to the conclusions.

2. Models

We employed the RMF model to construct the EOSs. In this framework, nucleons interact through the exchange of mesons, including scalar σ, vector ω, and isovector ρ mesons. The model also includes nonlinear terms for the σ and ω mesons, as well as an additional ω − ρ coupling term (Bao et al. 2014; Shen et al. 2020). The Lagrangian density of NS matter is written as

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{RMF}}& =\sum _{i=p,n}{\overline{\psi }}_i\{i{\gamma }_{\mu }{\partial }^{\mu }-(M_N+g_{\sigma }\sigma )\hfill \\ \hfill & -{\gamma }_{\mu }[g_{\omega }{\omega }^{\mu }+\frac{g_{\rho }}{2}{\tau }_a{\rho }^{a\mu }]\}{\psi }_i\hfill \\ \hfill & +\frac{1}{2}{\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -\frac{1}{2}{m}_{\sigma }^2{\sigma }^2-\frac{1}{3}{g}_2{\sigma }^3-\frac{1}{4}{g}_3{\sigma }^4\hfill \\ \hfill & -\frac{1}{4}{W}_{\mu \nu }{W}^{\mu \nu }+\frac{1}{2}{m}_{\omega }^2{\omega }_{\mu }{\omega }^{\mu }+\frac{1}{4}{c}_3{\left({\omega }_{\mu }{\omega }^{\mu }\right)}^2\hfill \\ \hfill & -\frac{1}{4}{R}_{\mu \nu }^a{R}^{a\mu \nu }+\frac{1}{2}{m}_{\rho }^2{\rho }_{\mu }^a{\rho }^{a\mu }\hfill \\ \hfill & +{\mathrm{\Lambda }}_{\mathrm{v}}\left(g_{\omega }^2{\omega }_{\mu }{\omega }^{\mu }\right)\left(g_{\rho }^2{\rho }_{\mu }^a{\rho }^{a\mu }\right)\hfill \\ \hfill & +\sum _{l=e,\mu }{\overline{\psi }}_l(i{\gamma }_{\mu }{\partial }^{\mu }-m_l){\psi }_l,\hfill \end{array}$$(1)

where W^μν and R^αμν are the antisymmetric field tensors corresponding to ω^μ and ρ^aμ, respectively.

In order to explore the influence of symmetry energy on rotating MSPs, we provide a brief description of the symmetry energy, E_sym, and its slope, L. The symmetry energy of nuclear matter is defined as

$$\begin{array}{c}\hfill E_{\mathrm{sym}}({\rho }_B)={\frac{1}{2}\frac{{\partial }^2\epsilon ({\rho }_B,\delta )/{\rho }_B}{\partial {\delta }^2}|}_{\delta =0},\end{array}$$(2)

where ρ_B is the baryon number density, and δ is the isospin asymmetry defined as δ = (ρ_n − ρ_p)/(ρ_n + ρ_p), with ρ_n and ρ_p being the neutron and proton number densities, respectively. ε(ρ_B, δ)/ρ_B is the energy per baryon. Under the RMF framework, the symmetry energy expression can be derived from Bao et al. (2014), Hu et al. (2020)

$$\begin{array}{c}\hfill E_{\mathrm{sym}}({\rho }_B)=\frac{k_F^2}{6\sqrt{M_N^{\ast 2}+k_F^2}}+\frac{g_{\rho }^2{\rho }_B}{8(m_{\rho }^2+2{\mathrm{\Lambda }}_V{g}_{\omega }^2{g}_{\rho }^2{\omega }^2)}.\end{array}$$(3)

k_F is the Fermi momentum of symmetric nuclear matter and M_N^* = M_N + g_σσ is the effective nucleon mass, with M_N being the free nucleon mass. The symmetry energy slope is given as

$$\begin{array}{c}\hfill L=3{\rho }_0\left(\frac{\partial E_{\mathrm{sym}}}{\partial {\rho }_B}\right),\end{array}$$(4)

where ρ₀ is the saturation density. We used a series of EOSs based on the TM1/IUFSU parameter set, which exhibit different E_sym and L at the saturation density, while keeping the same values for the other saturation properties. These EOSs feature a uniform symmetry energy at a subsaturation density of 0.11 fm⁻³ (E_sym(0.11 fm⁻³) = 28.05 MeV (TM1), E_sym(0.11 fm⁻³) = 26.78 MeV (IUFSU)) (Bao et al. 2014). Therefore, the symmetry energy coefficient at saturation density ρ₀, E_sym(ρ₀), is also affected by the change in the symmetry energy slope. Through this approach, we can focus on the effect of a single variable among the saturation parameters, the symmetry energy slope, L. In Table 1, we present the mass and coupling constants of the nucleon and mesons in the original TM1 and IUFSU model. Table 2 lists the saturation properties, while Tables 3 and 4 list the corresponding coupling constants for generated models. It is shown that the parameters in this TM1 series yield different symmetry energies at ρ₀, while a smaller L corresponds to a smaller E_sym(ρ₀). The configuration with a symmetry energy slope of L = 40 MeV is commonly referred to as TM1e, which is more compatible with the GW170817 constraint compared to other values of L (Ji et al. 2019). The detailed EOS table obtained with TM1e can be found in Shen et al. (2020). In this work, we utilize unified EOSs in which both the core and crust EOSs share the same value of symmetry energy slope, L. To investigate the effects of unified EOSs, we also compare with the results using non-unified EOSs that match the crust with different L values and maintain the same core EOS by using TM1e.

For nonrotating NSs, the equations of the hydrostatic equilibrium are the Tolman-Oppenheimer-Volkof (TOV) equations, which read

$$\begin{array}{cc}\hfill \frac{\mathrm{d}p(r)}{\mathrm{d}r}& =\frac{-GM(r)\epsilon (r)}{r^2{c}^2}[1+\frac{p(r)}{\epsilon (r)}][1+\frac{4\pi p(r)r^3}{M(r)c^2}][1-\frac{2GM(r)}{c^{2}r}{]}^{-1},\hfill \\ \hfill \frac{\mathrm{d}M(r)}{\mathrm{d}r}& =4\pi r^2\frac{\epsilon (r)}{c^2},\hfill \end{array}$$(5)

where p(r) and ε(r) represent the pressure and total energy density (including the rest-mass energy) at radius r, respectively. M(r) denotes the mass enclosed within a radius, r. In the following, we adopt the units G = c = 1.

The rotation of massive NSs would distort the spacetime metric, and produces the frame-dragging effect, which leads to an increase in the static mass. Since the observed MSPs have periods on the order of milliseconds rather than sub-milliseconds, we adopted the slow-rotation Hartle-Thorne approximation, where the spin is expanded to quadratic order (Hartle & Thorne 1968; Berti et al. 2005; Ordaz 2019; Lopes 2024). The parameter ω_r = ω_r(r) = dϕ/dt represents the angular velocity acquired by a free-falling observer due to the dragging of the inertial frame. The perturbed Schwarzschield metric for a slowly rotating star is written as

$$\begin{array}{cc}\hfill \mathrm{d}{s}^2& =-e^{2\phi }(1+2h(r,\theta ))\mathrm{d}{t}^2\hfill \\ \hfill & +[1+\frac{2m(r,\theta )}{r-2M(r)}][1-\frac{2M(r)}{r}{]}^{-1}\mathrm{d}{r}^2\hfill \\ \hfill & +r^2(1+2k(r,\theta ))[\mathrm{d}{\theta }^2+{sin}^2\theta {(\mathrm{d}\varphi -{\omega }_r\mathrm{d}t)}^2],\hfill \end{array}$$(6)

in which the metric function, φ, fulfills the equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}\phi }{\mathrm{d}r}=\frac{M(r)+4\pi r^{3}p(r)}{r(r-2M(r))}.\end{array}$$(7)

h(r, θ), m(r, θ), and k(r, θ) can be considered to be the perturbation terms. Expanding them up to the second order, one obtains:

$$\begin{array}{cc}\hfill h(r,\theta )& =h_0(r)+h_2(r)P_2(cos\theta ),\hfill \\ \hfill m(r,\theta )& =m_0(r)+m_2(r)P_2(cos\theta ),\hfill \\ \hfill k(r,\theta )& =k_0(r)+k_2(r)P_2(cos\theta ),\hfill \end{array}$$(8)

in which P₂(cos(θ) is the Legendre polynominal. The monopole pressure perturbation, p₀^*, and the quadrupole pressure perturbation, p₂^*, are defined in

$$\begin{array}{c}\hfill \mathrm{\Delta }p=(\epsilon +p)[p_0^{\ast }+p_2^{\ast }{P}_2(cos(\theta )].\end{array}$$(9)

Δp denotes the pressure change induced by rotation at a fixed point in the frame corotating with the star. The corresponding perturbed energy density, Δε, and perturbed number density, Δn, reads

$$\begin{array}{cc}\hfill \mathrm{\Delta }\epsilon & =(\epsilon +p)[p_0^{\ast }+p_2^{\ast }{P}_2(cos(\theta )](\mathrm{d}\epsilon /\mathrm{d}p),\hfill \\ \hfill \mathrm{\Delta }n& =(\epsilon +p)[p_0^{\ast }+p_2^{\ast }{P}_2(cos(\theta )](\mathrm{d}n/\mathrm{d}p).\hfill \end{array}$$(10)

The perturbative terms, m₀, p₀^*, and metric functions, ν₂  ≡  h₂ + k₂, h₂, were obtained by applying Einstein’s field equations:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{m}_0}{\mathrm{d}r}& =4\pi r^2\frac{\mathrm{d}\epsilon }{\mathrm{d}p}(\epsilon +p)p_0^{\ast }+\frac{1}{12}{j}^2{r}^4(\frac{d\overline{{\omega }_r}}{\mathrm{d}r}{)}^2-\frac{1}{3}{r}^3\frac{\mathrm{d}{j}^2}{\mathrm{d}r}{\overline{{\omega }_r}}^2,\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{p}_0^{\ast }}{\mathrm{d}r}& =-\frac{m_0(1+8\pi r^{2}p)}{{(r-2M(r))}^2}-\frac{p_0^{\ast }(\epsilon +p)4\pi r^2}{(r-2M(r))}\hfill \\ \hfill & +\frac{1}{12}\frac{r^4{j}^2}{(r-2M(r))}(\frac{\mathrm{d}\overline{{\omega }_r}}{\mathrm{d}r}{)}^2+\frac{1}{3}\frac{\mathrm{d}}{\mathrm{d}r}(\frac{r^3{j}^2{\overline{{\omega }_r}}^2}{r-2M(r)}),\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\nu }_2}{\mathrm{d}r}& =-2h_2(\frac{\mathrm{d}\phi }{\mathrm{d}r})+(\frac{1}{r}+\frac{\mathrm{d}\phi }{\mathrm{d}r})[-\frac{1}{3}{r}^3\frac{\mathrm{d}{j}^2}{\mathrm{d}r}{\overline{{\omega }_r}}^2+\frac{1}{6}{j}^2{r}^4(\frac{\mathrm{d}\overline{{\omega }_r}}{\mathrm{d}r}{)}^2],\hfill \\ \hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{h}_2}{\mathrm{d}r}& =-2h_2(\frac{\mathrm{d}\phi }{\mathrm{d}r})+\frac{r}{r-2M(r)}(2\frac{\mathrm{d}\phi }{\mathrm{d}r}{)}^{-1}[8\pi (\epsilon +p)-\frac{4M(r)}{r^3}]h_2\hfill \\ \hfill & -\frac{4{\nu }_2}{r(r-2M(r))}(2\frac{\mathrm{d}\phi }{\mathrm{d}r}{)}^{-1}\hfill \\ \hfill & +\frac{1}{6}[\frac{\mathrm{d}\phi }{\mathrm{d}r}r-\frac{1}{r-2M(r)}(2\frac{\mathrm{d}\phi }{\mathrm{d}r}{)}^{-1}]r^3{j}^2(\frac{\mathrm{d}\overline{{\omega }_r}}{\mathrm{d}r}{)}^2\hfill \\ \hfill & -\frac{1}{3}[\frac{\mathrm{d}\phi }{\mathrm{d}r}r-\frac{1}{r-2M(r)}(2\frac{\mathrm{d}\phi }{\mathrm{d}r}{)}^{-1}]r^2(\frac{\mathrm{d}{j}^2}{\mathrm{d}r}){\overline{{\omega }_r}}^2,\hfill \end{array}$$(14)

with the boundary conditions: m₀(0) = p₀^*(0) = h₂(0) = ν₂(0) = 0.

At the first order in spin, we could extract the angular momentum. The angular velocity of the star relative to the local inertial frame, $\overline{{\omega }_r}=\mathrm{\Omega }-{\omega }_r$, satisfies the equation (Hartle 1967)

$$\begin{array}{c}\hfill \frac{1}{r^4}\frac{\mathrm{d}}{\mathrm{d}r}(r^{4}j\frac{\mathrm{d}\overline{{\omega }_r}}{\mathrm{d}r})+\frac{4}{r}\frac{\mathrm{d}j}{\mathrm{d}r}\overline{{\omega }_r}=0,\end{array}$$(15)

with

$$\begin{array}{c}\hfill j=j(r)=e^{-\phi }[1-\frac{2M(r)}{r}{]}^{1/2}.\end{array}$$(16)

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

$$\begin{array}{c}\hfill J=\frac{1}{6}{R}^4(\frac{\mathrm{d}\overline{{\omega }_r}}{\mathrm{d}r}){|}_{r=R}.\end{array}$$(17)

At the second order in spin, the rotation causes the NSs to become ellipsoids. The eccentricity is defined as

$$\begin{array}{c}\hfill e=\sqrt{1-(\frac{R_{\mathit{po}}}{R_{\mathit{eq}}}{)}^2},\end{array}$$(18)

where R_po and R_eq are the polar and equator radii, respectively. Thus, e = 0 refers to the spherical limit and 0 < e < 1 corresponds to oblate sphere for rotating NSs. In the Hartle-Thorne approach, the polar and equatorial radii are obtained from the nonrotating star radius, R,

$$\begin{array}{cc}\hfill R_{\mathit{po}}& \approx R+{\xi }_0(R)+{\xi }_2(R),\hfill \\ \hfill R_{\mathit{eq}}& \approx R+{\xi }_0(R)-\frac{1}{2}{\xi }_2(R),\hfill \end{array}$$

in which ξ₀ and ξ₂ are referred to the spherical and quadrupole stretching functions,

$$\begin{array}{cc}\hfill {\xi }_0& =-p_0^{\ast }(\epsilon +p)(\frac{\mathrm{d}p}{\mathrm{d}r}{)}^{-1},\hfill \\ \hfill {\xi }_2& =-p_2^{\ast }(\epsilon +p)(\frac{\mathrm{d}p}{\mathrm{d}r}{)}^{-1}.\hfill \end{array}$$(19)

The quadrupole pressure perturbation, p₂^*, is given by

$$\begin{array}{c}\hfill p_2^{\ast }=-h_2-\frac{1}{3}{r}^2{e}^{-2\phi }{\overline{{\omega }_r}}^2.\end{array}$$(20)

3. Results

In Fig. 1, the pressures of pure neutron matter (PNM) and symmetric nuclear matter (SNM) are plotted applying, respectively, the TM1 and IUFSU parameter families. The family TM1 parameter sets with L = 40, 60, 80, 111 MeV are applied, while the IUFSU family adopt L = 47.2, 60, 80, 110 MeV. The ab initio calculations (Hebeler et al. 2013) and experimental constraints obtained from collective flow data in heavy-ion collisions (Danielewicz et al. 2002) are included in colored areas. The PNM constraint favors smaller values of L at higher densities (ρ_B > 0.06 fm⁻³), but at lower densities, smaller L in turn leads to higher pressure. In the right panel of SNM, the pressure, p, is independent of the symmetry energy slope, L. The SNM constraint area prefers the softer IUFSU family, while the TM1 result lies very close to the edge of the constraint. Both the TM1 and IUFSU parameter families are adjusted based on the properties of finite nuclei and fulfill the saturation property ranges. These parameter sets can help to understand the effect of symmetry energy on the properties of rotating NSs, though some of them exceed the constraint area.

Fig. 1. Pressure as a function of baryon density, ρB, for PNM and SNM. The unit of L is mega-electronvolts. Ab initio results from Hebeler et al. (2013) for PNM and the experimental flow constraint from Danielewicz et al. (2002) for SNM are shown.

In Fig. 2(a), the mass radius relations for different TM1 EOSs are displayed (these results can also be found in Ji et al. 2019; Hu et al. 2020). These results were obtained by using the unified EOSs with L = 40, 60, 80, 111 MeV and the non-unified ones with a core L = 40 MeV matched to crusts with L = 60, 80, 111 MeV. The constraints from recent observations are also shown, including those from PSR J0740+6620 (Riley et al. 2021), PSR J0030+0451 (Riley et al. 2019), PSR J1614-2230 (Arzoumanian et al. 2018), and PSR J0348+0432 (Antoniadis et al. 2013). We find that in the case of the maximum-mass NS, symmetry energy effects are nearly negligible, except for the original TM1 model, which has L = 111 MeV. However, the results with L = 111 MeV exceed the constraint derived from PSR J0030+0451, and therefore the unified EOS with L = 111 MeV is not included in the following discussion. Other EOSs in Fig. 2 (a) lead to very similar maximum masses: approximately M_max ≈ 2.125 M_⊙ for unified EOS with L = 40 MeV and all the non-unified EOSs, and M_max ≈ 2.120 M_⊙ for the unified EOSs with L = 60 and 80 MeV. On the other hand, significant differences can be observed in the radius, especially for NSs with smaller masses. It is shown that the differences between these EOSs increase as the mass of NSs decreases. For the unified EOSs, a smaller slope, L, corresponds to a smaller radius. In contrast, for the non-unified EOSs with the same core EOS (L = 40 MeV) but different crusts, a larger L for the crust EOS corresponds to a smaller radius. This behavior is discussed in Ji et al. (2019), where it is pointed out that a larger L for the crust corresponds to a softer EOS. For a unified EOS, if it has a stiff core, the crust is softer. This L-driven crust EOS stiffness effect on the NS radius is consistent with the results in Fortin et al. (2016), Suleiman et al. (2021), in which the crust EOS differences stem from the glue-density choice. In Fig. 2 (b), we performed a similar calculation using the family IUFSU parameter sets to assess the model dependence. The original IUFSU parameter gives L = 47.2 MeV. Whether the symmetry energy slope, L, is varied within unified or non-unified schemes, the maximum mass of static NSs and the corresponding radius are almost unchanged. The L dependence of the radius is consistent with that found for the TM1 parameter sets. Fig. 2 (c) and (d) show the relative radius differences with respect to the unified EOS L = 40 MeV for family TM1 parameter sets and unified EOS L = 47.2 MeV for family IUFSU parameter sets. A different L results in almost the same maximum mass and a less than 2% difference in radius for the maximum-mass NSs. Even for the 1.4 M_⊙ NSs, the radius relative difference is less than 5%.

Fig. 2. Mass-radius relations for different symmetry energy slopes, L (upper panel). Uncertainties on the radius due to the use of EOSs with different L (lower panel). The unit of L is mega-electronvolts. Results from unified (L (crust) = L (core)) and non-unified EOSs (L (crust)≠ L (core)) are shown.

We now concentrate our study on the bulk properties of rotating MSPs. We implemented numerical codes to solve the mass-radius relation equations by fixing the period, P, with various EOSs, in order to study the effects of the symmetry energy slope, L, on rotating MSPs. In Fig. 3, we present the relationship between the Kepler frequency (mass-shedding limit), $f_K=\frac{{\mathrm{\Omega }}_K}{2\pi }$, and the static mass (gravitational mass of nonrotating NS) of NSs, using an empirical formula for relativistic stars (Haensel & Zdunik 1989; Friedman et al. 1989; Lopes 2024),

Fig. 3. Kepler frequency as a function of the NSs static mass for different EOSs. The horizontal lines correspond to the periods of 1 ms, 2 ms, and 5 ms, and the fastest observed MSP frequency of 716 Hz. The vertical line indicates the canonical mass of 1.4 M⊙.

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_K=\chi \sqrt{\frac{M}{R^3}}.\end{array}$$(21)

Under the Newtonian limit, the value of χ is 1, while under the effect of general relativity, it satisfies the condition χ < 1 (Paschalidis & Stergioulas 2017). Haensel et al. (2009) discussed the value of the parameter χ, whose error holds within a few percent for rotating NSs. In this figure, χ = 0.67 is used, following Paschalidis & Stergioulas (2017). For a canonical 1.4 M_⊙ NS, the period corresponding to the limit frequency will not exceed 1 ms, except for the unified L = 80 MeV result in the right panel, and remain safe at f = 716 Hz. For NSs with static masses of around ∼2 M_⊙, the period can be less than 1 ms while still remaining stable; this is suitable for both TM1 and IUFSU families. A softer EOS, characterized by a smaller slope for unified EOSs or a larger crust slope for non-unified EOSs, results in a larger Kepler frequency. For static mass of M > 0.5 M_⊙, MSPs can be stable with a period of P = 2 ms. However, to achieve P = 1 ms, the stable MSP must have a static mass of M > 1.6 M_⊙ for L = 80 MeV (unified, TM1) and M > 1.5 M_⊙ for L = 40 MeV (unified, TM1). For the non-unified EOS with L = 111 + 40 MeV (crust+core, TM1), the stable static mass lower limit for MSPs is M > 1.4 M_⊙ to achieve P = 1 ms. Comparatively, family IUFSU could support a higher frequency limit with the same static mass and L. This indicates that the minimum stable mass for MSPs decreases with decreasing L or softer EOSs for unified EOSs. From this figure, we can conclude that softer EOSs correspond to a smaller minimum stable static mass for the same period, P. With varying L, the curves exhibit no intersection points for M > 0.15 M_⊙, which differs from the findings in Lopes (2024). However, the above relationship between EOS stiffness and minimum stable static mass for the same P remains consistent with Lopes (2024).

To examine the effects of symmetry energy and its slope on the bulk properties of MSPs, we solved the TOV equations and the equations in Hartle’s formalism. For a given period, a stable mass-radius curve extends from the minimum mass at which mass-shedding instability occurs to the maximum mass. Beyond this maximum mass, MSPs become unstable with respect to axi-symmetric perturbations (Bejger et al. 2007). In Fig. 4, we display the equatorial radius as a function of the gravitational mass for MSPs with fixed period P = 1 ∼ 5 ms (left panels) and the MSP mass as a function of the central density, ρ_c (right panels), using unified TM1 EOSs. In the right panel, the region of high densities beyond the maximum mass point is unstable. As is shown in Appendix A, the results of P = 1 ms are beyond the slow rotation limit. However, they can still be considered as a first estimate of the trend with the rotational period, and the actual effect is expected to be more significant. The purple crosses mark the rotational mass (gravitational mass of rotating NS) at the Kepler limit. The results for static NSs are also shown for comparison. Different values of the slope in unified EOSs, L = 40,  60,  80 MeV, are presented in the upper, middle, and lower panels, respectively. With the same period, P, and rotational mass, the effect of the symmetry energy slope, L, on the equatorial radius, R_eq, is similar to that on the static NS radius. The EOSs with larger L result in a slightly larger R_eq. This effect decreases with increasing rotational mass. Fig. 5 is the same as Fig. 4, but it uses non-unified TM1 EOSs, which share the same core with L = 40 MeV and different crusts with L = 111,  80,  60 MeV in the upper, middle, and lower panels. Since there are only differences in the crust that correspond to the sub-saturation density range, the effect of L on the M_rot − R_eq relation (with same P) is mainly reflected in low-mass MSPs, specifically those with M_rot ≲ 0.5 M_⊙, which could be stable with P ≥ 2 ms. Both Fig. 4 and Fig. 5 indicate that the rotational mass-central density relations (M_rot − ρ_c) show no clear L dependence for M_rot ≳ 1 M_⊙. It is shown that, as the spin period decreases, the curves shift toward a higher mass and larger radius. The results indicate that with faster rotation, the change in the mass with respect to the nonrotating configuration also grows faster. The relevant quantities (ΔM, ΔR) are summarized in Table B.1 and Table B.2. The rotation effect leads to an increase in both the maximum mass (M_max) and the equatorial radius at the maximum mass, R_eq(max), while lowering the central density, ρ_c. The radius of a 1.4 M_⊙ star, denoted as R_(1.4), is significantly influenced by the rotation effect. For a unified EOS with L = 40 MeV, the equatorial radius of the canonical 1.4 M_⊙ star increases from 13.157 km to 13.605 km at a period of P = 2 ms (note that P = 1 ms exceeds the Kepler limit). In the case of a 2 M_⊙ star, the equatorial radius increases from 12.704 km to 13.017 km for P = 2 ms. The change in maximum mass, ΔM, for different L cases shows little variation. The rotation lowers the central density of maximum mass MSPs, ρ_c(max), but the effect is minimal until P < 2 ms.

Fig. 4. Rotational mass-equatorial radius relations and central-density dependence. Left panels show rotational mass-equatorial radius sequences for unified TM1 EOSs with varying symmetry energy slope, L, and spin period, P. Right panels give the rotational mass, Mrot, as a function of central density, ρc. The purple crosses mark the rotational mass at the Kepler limit.

Fig. 5. Same as Fig. 4, but using non-unified TM1 EOSs with a different crust.

In Fig. 6, we show the relationship between the eccentricity, e, of the MSP and the gravitational mass under varying L and period, P. Panels (a), (b), and (c) show results for different periods with a given L = 40, 60, 80 MeV. Panel (d) provides a comparison for different EOSs with the same periods, P = 1, 2, 5 ms. As a measurement of the oblate shape of the star deformation, e = 0 indicates the sphere shape (P = 0). It can be observed that the eccentricity, e, increases as the rotational period, P, decreases for a given mass, and a heavier MSP exhibits a smaller e, indicating less deformation. This trend holds for all masses, but it should be noted that configurations beyond the purple cross-marks exceed the Kepler limit. In particular, for P = 1 ms, the Kepler limit is reached even for canonical NSs. As was mentioned in Fig. 2, the EOSs we used have nearly indistinguishable maximum masses, and there is very little difference among these EOSs on e for MSP with growing mass, especially for those greater than 1.5 M_⊙ (heavier than 2 M_⊙ for P = 1 ms). For low-mass MSPs, the curves in Fig. 6 are clearly separated with a given P, revealing the enhanced crustal influence on the induced eccentricity. For MSPs with masses less than 1.0 M_⊙, a trend emerges that softer EOSs (lower L for unified EOSs, and larger L for the crust) lead to smaller e. When P = 2 ms and M_rot = 1.0 M_⊙, the unified EOS with L = 40 MeV results in e = 0.38, while the one with L = 80 MeV yields e = 0.40. In Table B.1 and Table B.2, we summarize the eccentricities of 1.4 M_⊙ and 2.0 M_⊙ nonrotating NSs and MSPs (e_1.4 and e_2.0). Overall, both e_1.4 and e_2.0 slightly increase as the EOS stiffness grows (from top to bottom in Tables B.1 and B.2) with a fixed P. For example, compared to the results for e_2.0 with unified L = 80 MeV, those with L = 40 MeV and L = 60 MeV show slightly smaller values of e_2.0 for the same period, P, with differences of around 1%.

Fig. 6. Eccentricity, e, as a function of the gravitational mass of the MSPs for the family TM1 parameter set. Panels (a), (b), and (c) represent results from unified EOS with L = 40, 60, 80 MeV, respectively, while panel (d) shows a comparison of different EOSs.

Fig. 7 depicts the results of moment of inertia $I=\frac{J}{\mathrm{\Omega }}$ as a function of the rotational mass. Similarly to Fig. 6, panel (a) shows results for different periods with a given L = 40 MeV. Since the I results for P = 3, 4 ms are too close to those for P = 5 ms, panel (b) provides a comparison of different EOSs with the same periods, P = 5 ms. The moment of inertia, I, depends on both the mass and the square of the radius. At a fixed rotational mass, M_rot, a smaller period, P (faster spin), reduces the central density (see Figs. 4 and 5), making the star more dilute; because J increases sub-linearly with Ω, the moment of inertia I = J/Ω decreases (panel (a)). This also explains why the L effect on I is weakened at P = 5 ms (same trend for P = 4 ∼ 1 ms): the shifts in mass and radius produced by rotation are larger than those arising from varying the symmetry energy slope, L. A softer EOS leads to a smaller I mainly due to a smaller radius. These results are consistent with the mass-radius relations discussed above. The constraint $I_{1.4}=1.{51}_{-0.30}^{+0.20}\times {10}^{45}$ g cm² (Landry et al. 2020) was derived from measurements of pulsar masses, gravitational wave event GW170817 and GW190425, and X-ray data from NICER.

Fig. 7. Moment of inertia, I, as a function of the gravitational mass of the MSPs for the family TM1 parameter set. Panel (a) represents results from unified EOS L = 40 MeV, while panel (b) shows a comparison of different EOSs with P = 5 ms.

4. Conclusions

In this work, we have studied the effect of symmetry energy on the bulk properties of MSPs. We utilized the TM1 parameter set family within the RMF framework for our analysis. The EOSs with a different symmetry energy slope, L, lead to variations in the macroscopic properties of MSPs, such as the maximum mass, central density, Kepler frequency, equatorial radius, eccentricity, and moment of inertia. The IUFSU family exhibits a similar trend to that of the TM1 family on the L dependence. Within the models explored in this work, a soft EOS corresponding to smaller L for a unified EOS or larger L for the crust is favored by the recent observations of NS radius. Additionally, such softer EOSs yield slightly higher Kepler-limit frequencies.

Both the eccentricity, e, and the equatorial radius, R_eq, are sensitive to the symmetry energy slope, L, and the pulsar period, P, but are not directly coupled to each other. They can be treated as independent observables. Measuring P and compactness (from gravitational wave event) allows one to infer e via the universal relation. Together with observed R_eq, the two quantities could provide complementary constraints on L. Moreover, the effects of symmetry energy and its slope are particularly pronounced for less massive NSs. Future measurements of low-mass MSPs may help us constrain the symmetry energy slope more effectively.

With large-scale surveys such as FAST, radio telescopes are discovering an increasing number of pulsars. Measurements of period, mass, and eccentricity could potentially provide further constraints on the symmetry energy slope.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grants No. 12305148, No. 12175109, 12505159), Hebei Natural Science Foundation (No. A2023203055), and Fundamental Research Program of ShanxiProvince (No. 202303021211135), and Natural Science Foundation of Shandong Province under Grants No. ZR2023QA112, and Fundamental Research Funds for the Central Universities under Grant No. 22CX06006A.

References

Appendix A: Comparison with RNS

In appendix A, we compare rotational mass-radius results obtained with the slow rotation approximation to those computed with the RNS code¹ (Stergioulas & Friedman 1995; Stergioulas 2003) in Fig. A.1. The results are obtained using unified TM1 parameter sets with L = 40, 60, 80 MeV. The rotational mass-equatorial radius curves obtained using the slow rotation approximation (dash lines) are slightly lower than those from the RNS code (solid lines) with the same period P. The differences are minimal for P = 5 ms, distinguishable for P = 2 ms (Ω ≈ 3141.6 rad/s), and more pronounced when P = 1 ms. Comparatively, with larger gravitational masses, the differences become smaller, especially for the maximum rotational mass MSPs. These differences arise because the rotational-energy contribution to the gravitational mass is only approximate calculated in the Hartle-Thorne treatment. And the gap between approximation and RNS results becomes more pronounced with rapid rotation. This difference manifests primarily in the equatorial radius (and hence in the eccentricity and moment of inertia). Although such differences exist, the slow-rotation approximation has still proven useful for obtaining a quick estimate of rotational effects. In Fig. 4, purple crosses mark the Kepler limit; compared with those in Fig. A.1 (crosses between solid and Kepler lines), they share nearly the same rotational mass but exhibit a clear equatorial radius offset.

Fig. A.1. The rotational mass-equatorial radius relations derived from RNS code and Hartle-Thorne approximation (H&T).

Appendix B: Properties of rotating MSPs

The following two tables list the properties of nonrotating NSs and rotating MSPs for the TM1 EOSs.

All Tables

All Figures

	Fig. 1. Pressure as a function of baryon density, ρB, for PNM and SNM. The unit of L is mega-electronvolts. Ab initio results from Hebeler et al. (2013) for PNM and the experimental flow constraint from Danielewicz et al. (2002) for SNM are shown.
In the text

	Fig. 2. Mass-radius relations for different symmetry energy slopes, L (upper panel). Uncertainties on the radius due to the use of EOSs with different L (lower panel). The unit of L is mega-electronvolts. Results from unified (L (crust) = L (core)) and non-unified EOSs (L (crust)≠ L (core)) are shown.
In the text

	Fig. 3. Kepler frequency as a function of the NSs static mass for different EOSs. The horizontal lines correspond to the periods of 1 ms, 2 ms, and 5 ms, and the fastest observed MSP frequency of 716 Hz. The vertical line indicates the canonical mass of 1.4 M⊙.
In the text

	Fig. 4. Rotational mass-equatorial radius relations and central-density dependence. Left panels show rotational mass-equatorial radius sequences for unified TM1 EOSs with varying symmetry energy slope, L, and spin period, P. Right panels give the rotational mass, Mrot, as a function of central density, ρc. The purple crosses mark the rotational mass at the Kepler limit.
In the text

	Fig. 5. Same as Fig. 4, but using non-unified TM1 EOSs with a different crust.
In the text

	Fig. 6. Eccentricity, e, as a function of the gravitational mass of the MSPs for the family TM1 parameter set. Panels (a), (b), and (c) represent results from unified EOS with L = 40, 60, 80 MeV, respectively, while panel (d) shows a comparison of different EOSs.
In the text

	Fig. 7. Moment of inertia, I, as a function of the gravitational mass of the MSPs for the family TM1 parameter set. Panel (a) represents results from unified EOS L = 40 MeV, while panel (b) shows a comparison of different EOSs with P = 5 ms.
In the text

	Fig. A.1. The rotational mass-equatorial radius relations derived from RNS code and Hartle-Thorne approximation (H&T).
In the text