Radial Oscillations in Neutron Stars from Unified Hadronic and Quarkyonic Equation of States

Abstract

We study radial oscillations in non-rotating neutron stars by considering the unified equation of states (EoSs), which support the 2 M_⊙ star criterion. We solve the Sturm–Liouville problem to compute the 20 lowest radial oscillation modes and their eigenfunctions for a neutron star modeled with eight selected unified EoSs from distinct Skyrme–Hartree–Fock, relativistic mean field and quarkyonic models. We compare the behavior of the computed eigenfrequency for an NS modeled with hadronic to one with quarkyonic EoSs while varying the central densities. The lowest-order f-mode frequency varies substantially between the two classes of the EoS at 1.4 M_⊙ but vanishes at their respective maximum masses, consistent with the stability criterion $\partial M/\partial {\rho }_c>0$. Moreover, we also compute large frequency separation and discover that higher-order mode frequencies are significantly reduced by incorporating a crust in the EoS.

1. Introduction

Neutron stars (NSs) are the collapsing cores of previously massive stars that form after supernova explosions [1]. An NS’s central density is predicted to be around ten times that of the nuclear saturation density ($n_s\approx$ 0.16 fm⁻³). Some unusual phases, such as meson condensation [2,3,4] and quark deconfinement [5,6,7], can be achieved at such high densities. It is impossible to achieve these circumstances on Earth. As a result, studying NSs gives us unique insights into the physics of strongly interacting nuclear matter and phase transitions at ultra-high densities.

Understanding NS also necessitates knowledge of several scientific fields, such as nuclear physics, particle physics, astrophysics and gravitational physics. Even while we know a great deal about how an NS forms, we know relatively little about its internal composition. The one known fact regarding NSs is that characteristics, such as the mass and radius are strongly influenced by the equation of state (EoS) of dense matter. Our overall objective in this field is to identify the properties that may be observed so that we can utilize the generated observational data to construct an accurate EoS, thereby bringing us one step closer to completely comprehending the interior of an NS.

NSs are shown to pulse with different quasi-normal modes (QNMs) in which infinitesimal perturbations induce oscillations whose amplitude decays exponentially with time due to various damping mechanics. Dynamical instabilities, such as mass accretion and tidal forces from a nearby binary companion [8,9], star-quakes generated by fissures in the crust [10], and supernova explosions, are a few sources that might induce oscillations inside an NS [11].

QNMs are categorized into two main groups based on the motion of the pulsations: radial and non-radial oscillation, and each of them is further classified based on the restoring force that acts on the displaced mass element to bring the system back to equilibrium. These restoring forces may be gravity (g-mode), a pressure gradient (p-mode), Coriolis forces (r-mode), magnetic fields and centrifugal force in rotating NS, distinguished by their frequency range. Here, we focus on pressure being our restoring force.

As QNMs are sensitive to the internal composition and EoS of dense matter, we analyze the internal structure of the star and identify the thermodynamic parameters of the NS’s interior using asteroseismology for the study of stellar pulsations in the general relativistic frame. The study of frequencies can thus assist us in indirectly probing within an NS and discovering how strongly interacting nuclear matter behaves at such high densities, thereby further constraining the choice of EoS [12,13,14,15,16,17,18,19,20,21].

Although radial oscillations cannot directly emit gravitational waves (GWs), they can couple with non-radial oscillations, amplifying them and creating a stronger GW that might be detected [17,22]. During the formation of hyper-massive NS through a binary NS merger, an SGRB is emitted, which can be modulated by the radial oscillations [23]. As a result, it is not only valuable in understanding the physics of dense nuclear matter inside an NS but also has some application in GW physics.

Chandrasekhar investigated radial oscillations in NS, which is the simplest mode of oscillation, for the first time in 1964 [24]. Following that, it was researched by other authors, such as Harrison et al. [25] and Chanmugam [26], employing zero-temperature EoS and finite-temperature EoS for proto-NSs by Gondek et al. [27] and strange stars by Benvenuto and Horvath [28] and Gondek and Zdunik [29]. Glass and Lindblom performed the first major investigation of radial oscillations in 1983 [30].

Their numerical results were later rectified by Väth and Chanmugam in 1992 [21] and re-examined by Kokkotas and Ruoff [13] using two alternative numerical approaches that included six more zero-temperature EoS. Their results suggested that oscillations become unstable after the central density at which the NS reaches its maximum mass. This was due to the fact that they used the equilibrium adiabatic index in all of their equations. However, if different adiabatic indices connected to the physical circumstances inside the NS are employed [27] and the slowness of weakly interacting processes is taken into account, stability can be extended beyond that central density [26,31].

The radial modes of oscillation are calculated in this work using eight unified EoSs based on the relativistic mean field (RMF), Skyrme–Hartree–Fock (SHF) [32,33] and quarkyonic models [34]. We solve the Sturm–Liouville eigenvalue problem [35,36,37] with the assumption that our NS is non-rotating and has a zero magnetic field. We base our theory on the fact that the oscillations are small enough to use linear theory. Furthermore, they are adiabatic such that the damping time scale is much longer than the oscillation period [38,39].

This work is structured as follows: in the upcoming Section 2, we go over the theoretical formalism, starting with hydrostatic equilibrium and stellar structure equations in general relativity in Section 2.1, followed by radial oscillation equations in Section 2.2. We discuss our numerical approach in Section 2.2.1. In Section 3, we offer a quick summary of our chosen eight EoS and the reasons for selecting them. In Section 4, we provide our numerical results. Section 4.1 corresponds to the computation of the mass–radius. In Section 4.2, we present our calculated radial mode frequencies in order to analyze them, followed by their chances and methods of detection in Section 4.3. Finally, we summarize and conclude as well as suggest opportunities for further improving our current work in Section 5.

2. Theoretical Formalism

The massive gravitational forces of a non-rotating NS’s interior allow only slight deviations from spherical symmetry, resulting in a nearly perfect formed sphere in its equilibrium condition. As a result, our assumption that the star is spherically symmetric is a reasonable approximation. The gravitational field of such a body is itself spherically symmetric and is given by the Schwarzschild metric in the form of (we use the mostly positive signature $(-+++)$) [40]:

$$d{s}^2=-e^{2\nu }{c}^{2}d{t}^2+e^{2\lambda }d{r}^2+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2),$$

(1)

where $\lambda \equiv \lambda \left(r\right)$ and $\nu \equiv \nu \left(r\right)$ follow their own set of equations. Here, the energy-momentum tensor $T_{\mu \nu }$ has the form of a perfect fluid:

$$T_{\mu \nu }=(P+\epsilon )u_{\mu }{u}_{\nu }+P{g}_{\mu \nu },$$

(2)

where P is the pressure, $\epsilon$ represents the energy density and $u_{\mu }$ is the covariant velocity. Since we have spherical symmetry and are only going to consider motion along the radial direction, only the components $u_0$ and $u_1$ are non-zero.

2.1. Hydrostatic Equilibrium Equations

In a state of hydrostatic equilibrium, all quantities are time-independent. Therefore, even $u_1$ is zero. From Einstein’s field equations, using the Schwarzschild metric in Equation (1) in equilibrium and applying the boundary condition $\lambda (r=0)=0$, we find that:

$$e^{-2\lambda \left(r\right)}=\left(1-\frac{2Gm}{c^{2}r}\right),$$

(3)

where the mass m can be integrated using:

$$\frac{dm}{dr}=\frac{4\pi r^2\epsilon }{c^2}.$$

(4)

Similarly, using the law of conservation of momentum, we obtain [41]:

$$\frac{d\nu }{dr}=-\frac{1}{P+\epsilon }\frac{dP}{dr}.$$

(5)

Finally, using Equation (5) and Einstein’s field equations, we obtain:

$$\frac{dP}{dr}=-\frac{Gm}{c^2{r}^2}\frac{\left(P+\epsilon \right)\left(1+\frac{4\pi r^{3}P}{m{c}^2}\right)}{\left(1-\frac{2Gm}{c^{2}r}\right)}$$

(6)

Equations (3) and (5) define the behavior of the metric functions inside the NS where $r<R$. At the surface, i.e., at $r=R$, they satisfy the boundary condition,

$$e^{2\nu \left(R\right)}=e^{-2\lambda \left(R\right)}=\left(1-\frac{2GM}{c^{2}R}\right).$$

(7)

Equation (7) stays true even outside the star, where R should be replaced by r for $r>R$ as it attains the familiar form of the Schwarzschild solution.

Equations (4) and (6) are collectively known as the Tolman–Oppenheimer–Volkoff (TOV) equations [42,43]. These equations express the equilibrium at each step of the radius r, between the internal pressure of the overlying material against the gravitational force of attraction. These equations can be interpreted if we consider a shell of radius r and thickness $dr$, with a pressure difference of $dp$ in the exterior with respect to its interior and evaluate the net force on each side. The only thing needed to solve the structure equations of NSs is the Equation of State (EoS) of dense matter, i.e., the relation between the pressure and the energy density, which enters the TOV equations.

For a given EoS, the TOV equations can be integrated from the origin with the boundary conditions $m(r=0)=0$ and $P(r=0)=P_c$, where $P_c$ is the central pressure, until the pressure becomes zero. The point R, where the pressure vanishes, provides the circumferential radius of the star. The integration of Equation (4) from zero to R gives its total mass $m\left(R\right)=M$.

2.2. Radial Oscillation Equations

Keeping the spherical symmetry of the background equilibrium configuration, we perturb both fluid and spacetime variables. Assuming a harmonic time dependence for the radial displacement of the fluid element located at position r in the unperturbed model,

$$\delta r(r,t)=X\left(r\right)e^{i\omega t},$$

(8)

the linearized radial perturbation equations can be written as [13]

$$\begin{array}{cc}\hfill c_s^2{X}^{″}+\left({\left(c_s^2\right)}^{\prime }-Z+\frac{4\pi G}{c^4}r\gamma P{e}^{2\lambda }-{\nu }^{\prime }{c}^2\right)X^{\prime }& \\ \hfill +\left[2{\left({\nu }^{\prime }\right)}^2{c}^2+\frac{2Gm}{r^3}{e}^{2\lambda }-Z^{\prime }-\frac{4\pi G}{c^4}(P+\epsilon )Zr{e}^{2\lambda }+{\omega }^2{e}^{2\lambda -2\nu }\right]X& \hfill =0,\hfill \end{array}$$

(9)

where primes denote differentiation with respect to radial coordinate r and $c_s^2=\frac{dP}{d\epsilon }$ is the speed of sound squared in units of $c^2$. $\gamma$ is the adiabatic index, related to the speed of sound by

$$\gamma =\left(\frac{P+\epsilon }{P}\right)c_s^2,$$

(10)

and

$$Z\left(r\right)=c_s^2\left({\nu }^{\prime }-\frac{2}{r}\right).$$

(11)

Now, we re-define the displacement function as

$$\zeta =r^2{e}^{-\nu }X.$$

(12)

Equation (9) can be rewritten for $\zeta$ as

$$\frac{d}{dr}\left(H\frac{d\zeta }{dr}\right)+\left({\omega }^{2}W+Q\right)\zeta =0,$$

(13)

with

$$\begin{array}{cc}\hfill H& \hfill =r^{-2}(P+\epsilon )e^{\lambda +3\nu }{c}_s^2\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill W& \hfill =r^{-2}(P+\epsilon )e^{3\lambda +\nu }\hfill \end{array}$$

(14b)

$$\begin{array}{cc}\hfill Q& \hfill =r^{-2}(P+\epsilon )e^{\lambda +3\nu }\left({\left({\nu }^{\prime }\right)}^2+\frac{4}{r}{\nu }^{\prime }-\frac{8\pi G}{c^4}{e}^{2\lambda }P\right).\hfill \end{array}$$

(14c)

H, W and Q are functions of radial coordinate r and can be calculated using the unperturbed background configuration. Note, Equation (13) explicitly shows its self-adjoint nature. The Lagrangian variation of the pressure now takes the simple form

$$\Delta P=-r^{-2}{e}^{\nu }(P+\epsilon )c_s^2{\zeta }^{\prime }.$$

(15)

The boundary condition at the center is

$$X(r=0)=0,$$

(16)

because the fluid element there cannot be displaced for radial oscillations, and, at the stellar surface, the Lagrangian variation of pressure should vanish

$$\Delta P\left(R\right)=0.$$

(17)

The differential Equation (13) subject to boundary conditions Equations (16) and (17) is a Sturm–Liouville eigenvalue problem. The eigenvalues ${\omega }_n^2$ are real and form an infinite discrete sequence with

$${\omega }_0^2<{\omega }_1^2<{\omega }_2^2<...$$

The eigenfunction of the n-th mode has precisely n number of zeros between the center and the surface of the star. $\omega$ is real for ${\omega }^2>0$; thus, the solution is purely oscillatory. However, for ${\omega }^2<0$, we have an imaginary frequency, which corresponds to an exponentially growing solution. Since the general solution is always a superposition of all such solutions ${\omega }_{n}s$, the presence of an exponentially growing solution corresponds to instability in the radial oscillations. For NSs, the fundamental mode ${\omega }_0$ becomes unstable exactly at the central density ${\rho }_c$ greater than the ${\rho }_{\mathrm{critical}}$ corresponding to the maximum mass configuration. The star will eventually collapse into a black hole in that case.

2.2.1. Numerical Method

For numerical integration, we write Equation (13) as a system of two first-order differential equations in $\zeta$ and $\eta =H{\zeta }^{\prime }$:

$$\frac{d\zeta }{dr}=\frac{\eta }{H}$$

(18)

and

$$\frac{d\eta }{dr}=-\left({\omega }^{2}W+Q\right)\zeta .$$

(19)

Expanding $\zeta$ and $\eta$ close to the origin and comparing the leading order coefficients gives ${\eta }_0=3{\zeta }_0{H}_0$ [13]. Here, ${\eta }_0$, ${\zeta }_0$ and $H_0$ are their corresponding values at $r=r_{\min }$, where $r_{\min }$ is the smallest radial coordinate considered in the integration of the TOV equations. Choosing ${\eta }_0=1$, we find ${\zeta }_0=1/\left(3H_0\right)$. At the surface, boundary condition Equation (17) implies

$$\eta \left(R\right)=0.$$

(20)

For an arbitrary value of $\omega$, the integration can be done from the stellar center to its surface, where $\eta \left(R\right)$ is obtained. Those $\omega$s for which Equation (20) is satisfied are the eigenfrequencies of the radial oscillation.

3. Equation of State

To solve the radial oscillation of the NS, we chose five unified hadronic EoSs based on the RMF and SHF models with different parameterization [32,33] and three quarkyonic EoSs [34,44], viz:

• NL3: The famous NL3 is based on a non-linear interaction, where only the $\sigma$-meson self-coupling term is included, while the cross-coupling terms are not considered [45].
• IOPB: Interaction with higher-order couplings, including self-coupling of $\rho$-mesons and $\omega$–$\rho$ cross-coupling terms [46].
• DD2: A density-dependent interaction with experimental values of proton and neutron mass—$m_p$ and $m_n$. This model can provide an accurate description of the composition and the thermodynamic quantities over a large range of densities [47].
• DDME2: Another effective mean-field interaction with density-dependent meson–nucleon couplings [48].
• SLy4: Based on the Skyrme-Lyon model, this interaction is suitable for calculating the properties of neutron-rich matter. This model can describe both the NS crust as well as the liquid core [49].
• Q1–3: A quark-to-hadron crossover transition model with leptons and nucleons that coexist in quarkyonic phase and degenerate in momentum space. The chosen sets of parameters are as follows [34]:
  • Q1: L = 30 MeV, $\Lambda$= 1400 MeV, $n_t$ = 0.3 fm⁻³.
  • Q2: L = 30 MeV, $\Lambda$= 800 MeV, $n_t$ = 0.3 fm⁻³.
  • Q3: L = 50 MeV, $\Lambda$= 1400 MeV, $n_t$ = 0.4 fm⁻³.
  Here, L is the slope parameter, $\Lambda$ is the baryonic shell thickness parameter, and $n_t$ is the transition density between the nucleonic and quarkyonic phases.

Figure 1 shows the pressure calculated with five unified hadronic EoSs: NL3 [45], IOPB [46], DD2 [47], DDME2 [48] and SLy4 [49] using dashed lines and with three quarkyonic EoSs Q1–3 [34] using solid lines.

Among these chosen EoSs, NL3 yields the stiffer EoS, followed by DDME2 and DD2. SLy4 is initially the softest at lower densities but becomes stiffer than the remaining two EoSs as the density rises. On the contrary, IOPB begins as stiffer than SLy4 but becomes softer with increasing density. At lower densities, Q1 and Q2 begin as softer EoSs, whereas Q3 is stiffer. These quarkyonic EoSs exhibit a crossover transition at their respective transition densities; however, unlike the usual phase transition, which generally softens EoSs, [50], they abruptly stiffen up when quarks flow off of nucleons and fill the lower momenta states, causing a spike in the pressure. Soft-to-stiff transition is a unique feature of crossover transitions; see [51,52], which can be dynamically constructed in the quarkyonic model. We use the above-mentioned EoSs for the calculation of radial oscillations and NS properties.

4. Results and Discussion

In the relativistic frame, the radial modes and global properties for eight separate sets, each with its own EoS, are computed concurrently. The Runge–Kutta method of the fourth order is employed to solve the set of coupled first-order differential Equations (4)–(6) employing a logarithmic step size ranging from 6 m near the center to about 1 m near the surface. This unequal variation aids in recording the rapid change in $\gamma$ near the surface. The solution is iterated until the boundary conditions are fulfilled. Using these solutions, we solve another set of coupled first-order differential Equations (18) and (19) using Equation (14).

Here, the boundary condition, Equation (20) (which is associated with a change of sign) can, thus, be calculated numerically by root-finding methods, such as the bisection method or the shooting method. We implemented the bisection method in our calculation. In all the EoSs used here, at least 20 nodes lie well within $f<50$ kHz. The computational findings will be discussed in detail in the following subsections. It is worth noting that the eigenfrequencies in Section 2.2 were angular frequencies ($\omega$), but we consider linear frequencies (f) where $\omega =2\pi f$ in our results since it is more prevalent in asteroseismology.

Table 1. Frequencies $f_n$ (in kHz) of the 20 lowest radial oscillation modes for our choice of five hadronic and three quarkyonic unified EoSs taken at 1.4 M_⊙.

Order	Hadronic EoSs					Quarkyonic EoSs
n	DD2	DDME2	IOPB	NL3	SLy4	Q1	Q2	Q3
0	3.11	3.25	2.87	2.59	3.01	5.61	4.95	3.45
1	6.43	6.52	6.15	5.56	6.96	8.03	8.11	7.45
2	8.65	8.56	8.28	7.28	9.60	9.63	10.21	9.18
3	9.36	9.31	8.99	8.30	10.81	11.75	11.97	10.71
4	11.48	11.47	11.30	10.57	12.33	14.35	14.21	12.64
5	13.18	13.13	12.74	11.17	14.59	15.88	15.95	14.38
6	14.43	14.42	13.95	13.15	16.54	17.51	17.86	15.79
7	16.18	16.14	15.78	14.35	17.89	19.98	20.51	18.03
8	18.14	18.13	17.51	15.77	19.97	21.95	22.15	19.77
9	19.28	19.22	18.72	16.97	22.08	24.39	23.99	21.46
10	21.12	21.15	20.64	18.83	23.61	25.88	26.47	23.59
11	23.17	23.14	22.17	19.93	25.4	27.71	28.49	25.50
12	24.40	24.31	23.62	21.48	27.59	30.25	30.23	27.19
13	26.03	26.06	25.22	23.38	29.41	32.53	32.52	29.03
14	28.15	28.24	27.27	24.38	30.97	34.18	34.66	30.79
15	29.65	29.53	28.30	26.02	32.99	35.99	36.69	32.67
16	31.07	31.09	30.38	27.46	35.15	38.34	38.67	34.53
17	33.08	33.13	31.89	28.93	36.72	40.75	40.84	36.36
18	34.69	34.54	33.44	30.29	38.55	42.44	42.93	38.32
19	36.15	36.24	34.98	31.92	40.63	44.34	44.88	40.21

lists the frequencies of the 20 lowest radial oscillation nodes, measured at 1.4 M_⊙, for our selection of EoSs.

4.1. M–R Relation

For an EoS, computations begin at the star’s center with the initial parameter central density ${\rho }_c$ as input and continue solving the TOV Equations (4)) and (6) until $P<0$ to obtain the radius R and mass M of an NS. We repeat the procedure by changing ${\rho }_c$ in equal increments and performing the same calculations to obtain the mass–radius (M–R) plot. Figure 2 shows the mass–radius plot for the eight EoSs considered in this work.

NICER helps in the study of exotic matter and NS composition. In 2019, the analysis of the NICER data of PSR J0030+0451 by Miller et al. [53] and Riley et al. [54] gave the measurements of mass M and radius R as $M=1.{44}_{-0.14}^{+0.15}$ M_⊙ and $R=13.{02}_{-1.06}^{+1.24}$ km as well as $M=1.{34}_{-0.16}^{+0.15}$ M_⊙ and $R=12.{71}_{-1.19}^{+1.14}$ km, respectively. These are shown as the lower pink and magenta boxes in Figure 2. Cromartie et al. [57] and Antoniadis et al. [58], through data from radio observations from PSR J0740 + 6620, calculated the NS mass as $M=2.{14}_{-0.09}^{+0.1}$ M_⊙ and $M=2.{01}_{-0.04}^{+0.04}$ M_⊙, respectively. Fonseca et al. [59] calculated the mass of this pulsar using a model-averaged estimation.

This gives the lower limit for the maximum mass of the NS as $M=2.{08}_{-0.07}^{+0.07}$ M_⊙. The NICER data and XMM Newton data of the millisecond pulsar PSR J0740+6620 in 2021 by Miller et al. gave the radius to be $R=13.7_{-1.5}^{+2.6}$ km [55], while Riley et al. calculated it to be $R=12.{39}_{-0.98}^{+1.30}$ km [60]. These are represented by the upper dark and light blue horizontal boxes in the figure. We also have the data from the GW detection of a black hole and a compact object merger, GW190814 [56], which tells us the secondary mass of the compact object is $M=2.{59}_{-0.09}^{+0.08}$ M_⊙ as depicted by the green horizontal band.

Among the quarkyonic EoSs, for a fixed L and $n_t$ values, one having higher $\Lambda$ has greater maximum mass and a more significant bending towards large radii, due to a higher peak in the sound speed. A smaller transition density $n_t$ and a larger L also lead to higher maximum mass and larger radii. As already verified by the authors in refs. [34,45,46,47,48,49], these EoSs are consistent with the astrophysical data from pulsars and GWs and thereby are a good choice for our radial oscillation calculations.

4.2. Radial Oscillation Modes

In Figure 3, we present the radial profile of ${\eta }_n$ and ${\zeta }_n$ for the f-mode (n = 0) as well as excited p-modes (n = 1, 2, 10, 15, 19) in 1.4 M_⊙ NS of SLy4 and Q1 EoS. Following the property of a Sturm–Liouville system, the eigenfunction ${\zeta }_n$ has exactly n zeros (nodes) in the region $0<r<R$. While $\zeta$ is related to the radial displacement function by Equation (12) and using Equations (14), (15) and (18) together, we find that $\eta$ is associated with the Lagrangian variation of pressure by $\eta =-\Delta P{e}^{\lambda +2\nu }$. Similar to ${\zeta }_n$, ${\eta }_n$ also has precisely n zeros (nodes) between the center and the stellar surface.

We can observe the growing amplitude of ${\zeta }_n$ as the radial coordinate increases, whereas ${\eta }_n$ oscillates with a decaying amplitude and vanishes finally at the surface. Since ${\eta }_n$ is continuous and so is $\Delta P$, it follows that the system always oscillates close to equilibrium [61]. The radial functions ${\zeta }_n$ and ${\eta }_n$ are smooth for hadronic SLy4 EoS; however, close to $0.8R$, we observe abrupt changes for quarkyonic Q1 EoS. This abrupt behavior is attributed to the jump in the adiabatic index $\gamma$ due to the crossover transition at $0.3$ fm⁻³. For higher-order modes, some of the nodes move across the core–crust transition and lie in the crust $(0.9R\lesssim r\le R)$, where ${\zeta }_n$ changes signs rapidly with a large amplitude and appears as parallel vertical lines; see the higher modes (n = 10, 15, 19) in Figure 3. ${\eta }_n$, on the other hand, possesses a small amplitude in the crust.

In Figure 4, we investigate the dependence of eigenfrequencies on the central densities ${\rho }_c$ for our chosen EoSs, comparing a set of two EoSs at a time. We plot the first four radial modes and find that when the density increases, regardless of EoS, we approach our stability limit. The instability point is defined by the presence of a zero eigenvalue for f-mode and corresponds to the central density ${\rho }_{\mathrm{critical}}$ where the star approaches its maximum mass $M_{\max }$ [13].

Moreover, with a lower central density, an NS can be modeled as a homogeneous, non-relativistic body [21,62,63] with an angular frequency that follows ${\omega }^2\propto \rho (4\gamma -3)$ [64]. This implies that the frequency is only affected by density because $\gamma$ is fairly constant at lower densities. Thus, regardless of EoS, the frequency approaches zero when the star’s central density is sufficiently low, as illustrated in Figure 4.

Furthermore, we can see a series of ‘avoided crossings’ between separate modes where the frequencies of two consecutive modes from different families repel each other when they approach one another [13,29]. We have two independent families of radial oscillation modes, one at the high density core and the other at the low density envelope, partitioned at the neutron drip density. This neutron drip is linked with any realistic EoSs; hence, the avoided crossings phenomenon, as shown in the figure, is likewise a feature of such EoSs. At the point of the avoided crossings, the solution to the eigenvalue problem shifts from being a standing wave localized primarily in the envelope to one that is localized mostly in the core [29].

We also notice that quarkyonic EoSs have multiple peaks for higher-order modes; see, for example, $n=3$ in Figure 4. As quarkyonic EoSs have a narrow peak of adiabatic index above the transition density, this local peak is a fine interior structure sensitive to higher-order modes. This unique feature in higher-order modes could be a strong indication of a crossover transition, if observed.

In Figure 5, we display the NS’s f-mode frequency vs. mass M to investigate the relationship between radial oscillation and stability more thoroughly. The curves for the quarkyonic EoSs follow the same trend as that of hadronic EoSs, reaching exactly zero at $M_{\max }$ (from the M–R curve). However, in the case of quarkyonic EoSs, the f-mode frequency increases rapidly near 1.4 M_⊙, at which point, quarks begin to drip out from nucleons; however, in the case of hadronic EoSs, it remains practically constant until it falls at $M_{\max }$ [65]. Beyond this limit, ${\omega }^2$ becomes negative; therefore, with an imaginary frequency, the star can no longer recover from minor radial perturbations and ultimately collapses into a black hole. Thus, this result is consistent with the stability condition $\partial M/\partial {\rho }_c>0$ [25,62].

Figure 6 depicts the variation of the so-called large separation, the difference between consecutive nodes $\Delta f_n=f_{n+1}-f_n$, with $f_n$ calculated at 1.4 M_⊙, which aids in understanding the physics of the star interiors and is a commonly used quantity in asteroseismology. When an EoS is considered without a crust, the variation is smooth and consistent with Sagun et al. [66], and the first large separation $\Delta f_0$ is larger than the rest. Furthermore, as we decrease the central density from its critical point, where the star reaches its maximum mass, $\Delta f$ falls as well; see the difference between subsequent nodes in Figure 4. However, when we consider unified EoSs, such as IOPB and Q1, we obtain uneven fluctuations.

This is because the inner crust, the region in which the behavior of the adiabatic index is no longer monotonic, is believed to be characterized by complex structures collectively known as nuclear pasta. Outside nuclear pasta, the adiabatic index is about $\gamma =4/3$ as determined by relativistic electron gas, while, inside the pasta, $\gamma \gtrsim 2$. For the lowest-order mode ($n=0$), the crust does not play a large role in radial oscillation, since it typically accounts for less than 10 percent of the stellar radius, and the oscillation nodes lie deeply in the NS core; see Figure 3. Some of the oscillation nodes are in the crust for higher-order modes, which is equivalent to lower-order modes of NS without a crust.

In other words, in order to reach identical oscillation mode in an NS without crust, NSs with a crust need to have additional nodes in the crust. We also note that, whenever a new node appears in the pasta region, $\Delta f$ shows a peak, meaning that the k-th peak will have k crustal nodes and ($n-k$) nodes inside the core, where n is the total number of nodes. Furthermore, the points between the k-th and ($k+1$)-th peaks will also have k and ($n-k$) nodes in the crust and the core, respectively. As a result, the crust modulates the eigenfrequency significantly.

4.3. Detectability

We may be able to learn about radial oscillation frequencies by examining the emission mechanism of SGRBs, as such oscillations influence the emission from a hypermassive NS generated following a binary NS merger [23]. These short hard gamma ray bursts (SGRBs) that accompany a compact binary system merger should not be confused with soft gamma repeater (SGR) bursts, which are thought to be magneto-emissions from the crust of magnetars [67]. Radial oscillation can also couple with non-radial oscillation and amplify GWs [17,22]. Nonetheless, current GW detectors, such as Advanced LIGO, Advanced Virgo and KAGRA, are projected to have a sensitivity of $\sim 2\times {10}^{-22}$–$4\times {10}^{-24}$ strain/$\sqrt{\mathrm{Hz}}$ at ∼20 Hz–4 kHz [68,69].

Even third-generation ground-based detectors, such as the 40 km long Cosmic Explorer, may have a sensitivity under ${10}^{-25}$ strain/$\sqrt{\mathrm{Hz}}$ above a few kHz [70], while the underground 10 km long Einstein Telescope is expected to have a sensitivity $>3\times {10}^{-25}$ strain/$\sqrt{\mathrm{Hz}}$ at 100 Hz and $\sim 6\times {10}^{-24}$ strain/$\sqrt{\mathrm{Hz}}$ at $\sim 10$ kHz [71]. As a result, neither the current nor the next generation of GW detectors could achieve the requisite sensitivities at the full frequency range in our study. However, the detector sensitivities can be considerably increased by optical reconfiguration or with the use of advanced quantum techniques [68,72]. This could help us reach our target level, or we may have to wait for the fourth generation [73].

5. Conclusions

In this study, we investigated radial oscillations of NSs while considering eight realistic EoSs based on the hadronic RMF, SHF and quarkyonic models [32,33,34]. Quarkyonic EoSs differ from hadronic EoSs by a peak in the speed of sound at which quarks drip out of nucleons. We computed the mass and radius from the TOV equations for each EoS and verified that they are compatible with astrophysical observational evidence from pulsars and GWs [53,54,55,56,57,58]. The radial oscillation equations were then solved considering infinitesimal adiabatic perturbations to calculate the modes of oscillation.

Through our analysis of highly excited modes, we observed that the functions ${\eta }_n\left(r\right)$ and ${\zeta }_n\left(r\right)$ undergo sign changes n times within the star. We compared the behavior of these functions in both hadronic and quarkyonic stars. Notably, in the case of quarkyonic stars, we observed an abrupt change in these functions near a radius of $0.8R$. This is due to the sudden jump in the adiabatic index that occurs during the crossover transition from hadronic matter to a quark–hadron mixture. We also investigated the behavior of the first four radial modes as their central densities change.

Our results for realistic EoSs resemble those of Kokkotas and Ruoff [13], and we also found that quarkyonic EoSs show multiple peaks at higher order. We then examined the nature of only the f-mode frequency with NS mass and observed that $f_0$ falls to zero at the maximum mass corresponding to its critical density. This outcome, regardless of EoS, is compatible with the static stability criterion $\partial M/\partial {\rho }_c>0$. In the M–R relation, quarkyonic EoSs have a significant bend towards a higher radius at higher mass. This unique feature is also shown in the f-mode frequency vs. mass curves. Hadronic EoSs show a flat f-mode frequency; however, the f-mode frequency of quarkyonic EoSs increases after the appearance of quarks, forming a peak.

Furthermore, this peak frequency is higher than the f-mode frequency of any hadronic EoS. If we observe an f-mode frequency of $f>3.5$ kHz, it is a hint of a strong crossover transition. This offers a chance to break the degeneracy of EoSs producing a similar radius, e.g., SLy4 has an identical radius but different f-mode frequency with Q1–3 for canonical mass NSs. In addition, we studied how the nature of the large separation $\Delta f_n$ is smooth and evenly spaced for an EoS without a crust. However, its nature changes dramatically and becomes uneven with eigenfrequencies being squeezed within a shorter range when we include the crust. Finally, we discussed the possibility of detecting such radial frequencies when coupled with GWs through using our current and third-generation ground-based GW detectors.

Our comprehension of the internal structure of an NS will surely be further deepened by a reasonable application of this study in more realistic environments considering rotation, non-zero temperature and magnetic fields [74]. Additionally, research on the potential for radial oscillation frequency detection and efforts to increase detector sensitivity will provide a chance for observational validation of our theory. We hope to explore these options in the near future.