Radial oscillations of neutron stars within density-dependent relativistic-mean field models

Neutron stars (NSs) are stellar remnants formed from the supernova explosions of massive stars, with central densities several times greater than the nuclear saturation density. They serve as natural laboratories for investigating the equation of state (EOS) of dense matter under extreme conditions, providing insights into quantum chromodynamics (QCD), and testing the validity of modern nuclear many-body frameworks [1].

Astronomical observation technologies for NSs have undergone significant improvements over the past decades. The successful operation of gravitational wave (GW) detectors by the LIGO Scientific Collaboration and Virgo Collaboration (LVC) in 2015 marked the advent of multi-messenger astronomy with GWs. Subsequently, the detection of the GW170817 event—the first binary neutron star (BNS) merger—by LVC provided an estimate of the tidal deformability for NSs of approximately $ 1.4M_\odot $, introducing a new constraint on the EOS of dense matter.

Most observed NS masses are approximately $ 1.4M_\odot $, although some exceed $ 2M_\odot $. Precise mass measurements of massive pulsars, such as PSR J1614-2230 [2], PSR J0348+0432 [3], and PSR J0740+6620 [4], using the Shapiro effect, require the predicted maximum NS mass to surpass $ 2M_\odot $. Recent simultaneous measurements of the masses and radii of PSR J0030+0451 [5, 6], PSR J0740+6620 [7, 8], PSR J0437-4715 [9, 10], and PSR J1231-1411 [11] by the Neutron Star Interior Composition Explorer (NICER) have provided further constraints on EOS behavior.

Understanding the internal composition of NSs is crucial for constructing their EOS. Oscillations of NSs offer a valuable probe for exploring their internal structures, akin to the role of asteroseismology in studying ordinary stars. These oscillations are categorized into radial and non-radial types, with the latter further divided into modes based on restoring forces [12]. Radial oscillations, the simplest type, involve expansions and contractions around the star's equilibrium structure while preserving its spherical symmetry [13−20].

Although radial oscillations do not directly emit gravitational waves, they can interact with non-radial oscillations, amplifying gravitational wave signals [16]. These interactions may be detectable with next-generation detectors such as the Cosmic Explorer [21] and Einstein Telescope [22]. Furthermore, short gamma-ray bursts (SGRBs) generated by binary NS mergers are influenced by radial oscillations [23]. Consequently, recent studies have examined radial oscillations in NSs considering various factors such as strangeness, quark-hadron phase transitions, dark matter, and proto-neutron stars [18, 19, 24−33].

The EOS of NSs can be derived from nuclear density functional theories, where nucleon–nucleon (NN) interactions are determined by fitting ground-state properties of finite nuclei or empirical saturation properties of infinite nuclear matter [34−39]. The relativistic mean-field (RMF) model, introduced by Walecka using the Hartree approximation [40], has achieved significant success in nuclear physics and astrophysical studies. Enhancements through nonlinear and coupled terms enable the RMF model to accurately describe the ground-state properties of most nuclei [41−44]. The density-dependent RMF (DDRMF) and relativistic Hartree-Fock (DDRHF) models incorporate density-dependent meson-nucleon couplings to account for medium effects [45], making DDRMF parameterizations widely applicable for NS studies. Certain parameter sets accurately describe both finite nuclei and extreme NS configurations [46, 47].

In this study, we systematically investigate the radial oscillations of NSs using EOSs derived from DDRMF parameterizations, which provide accurate ground-state properties for finite nuclei. Most existing numerical methods for solving radial oscillation equations rely on the shooting method, with limited applications of the finite difference method (FDM) [13, 14, 24, 48]. However, if the conservation form of the differential equation is not considered, the FDM fails to handle discontinuities in the adiabatic index of the EOS. Therefore, we discuss this problem and introduce the finite volume method (FVM), which satisfies the conservation property. The radial oscillation equations, FVM, and DDRMF models are introduced in Section II. The results and discussions are presented in Section III. The conclusions are given in Section IV.

A static NS can be regarded as a sphere because of the massive gravitational field, described by the Schwarzschild metric [49]

where $ e^{\lambda} $ and $ e^{\nu} $ are the metric functions.

The Tolman-Oppenheimer-Volkoff(TOV) equations describing the mass and radius of NS [50, 51] can be obtained by solving the Einstein field equation with the above-defined metric as

where the corresponding metric functions, $ \lambda(r) $ and $ \nu(r) $ express as

and they are in the surface of NS related with the total mass, M and radius R

After giving a EOS of NS matter, i.e., $ P-\varepsilon $ relation, the TOV equations as first-order ordinary differential equations can be solved with initial conditions $ m(r=0)=0 $, $ P(r=0)=P_{c} $, and the boundary condition $ P(r=R)=0 $ at surface, R, where $ P_{c} $ is the central pressure of NS.

Considering a spherically symmetric system only with a small perturbative radial movement,

the metric Eq. (1) is no longer time independent. The radial oscillation properties can be obtained from the Einstein field equation based on the static equilibrium structure [52]. One can define the small perturbations of the dimensionless quantities $ \xi=\Delta r/r $ and $ \eta=\Delta P/P $, where $ \Delta r $ and $ \Delta P $ represent the radial displacement and the pressure perturbation. They are governed by [49, 53, 54]

where ω is the eigenfrequency of radial oscillation and Γ is the adiabatic index of NS matter

where $ c_{s} $ is the speed of sound.

In addition, ξ and η in Eq. (6) satisfy two boundary conditions. At the center region, $ r=0 $

and at the surface of NS, $ r=R $

A new quantity can be defined to present the perturbative movement as, $ \zeta=r^{2}e^{-\nu}\Delta r $. The Eq. (6) can be transferred to a linearized radial perturbation equations [14] as

with

The boundary conditions are

where $ \Delta P $ is Lagrangian variation of the pressure defined as

The eigenvalues $ \omega^{2}_{n} $ in Eq. (10) are discrete real numbers, which follow

For the imaginary solution, according to Eq. (5), the star will become unstable [55], while oscillations will be harmonic and stable for real ω. When the central density exceeds the critical density corresponding to the maximum mass, $ \omega_{0} $ will become imaginary. The radial oscillation frequency is defined as

A static NS can be considered a sphere because of the massive gravitational field, described by the Schwarzschild metric [49]

where $ e^{\lambda} $ and $ e^{\nu} $ are the metric functions.

The Tolman-Oppenheimer-Volkoff (TOV) equations describing the mass and radius of NSs [50, 51] can be obtained by solving the Einstein field equation with the above-defined metric as

where the corresponding metric functions, $ \lambda(r) $ and $ \nu(r) $, are expressed as

and they are on the surface of NSs related with the total mass M and radius R:

After obtaining an EOS of NS matter, i.e., the $ P-\varepsilon $ relation, the TOV equations as first-order ordinary differential equations can be solved with initial conditions $ m(r=0)=0 $ and $ P(r=0)=P_{c} $ and the boundary condition $ P(r=R)=0 $ at the surface (R), where $ P_{c} $ is the central pressure of the NS.

Considering a spherically symmetric system only with a small perturbative radial movement,

the metric Eq. (1) is no longer time independent. The radial oscillation properties can be obtained from the Einstein field equation based on the static equilibrium structure [52]. We can define the small perturbations of the dimensionless quantities $ \xi=\Delta r/r $ and $ \eta=\Delta P/P $, where $ \Delta r $ and $ \Delta P $ represent the radial displacement and pressure perturbation, respectively. They are governed by [49, 53, 54]

where ω is the eigenfrequency of radial oscillation, and Γ is the adiabatic index of NS matter:

where $ c_{s} $ is the speed of sound.

In addition, ξ and η in Eq. (6) satisfy two boundary conditions. At the center region, $ r=0 $

and at the surface of NS, $ r=R $

A new quantity can be defined to present the perturbative movement: $ \zeta=r^{2}{\rm e}^{-\nu}\Delta r $. Eq. (6) can be transferred to a linearized radial perturbation equations [14] as

with

The boundary conditions are

where $ \Delta P $ is the Lagrangian variation of the pressure, defined as

The eigenvalues $ \omega^{2}_{n} $ in Eq. (10) are discrete real numbers, which follow

For the imaginary solution, according to Eq. (5), the star will become unstable [55], whereas oscillations will be harmonic and stable for real ω values. When the central density exceeds the critical density corresponding to the maximum mass, $ \omega_{0} $ will become imaginary. The radial oscillation frequency is defined as

The FDM is a discretized scheme to treat the Sturm-Liouville eigenvalue equation such as Eq. (10), which requires the eigenfunction changes smoothly. However, due to crust-core phase transition, the adiabatic index, Γ has a discontinuity at crust region. Once the crust is considered, the eigenvalues obtained by FDM have the obvious differences with those from shooting method. It was found that the finite volume method (FVM), which is applied to solve the heat transfer equation [56] is a good tool to deal with such problem. We found that in the presence of a discontinuous adiabatic index, it is essential to carefully consider the conservative form of the differential equation—an issue that was not addressed in previous studies [13, 14, 24, 48]. Similar to the FDM, a mesh is constructed in FVM. The differential equations form an equilibrium equation at each element of the grid. Using the integral form of the equilibrium equation and the midpoint theorem, a discrete equation in the coordinate space can be obtained.

After solving the TOV equation, the radius in the eigenvalue equation, Eq. (10), is divided into equal lattices N to get the step, h

Other quantities with the lower corner label $ r_{i} $ represent their corresponding values at $ r_{i} $. One can integrate Eq. (10) over the discrete cells between $ [r_{i-1/2},\; r_{i+1/2}] $,

where the first part can be integrated analytically,

$ H \frac{\text{d} \zeta}{\text{d} r} $ is rewritten as,

We integrate it in the interval $ \left [r_{i-1}, r_{i} \right ] $, and take an approximation,

so

where $ \int_{r_{i-1}}^{r_{i}} \dfrac{1}{H} \text{d}r $ can be written approximately to $ \dfrac{r_{i}-r_{i-1}}{H_{r_{i-1/2}}} $. In the same way, we can obtain $ \left( H \frac{\text{d} \zeta}{\text{d} r}\right)_{r_{i+1/2}} $. Meanwhile the other parts can be expressed as

For a uniform grid, $ r_{i}-r_{i-1}=h $, by combining above equations, one can get a discretized equation in lattice space,

which can also be obtained by calculating the derivative of the term $ H\zeta' $ in a staggered grid in FDM.

There are two boundary conditions at $ r=0 $ and $ r=R $,

A matrix form can be generated with Eq. (23),

where A is the tridiagonal matrix with the coefficients,

where

and $ H_{r_{N+1/2}} $ is approximated to $ H_{r_{N}} $.

The FDM obtains the discretized equation by approximating the first and second derivatives of ζ using central differences. One advantage of Eq. (23) is that it preserves the conservation form of the differential equation while allowing discontinuities in the coefficients. Additionally, a non-uniform grid can be adopted, and when integrating the coefficients, one may employ the rectangle rule or Simpson’s rule. These strategies enhance the capability of coarse grids to handle oscillatory equations.

The FDM is a discretized scheme to treat the Sturm-Liouville eigenvalue equation such as Eq. (10), which requires that the eigenfunction changes smoothly. However, owing to crust-core phase transition, the adiabatic index (Γ) has a discontinuity in the crust region. When the crust is considered, the eigenvalues obtained using the FDM have apparent differences from those obtained using the shooting method. The FVM, which is applied to solve the heat transfer equation [56], is an effective tool for solving such problems. We found that, in the presence of a discontinuous adiabatic index, the conservative form of the differential equation must be carefully considered—a challenge that was not addressed in previous studies [13, 14, 24, 48]. Similar to that in the FDM, a mesh is constructed in the FVM. The differential equations form an equilibrium equation at each element of the grid. A discrete equation in the coordinate space can be obtained using the integral form of the equilibrium equation and the midpoint theorem.

After the TOV equation is solved, the radius in the eigenvalue equation, Eq. (10), is divided into equal lattices N to obtain the step, h:

Other quantities with the lower corner label $ r_{i} $ represent their corresponding values at $ r_{i} $. We can integrate Eq. (10) over the discrete cells in the interval $ [r_{i-1/2},\; r_{i+1/2}] $,

where the first part can be integrated analytically,

$ H \dfrac{\text{d} \zeta}{\text{d} r} $ is rewritten as

We integrate it in the interval $ \left [r_{i-1}, r_{i} \right ] $ and take an approximation

so

where $ \int_{r_{i-1}}^{r_{i}} \dfrac{1}{H} \text{d}r $can be expressed approximately as $ \dfrac{r_{i}-r_{i-1}}{H_{r_{i-1/2}}} $. Similarly, we can obtain $ \left( H \frac{\text{d} \zeta}{\text{d} r}\right)_{r_{i+1/2}} $. Meanwhile the other parts can be expressed as

For a uniform grid, $ r_{i}-r_{i-1}=h $, by combining above equations, we can obtain a discretized equation in lattice space:

which can also be obtained by calculating the derivative of the term $ H\zeta' $ in a staggered grid in the FDM.

Two boundary conditions exist at $ r=0 $ and $ r=R $,

A matrix form can be generated using Eq. (23):

where A is the tridiagonal matrix with the coefficients,

where

and $ H_{r_{N+1/2}} $ is approximated to $ H_{r_{N}} $.

The FDM obtains the discretized equation by approximating the first and second derivatives of ζ using central differences. An advantage of Eq. (23) is that it preserves the conservation form of the differential equation while enabling discontinuities in the coefficients. Additionally, a non-uniform grid can be adopted, and when integrating the coefficients, we may employ the rectangle rule or Simpson’s rule. These strategies enhance the capability of coarse grids to solve oscillatory equations.

In our work, the EOSs to solve the radial oscillation equation of NS are generated by the DDRMF model. It is an effective field theory of interacting mesons and nucleons, which takes into account the nucleon-nucleon interaction in the dense matter affected by the nuclear medium. We adopted $ 14 $ groups of parameters of the DDRMF model, DDLZ1[57], DDMEX[58], DDMEX1[59], DDMEX2[59], DDMEXY[59], DDME2[60], DDME1[61], DD2[62], PKDD[63], TW99[64], DDV[65], DDVT[65], DDVTD[65], DDT[66], which can precisely describe the ground-state properties of finite nuclei. The maximum NS mass and corresponding radius from most of these parameters have been confirmed to meet the astronomical observation data [46]. The theoretical framework of DDRMF model was given in detail at Ref. [67]. The Lagrangian of the DDRMF model can be expressed as,

where $ W_{\mu\nu}=\partial_\mu\omega_\nu-\partial_\nu\omega_\mu $ and $ \vec{R}_{\mu\nu}=\partial_\mu\vec\rho_\nu-\partial_\nu\vec\rho_\mu $ are the asymmetric tensor fields. With the mean-field and no-sea approximations, the energy density, ε, and pressure, P, of infinite nuclear matter in DDRMF model can be written as

Here, the contributions from the kinetic energy are

$ \gamma=2 $ is the spin degeneracy factor.

The binding energy per nucleon of the isospin asymmetric nuclear matter, $ E/A=\varepsilon/\rho_N-M_N $, can be expanded in terms of the isospin asymmetry, $ \delta=\left(\rho_n-\rho_p\right) /\rho_N $, as

where $ \rho_N=\rho_n+\rho_p $ is the baryon density with $ \rho_n $ and $ \rho_p $ denoting the neutron and proton densities, respectively. The first term $ E_{\rm SNM}/A $ is the binding energy per nucleon of the symmetric nuclear matter. The symmetry energy $ E_{\text {sym}}(\rho_N) $ can be expressed as

The odd-order terms in δ is attributed to the exchange symmetry between protons and neutrons in nuclear matter. The high-order terms are usually very small so we neglect them.

Around the normal nuclear matter saturation density $ \rho_0 $, $ E_{\rm SNM}/A $ can be expanded, for example, up to fourth-order in terms of $ \chi=\dfrac{\rho_N-\rho_0}{3\rho_0} $,

The first term $ E/A $ is the binding energy per nucleon at the saturation point. The curvature parameter K is the imcompressibility and Q is the skewness coefficients at $ \rho_0 $. They can be written as

The NS is a extreme isospin asymmetry system, where the symmetry energy and its density-dependence play important roles. Similarly, the symmetry energy $ E_{\text {sym}}(\rho_N) $ can also be expanded around $ \rho_0 $ as

where the $ L,\; K_{\text {sym }},\; Q_{\text {sym }} $ are the slope parameter, curvature parameter, and skewness coefficients of the $ E_{\rm sym}(\rho) $ at $ \rho_0 $, respectively, whose definitions are similar to Eqs. 33,

The saturation properties of symmetric nuclear matter with different DDRMF effective interactions are calculated to investigate the relationship between saturation properties and radial oscillation frequencies, which are listed in Table 1, i.e., the saturation density, $ \rho_0 $, the binding energy per nucleon, $ E/A $, incompressibility, K, skewness of the binding energy per nucleon, Q, the symmetry energy, $ E_\text{sym} $, the slope of symmetry energy, L, the curvature and skewness of the symmetry energy, $ K_{\rm sym},\; Q_{\rm sym} $, and the effective neutron and proton masses, $ M^*_n $ and $ M^*_p $. Their binding energy per nucleon, saturation density, and incompressibility are almost same. The symmetry energies and their slopes are consistent with each other except those from DDMEX2 and PKDD sets, whose $ Q_{\rm sym} $ are very small. The effective nucleon masses from DDVT, DDVTD, and DDT are obviously larger than another ones. The tensor coupling between ω meson and nucleon are considered in these three sets.

	$ \rho_{0}[\rm fm^{-3}] $	$ E/A[\rm MeV] $	$ K[\rm MeV] $	$ Q[\rm MeV] $	$ E_{\rm sym}[\rm MeV] $	$ L[\rm MeV] $	$ K_{\rm sym}[\rm MeV] $	$ Q_{\rm sym}[\rm MeV] $	$ M_n^{*}/M_N $	$ M_p^{*}/M_N $
DDLZ1	0.1581	−16.0584	230.3336	1331.7016	32.1895	42.1091	−21.0259	919.4752	0.5581	0.5581
DDMEX	0.1516	−16.1113	267.1415	859.7477	32.3639	49.9555	−73.6292	793.8295	0.5554	0.5554
DDMEX1	0.1506	−16.0366	291.8938	938.0250	31.8400	53.4410	−66.8187	707.2811	0.5709	0.5709
DDMEX2	0.1521	−16.0376	255.7117	631.8805	35.3064	86.8583	−55.3598	99.1595	0.5780	0.5780
DDMEXY	0.1536	−16.0241	268.6004	783.4160	32.0443	53.2243	−74.2635	736.0733	0.5811	0.5811
DDME2	0.1521	−16.1417	251.8909	480.5739	32.3178	51.2763	−87.1828	776.7706	0.5718	0.5718
DDME1	0.1521	−16.2327	244.4246	317.6665	33.0687	55.4385	−101.0120	706.4165	0.5776	0.5776
DD2	0.1491	−16.6679	241.4682	169.3215	31.6642	55.0275	−93.0319	599.0228	0.5627	0.5614
PKDD	0.1496	−16.9145	261.3327	−118.5887	36.7950	90.2315	−80.4556	24.7448	0.5713	0.5699
TW99	0.1531	−16.2471	240.5921	−539.9163	32.7742	55.3160	−124.7304	538.8057	0.5549	0.5549
DDV	0.1516	−16.9278	241.1875	−612.3358	33.6850	69.8717	−97.7363	259.8839	0.5869	0.5852
DDVT	0.1536	−16.9155	239.1848	−743.0882	31.5692	42.3958	−118.8231	896.2710	0.6670	0.6657
DDVTD	0.1535	−16.9165	239.3400	−762.4664	31.8052	42.5800	−117.3942	873.0518	0.6673	0.6660
DDT	0.1551	−16.8727	229.2688	90.7116	31.2197	34.2543	−69.2275	1208.8295	0.6574	0.6560

Table 1.  Nuclear matter properties at saturation density generated by DDRMF parameterizations.

In our study, the EOSs to solve the radial oscillation equation of NSs are generated using the DDRMF model. It is an effective field theory of interacting mesons and nucleons that considers the NN interaction in the dense matter affected by a nuclear medium. We adopt $ 14 $ groups of parameters of the DDRMF model, DDLZ1 [57], DDMEX [58], DDMEX1 [59], DDMEX2 [59], DDMEXY [59], DDME2 [60], DDME1 [61], DD2 [62], PKDD [63], TW99 [64], DDV [65], DDVT [65], DDVTD [65], and DDT [66], which can precisely describe the ground-state properties of finite nuclei. The maximum NS mass and corresponding radius from most of these parameters have been confirmed to satisfy the astronomical observation data [46]. The theoretical framework of the DDRMF model is given in detail in Ref. [67]. The Lagrangian of the DDRMF model can be expressed as

where $ W_{\mu\nu}=\partial_\mu\omega_\nu-\partial_\nu\omega_\mu $ and $ \boldsymbol{R}_{\mu\nu}=\partial_\mu\boldsymbol\rho_\nu-\partial_\nu\boldsymbol\rho_\mu $ are the asymmetric tensor fields. With the mean-field and no-sea approximations, the energy density (ε) and pressure (P) of infinite nuclear matter in DDRMF model can be expressed as

Here, the contributions from the kinetic energy are

$ \gamma=2 $ is the spin degeneracy factor.

The binding energy per nucleon of the isospin asymmetric nuclear matter, $ E/A=\varepsilon/\rho_N-M_N $, can be expanded in terms of the isospin asymmetry, $ \delta=\left(\rho_n-\rho_p\right) /\rho_N $, as

where $ \rho_N=\rho_n+\rho_p $ is the baryon density, with $ \rho_n $ and $ \rho_p $ denoting the neutron and proton densities, respectively. The first term $ E_{\rm SNM}/A $ is the binding energy per nucleon of the symmetric nuclear matter. The symmetry energy $ E_{\text {sym}}(\rho_N) $ can be expressed as

The odd-order terms in δ is attributed to the exchange symmetry between protons and neutrons in nuclear matter. The high-order terms are often very small; therefore, we neglect them.

Around the normal nuclear matter saturation density $ \rho_0 $, $ E_{\rm SNM}/A $ can be expanded, for example, up to fourth-order in terms of $ \chi=\dfrac{\rho_N-\rho_0}{3\rho_0} $:

The first term $ E/A $ is the binding energy per nucleon at the saturation point. The curvature parameter K is the imcompressibility, and Q is the skewness coefficients at $ \rho_0 $. They can be expressed as

The NS is an extreme isospin asymmetry system, where the symmetry energy and its density-dependence play important roles. Similarly, the symmetry energy $ E_{\text {sym}}(\rho_N) $ can also be expanded around $ \rho_0 $ as

where $ L,\; K_{\text {sym }},\; Q_{\text {sym }} $ are the slope parameter, curvature parameter, and skewness coefficients of $ E_{\rm sym}(\rho) $ at $ \rho_0 $, respectively, whose definitions are similar to those in Eqs. (33):

The saturation properties of symmetric nuclear matter with different DDRMF effective interactions are calculated to investigate the relationship between saturation properties and radial oscillation frequencies, which are listed in Table 1, i.e., the saturation density ($ \rho_0 $), binding energy per nucleon ($ E/A $), incompressibility (K), skewness of the binding energy per nucleon (Q), symmetry energy ($ E_\text{sym} $), slope of symmetry energy (L), curvature and skewness of the symmetry energy ($ K_{\rm sym},\; Q_{\rm sym} $), and effective neutron and proton masses ($ M^*_n $ and $ M^*_p $). Their binding energies per nucleon, saturation densities, and incompressibilities are almost the same. The symmetry energies and their slopes are consistent with each other except those from the DDMEX2 and PKDD sets, whose $ Q_{\rm sym} $ values are very small. The effective nucleon masses from DDVT, DDVTD, and DDT are evidently larger than the others. The tensor coupling between the ω meson and nucleon are considered in these three sets.

	$ \rho_{0}/\rm {fm^{-3}} $	$ (E/A) /\rm MeV $	$ K/\rm MeV $	$ Q / \rm MeV $	$ E_{\rm sym}/\rm MeV $	$ L /\rm MeV $	$ K_{\rm sym}/ \rm MeV $	$ Q_{\rm sym}/ \rm MeV $	$ M_n^{*}/M_N $	$ M_p^{*}/M_N $
DDLZ1	0.1581	−16.0584	230.3336	1331.7016	32.1895	42.1091	−21.0259	919.4752	0.5581	0.5581
DDMEX	0.1516	−16.1113	267.1415	859.7477	32.3639	49.9555	−73.6292	793.8295	0.5554	0.5554
DDMEX1	0.1506	−16.0366	291.8938	938.0250	31.8400	53.4410	−66.8187	707.2811	0.5709	0.5709
DDMEX2	0.1521	−16.0376	255.7117	631.8805	35.3064	86.8583	−55.3598	99.1595	0.5780	0.5780
DDMEXY	0.1536	−16.0241	268.6004	783.4160	32.0443	53.2243	−74.2635	736.0733	0.5811	0.5811
DDME2	0.1521	−16.1417	251.8909	480.5739	32.3178	51.2763	−87.1828	776.7706	0.5718	0.5718
DDME1	0.1521	−16.2327	244.4246	317.6665	33.0687	55.4385	−101.0120	706.4165	0.5776	0.5776
DD2	0.1491	−16.6679	241.4682	169.3215	31.6642	55.0275	−93.0319	599.0228	0.5627	0.5614
PKDD	0.1496	−16.9145	261.3327	−118.5887	36.7950	90.2315	−80.4556	24.7448	0.5713	0.5699
TW99	0.1531	−16.2471	240.5921	−539.9163	32.7742	55.3160	−124.7304	538.8057	0.5549	0.5549
DDV	0.1516	−16.9278	241.1875	−612.3358	33.6850	69.8717	−97.7363	259.8839	0.5869	0.5852
DDVT	0.1536	−16.9155	239.1848	−743.0882	31.5692	42.3958	−118.8231	896.2710	0.6670	0.6657
DDVTD	0.1535	−16.9165	239.3400	−762.4664	31.8052	42.5800	−117.3942	873.0518	0.6673	0.6660
DDT	0.1551	−16.8727	229.2688	90.7116	31.2197	34.2543	−69.2275	1208.8295	0.6574	0.6560

Table 1.  Nuclear matter properties at saturation density generated using DDRMF parameterizations.

Initially, the mass-radius (M-R) relations of NSs based on EOSs with 14 DDRMF parameterizations are illustrated in Fig. 1, where the crust EOS is derived from the original TM1 parameter set. The maximum masses of NSs predicted by these EOSs range from approximately $ 1.9M_\odot $ to $ 2.6M_\odot $, consistent with recent astronomical observations [68]. These EOSs can be classified into three types:

Figure 1.  (color online) The mass-radius relations of NS from 14 DDRMF parameter sets.

1. Low-Mass Group: Represented by dashed lines, this group includes EOSs with maximum masses below $ 2.0M_\odot $ and smaller radii at $ 1.4M_\odot $. The skewness coefficient of the binding energy per nucleon, Q, is notably low in this group.

2. Intermediate-Mass Group: EOSs in this category, shown with solid lines, feature maximum masses exceeding $ 2.0M_\odot $ and radii at $ 1.4M_\odot $ within an intermediate range.

3. High-Radius Group: Indicated by dotted lines, this group includes EOSs with the largest radii at $ 1.4M_\odot $. These EOSs exhibit significantly higher values of symmetry energy and its slope compared to the other parameter sets, where the $ Q_{\rm sym} $ are also smaller.

In Fig. 2, the first four eigenfrequencies of NS radial oscillations are presented for the EOS derived from the DDME2 parameterization, solved by the shooting method, FDM, and FVM. The results from the shooting method and FDM are in complete agreement; however, notable discrepancies arise between these methods and FVM, particularly for higher-order modes at larger NS masses. These differences stem from the presence of the crust, where the nodes of higher-order oscillations are primarily concentrated. The discontinuities of the speed of sound within the crust region significantly affect the accuracy of the FDM, which leads to the observed fluctuations in oscillation frequencies, since the speed of sound, $ c^2_s $ in H function in Eq. (10) must be differentiated. Additionally, as the frequency mode increases, the influence of the crust becomes more prominent, and the fluctuations in the results of the FDM become greater.

Figure 2.  (color online) The first four eigenfrequencies of the radial oscillation as functions of NS masses with the EOS obtained by DDME2 parameter set from shooting method, FDM, and FVM.

For the first four eigenfrequencies, the error between the FVM and the shooting method is less than 0.1% when using $ N=3000 $ grid points, and the differences for the first 15 eigenfrequencies remain below 1%. Even with a reduced grid number, $ N=100 $, the fundamental mode ($ \nu_0 $) from FVM remains consistent with the shooting method. The frequency behavior across these modes as a function of NS mass aligns well with previously reported results, such as those in Refs. [30, 69].

The frequencies of the first four modes as functions of NS mass, based on the EOSs from 14 DDRMF parameterizations, are shown in Fig. 3. These frequencies generally increase with NS masses in the lower mass region. For the fundamental mode ($ \nu_0 $) at panel (a), the frequency reaches a maximum value below $ 1.0M_\odot $ and then decreases in the higher mass region for EOSs with smaller symmetry energy. In contrast, for EOSs with larger symmetry energy ($ E_{\rm sym} $) and slope of symmetry energy (L), such as DDMEX2 and PKDD, the maximum value of $ \nu_0 $ occurs around $ 1.5M_\odot $. EOSs with relatively large masses exhibit different variations under different frequency modes. At the maximum NS mass, $ \nu_0 $ becomes zero, satisfying the stability condition of the NS. In contrast, the frequencies of higher modes ($ \nu_1 $, $ \nu_2 $, and $ \nu_3 $) continue to increase with NS mass. These frequencies exhibit a strong correlation with the density-dependent behavior of the symmetry energy in the EOSs, a topic that will be discussed in more detail later.

Figure 3.  (color online) The first four frequencies of the radial oscillation as functions of NS masses obtained with the EOSs from different DDRMF parameterizations.

The first four lowest frequencies at $ 1.4M_{\odot} $ and $ 2M_{\odot} $ for the current DDRMF EOSs are shown at panel (a) and panel (b) in Fig. 4, respectively. We begin by listing the fundamental frequencies from smallest to largest across the different DDRMF parameterizations. However, the higher-order oscillation modes do not display the consistent growth trend as the fundamental mode. Notably, for the higher frequencies at $ 1.4M_\odot $, the softer EOSs tend to have larger frequency differences. Meanwhile DDMEX2, PKDD and DDT, which have a very special skewness in the symmetry energy, also exhibit a staggering behavior. Furthermore, the frequency evolutions for $ \nu_1 $, $ \nu_2 $, and $ \nu_3 $ at $ 2M_{\odot} $ are consistent with each other, and these frequencies do not need to approach zero at the maximum NS mass.

Figure 4.  (color online) The first four frequencies from the DDRMF EOSs at $1.4M_{\odot }$ ((a) panel) and $2M_{\odot }$ are shown ((b) panel).

The frequency difference, $ \Delta \nu_n = \nu_{n+1} - \nu_n $, at $ 1.4M_{\odot} $ as a function of frequency $ \nu_n $ for partial DDRMF parameterizations is shown in Fig. 5. Compared with the EOSs without crust, the EOSs with crust show that $ \Delta\nu_n $ varies unevenly with $ \nu_n $. This discrepancy arises because some oscillation nodes are located in the NS crust at higher-order modes. The non-uniform pasta structure in the crust affects the oscillation frequency due to the non-monotonic behaviour of the adiabatic index. The distribution of the first large separation $ \Delta\nu_0 $ is roughly opposite to the corresponding radius distribution, and the third large separation $ \Delta\nu_2 $ is generally the smallest among them. The formation of a new node in the crust region typically corresponds to a peak in $ \Delta \nu_n $. In general, the frequency difference decreases for higher-order oscillation modes, as there are many nodes in the crust region that have a weaker influence on the frequencies.

Figure 5.  (color online) The relationships between oscillation frequencies $\nu_n$ and their differences $\Delta \nu_n=\nu_{n+1}-\nu_n$ at $1.4M_{\odot }$.

The radial displacement and pressure perturbation, $ \xi(r) $ and $ \eta(r) $ for the $ \nu_3 $ mode as functions of radius at $ 1.4M_{\odot} $, are plotted at panel (a) and panel (b) in Fig. 6, respectively. The cross symbols represent the NS radii at the crust-core transition. They are the eigenfunctions of Eq. (6), with the nodes corresponding to the mode number. Both $ \xi(r) $ and $ \eta(r) $ exhibit more pronounced oscillations for higher modes. Their magnitudes undergo rapid changes between $ 0.9R $ and R, particularly after the crust-core phase transition, with this effect being more prominent for harder EOSs. The last nodes are located in the NS crust region, as discussed in previous studies [69, 70]. Additionally, at $ r \sim 0.7R $, ξ and at $ r \sim 0.55R $, η from different EOSs have similar values.

Figure 6.  (color online) The radial displacement ((a) panel) and pressure perturbation ((b) panel) as functions of NS radius for NS with $1.4M_{\odot }$.

Finally, the correlations between the nuclear saturation properties and oscillation frequencies are investigated using the EOSs from the DDRMF parameter sets. Meanwhile, Skyrme-like models (SLy4, SLy9[71−73], BSK22, BSK24[74−82]) and non-linear relativistic mean-field (NL-RMF) models (TM1[35, 43, 83] and IUFSU[84]) were added for comparison. It is found that the slope of symmetry energy, L, and the skewness coefficient of symmetry energy, $ Q_{\rm sym} $, exhibit a strong linear correlation with the oscillation frequencies of the first excitation mode, $ \nu_1 $, at both $ 1.4 M_\odot $ and $ 2.0 M_\odot $, as shown in Fig. 7. For the slope of symmetry energy, the correlation coefficients with the frequencies are negative, with $ R_{1.4} = -0.9681 $ and $ R_{2.0} = -0.9829 $, respectively. The EOSs with larger values of L from the DDMEX2 and PKDD sets yield lower frequencies around $ 5.8 $ kHz, while for the smallest L from the DDT set, the frequency is around $ 7.0 $ kHz. This suggests that oscillation frequencies could serve as a new probe to constrain the slope of symmetry energy through NS observations. On the other hand, $ Q_{\rm sym} $ has a positive correlation with $ \nu_1 $, with correlation coefficients of $ R_{1.4} = 0.8917 $ and $ R_{2.0} = 0.9335 $. We also notice that the results of SLy4 and IUFSU models are a little deviation from the inverse linear relation. $ Q_{\rm sym} $ corresponds to the third term of the density dependence in symmetry energy, which becomes dominant in the high-density region, but cannot be efficiently measured with current nuclear experiments.

Figure 7.  (color online) The correlations between the frequencies $\nu_{1}$ and the slope L and skewness coefficient $Q_{\rm{sym}}$ of symmetric energy at $1.4M_{\odot }$ ((a, b) panels) and $2M_{\odot }$ ((c, d) panels).

Furthermore, the NS radius exhibits a strong negative linear dependence on the oscillation frequency of the third excitation mode, $ \nu_3 $, at $ 1.4M_\odot $ for the DDRMF models, with a correlation coefficient of $ R = -0.9775 $ as given in Fig. 8. Among all the DDRMF parameters, DDMEX2 set generates the largest NS radius at $ 1.4M_\odot $, $ 13.8 $ km, and the lowest oscillation frequency, $ \nu_3 = 8.7 $ kHz, while the EOS from the DDVTD set produces the smallest NS radius, $ 11.5 $ km, and the highest oscillation frequency, $ \nu_3 = 11.8 $ kHz. This strong correlation suggests that it is possible to determine the radius of a NS by detecting its radial oscillation frequency at higher modes.

Figure 8.  (color online) The relationship between NS radius and oscillation frequency $\nu_{3}$ at $1.4M_\odot$.

The mass-radius (M-R) relations of NSs based on EOSs with 14 DDRMF parameterizations are illustrated in Fig. 1, where the crust EOS is derived from the original TM1 parameter set. The maximum masses of NSs predicted by these EOSs range from approximately $ 1.9M_\odot $ to $ 2.6M_\odot $, consistent with recent astronomical observations [68]. These EOSs can be classified into three types:

Figure 1.  (color online) Mass-radius relations of NSs from 14 DDRMF parameter sets.

1. Low-Mass Group: Represented by dashed lines, this group includes EOSs with maximum masses below $ 2.0M_\odot $ and smaller radii at $ 1.4M_\odot $. The skewness coefficient of the binding energy per nucleon (Q) is notably low in this group.

2. Intermediate-Mass Group: EOSs in this category, shown as solid lines, feature maximum masses exceeding $ 2.0M_\odot $ and radii at $ 1.4M_\odot $ within an intermediate range.

3. High-Radius Group: Indicated by dotted lines, this group includes EOSs with the largest radii at $ 1.4M_\odot $. These EOSs exhibit significantly higher values of the symmetry energy and its slope than the other parameter sets, where the $ Q_{\rm sym} $ values are also smaller.

In Fig. 2, the first four eigenfrequencies of NS radial oscillations are presented for the EOS derived from the DDME2 parameterization, solved using the shooting method, FDM, and FVM. The results from the shooting method and FDM are in complete agreement; however, notable discrepancies occur between these methods and the FVM, particularly for higher-order modes at larger NS masses. These differences result from the presence of the crust, where the nodes of higher-order oscillations are primarily concentrated. The discontinuities of the speed of sound within the crust region significantly affect the accuracy of the FDM, which leads to the observed fluctuations in oscillation frequencies, because the speed of sound ($ c^2_s $) in the H function in Eq. (10) must be differentiated. Additionally, as the frequency mode increases, the influence of the crust becomes more prominent, and the fluctuations in the results of the FDM become greater.

Figure 2.  (color online) First four eigenfrequencies of the radial oscillation as functions of NS masses with the EOS obtained using the DDME2 parameter set via the shooting method, FDM, and FVM.

For the first four eigenfrequencies, the error between the FVM and shooting method is less than 0.1% when using $ N=3000 $ grid points, and the differences for the first 15 eigenfrequencies remain below 1%. Even with a reduced grid number, $ N=100 $, the fundamental mode ($ \nu_0 $) from the FVM remains consistent with that of the shooting method. The frequency behavior across these modes as a function of NS mass corresponds well with previously reported results, such as those in Refs. [30, 69].

The frequencies of the first four modes as functions of NS mass, based on the EOSs from 14 DDRMF parameterizations, are shown in Fig. 3. These frequencies generally increase with NS masses in the lower mass region. For the fundamental mode ($ \nu_0 $) in panel (a), the frequency reaches a maximum value below $ 1.0M_\odot $ and then decreases in the higher mass region for EOSs with smaller symmetry energies. In contrast, for EOSs with larger symmetry energies ($ E_{\rm sym} $) and slopes of symmetry energy (L), such as those from DDMEX2 and PKDD, the maximum value of $ \nu_0 $ occurs at approximately $ 1.5M_\odot $. EOSs with relatively large masses exhibit different variations under different frequency modes. At the maximum NS mass, $ \nu_0 $ becomes zero, satisfying the stability condition of the NS. In contrast, the frequencies of higher modes ($ \nu_1 $, $ \nu_2 $, and $ \nu_3 $) continue to increase with NS mass. These frequencies exhibit a strong correlation with the density-dependent behavior of the symmetry energy in the EOSs, a topic that will be discussed in detail later.

Figure 3.  (color online) First four frequencies of the radial oscillation as functions of NS masses obtained with the EOSs from different DDRMF parameterizations.

The first four lowest frequencies at $ 1.4M_{\odot} $ and $ 2M_{\odot} $ for the current DDRMF EOSs are shown in panels (a) and (b) in Fig. 4, respectively. We begin by listing the fundamental frequencies from smallest to largest across the different DDRMF parameterizations. However, the higher-order oscillation modes do not display the consistent growth trend as the fundamental mode. Notably, for the higher frequencies at $ 1.4M_\odot $, the softer EOSs tend to have larger frequency differences. Meanwhile DDMEX2, PKDD, and DDT, which have a very special skewness in the symmetry energy, also exhibit a staggering behavior. Furthermore, the frequency evolutions for $ \nu_1 $, $ \nu_2 $, and $ \nu_3 $ at $ 2M_{\odot} $ are consistent with each other, and these frequencies do not need to approach zero at the maximum NS mass.

Figure 4.  (color online) First four frequencies from the DDRMF EOSs at (a) $1.4M_{\odot }$ and (b) $2M_{\odot }$.

The frequency difference ($ \Delta \nu_n = \nu_{n+1} - \nu_n $) at $ 1.4M_{\odot} $ as a function of frequency $ \nu_n $ for partial DDRMF parameterizations is shown in Fig. 5. Compared with the EOSs without the crust, the EOSs with the crust show that $ \Delta\nu_n $ varies unevenly with $ \nu_n $. This discrepancy occurs because some oscillation nodes are located in the NS crust at higher-order modes. The non-uniform pasta structure in the crust affects the oscillation frequency owing to the non-monotonic behaviour of the adiabatic index. The distribution of the first large separation $ \Delta\nu_0 $ is generally opposite to the corresponding radius distribution, and the third large separation $ \Delta\nu_2 $ is generally the smallest among them. The formation of a new node in the crust region typically corresponds to a peak in $ \Delta \nu_n $. In general, the frequency difference decreases for higher-order oscillation modes, as the crust region has many nodes that have a weaker influence on the frequencies.

Figure 5.  (color online) Relationships between oscillation frequencies $\nu_n$ and their differences $\Delta \nu_n=\nu_{n+1}-\nu_n$ at $1.4M_{\odot }$.

The radial displacement and pressure perturbation, $ \xi(r) $ and $ \eta(r) $, for the $ \nu_3 $ mode as functions of radius at $ 1.4M_{\odot} $, are plotted in panels (a) and (b) of Fig. 6, respectively. The cross symbols represent the NS radii at the crust-core transition. They are the eigenfunctions of Eq. (6), with the nodes corresponding to the mode number. Both $ \xi(r) $ and $ \eta(r) $ exhibit more pronounced oscillations for higher modes. Their magnitudes undergo rapid changes between $ 0.9R $ and R, particularly after the crust-core phase transition, with this effect being more prominent for harder EOSs. The last nodes are located in the NS crust region, as discussed in previous studies [69, 70]. Additionally, at $ r \sim 0.7R $, ξ values from different EOSs are similar, and at $ r \sim 0.55R $, η values from different EOSs are similar.

Figure 6.  (color online) (a) Radial displacement and (b) pressure perturbation as functions of NS radius for an NS with $1.4M_{\odot }$.

Finally, the correlations between the nuclear saturation properties and oscillation frequencies are investigated using the EOSs from the DDRMF parameter sets. Meanwhile, Skyrme-like models (SLy4, SLy9 [71−73], BSK22, and BSK24 [74−82]) and non-linear relativistic mean-field (NL-RMF) models (TM1 [35, 43, 83] and IUFSU [84]) are included for comparison. We find that the slope of symmetry energy (L) and skewness coefficient of the symmetry energy ($ Q_{\rm sym} $) exhibit a strong linear correlation with the oscillation frequencies of the first excitation mode ($ \nu_1 $) at both $ 1.4 M_\odot $ and $ 2.0 M_\odot $, as shown in Fig. 7. For the slope of symmetry energy, the correlation coefficients with the frequencies are negative, with $ R_{1.4} = -0.9681 $ and $ R_{2.0} = -0.9829 $, respectively. The EOSs with larger values of L from the DDMEX2 and PKDD sets yield lower frequencies at approximately $ 5.8 $ kHz, whereas for the smallest L from the DDT set, the frequency is approximately $ 7.5 $ kHz. This suggests that oscillation frequencies could serve as a new probe to constrain the slope of symmetry energy through NS observations. In contrast, $ Q_{\rm sym} $ has a positive correlation with $ \nu_1 $, with correlation coefficients of $ R_{1.4} = 0.8917 $ and $ R_{2.0} = 0.9335 $. We also notice that the results of SLy4 and IUFSU models deviate slightly from the inverse linear relation. $ Q_{\rm sym} $ corresponds to the third term of the density dependence in the symmetry energy, which becomes dominant in the high-density region but cannot be efficiently measured with current nuclear experiments.

Figure 7.  (color online) Correlations between the frequencies $\nu_{1}$and the slope L and skewness coefficient $Q_{\rm{sym}}$ of symmetric energy at (a, b) $1.4M_{\odot }$ and (c, d) $2M_{\odot }$.

Furthermore, the NS radius exhibits a strong negative linear dependence on the oscillation frequency of the third excitation mode ($ \nu_3 $) at $ 1.4M_\odot $ for the DDRMF models, with a correlation coefficient of $ R = -0.9775 $ as given in Fig. 8. Among all the DDRMF parameters, the DDMEX2 set generates the largest NS radius at $ 1.4M_\odot $ ($ 13.8 $ km) and the lowest oscillation frequency ($ \nu_3 = 8.7 $ kHz), whereas the EOS from the DDVTD set produces the smallest NS radius ($ 11.5 $ km) and the highest oscillation frequency ($ \nu_3 = 11.8 $ kHz). This strong correlation suggests that the radius of an NS can be determined by detecting its radial oscillation frequency at higher modes.

Figure 8.  (color online) Relationship between NS radius and oscillation frequency $\nu_{3}$ at $1.4M_\odot$.

The radial oscillation of NSs was systematically investigated using EOSs from DDRMF parameter sets, which accurately describe the ground-state properties of finite nuclei. We examined the impact of the conservative form on the radial oscillation problem of NS with discontinuous coefficients and introduced the FVM. This approach preserves the conservation properties of the difference equations and yields accurate frequency calculations.

The EOSs from the 14 DDRMF parameterizations were used to calculate NS oscillation frequencies and the radial displacement and pressure perturbation, ξ and η, for various modes, from the fundamental mode to high-order modes. The variation tendencies for the high-order modes were consistent across the EOSs. Frequency differences varied unevenly with oscillation frequencies, and rapid changes in ξ and η were observed in the crust region for higher modes, where nodes of the eigenfunctions appear in the crust.

Correlations between the saturation properties of symmetric nuclear matter and oscillation frequencies were examined for the DDRMF EOSs, as well as six EOSs including Skyrme-like and NL-RMF models. Strong linear correlations were found between the slope of symmetry energy, L and the oscillation frequency of the first excitation mode, $ \nu_1 $, at $ 1.4M_\odot $ and $ 2.0M_\odot $, respectively. Larger values of L correspond to smaller values of $ \nu_1 $. Similarly, there is an inverse linear relationship between the skewness coefficient of symmetry energy, $ Q_{sym} $, and $ \nu_1 $. Additionally, the NS radius at $ 1.4M_\odot $ showed a strong negative correlation with the oscillation frequency of the third excitation mode, $ \nu_3 $. The EOSs of Skyrme-like and NL-RMF models exhibit the same properties as those of the DDRMF model, which indicates the universality of these relationships. These correlations provide new probes for constraining the density-dependent behavior of symmetry energy and suggest a method for measuring the NS radius by detecting its high-order radial oscillation frequencies. Radial oscillations of NSs alone cannot generate GWs. However, they can couple with non-radial oscillations of NSs, amplifying GWs and increasing their detectability. Although current GW detectors are unable to observe GWs in this frequency range, third-generation detectors, such as the Einstein Telescope, are designed to achieve a detection bandwidth spanning from 1 Hz to 10 kHz [22, 85]. Additionally, the Cosmic Explorer is expected to bring significant advancements in GW detection. Radial oscillations of NSs can also modulate short gamma-ray bursts (SGRBs), which may originate from the merger of two NSs. Therefore, detecting SGRBs could provide a means of confirming the oscillation frequencies. A promising detection approach, which utilizes the modulation of SGRBs and requires lower detector sensitivity, has been proposed [23]. However, it is still likely that third-generation GW detectors will be necessary to observe the target frequencies. Thus, it is reasonable to expect that third-generation GW detectors will validate these conclusions.

Future work will focus on non-radial oscillations. However, when addressing non-radial oscillation problems, the unique structure of the equation prevents isolating the frequency on one side and reformulating it as an eigenvalue problem in matrix form. Therefore, we are actively exploring potential methods to fully leverage the capabilities of FVM.

The radial oscillation of NSs was systematically investigated using EOSs from DDRMF parameter sets, which accurately describe the ground-state properties of finite nuclei. We examined the impact of the conservative form on the radial oscillation problem of NSs with discontinuous coefficients and introduced the FVM. This approach preserves the conservation properties of the difference equations and yields accurate frequency calculations.

The EOSs from 14 DDRMF parameterizations were used to calculate NS oscillation frequencies as well as the radial displacement (ξ) and pressure perturbation (η) for various modes, from the fundamental mode to high-order modes. The variation tendencies for the high-order modes were consistent across the EOSs. Frequency differences varied unevenly with oscillation frequencies, and rapid changes in ξ and η were observed in the crust region for higher modes, where nodes of the eigenfunctions appeared in the crust.

Correlations between the saturation properties of symmetric nuclear matter and oscillation frequencies were examined for the DDRMF EOSs, as well as six EOSs including Skyrme-like and NL-RMF models. Strong linear correlations were observed between the slope of symmetry energy (L) and oscillation frequency of the first excitation mode ($ \nu_1 $) at $ 1.4M_\odot $ and $ 2.0M_\odot $, respectively. Larger values of L correspond to smaller values of $ \nu_1 $. Similarly, an linear relationship was observed between the skewness coefficient of symmetry energy ($ Q_{\rm sym} $) and $ \nu_1 $. Additionally, the NS radius at $ 1.4M_\odot $ exhibited a strong negative correlation with the oscillation frequency of the third excitation mode ($ \nu_3 $). The EOSs of Skyrme-like and NL-RMF models exhibit the same properties as those of the DDRMF model, which indicates the universality of these relationships. These correlations provide new probes for constraining the density-dependent behavior of symmetry energy and suggest a method for measuring the NS radius by detecting its high-order radial oscillation frequencies. Radial oscillations of NSs alone cannot generate GWs. However, they can couple with non-radial oscillations of NSs, amplifying GWs and increasing their detectability. Although current GW detectors cannot observe GWs in this frequency range, third-generation detectors, such as the Einstein Telescope, are designed to achieve a detection bandwidth spanning from 1 Hz to 10 kHz [22, 85]. Additionally, the Cosmic Explorer is expected to introduce significant advancements in GW detection. Radial oscillations of NSs can also modulate short gamma-ray bursts (SGRBs), which may originate from the merger of two NSs. Therefore, detecting SGRBs could provide a method of confirming the oscillation frequencies. A promising detection approach, which utilizes the modulation of SGRBs and requires lower detector sensitivity, has been proposed [23]. However, third-generation GW detectors will still likely be necessary to observe the target frequencies. Thus, we reasonably expect that third-generation GW detectors will validate these conclusions.

Future work will focus on non-radial oscillations. However, when addressing non-radial oscillation problems, the unique structure of the equation prevents isolating the frequency on one side and reformulating it as an eigenvalue problem in matrix form. Therefore, we are actively exploring potential methods to fully leverage the capabilities of the FVM.