Analysis of Neutron Star f-mode Oscillations in General Relativity with Spectral Representation of Nuclear Equations of State

Neutron stars (NS) are among the densest objects existing in the present Universe with a measured and estimated mass in the range of 1–2 M_⊙, confined in a sphere of radius ∼ 10 km. These tiny but massive stars harbor matter under extreme physical conditions that remain inaccessible to direct observation. Hence, our knowledge of an NS's interior, including the behavior of matter at such extreme densities, remains theoretical and is the subject of ongoing research and modeling. However, the new generations of telescopes have provided valuable observation data in the recent past and helped refine the theoretical description—called the equation of state (EOS). The EOS of matter inside an NS plays a critical role in determining various properties of the NS, such as its mass, radius, moment of inertia, tidal deformability, oscillation frequencies, and the structure of the internal layers. These are essential for understanding the observational characteristics of NSs, including their gravitational wave (GW) signals, electromagnetic radiation, and thermal emissions. The maximum observed mass of an NS, ∼2 M_⊙, ruled out the softer EOS (Demorest et al. 2010; Antoniadis et al. 2013; Cromartie et al. 2019; Fonseca et al. 2021). The International Space Station–based X-ray telescope and instrument package, the Neutron Star Interior Composition Explorer (NICER; Gendreau & Arzoumanian 2017) measures the mass and radius of a pulsar simultaneously for the first time. Electromagnetic radiation from the hot spots at the magnetic poles of a rotating NS provide us with a direct route of observation when the radiation crosses our line of sight (Miller et al. 2021; Riley et al. 2021). The detection of GWs from NS mergers events have provided stringent limits on tidal deformability for an NS with a canonical mass of 1.4 M_⊙, adding valuable insights into the EOS (Abbott et al. 2018). The nuclear physics data at the moderate densities combined with recent NICER observations can further improve existing constraints. Additionally, new X-ray observatories like the European Space Agency's Athena mission (Nandra et al. 2013) and NASA's Lynx mission (Kouveliotou et al. 2014) are expected to offer even higher-resolution X-ray observations of NSs. The third generation Einstein Telescope (Branchesi et al. 2023) will have a sensitivity 10 times better compared to second-generation detectors. Advanced LIGO and Virgo could be optimized for the detection of gravitational signals at low frequency. Ongoing and future surveys, such as the Square Kilometre Array (Kramer & Stappers 2015; Watts et al. 2015), are poised to discover more pulsars and provide a better understanding of the population of NSs and their properties. The coordinated observations of NS mergers at multiple wavelengths, such as X-rays, gamma-rays, and radio waves will definitely help enhance our understanding of the EOS of dense matter.

The NS core is believed to be made up of neutrons, along with an admixture of protons and leptons. There are several EOS models, which can be broadly categorized into nonrelativistic and relativistic models. Among them, the relativistic mean-field (RMF) model varieties have been successfully applied to study NS properties as well as finite nuclei properties. Over the decades additional mesons and nonlinear self-interaction as well as cross-coupling terms for all the mesons have been added to the original RMF, i.e., the Walecka (1974) Model. Various sets of tabulated RMF EOSs, compatible with recent astrophysical observations (Nandi et al. 2019), are available, ³ depending upon the form of interactions in the Lagrangian density. The associated couplings for the interactions are fitted to various nuclear matter parameters such as binding energy per particle, incompressibility, the symmetry energy coefficients and their derivatives, which are determined from terrestrial experiments at nuclear saturation density. Nevertheless the EOS provides a reasonable description of nuclear matter at suprasaturation densities. The nuclear matter parameters define the NS EOS and are correlated to NS observables (Agrawal et al. 2021; Ghosh et al. 2022; Pradhan et al. 2023a; Malik et al. 2023). At these high densities and extreme pressures, different models consider the possibility of strange matter like hyperons, quarks, and Bose–Einstein condensates (Banik et al. 2014; Char et al. 2015; Malik et al. 2021). It is a challenge to unravel the NS core composition and precisely determine its EOS owing to our limited knowledge of matter under extreme physical conditions.

In the absence of any universally acceptable model EOS, subject to limited observational data, several recent works have used spectral representations of realistic EOS as an accurate and efficient way of parameterizing the high-density sections of the NS EOS. (Abbott et al. 2018; Miller et al. 2019; Raaijmakers et al. 2020). This is a powerful theoretical tool, used to describe the EOS in terms of spectral functions, incorporating our current knowledge of fundamental physics. The tabulated version of the realistic EOS often infuses numerical errors in the calculation of the relevant thermodynamic quantities (Servignat et al. 2023). A robust analytic representation of the EOS like spectral representation (Lindblom 2010, 2018, 2022), generalized piecewise polytropic representation (O'Boyle et al. 2020), or speed of sound parameterization (Tews et al. 2018; Greif et al. 2019) becomes quite useful in scenarios where we need numerical determination of thermodynamic quantities associated with a cold β-equilibrated EOS at any arbitrary energy density inside the NS.

The structure parameters like mass, radius, and tidal deformability have taken crucial roles to address the uncertainties in NS EOS. Quasiperiodic oscillations are also valuable tools for studying compact objects and the processes occurring in their binary systems. NSs can exhibit a variety of oscillation modes mainly observed in the X-ray light curves of NSs in various types of binary systems, including low-mass X-ray binaries and NS X-ray binaries. A general relativistic description of oscillations due to nonradial perturbations involves what are called the quasi-normal modes (QNMs). Einstein equations govern the dynamics both inside and outside the star. QNMs are different than the normal modes describing Newtonian oscillations of matter inside the star. These modes are damped due to the emission of GWs from the star. The emitted GWs, given by the perturbations of the spacetime metric tensor, carry information about the internal structure and composition of the NS. The real part of the complex-valued QNM frequency represents the actual frequency of the oscillations, and the imaginary part corresponds to an exponential damping of the wave. Outside the star, determination of the frequency and damping time of the QNMs becomes an eigenfrequency problem subject to appropriate boundary conditions both at the surface of the NS and at spatial infinity. The internal structure and dynamics of the matter inside it influence the frequencies and damping times of the QNMs via the boundary condition at the surface of the NS. Inside the star, the oscillations of the matter and the GWs are coupled to each other via the Einstein equations.

The study of QNMs of oscillating astrophysical objects started with spherically symmetric black hole spacetimes (see Regge & Wheeler 1957; Vishveshwara 1970; Zerilli 1970) and was later extended to compact astrophysical stars (Thorne & Campolattaro 1967) including NSs (see the review by Kokkotas & Schmidt 1999 for references). The modes of oscillation inside the NS are characterized based on the physics behind the restoring force. For example, the pressure of the fluid (which is generally considered incompressible) in the NS interior acts as a restoring force for p-modes. Buoyancy restores the so-called gravity g-modes, which originate from internal temperature or composition gradients. The most fundamental mode of oscillation is the f-modes, which arise out of density perturbations. They typically have a frequency of around 2.4 kHz (Kokkotas et al. 2001; Rezzolla 2003) and have the strongest coupling with the emission of GWs. Determination of the eigenfrequencies while ignoring the metric perturbations during fluid oscillations (relativistic Cowling approximation) to simplify the calculations overestimates the frequency of certain fundamental modes (see the discussion on the accuracy of this approximation in Sotani & Takiwaki 2020). We can improve the accuracy by introducing the emission of GWs and using a post-Newtonian approximation to estimate the damping rate of the f-modes as a result. This way, one can obtain scaling relations between the bulk properties of the star, e.g., compactness, mass, moment of inertia, and mode frequency, and the damping time (Detweiler 1975; Andersson & Kokkotas 1996, 1998; Kokkotas et al. 2001; Benhar et al. 2004; Andersson 2019). When the problem is treated general relativistically, we get f-modes for the QNM.

Now, several works in the literature have demonstrated that for f-modes, the relativistic analysis also leads to very similar scaling relations. Such scaling relations can thus be taken as universal (see, for example, Andersson & Kokkotas 1996, 1998; Lattimer & Prakash 2001; Tsui & Leung 2005; Lau et al. 2010; Yagi & Yunes 2013). This feature of universality can be used as a landmark. One can check how well the calculated frequencies and damping times of the QNMs for a given EOS fit the scaling relations as a sort of consistency requirement (for details, see Andersson 2021). This is a strategy we consistently follow throughout this paper. Another feature of the NS, tidal deformability, also encodes important information about the shape of the star. An external tidal gravitational field acting on the NS contributes to the quadrupole and higher moments of the matter distribution inside (for more details, see Chan et al. 2014; Poisson & Will 2014; Andersson & Pnigouras 2020, 2021; Pratten et al. 2020).

Apart from the f-modes, there are also other modes that typically correspond to one or more additional physical processes. These modes mostly have frequencies higher than those of f-modes. For example, the frequency of the g-modes, which appear when the temperature and/or the pressure of the NS vary from one point inside the star to another, can lie in a range starting from the order of 100 kHz (Reisenegger & Goldreich 1992). For this reason, one expects g-modes to appear when a phase transition occurs inside the NS (Andersson 2019). Apart from these, p-modes and r-modes can also occur due to pressure and spin-related instabilities in the NS, respectively (for details, see the review by Andersson 2019 and the references therein). However, in this work, we focus only on f-modes for the NS, determining their frequencies and damping times from fully general relativistic calculations.

Recent publications (Sotani 2021; Kunjipurayil et al. 2022; Pradhan et al. 2022, 2023b) show efforts to estimate nuclear matter parameters using stellar oscillation modes. Pradhan et al. (2022) reported how the effective nucleon mass (m^⋆) at saturation is correlated with the f-mode frequencies for a set of EOSs generated within the RMF formalism. A strong correlation of the stellar oscillation frequencies with symmetry energy parameters like J_sym,0, K_sym,0, L_sym,0, and so on would allow us to infer the composition of matter at very high densities in the NS interior (Patra et al. 2023).

In this work we have prepared a set of ∼9000 nucleonic EOSs though a Bayesian analysis of RMF model parameters. Constraints of recent interest like the chiral effective field theory (χEFT) and perturbative quantum chromodynamics (pQCD) at relevant energy densities are imposed to infer the parameter set. Then we determine a highly accurate spectral representation of the generated set and numerically solve the fully general relativistic perturbation equations for a set of equilibrium configuration NSs with varying mass for each EOS. We report the corresponding f-mode frequencies and the GW damping times over the EOS set and finally explore correlations of the frequencies and damping times with nuclear matter parameters and stellar observables.

The paper is organized in the following manner. Section 2 provides a brief overview of the RMF formalism for the nuclear EOS at zero temperature followed by a brief description framework for spectral fitting of the EOS, and the formalism for nonradial quadrupolar oscillations of a nonrotating NS. Section 3 contains the results for a set of parameterized spectral versions of the entire EOS set. We discuss the results and conclude our work in Section 4.

2.1. Relativistic Mean-field Equation of State

We generate a set of EOSs within a description of nuclear matter based on a relativistic field theoretical approach where the nuclear interaction between nucleons is introduced through the exchange of the scalar–isoscalar meson σ, the vector–isoscalar meson ω, and the vector–isovector meson ϱ.

The Lagrangian density is given by (Dutra et al. 2014; Malik et al. 2023)

$$\begin{eqnarray}&&{ \mathcal L }={{ \mathcal L }}_{N}+{{ \mathcal L }}_{M}+{{ \mathcal L }}_{\mathrm{NL}}+{{ \mathcal L }}_{\mathrm{leptons}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray*}&&{{ \mathcal L }}_{N}=\bar{{\rm{\Psi }}}\left[{\gamma }^{\mu }\left(i{\partial }_{\mu }-{g}_{\omega }{\omega }_{\mu }-\displaystyle \frac{1}{2}{g}_{\varrho }{\boldsymbol{t}}\cdot {{\boldsymbol{\varrho }}}_{\mu }\right)-\left({m}_{N}-{g}_{\sigma }\sigma \right)\right]{\rm{\Psi }},\end{eqnarray*}$$

denotes the Dirac equation for the nucleon doublet (neutron and proton) with bare mass m_N , Ψ is a Dirac spinor, γ^μ are the Dirac matrices, and t is the isospin operator. ${{ \mathcal L }}_{M}$ represents the mesons, which is given by

$$\begin{eqnarray*}\begin{array}{rcl}{{ \mathcal L }}_{M} & = & \displaystyle \frac{1}{2}\left[{\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -{m}_{\sigma }^{2}{\sigma }^{2}\right]-\displaystyle \frac{1}{4}{F}_{\mu \nu }^{(\omega )}{F}^{(\omega )\mu \nu }+\displaystyle \frac{1}{2}{m}_{\omega }^{2}{\omega }_{\mu }{\omega }^{\mu }\\ & & -\ \displaystyle \frac{1}{4}{{\boldsymbol{F}}}_{\mu \nu }^{(\varrho )}\cdot {{\boldsymbol{F}}}^{(\varrho )\mu \nu }+\displaystyle \frac{1}{2}{m}_{\varrho }^{2}{{\boldsymbol{\varrho }}}_{\mu }\cdot {{\boldsymbol{\varrho }}}^{\mu }{\boldsymbol{,}}\end{array}\end{eqnarray*}$$

where F^{(ω,ϱ)μ ν} = ∂^μ A^(ω,ϱ)ν − ∂^ν A^(ω,ϱ)μ are the vector meson tensors.

$$\begin{eqnarray*}\begin{array}{rcl}{{ \mathcal L }}_{\mathrm{NL}} & = & -\displaystyle \frac{1}{3}b\,{m}_{N}\,{g}_{\sigma }^{3}{\left(\sigma \right)}^{3}-\displaystyle \frac{1}{4}c{\left({g}_{\sigma }\sigma \right)}^{4}+\displaystyle \frac{\xi }{4!}{g}_{\omega }^{4}{\left({\omega }_{\mu }{\omega }^{\mu }\right)}^{2}\\ & & +\,{{\rm{\Lambda }}}_{\omega }{g}_{\varrho }^{2}{{\boldsymbol{\varrho }}}_{\mu }\cdot {{\boldsymbol{\varrho }}}^{\mu }{g}_{\omega }^{2}{\omega }_{\mu }{\omega }^{\mu },\end{array}\end{eqnarray*}$$

contains nonlinear terms with parameters b, c, ξ, and Λ_ω to take care of the high-density behavior of the matter. The g_i values are the couplings of the nucleons to the meson fields i = σ, ω, and ϱ, with masses m_i . Finally, the Lagrangian density for the leptons is given as ${{ \mathcal L }}_{\mathrm{leptons}}=\bar{{{\rm{\Psi }}}_{l}}\left[{\gamma }^{\mu }\left(i{\partial }_{\mu }-{m}_{l}\right){{\rm{\Psi }}}_{l}\right]$, where Ψ_l (l = e⁻, μ⁻) denotes the lepton spinor; leptons are considered noninteracting. The equations of motion for the mesons fields are obtained from Euler–Lagrange equations

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{{g}_{\sigma }}{{m}_{\sigma ,\mathrm{eff}}^{2}}\displaystyle \sum _{i}{\rho }_{i}^{s},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\omega =\displaystyle \frac{{g}_{\omega }}{{m}_{\omega ,\mathrm{eff}}^{2}}\displaystyle \sum _{i}{\rho }_{i},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\varrho =\displaystyle \frac{{g}_{\varrho }}{{m}_{\varrho ,\mathrm{eff}}^{2}}\displaystyle \sum _{i}{t}_{3}{\rho }_{i},\end{eqnarray} \tag{ 4 }$$

where ${\rho }_{i}^{s}$ and ρ_i are, respectively, the scalar density and the number density of nucleon i. The effect of the nonlinear terms on the magnitude of the meson fields is clearly shown through the following terms

$$\begin{eqnarray}&&{m}_{\sigma ,\mathrm{eff}}^{2}={m}_{\sigma }^{2}+{{bg}}_{\sigma }^{3}\sigma +{{cg}}_{\sigma }^{4}{\sigma }^{2},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{m}_{\omega ,\mathrm{eff}}^{2}={m}_{\omega }^{2}+\displaystyle \frac{\xi }{3!}{g}_{\omega }^{4}{\omega }^{2}+2{{\rm{\Lambda }}}_{\omega }{g}_{\varrho }^{2}{g}_{\omega }^{2}{\varrho }^{2},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{m}_{\varrho ,\mathrm{eff}}^{2}={m}_{\varrho }^{2}+2{{\rm{\Lambda }}}_{\omega }{g}_{\omega }^{2}{g}_{\varrho }^{2}{\omega }^{2},\end{eqnarray} \tag{ 7 }$$

and are solved using the RMF approximation with the constraint of β equilibrium, i.e., μ_n − μ_p = μ_e , where the μ values are the chemical potential of the respective particles. Also, total baryon number (n_p + n_n ) conservation and the charge neutrality condition n_p = n_e are enforced.

The pressure and energy density of the baryons and leptons are given by the following expressions

$$\begin{eqnarray}\begin{array}{rcl}\epsilon & = & \displaystyle \sum _{i=n,p,e,\mu }\displaystyle \frac{1}{{\pi }^{2}}{\displaystyle \int }_{0}^{{k}_{\mathrm{Fi}}}\sqrt{{k}^{2}+{m}_{i}^{* }}\,{k}^{2}\,{dk}\\ & & +\ \displaystyle \frac{1}{2}{m}_{\sigma }^{2}{\sigma }^{2}+\displaystyle \frac{1}{2}{m}_{\omega }^{2}{\omega }^{2}+\displaystyle \frac{1}{2}{m}_{\varrho }^{2}{\varrho }^{2}\\ & & +\ \displaystyle \frac{b}{3}{\left({g}_{\sigma }\sigma \right)}^{3}+\displaystyle \frac{c}{4}{\left({g}_{\sigma }\sigma \right)}^{4}+\displaystyle \frac{\xi }{8}{\left({g}_{\omega }\omega \right)}^{4}+{{\rm{\Lambda }}}_{\omega }{\left({g}_{\varrho }{g}_{\omega }\varrho \omega \right)}^{2},\end{array}\end{eqnarray} \tag{ 8 }$$

where ${m}_{i}^{* }={m}_{i}-{g}_{\sigma }\sigma$ for protons and neutrons and ${m}_{i}^{* }={m}_{i}$ for electrons and muons, and k_Fi is the Fermi moment of particle i. The pressure is determined from the thermodynamic relation for every component

$$\begin{eqnarray}&&p=\displaystyle \sum _{i}{\mu }_{i}{\rho }_{i}-\epsilon .\end{eqnarray} \tag{ 9 }$$

We generate our EOS set using a Bayesian setup with minimal constraints as done in Malik et al. (2023). These constraints include nuclear saturation properties such as a saturation density of ρ₀ = 0.153 ± 0.005 fm⁻³, a binding energy per nucleon of e₀ = −16.1 ± 0.2 MeV, a nuclear incompressibility of K₀ = 230 ± 40 MeV, and a symmetry energy of J_sym,0 = 32.5 ± 1.8 MeV at the saturation density. Furthermore, we impose a low-density constraint for the pure neutron matter pressure from EFT (Hebeler et al. 2010) with twice the uncertainty and demand that the maximum NS mass corresponding to the EOS is above the observed 2 M_⊙. Additionally, we employ a pQCD-derived constraint on the pressure at 8 times the baryon density (n_s = 0.16 fm⁻³) for a pQCD scale X = 4 (Komoltsev & Kurkela 2022). Finally, we choose ∼9000 nucleonic EOSs for the present calculation from the posterior.

2.2. Spectral Representation of the Equation of State

We find the spectral parameters for our set of EOSs by fitting the tabular RMF EOSs with a spectral version. The methodology in Lindblom (2010) leads us to faithful representations of realistic EOSs for nuclear matter inside the NS by constructing spectral expansions of key thermodynamic quantities. Linear combinations of a complete set of functions like the Fourier basis functions allow us to expand a thermodynamic quantity like the adiabatic index, defined as

$$\begin{eqnarray}&&{\rm{\Gamma }}(p)=\displaystyle \frac{\epsilon +p}{p}\displaystyle \frac{{dp}}{d\epsilon }.\end{eqnarray} \tag{ 10 }$$

The coefficients of such basis functions in the expansion uniquely determine any given physically valid EOS. One can find that such a representation satisfies the basic necessary thermodynamic conditions. To obtain a pressure-based form, we expand the adiabatic index Γ(p) in terms of the basis functions Φ_k (p) as

$$\begin{eqnarray}&&{\rm{\Gamma }}(p)=\exp \left[\displaystyle \sum _{k}{\gamma }_{k}{{\rm{\Phi }}}_{k}(p)\right].\end{eqnarray} \tag{ 11 }$$

We choose Φ_k (p) = Φ_k (x), where x is a dimensionless logarithmic pressure variable $x=\mathrm{log}(p/{p}_{0})$ as in Lindblom (2010), and p₀ is the minimum pressure where the high-density part of the EOS is joined to the low-density part of the Baym–Pethick–Sutherland (Baym et al. 1971) EOS. The γ_k values are the spectral coefficients.

Numerical integration of the following gives us (p) for the obtained expansion of Γ(p)

$$\begin{eqnarray}&&\epsilon (p)=\displaystyle \frac{{\epsilon }_{0}}{\mu (p)}+\displaystyle \frac{1}{\mu (p)}{\int }_{{p}_{0}}^{p}\displaystyle \frac{\mu (p^{\prime} )}{{\rm{\Gamma }}(p^{\prime} )}{dp}^{\prime} ,\end{eqnarray} \tag{ 12 }$$

where μ(p) is the chemical potential defined as

$$\begin{eqnarray}&&\mu (p)=\exp \left[-{\int }_{{p}_{0}}^{p}\displaystyle \frac{{dp}^{\prime} }{p^{\prime} {\rm{\Gamma }}(p^{\prime} )}\right].\end{eqnarray} \tag{ 13 }$$

The energy density ₀ at the minimum pressure p₀ is the constant of integration. We match the low-density EOS at this pressure generally occurring at a number density slightly below the nuclear saturation density such that there is no unexpected phase transition at the junction. One can determine an empirical fit to any given EOS with arbitrary accuracy by increasing the number of basis functions used. There are four spectral coefficients corresponding to four basis functions as in Lindblom (2010) for fitting the set of RMF EOSs we are using. We plot the distribution of the spectral parameters γ₀, γ₁, γ₂, γ₃, and p₀ in Figure 2. Apart from the coefficients, the parameters ₀ and x_max have fixed values. We optimize the choice of spectral coefficients, γ_k , by minimizing the differences between _fit(x_i , γ_k ) and the energy density _i = (x_i ) values from the tabular version of the EOSs for a specific set of x values (or p values). The following residual quantifies the faithfulness of the fitting

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}({\gamma }_{k})=\displaystyle \sum _{i=1}^{N}\displaystyle \frac{1}{N}{\left\{\mathrm{log}\left[\displaystyle \frac{{\epsilon }_{\mathrm{fit}}({x}_{i},{\gamma }_{k})}{{\epsilon }_{i}}\right]\right\}}^{2},\end{eqnarray} \tag{ 14 }$$

where we sum over the chosen pressure values between ${p}_{0}\leqslant {p}_{i}\leqslant {p}_{\max };$ p_max is the pressure at the center of the nonrotating maximum-mass NS for the given EOS. The spectral version of the EOS serves as input for the Tolman–Oppenheimer–Volkoff equations, which give us the equilibrium configuration of the star.

2.3. Metric with Perturbations

To determine the frequencies of the f-modes, we solve the Einstein equations G_{μ ν} = 8π GT_{μ ν} . We assume here that the GWs are perturbations of the static background spacetime metric for the nonrotating NS. The perturbed metric has the form

$$\begin{eqnarray}&&{g}_{\mu \nu }={g}_{\mu \nu }^{0}+{h}_{\mu \nu }.\end{eqnarray} \tag{ 15 }$$

The time dependence of the perturbed metric components can be expressed by the factor e^{i ω t} for a wave mode. Here ω is complex as the waves will decay due to the imposed open boundary conditions, which will be discussed shortly. The real part of ω is the oscillation frequency, while the imaginary part gives the inverse of the GW damping time (if it is positive) of the wave mode. These oscillations are studied considering the linearized Einstein equations, coupled with the equations of hydrodynamics, for our set of EOS, with suitably posed boundary conditions. A small perturbation h_{μ ν} is assumed on a static spherically symmetric background metric given by

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & -{e}^{\nu (r)}\left[1+{r}^{l}{H}_{0}(r){e}^{i\omega t}{Y}_{{lm}}(\phi ,\theta )\right]{c}^{2}{{dt}}^{2}\\ & & +\,{e}^{\lambda (r)}\left[1-{r}^{l}{H}_{0}(r){e}^{i\omega t}{Y}_{{lm}}(\phi ,\theta )\right]{{dr}}^{2}\\ & & +\,\left[1-{r}^{l}K(r){e}^{i\omega t}{Y}_{{lm}}(\phi ,\theta )\right]{r}^{2}d{\omega }^{2}\\ & & -\,2i\omega {r}^{l+1}{H}_{1}(r){e}^{i\omega t}{Y}_{{lm}}(\phi ,\theta ){dtdr}{.}\end{array}\end{eqnarray} \tag{ 16 }$$

Here ${e}^{\lambda (r)}=\tfrac{1}{1-\tfrac{2{Gm}(r)}{{c}^{2}r}}$ and ${e}^{\nu (r)}\,=\exp \left(\displaystyle \frac{2G}{{c}^{2}}{\int }_{0}^{r}\left\{\displaystyle \frac{m(r^{\prime} )+\tfrac{4\pi p(r^{\prime} )r{{\prime} }^{3}}{{c}^{2}}}{r^{\prime} \left[r^{\prime} -\tfrac{2m(r^{\prime} )G}{{c}^{2}}\right]}\right\}{dr}^{\prime} \right){e}^{{\nu }_{0}};$ $m(r^{\prime} )$ with $p(r^{\prime} )$ being the enclosed mass and pressure of the star at $r^{\prime}$. Here H₀, H₁, and K are radial perturbations of the metric and the angular part is contained in the spherical harmonics ${Y}_{m}^{l}$ with l denoting the orbital angular momentum number and m is the azimuthal number.

2.4. Lindblom–Detweiler Equations

The perturbations of the energy–momentum tensor of the fluid also need to be taken into account in the Einstein equation. The components of the Lagrangian displacement vector ξ^a (r, θ, ϕ) describe the perturbations of the fluid inside the star

$$\begin{eqnarray}&&{\xi }^{r}={r}^{l-1}{e}^{-\lambda /2}{{WY}}_{m}^{l}{e}^{i\omega t},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\xi }^{\theta }=-{r}^{l-2}V{\partial }_{\theta }{Y}_{m}^{l}{e}^{i\omega t},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\xi }^{\phi }=-\displaystyle \frac{{r}^{l-2}}{{\sin }^{2}\theta }V{\partial }_{\theta }{Y}_{m}^{l}{e}^{i\omega t},\end{eqnarray} \tag{ 19 }$$

where W and V are functions with respect to r that describe fluid perturbations. Fluid perturbations exist only inside the star. The introduction of a new variable X leads to simplification of the perturbation equations for the NS interior. X is related to the rest of the perturbation functions through the following equations

$$\begin{eqnarray}\begin{array}{rcl}{H}_{0} & = & \left\{8\pi {r}^{2}{e}^{-\nu /2}X-\left[(n+1)Q\right]\right.\\ & & -\,{\omega }^{2}{r}^{2}{e}^{-(\nu +\lambda )}{H}_{1}+\left[n-{\omega }^{2}{r}^{2}{e}^{-\nu }\right.\\ & & -\,\left.\left.{e}^{\lambda }Q(Q-{e}^{-\lambda })\right]K\right\}{\left(2b+n+Q\right)}^{-1},\end{array}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&V=\left[\displaystyle \frac{X}{\epsilon +p}-\displaystyle \frac{Q}{{r}^{2}}{e}^{(\nu +\lambda )/2}W-{e}^{\nu /2}\displaystyle \frac{{H}_{0}}{2}\right]\displaystyle \frac{{e}^{\nu /2}}{{\omega }^{2}},\end{eqnarray} \tag{ 21 }$$

where n = (l − 1)(l + 2)/2, b = Gm/rc², Q = b + 4π Gr² p/c⁴, and is the local energy density. The GW equations can be written as a set of four coupled linear differential equations for the four perturbation functions H₁, K, W, and X, which do not diverge inside the star for any given value of ω (Lindblom & Detweiler 1983)

$$\begin{eqnarray}\begin{array}{rcl}r\displaystyle \frac{{{dH}}_{1}}{{dr}} & = & -\left[l+1+2{{be}}^{\lambda }4\pi {r}^{2}{e}^{\lambda }(p-\epsilon )\right]{H}_{1}\\ & & +\ {e}^{\lambda }\left[{H}_{0}+K-16\pi (\epsilon +p)V\right],\end{array}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}\begin{array}{rcl}r\displaystyle \frac{{dK}}{{dr}} & = & {H}_{0}+(n+1){H}_{1}+\left[{e}^{\lambda }Q-l-1\right]K\\ & & -\ 8\pi (\epsilon +p){e}^{(\lambda /2)}W,\end{array}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\begin{array}{rcl}r\displaystyle \frac{{dW}}{{dr}} & = & -(l+1)\left[W+{{le}}^{\lambda /2}V\right]\\ & & +\ {r}^{2}{e}^{\lambda /2}\left[\displaystyle \frac{{e}^{-\nu /2}X}{(\epsilon +p){c}_{\mathrm{eq}}^{2}}+\displaystyle \frac{{H}_{0}}{2}+K\right],\\ r\displaystyle \frac{{dX}}{{dr}} & = & -{lX}+\displaystyle \frac{(\epsilon +p){e}^{(\nu /2)}}{2}\left\{\Space{0ex}{3.28ex}{0ex}(3{e}^{\lambda }Q-1)K\right.\\ & & -\ \displaystyle \frac{4(n+1){e}^{\lambda }Q}{{r}^{2}}V+(1-{e}^{\lambda }Q){H}_{0}\\ & & +\ ({r}^{2}{\omega }^{2}{e}^{-\nu }+n+1){H}_{1}\\ & & -\ \left[\Space{0ex}{3.25ex}{0ex}8\pi (\epsilon +p){e}^{\lambda /2}+2{\omega }^{2}{e}^{\lambda /2-\nu }{r}^{2}\right.\\ & & \left.\left.\times \ \displaystyle \frac{d}{{dr}}\left(\displaystyle \frac{{e}^{-\lambda /2}}{{r}^{2}}\displaystyle \frac{d\nu }{{dr}}\right)\right]W\right\},\end{array}\end{eqnarray} \tag{ 24 }$$

where ${c}_{\mathrm{eq}}^{2}={dp}/d\epsilon$ is the equilibrium speed of sound of NS matter undergoing oscillations.

2.5. Boundary Conditions for the Perturbation Equations

We discuss the boundary conditions imposed here to solve the Einstein equations inside the star. For that, one needs to integrate the perturbation equations from the center of the star to a certain radius inside the star. The boundary condition for the perturbation functions at the center of the star r = 0 are

$$\begin{eqnarray}\begin{array}{rcl}W(0) & = & 1,\\ X(0) & = & ({\epsilon }_{0}+{p}_{0}){e}^{{\nu }_{0}/2}\\ & & \times \,\left\{\left[\displaystyle \frac{4\pi }{3}({\epsilon }_{0}+3{p}_{0})-\displaystyle \frac{{\omega }^{2}}{l}{e}^{-{\nu }_{0}}\right]W(0)+\displaystyle \frac{K(0)}{2}\right\},\\ {H}_{1}(0) & = & \displaystyle \frac{{lK}(0)+8\pi ({\epsilon }_{0}+{p}_{0})W(0)}{n+1}.\end{array}\end{eqnarray} \tag{ 25 }$$

The condition Δp = 0 at r = R, the star's surface, requires the perturbation function X to be equal to 0 at the surface of the star. We assign some small arbitrary values to the functions H₁, K, and W at the star's surface and integrate backward to reach the point at which the integration from the center of the star ends and join the two forward and backward solutions. To find the QNM frequency for a given star, we solve the Zerilli equation given as (Fackerell 1971; Lindblom & Detweiler 1983)

$$\begin{eqnarray}&&{d}^{2}Z/{{dr}}^{\star 2}=\left({V}_{Z}(r)-{\omega }^{2}\right)Z({r}^{\star }),\end{eqnarray} \tag{ 26 }$$

where ${r}^{\star }=r+2{Mlog}\left(r/2M-1\right)$ and V_Z (r) is the effective potential outside the star(r ≥ R). The value of the Zerilli function Z(r^⋆), as given in Equation (20) of Kunjipurayil et al. (2022), depends only on H₁ and K as the fluid perturbations W, V, and X become nonexistent outside the star. We calculate the value of Z(r*) at the surface of the star from the values of H₁ and K at the surface. Outside the star Equation (26) is numerically integrated starting from the surface as done in Kunjipurayil et al. (2022) and Zhao & Lattimer (2022) until a distance equivalent to r = 25ω⁻¹. We match this value of Z at r = 25ω⁻¹ with that obtained from an asymptotic expansion of Z valid far away from the surface of the NS. We impose the outgoing boundary condition, which requires the amplitude of incoming GWs far away from the surface of the star to be zero (Lindblom & Detweiler 1983) and solve for ω as done in Kunjipurayil et al. (2022) and Zhao & Lattimer (2022).

Our primary aim is to comprehensively examine the f-mode oscillations, one of the nonradial oscillatory modes of NSs over a broad range of realistic EOSs. To achieve this, we employ a linearized and fully general relativistic approach. For this work, we prepare a set of realistic nuclear EOSs that only have nucleonic composition, obtained with the RMF approach. We utilize around 9000 nucleonic EOSs that were acquired through Bayesian inference with minimal restrictions imposed on them, as previously mentioned in Section 2.1. The ranges of nuclear saturation properties, including the binding energy e₀, nuclear incompressibility K₀, skewness parameter Q₀, symmetry energy J_sym,0, its slope L_sym,0, and its curvature K_sym,0 at the saturation number density n₀, are presented in Table 1.

Table 1. Nuclear Matter Parameters Calculated at Saturation Density for the 9000 Equations of State Used in This Work, along with the Median and 90% Confidence Interval

Model	e0	n0	K0	Q0	Jsym,0	Lsym,0	Ksym,0	m⋆
	(MeV)	(fm−3)	(MeV)	(MeV)	(MeV)	(MeV)	(MeV)
NL	$-{16.10}_{-0.33}^{+0.33}$	${0.152}_{-0.005}^{+0.008}$	${255}_{-40}^{+42}$	$-{451}_{-71}^{+154}$	${32.07}_{-2.77}^{+2.52}$	${32}_{-21}^{+13}$	$-{96}_{-54}^{+96}$	${0.74}_{-0.10}^{+0.03}$

Download table as:  ASCIITypeset image

Our tabulated set of EOSs is given a functional representation by the spectral decomposition method (see Section 2.2). This is a reliable approach which is often used in the context of NS interior structure calculations. Using a minimal number of spectral coefficients, we achieved a fit of all our ∼9000 EOS with a precision of less than 0.5% using the Levenberg–Marquardt method for optimization of the residual errors in Equation (14) as done in Lindblom (2010, 2018).

In Figure 1, we plot the entire EOS set employed in this work and compare with constraints derived from pQCD at baryon densities of 5n_s and 8n_s , with n_s = 0.16 fm⁻³. The plot is divided into two panels, representing two different renormalization scales, (X = 1) and (X = 4), for the pQCD constraints. The left panel focuses on the renormalization scale (X = 1) of the pQCD constraints, with the green patch indicating the energy densities and pressure points for each EOS at the respective baryon densities that are in agreement with the pQCD constraints. All the EOSs in our set meet the requirements set by the pQCD results up to 8n_s and a renormalization scale of four. However, the central density for maximum-mass NS for all the EOS in our set is below 7n_s .

Zoom In Zoom Out Reset image size

In Figure 2, we present a corner plot showing the distributions of the spectral parameters for our set of EOSs. This visual representation is essential to understand the statistical properties and correlations among the parameters. The diagonal plots contain representations of the marginalized 1D distributions of the individual parameters, while the off-diagonal plots illustrate the 2D distributions for pairs of parameters, giving insight into their mutual relationships. The median value of each spectral parameter along with error are shown as a label on the corresponding diagonal plot and the dashed vertical lines mark the limits of the 68% confidence interval (CI) of the parameter. The color gradient, which changes from dark to light blue in the 2D graphs, symbolizes the CIs 1σ, 2σ, and 3σ. The elliptical shapes observed in the 2D plots indicate strong correlations among the parameters. If the shape approaches a circle, it signifies a weak correlation between the given parameter pair.

Zoom In Zoom Out Reset image size

An interesting observation from the figure is that most of the γ_i parameters show strong correlations with their subsequent γ_i+1 parameters. This trend is expected; to ensure a smooth EOS, an increase in a particular γ must be related to the previous one. This aspect of parameter correlation is essential for maintaining a coherent and smooth representation of the EOS, thus confirming the effectiveness of the spectral decomposition format in capturing the essential characteristics of the dense matter EOS.

In Figure 3, we plot the mass–radius (MR) relation of the NS (left) and the mass–tidal deformability (right) of all EOSs considered using the spectral representation of the EOSs (see Section 2.1). The gray zones in the left panel represent the 90% (dark) and 50% (light) CIs for the binary components of the GW170817 event (Abbott et al. 2019). The predictions of NS radius measurements from NICER X-ray data for the PSR J0030+0451 and PSR J0740+6620 pulsars are also included, with the 1σ (68%) confidence zone of the 2D posterior distribution in the MR domain of the millisecond pulsar PSR J0030+0451 shown in brown and yellow from Riley et al. (2019) and Miller et al. (2019), respectively, and PSR J0740+6620 is shown in magenta (Miller et al. 2021; Riley et al. 2021). Finally, the blue bar in the left panel represent the 90% CI radius of a 2.08 M_⊙ NS as reported in Miller et al. (2021), who combined NICER data and X-ray observations. The blue bar in the right panel shows the bounds on the tidal deformability of a 1.36 M_⊙ NS as given in Abbott et al. (2018). All MR curves have a maximum mass greater than 2 M_⊙, since this was imposed as a constraint during Bayesian sampling. Despite the dearth of constraints on radius and tidal deformability during sampling, the predictions of NS mass and radius are well within the observational bounds over the entire EOS set with the minimal requirements that were set. The obtained radius of a 1.4 M_⊙ NS in our set is 12.03–13.20 km, and the deformability of the tidal channel for a 1.36 M_⊙ star is 431–659 within the 90% CI. The maximum mass of the NS ranges from 2 to 2.25 M_⊙, which is consistent with the results reported in Rezzolla et al. (2018). It is worth noting that we applied pQCD constraints to our EOS set, which limited the maximum mass of the NS to a certain range and eliminated some of the stiffer EOSs.

Zoom In Zoom Out Reset image size

We display the f-mode frequencies in Figure 4 (left) versus NS mass. They increase from 1.6 to 2.6 kHz for variations of the NS mass from 1.2 to 2.25 M_⊙. It is worth noting that the f-mode frequency within the Cowling approximation for NSs with a mass greater than 2 M_⊙ is in the range of 2.1–2.7 kHz and 2.3–2.65 kHz for the two types of EOSs used in Kumar et al. (2023), respectively. It is a widely known fact that the calculation of the f-mode frequencies using the Cowling approximation overestimates them by up to 10% to 30% for NSs with masses ranging from 1.0 to 2.5 M_⊙, when compared to the frequency obtained from the linearized general relativistic formalism (Benhar et al. 2004; Doneva et al. 2013; Pradhan et al. 2022). In the right panel, we plot the distribution of f-mode frequencies for NS masses ranging from 1.4 to 2.0 M_⊙. The distribution for a given mass is accompanied by its median value and the 90% CI bounds. This variation in the f-mode frequencies across various NS masses is indicative of the current understanding of the EOSs for NSs. The measurements of f-mode frequencies (with precision beyond the current uncertainty of approximately 0.2 kHz) by next-generation detectors could potentially lead to a more refined understanding of the EOS domain for NS interiors.

Zoom In Zoom Out Reset image size

Next we calculate the GW damping times corresponding to the f-mode frequencies for our EOS set. The inverse of the imaginary part of ω in the metric perturbations gives us the damping time as discussed in Section 3. In Figure 5 we plot the variation of damping times with stellar mass. We obtain a monotonically decreasing trend for the damping times with an increase in NS mass until the maximum mass for all EOSs. The damping time of a 1.4 M_⊙(2.0 M_⊙) NS lies between 0.20 and 0.25 s (0.13–0.15 s). The range of damping times in this work, for the two NS masses mentioned, is in good agreement with the range of damping times reported in Pradhan et al. (2022) for their set of RMF EOSs.

Zoom In Zoom Out Reset image size

In Figure 6 we explore the linear dependence between various spectral parameters, nuclear saturation properties, the pressure of NS matter at 2 and 3 times the fixed baryon density n_s , and other NS properties. We also examine how the frequency of the f-mode oscillation for different masses relates to the individual components of our nuclear EOS set. We calculate the Pearson correlation coefficients between these quantities using the entire ∼9000 EOS set. The heat map's color spectrum, ranging from green to red, indicates the strength of the correlation, from negative to positive, respectively. A deeper red (green) hue indicates a stronger positive (negative) correlation, with a magnitude of coefficients above 0.8 indicating a particularly strong linear relationship. The correlation between the spectral parameters γ₁, γ₂, and γ₃ is exceptionally high, clearly indicating a linear relationship between them. This is also true for the pair m^⋆ and the γ_i , suggesting that the effective mass of the nucleon has a strong influence on these spectral parameters. Another finding is the strong correlation, greater than 0.9, between P(2n_s ) and P(3n_s ), which is the pressure β at 2 and 3 times n_s and various properties of the NS such as radius, deformability of the tidal, frequency of the f-mode, and damping time for NS masses of 1.4 M_⊙ and 2.0 M_⊙. This suggests a considerable degree of linear dependence between the properties of NSs of different masses and the overall EOS, rather than individual nuclear saturation properties. It is worth emphasizing the strong correlation previously reported in Pradhan et al. (2022) between the fundamental f-mode frequency of different NS masses and the effective nucleon mass at saturation m^⋆ appears to be less pronounced in our study using the chosen set of EOSs. This difference can be mainly attributed to certain choices regarding the set of EOSs. The set employed in our research was generated through Bayesian inference subject to constraints. Specifically, constraints from χEFT for our low-density pure neutron matter and pQCD at 8 times n_s were included. These constraints were not considered by the other authors. An additional notable distinction is the inclusion of an interaction ω⁴ in the Lagrangian that we studied. It is important to note that this ω⁴ term tends to soften the EOS, which is a characteristic favored by pQCD analysis (Malik et al. 2023). The heat map suggests that the maximum radius of the NS is closely related to its R_1.4, R_2.0, and f_2.0 values. Furthermore, the f-mode frequencies of NSs with 1.4 M_⊙ and 2.0 M_⊙ are strongly correlated with their tidal deformabilities Λ_1.4 and Λ_2.0, respectively, providing insight into how the internal structure of the NS responds to tidal forces.

Asteroseismology, the study of stellar oscillations, allows us to infer global NS properties such as mass, radius, and average stellar density from f-mode frequencies and the GW damping time through certain empirical relations. We describe the results for some selected relations in the following paragraphs.

Zoom In Zoom Out Reset image size

In Figure 7 we plot the data with a straight line fit using

$$\begin{eqnarray}&&{\nu }_{f}=a\sqrt{M/{R}^{3}}+b.\end{eqnarray} \tag{ 27 }$$

This empirical relation is well known since the publication of Andersson & Kokkotas (1998). We find the slope a to be 35.94 ± 0.0113 kHz km and the intercept b to be 0.6260 ± 0.0004 kHz. We can compare these values with the values reported in Pradhan et al. (2022; see Table 2). They calculated the f-mode frequencies with a different set of RMF EOSs and found the slope equal to 36.20 kHz km and the intercept equal to 0.535 kHz. Kumar et al. (2023) report the same coefficients with values a = 22.27 and 26.76 and b = 1.520 and 1.348 corresponding to determination of f-modes with the Cowling approximation for a set of hybrid and a set of purely nucleonic EOSs with density-dependent couplings. The values for a and b thus appears to be dependent on the set of EOSs considered, which makes the empirical Equation (27) somewhat EOS dependent. The relative percentage error for the aforementioned fit in this work is 2.84%.

Zoom In Zoom Out Reset image size

Table 2. Coefficients of the Fitting Equation for the Empirical Relation in Equation (27)

		[Full General Relativity]		[Full General Relativity]	[Cowling]
Empirical Relation	Fitting Function	Coefficients of the Current Work	Relative Error	Coefficients of Pradhan et al. (2022)	Coefficients of Kumar et al. (2023)
$f\mbox{--}\sqrt{M/{R}^{3}}$	ax + b	a = 35.94 ± 0.0113	2.84%	a = 36.20	a = 22.27, 26.76
		b = 0.6260 ± 0.0004		b = 0.535	b = 1.520, 1.348

Download table as:  ASCIITypeset image

We display another empirical relation

$$\begin{eqnarray}&&\displaystyle \frac{{R}^{4}}{{M}^{3}{\tau }_{f}}=a\displaystyle \frac{M}{R}+b,\end{eqnarray} \tag{ 28 }$$

in the right panel of Figure 7. τ_f is the GW damping time, equal to the inverse of the imaginary part of ω. Andersson & Kokkotas (1998) mentioned the existence of this scaling between a dimensionless damping time and compactness for the first time. Both the coefficients a and b are dimensionless and we find the following values (see Table 3) of a = −0.3051 ± 0.0001 and b = 0.0939 ± 3.0 × 10⁻⁵ with a relative error of 0.29%. Pradhan et al. (2022) reported the values of a = −0.245 and b = 0.080. These empirical relations have poor validity beyond M/R > 0.25 (Tsui & Leung 2005).

Table 3. Coefficients of the Fitting Equation for the Empirical Relation in Equation (28)

		[Full General Relativity]		[Full General Relativity]
Empirical Relation	Fitting Function	Coefficients of the Current Work	Relative Error	Coefficients of Pradhan et al. (2022)
$\displaystyle \frac{{R}^{4}}{{M}^{3}{\tau }_{f}}\mbox{--}M/R$	ax + b	a = −0.3051 ± 0.0001	0.29%	a = −0.245
		b = 0.0939 ± 3.0 × 10−5		b = 0.080

Download table as:  ASCIITypeset image

Kumar et al. (2023) have reported a new scaling between the f-mode frequency and the radius for an NS of given stellar mass

$$\begin{eqnarray}&&{f}_{x}={{aR}}_{x}+b,\end{eqnarray} \tag{ 29 }$$

where x is the mass of the NS in M_⊙ units. We plot the relation with our results of f-mode frequencies for 1.4 M_⊙ and 2.0 M_⊙ in Figure 8.

Zoom In Zoom Out Reset image size

We find values (see Table 4) of a = −0.1933 ± 0.0008 kHz km⁻¹ and b = 4.195 ± 0.010 kHz for 1.4 M_⊙ and a =−0.2455 ± 0.0003 kHz km⁻¹ and b = 5.027 ± 0.003 kHz for 2.0 M_⊙. The maximum relative errors are 1.86% and 1.18%, respectively. The slope reported in Kumar et al. (2023) for masses over the range 1.6–2.4 M_⊙ is around −0.22 kHz km⁻¹ and the intercept is around 5.1 khz with a 1.5% relative error. A general relativistic treatment of the stellar oscillations in this work presents us with two different slopes for the two cases. This difference is due to the overestimation of f-mode frequencies that occurs when one employs Cowling approximation, as done in Kumar et al. (2023). From a specific observed value of the f-mode frequency for 1.4 M_⊙, we can infer a smaller value of the radius from the fit in this work as compared to the fit in Kumar et al. (2023).

Motivated by the study in Lin & Steiner (2023), we also look into the correlation between the radii as well as f-mode frequencies of two extreme-mass NSs, utilizing our nucleonic EOS set. In Figure 9, we show a scatter plot of the radii (f-mode frequencies) of NSs with masses of 1.34 M_⊙ and 2.07 M_⊙ on the left (right). On both the panels, a strong correlation is noticed, with a Pearson's correlation coefficient of approximately 0.9 for our EOS set. In the left panel, the red dashed line indicates the most probable region of the joint posterior distribution of the NICER radii measurements of PSR J0740+6620 and PSR 0030+0451 calculated in Lin & Steiner (2023). The authors have claimed a 48% probability (with a 5% false alarm rate) of identifying a strong and sharp phase transition by applying their method to NICER observations. We note how our entire nucleonic set overlaps with the most probable region within the red dashed line. In the same panel, the area marked with the blue dashed line represents the constraints after we impose a strong correlation of 0.9 among the radius data in the NICER observations of the pulsars PSR J0740+6620 and PSR 0030+0451. NICER's 1σ posterior distributions are used to obtain this constraint. This method is based on data reported in Riley et al. (2019, 2021). To be more specific, we produce 2000 random samples of R_1.34 and R_2.07 within the 1σ range of the NICER observed values marginalized over NS masses and implemented a correlation of approximately 0.9 between the two radii. The plot shows a noticeable difference in the slope of the data points corresponding to a nucleonic EOS and those derived from NICER observations.

Zoom In Zoom Out Reset image size

In the right panel, we translate the NICER data from Riley et al. (2019, 2021) as well as those from Lin & Steiner (2023) in the radius domain using Equation (29) to obtain the corresponding constraints on the f-mode frequency domain. Note that the universal relation presented in Table 4 applies to NS masses of 1.4 M_⊙ and 2.0 M_⊙. For Figure 9 we calculate the relations for NS masses of 1.34 M_⊙ and 2.07 M_⊙. The slope and the intercept values are very close to those presented in Table 4: −0.1794 (−0.2360) kHz km⁻¹ and 3.9891 (4.9283) kHz for an NS mass of 1.34 (2.07) M_⊙. The nucleonic EOS set exhibits a similar overlap between the radius and frequency domains with the NICER-inferred constraint from the data in Riley et al. (2019, 2021). However, the partial overlap between the nucleonic EOS set and this NICER bound suggests that only nucleonic degrees of freedom for larger NS masses may not be sufficient enough to explain the observations.

In Figure 10, we present a similar analysis as in Figure 9 for NS masses 1.44 M_⊙ and 2.07 M_⊙ but determine the NICER constraints from the data published by Miller et al. (2019, 2021). The slope and the intercept values for a 1.44 M_⊙ are −0.1869 kHz km⁻¹ and 4.1215 kHz, respectively. We observe a clear distinction between the NICER bound and the nucleonic EOS set in both the radius and frequency domains. This finding further strengthens the conclusion drawn from the previous figure.

Zoom In Zoom Out Reset image size

To investigate further, we utilize 12 EOS models with a maximum NS mass above 2.0 M_⊙ by simulating a hadron–quark phase transition inside the NS core from the CompOSE database, such as VQCD (Jokela et al. 2021), QHC21-AT (Kojo et al. 2022), QHC19-C (Baym et al. 2019), QHC21-BT, QHC21-D, QHC19-D, QHC21-CT, QHC21-B, QHC21-C, QHC18 (Baym et al. 2018), QHC19-B, and QHC21-DT. The red stars in both Figures 9 and 10 with the label "Hyb-EOS" in both panels are the calculated values for these EOSs in the corresponding domains. It is clear that all the red stars fit comfortably within our NICER-derived constraints, as determined by both Riley et al. (2019, 2021) and Miller et al. (2019, 2021), over the radius and the f-mode frequency domains. The possibility of a hadron–quark phase transition or other exotic degrees of freedom in the NS EOS appears to be promising. Please be aware that the blue dashed contours derived from the NICER data for the pulsars PSR J0030+0451 and PSR J0740+6620, whether in the radius domain or in the f-mode frequency domain, are calculated using a correlation of 0.9 between the radii of NSs with these two masses, and are entirely dependent on the nucleonic degrees of freedom. Examining the differences in these correlations with other exotic degrees of freedom is a subject for future research and is not within the purview of this study.

The universal relation in Equation (29) acts as a tool to compare data from two different channels: X-ray channel and GW channel. Observations of thermal emission from hot spots on pulsars provide an effective way of measuring the NS radius. Parameter estimation of parameterized models of NS EOSs using such radius measurements (Raaijmakers et al. 2021) allows us to narrow down viable EOSs by putting constraints on the relevant parameters. Radius measurements from X-ray data are essentially based on X-ray pulse profile modeling (Bogdanov et al. 2021). This involves making assumptions about the stellar atmosphere and the EOS of matter inside the star itself. In short, such measurements are model dependent. f-mode observations offer a more direct path to probe the NS EOS by measuring the bulk properties of the star through its response to gravity. Thus, X-ray observations, while valuable, require additional assumptions about the star's surface and internal structure, making them less direct for studying the deep core where the EOS is truly uncertain as yet. These observations work best together (Yunes et al. 2022), with X-rays providing complementary data points like mass and radius to refine the interpretation of the NS EOS extracted from GWs (Psaltis et al. 2014).

We have carefully calculated the QNM frequency, specifically focusing on the f-mode oscillations of NSs. This analysis was conducted using a linearized approach to full general relativity, which allowed for a detailed examination of these oscillations. In our study, we utilized a wide range of nuclear EOSs (9000), derived through Bayesian inference. These EOSs were subjected to a few constraints. First, they conform to the minimal constraints on nuclear saturation properties. Second, they incorporate constraints from χEFT relevant to low-density pure neutron matter. Additionally, constraints derived from pQCD at densities pertinent to the cores of NSs were also applied. Finally, our selection of EOSs ensured that the NS maximum mass was above 2 M_⊙, aligning with current astrophysical observations and theoretical predictions.

We then determine the spectral representation of our EOS set, which is one of the most widely accepted and optimal way to adapt any realistic EOS table into a functional form. This helps to reduce the computation time when calculating NS properties and is convenient for other astrophysical calculations. The MR and mass–tidal deformability relationships of NSs were calculated, and a detailed examination of the correlation among spectral parameters, various nuclear saturation properties, and NS properties was conducted. Empirical relations for f-mode frequencies and GW damping times were also explored, providing valuable insights into the composition of these compact objects.

One of the key findings is that the f-mode frequencies do not show a strong relationship with any single nuclear saturation property of the EOS when the 9000 EOSs are used. In fact, we find that some of the strong correlations reported in recent publications have become weak for our EOS set. This observation suggests that a more nuanced approach, possibly involving a multiparameter correlation analysis using machine learning techniques, might be necessary to unravel these complex relationships.

Interestingly, for NS with masses ∼ 1.4 M_⊙ and ∼2.0 M_⊙, we notice a strong correlation between their radii as well as f-mode frequencies. Using this correlation alongside NICER observations of PSR J0740+6620 and PSR 0030+0451, we find partial and minimal overlap with the observational data from Riley et al. (2019, 2021) and Miller et al. (2019, 2021), respectively, with our nucleonic EOS data set. This motivates us to further investigate these results for the presence of nonnucleonic degrees of freedom in the EOS. Hence, we employ a selection of 12 hybrid EOS from the CompOSE database to calculate radius and the f-mode frequency data, which appear to align more closely with the abovementioned NICER constraints particularly in regions where the nucleonic EOS shows no overlap. Note, that the NICER constraints were derived on the basis of a correlation of 0.9 among the radii of 1.4 and 2.07 M_⊙ NSs calculated with our nucleonic set. Thus our findings emphasize the need to systematically check for a signature of a hadron–quark phase transition in the NS interior.

A recent study by Lin & Steiner (2023) also demonstrates the correlation between the radii of NSs for the observed masses of the same pulsars, as detected by NICER. Their analysis suggests the existence of a sharp and strong phase transition with 48% probability. Although our methodology differs from that of Lin & Steiner (2023), our conclusions are similarly aligned. Our study provides a clearer distinction between nucleonic and nonnucleonic EOSs across two domains over NICER data from two complimentary analyses. By translating NICER radius measurements into the f-mode frequency domain using a universal relationship, we highlight an alternative avenue through GW astronomy to compare various EOS models and check for the presence of exotic phases through the GW channel.

For the first time, our analysis shows that the NICER constraints on radius prefer hybrid EOSs over purely nucleonic ones, which also holds true in the f-mode frequency domain. It is also expected that NSs with higher mass, for example 2.0 M_⊙, will have a different composition relative to 1.4 M_⊙ NSs. The conditions inside a higher-mass NS with a higher central density than that inside a 1.4 M_⊙ NS are more favorable for a deconfined quark phase. We expect future precise measurements of the f-mode frequency, particularly those for NSs with masses near the two extremes of viable NS mass values, will provide decisive evidence regarding the internal composition of NSs. Such measurements will be able to differentiate between purely nucleonic matter and other exotic phases, shedding light on the state of matter at extreme densities. We would do a thorough and organized study along similar lines taking into account a wide range of hybrid EOSs generated using multiple phenomenological models and report our results in a future communication.

T.M. would like to acknowledge the support from national funds from FCT (Fundaço para a Ciência e a Tecnologia, IEP, Portugal) under Projects No. UIDP/-04564/-2020, No. UIDB/-04564/-2020, with DOI identifiers 10.54499/UIDB/04564/2020 and 10.54499/UIDP/04564/2020, and 2022.06460.PTDC with the associated DOI identifier 10.54499/2022.06460.PTDC. T.M. is also grateful for the support of EURO-LABS "EUROpean Laboratories for Accelerator Based Science," which was funded by the Horizon Europe research and innovation program with grant agreement No. 101057511. D.G.R., Sw.B., and Sa.B. would like to acknowledge the financial support by DST-SERB, Govt. of India through the Core Research Grant (CRG/2020/003899) for the project entitled "Investigating the Equation of State of Neutron Stars through Gravitational Wave Emission." Finally, D.G.R. would like to thank colleagues A. Venneti and K. Nobleson for numerous suggestions regarding the work done in this publication.