Fundamental and higher-order excited modes of radial oscillation of neutron stars for various types of cold nucleonic and hyperonic matter

TL;DR

This study computes fundamental and excited radial oscillation modes of neutron stars using various realistic equations of state, revealing correlations that aid in constraining models and advancing neutron star asteroseismology.

Contribution

It provides detailed calculations of oscillation modes for different EoS models, enhancing understanding of neutron star internal structure and dynamics.

Findings

Correlation between oscillation frequency and adiabatic index

Constraints on realistic EoS models from oscillation data

Implications for neutron star asteroseismology

Abstract

This research paper complements our earlier qualitative study of the effect of viscosity and thermal conductivity on the radial oscillation and relaxation of non-rotating neutron stars. The fundamental and first two lowest-frequency excited modes of radial oscillation have been computed in the high nuclear density regime for a set of seven realistic equations of state (EoS) as functions of central energy density. Various types of zero-temperature EoS of cold nucleonic and hybrid nucleon-hyperon-quark matter models are used in the inner core to determine the internal structure in and around the hydrostatic equilibrium states and investigate the influence of each EoS on the dynamical behaviour of non-rotating neutron stars. We confirm the principal results of earlier, related studies that suggest an underlying correlation between the frequency spectrum of the fundamental oscillation mode and the variation of the adiabatic index over the high nuclear-density regime. We provide valuable information to impose further constraints on the plausible set of realistic EoS models, in addition to the practical applications for the rapidly evolving field of asteroseismology of compact objects.

Tables

Table 1. Table 1: Nucleonic and hybrid nucleon–hyperon–quark matter models of NS cores. The upper four models correspond to matter with nucleonic degrees of freedom only, whereas the lower three models involves non-nucleonic states of nuclear matter, such as kaon condensates or hyperons. We present some characteristic parameters of the maximal-mass-mass configuration for each EoS (marked in Fig. 6a by the symbols ⋄ ⋄ \diamond ): M max subscript 𝑀 max M_{\text{max}} is the maximal mass in units of M ⊙ subscript 𝑀 direct-product M_{\odot} , R 𝑅 R is the circumferential radius in units of km, r s / R subscript 𝑟 𝑠 𝑅 r_{s}/R is the compactness, respectively.

EoS	Constituents	${𝑴}_{\text{max}}/{𝑴}_{\odot }$	$𝑹$/km	${𝒓}_{𝒔}/𝑹$	Ref.
APR4	Nucleons	2.22	9.94	0.66	[41]
MPA1	Nucleons	2.47	11.25	0.65	[41]
MS1	Nucleons	2.78	13.27	0.62	[42]
SLy4	Nucleons	2.06	9.9	0.61	[43]
SQM1	Nuc., s quarks	1.56	8.54	0.54	[44]
H4	Nuc., $\Sigma ,\Lambda$	2.04	11.81	0.51	[38]
ALF1	u, d, s quarks	1.50	9.02	0.49	[39]

Table 2. Table 2: Possible mechanisms of neutron star cooling by various neutrino-emission processes due to nucleon–nucleon collisions assumed to take part in the core. The modified Urca process has the neutron and the proton branch, each including a direct and an inverse, where N = n 𝑁 𝑛 N=n or p 𝑝 p , respectively. [ 51 ]

Process	Reaction	Heat flux ${\text{[erg }\text{cm}}^{-3}{\text{s}}^{-1}]$	Local luminosity $\left[{\text{erg }\text{s}}^{-1}\right]$
Direct Urca	$n\to p+e+{\overline{\nu }}_e$ $p+e\to n+{\nu }_e$	$Q\sim 3\times {10}^{27}{T}_9^8$	$L_{\nu }\sim {10}^{46}{T}_9^8$
Modified Urca	$n+N\to p+e+N+{\overline{\nu }}_e$ $p+e+N\to n+e+N+{\nu }_e$	$Q\sim {10}^{20-22}{T}_9^8$	$L_{\nu }\sim {10}^{38-40}{T}_9^8$
Bresmmstrahlung	$N+N\to N+N+\nu +\overline{\nu }$	$Q\sim {10}^{18-20}{T}_9^8$	$L_{\nu }\sim {10}^{36-38}{T}_9^8$

Table 3. (a) APR4

${ϵ}_c$	$R$	$M$	${\nu }_0$	${\nu }_1$	${\nu }_2$
$\left({\text{GeV fm}}^{-3}\right)$	(km)	$\left(M_{\odot }\right)$	(kHz)	(kHz)	(kHz)
1.518	9.940	2.219	${0.398}^{\ast }$	6.012	9.142
1.420	10.051	2.217	0.684	6.131	9.261
1.320	10.172	2.210	1.145	6.279	9.409
1.220	10.304	2.195	1.470	6.418	9.542
1.120	10.441	2.171	1.799	6.576	9.698
1.030	10.589	2.131	2.143	6.738	9.868
0.930	10.739	2.068	2.471	6.883	10.005
0.830	10.886	1.974	2.797	6.965	10.196
0.730	11.008	1.832	3.108	6.937	10.254
0.630	11.080	1.625	3.278	6.769	8.984
0.530	11.038	1.336	3.063	6.254	6.645
0.430	10.836	0.977	2.269	4.191	4.156
1.330	10.474	0.600	1.001	0.984	0.998

Table 4. (a) APR4

${ϵ}_c$	$R$	$M$	${\nu }_0$	${\nu }_1$	${\nu }_2$
$\left({\text{GeV fm}}^{-3}\right)$	(km)	$\left(M_{\odot }\right)$	(kHz)	(kHz)	(kHz)
1.518	9.940	2.219	${0.398}^{\ast }$	6.012	9.142
1.420	10.051	2.217	0.684	6.131	9.261
1.320	10.172	2.210	1.145	6.279	9.409
1.220	10.304	2.195	1.470	6.418	9.542
1.120	10.441	2.171	1.799	6.576	9.698
1.030	10.589	2.131	2.143	6.738	9.868
0.930	10.739	2.068	2.471	6.883	10.005
0.830	10.886	1.974	2.797	6.965	10.196
0.730	11.008	1.832	3.108	6.937	10.254
0.630	11.080	1.625	3.278	6.769	8.984
0.530	11.038	1.336	3.063	6.254	6.645
0.430	10.836	0.977	2.269	4.191	4.156
1.330	10.474	0.600	1.001	0.984	0.998

Table 5. (b) H4

${ϵ}_c$	$R$	$M$	${\nu }_0$	${\nu }_1$	${\nu }_2$
$\left({\text{GeV fm}}^{-3}\right)$	(km)	$\left(M_{\odot }\right)$	(kHz)	(kHz)	(kHz)
1.179	11.811	2.039	0.777*	1.800	2.365
1.000	12.15	2.031	0.880	1.821	2.430
0.950	12.255	2.025	0.892	1.818	2.437
0.900	12.375	2.016	0.903	1.814	2.444
0.800	12.488	2.004	0.912	1.812	2.450
0.700	12.625	1.988	0.919	1.809	2.453
0.600	12.900	1.940	0.927	1.783	2.427
0.500	13.206	1.867	0.956	1.820	2.460
0.450	13.495	1.753	1.018	1.927	2.572
0.350	13.608	1.673	1.045	1.960	2.600
0.550	13.809	1.430	0.970	1.733	2.130
0.250	13.704	0.980	0.793	1.257	1.262
0.190	13.563	0.624	0.352	0.374	0.341

Table 6. (c) MPA1

${ϵ}_c$	$R$	$M$	${\nu }_0$	${\nu }_1$	${\nu }_2$
$\left({\text{GeV fm}}^{-3}\right)$	(km)	$\left(M_{\odot }\right)$	(kHz)	(kHz)	(kHz)
1.174	11.250	2.470	0.390*	7.820	11.960
1.100	11.353	2.469	0.985	7.914	12.088
1.000	11.507	2.458	1.608	7.997	12.166
0.950	11.593	2.448	1.868	8.043	12.221
0.850	11.772	2.411	2.320	8.111	12.263
0.800	11.864	2.382	2.526	8.132	12.237
0.750	11.955	2.343	2.725	8.151	12.194
0.700	12.042	2.291	2.907	8.160	12.094
0.600	12.195	2.129	3.188	8.082	11.617
0.500	12.246	1.846	3.290	7.612	10.550
0.400	12.034	1.392	3.250	6.479	7.720
0.350	11.765	1.106	2.907	5.181	5.447
0.280	11.189	0.703	0.501	0.512	0.608

Table 7. (d) MS1

${ϵ}_c$	$R$	$M$	${\nu }_0$	${\nu }_1$	${\nu }_2$
$\left({\text{GeV fm}}^{-3}\right)$	(km)	$\left(M_{\odot }\right)$	(kHz)	(kHz)	(kHz)
0.863	13.265	2.782	${3.270}^{\ast }$	7.873	11.117
0.800	13.410	2.780	3.284	7.741	10.830
0.750	13.535	2.773	3.290	7.609	10.544
0.680	13.710	2.752	3.288	7.385	10.054
0.630	13.863	2.725	3.280	7.198	9.635
0.560	14.070	2.663	3.268	6.906	8.880
0.500	14.252	2.573	3.257	6.625	8.070
0.450	14.394	2.459	3.243	6.280	7.315
0.390	14.523	2.247	3.118	5.702	6.337
0.330	14.555	1.923	2.844	4.996	5.232
0.270	14.392	1.484	2.416	4.175	3.914
0.210	13.965	1.012	1.665	3.237	2.380
0.150	13.184	0.560	0.502	2.182	0.631

Table 8. (e) Sly4

${ϵ}_c$	$R$	$M$	${\nu }_0$	${\nu }_1$	${\nu }_2$
$\left({\text{GeV fm}}^{-3}\right)$	(km)	$\left(M_{\odot }\right)$	(kHz)	(kHz)	(kHz)
2.056	9.238	2.025	1.818*	3.018	3.968
1.900	9.496	2.042	1.814	2.983	4.029
1.800	9.641	2.050	1.805	2.956	4.036
1.600	9.887	2.056	1.773	2.887	3.986
1.450	10.065	2.053	1.738	2.825	3.886
1.300	10.262	2.040	1.691	2.753	3.734
1.150	10.486	2.011	1.635	2.671	3.530
1.000	10.727	1.956	1.568	2.579	3.274
0.850	10.982	1.857	1.490	2.478	2.965
0.700	11.218	1.688	1.415	2.291	2.541
0.550	11.360	1.408	1.358	1.750	1.791
0.400	11.251	0.979	0.916	0.886	0.891
0.280	10.831	0.553	0.296	0.299	0.300

Table 9. (f) SQM1

${ϵ}_c$	$R$	$M$	${\nu }_0$	${\nu }_1$	${\nu }_2$
$\left({\text{GeV fm}}^{-3}\right)$	(km)	$\left(M_{\odot }\right)$	(kHz)	(kHz)	(kHz)
0.823	13.357	2.782	${3.165}^{\ast }$	8.109	11.703
0.800	13.410	2.780	3.207	8.057	11.544
0.770	13.484	2.776	3.251	7.956	11.297
0.740	13.560	2.770	3.278	7.815	10.990
0.710	13.637	2.762	3.290	7.612	10.550
0.680	13.720	2.752	3.286	7.345	9.967
0.660	13.775	2.742	3.278	7.134	9.481
0.630	13.863	2.725	3.262	6.778	8.511
0.600	13.951	2.703	3.240	6.272	7.301
0.570	14.040	2.674	3.035	5.474	5.954
0.550	14.100	2.651	2.750	4.820	4.894
0.520	14.190	2.608	2.005	3.120	2.958
0.490	14.280	2.554	0.512	0.524	0.608

Equations

$$p=p(\rho ),ϵ=ϵ(\rho ),$$

$${\rho }_{\text{nuc}}\simeq 2.3\times 10^{14}{\text{ g/cm}}^3.$$

$$p=K{\rho }^{\Gamma },$$

$$ϵ=\rho +\frac{1}{\Gamma -1}K{\rho }^{\Gamma }=\rho +\frac{1}{\Gamma -1}p,$$

$${\Gamma }_1=\frac{d\log p}{d\log \rho }=\frac{ϵ+p}{p}\frac{dp}{dϵ},$$

$$c_{s}^{2}=\frac{dp}{d\epsilon}\bigg{\rvert}_{s}\leq 1.$$

$${\bar{\Gamma }}_1=\frac{{\int }_0^R{\Gamma }_{1}p(4\pi r^2)dr}{{\int }_0^{R}p(4\pi r^2)dr},$$

$${\bar{\Gamma }}_1=\frac{ϵ+p}{p}{c}_s^2\ge \frac{4}{3}$$

$${\bar{\Gamma }}_1<\frac{4}{3}+K\frac{4M}{R},$$

$$p(\rho )=K_i{\rho }^{{\Gamma }_i}\text{for}{\rho }_i\le \rho <{\rho }_{i+1}(0\le i\le n),$$

$$p=K{\rho }^{\Gamma }+({\Gamma }_{\text{th}}-1)\rho {ϵ}_{\text{th}}$$

$$ϵ={ϵ}_{\text{th}}+{ϵ}_{\text{cold}}.$$

$$p=K\frac{\Gamma -{\Gamma }_{\text{th}}}{\Gamma -1}{\rho }^{\Gamma }+({\Gamma }_{\text{th}}-1)\rho {ϵ}_{\text{th}}.$$

$${\partial }_{\beta }{T}^{\alpha \beta }=-Q_{\nu }{u}^{\alpha },$$

$$T^{\alpha \beta }=T_{\text{pf}}^{\alpha \beta }+T_{\text{visc}}^{\alpha \beta }+T_{\text{heat}}^{\alpha \beta },$$

$$\begin{array}[]{l}T^{\alpha\beta}_{\text{pf}}=(\epsilon+p)u^{\alpha}u^{\beta}+ph^{\alpha\beta},\\[10.0pt] T^{\alpha\beta}_{\text{heat}}=Q^{\alpha}u^{\beta}+Q^{\beta}u^{\alpha}\\[10.0pt] T^{\alpha\beta}_{\text{visc}}=-\zeta\Theta h^{\alpha\beta}-2\eta\sigma^{\alpha\beta}\end{array}$$

$$h^{\alpha \beta }=g^{\alpha \beta }+u^{\alpha }{u}^{\beta }$$

$${\sigma }^{\alpha \beta }=\frac{1}{2}\left(u_{;\mu }^{\alpha }{h}^{\mu \beta }+u_{;\mu }^{\beta }{h}^{\mu \alpha }\right)-\frac{1}{3}\Theta h^{\alpha \beta }$$

$$\Theta =u_{;\alpha }^{\alpha }$$

$$Q^{\alpha }=-\kappa h^{\alpha \beta }\left(T_{,\beta }-T{a}_{\beta }\right)$$

$$d{s}^2=e^{\nu }d{t}^2-e^{\lambda }d{r}^2-r^{2}d{\Omega }^2,$$

$$\nu =\nu (t,r),\lambda =\lambda (t,r)$$

$$e^{\lambda }={\left(1-2m/r\right)}^{-1},$$

$$m(r)=4\pi {\int }_0^r\rho (r^{\prime })r^{\prime 2}d{r}^{\prime }.$$

$$e^{\lambda }=(1-r_s/r{)}^{-1}.$$

$$r_s=2GM/c^2,$$

$$z(r)={\int }_r^R{e}^{\lambda (r^{\prime })}d{r}^{\prime },$$

$$\overline{ϵ}=ϵ+(T_{\text{visc}}+T_{\text{heat}}{)}_0^0,\bar{p}=p-(T_{\text{visc}}+T_{\text{heat}}{)}_1^1$$

$$\frac{dm}{dr}=4\pi r^2ϵ$$

$$\frac{dp}{dr}=-\frac{(ϵ+p)(m+4\pi r^{3}p)}{r(r-2m)}+\frac{2\ddot{\lambda }+\dot{\lambda }(\dot{\lambda }-\dot{\nu })}{16\pi }{e}^{-\nu }$$

$$\frac{d\nu }{dr}=-\frac{2}{ϵ+p}\frac{dp}{dr},$$

Full text

Relativistic stars: structure, stability, and oscillations – Neutron stars – Equations of state of neutron-star matter

††preprint: APS/123-QED

I Introduction

Neutron stars (NSs) provide us with unique insights into the physics of the extremely dense and cold nuclear matter, which cannot be reached in terrestrial experiments. The central density of NSs is expected to reach up to several times of the nuclear saturation density ($n_{0}\simeq 0.16\text{ fm}^{-3}$). The constraints on the attributes of matter at supernuclear densities (see, e.g. [1] and references therein) and relatively low temperatures (compared to particle collisions producing comparable energy densities) rely heavily on observations of macroscopic equilibrium parameters of NSs. Such macroscopic parameters as the total gravitational mass $(M)$ and the circumferential radius $(R)$ of NSs in binary systems are determined by high-precision timing observations of radio and X-ray pulsars. [2, 3, 4] The recent discovery of high-mass NSs of masses up to around $M\approx 2.01M_{\odot}$ (e.g., PSR J1614-2230 has yielded a mass of $1.927\pm 0.017M_{\odot}$ [5] and the mass of PSR J0348+0432 has been measured at around $2.01\pm 0.04M_{\odot}$ [6]) has ruled out several equations of state (EoS) models (see 6a). Due to the tremendous advances in the measurements, precise masses for $\sim 35$ currently known NSs range from $1.17$ to $2.01$ $M_{\odot}$. Also more than a dozen radii are known in the range of $9.9$ to $11.2$ kms, but current estimates for radii are still dominated by systematic errors. [7] More strict constraints on the plausible EoS models will be imposed by yielding masses and radii of a few stars to $\sim 5\%$ precision from the recently lunched NICER mission [8] and from the upcoming LOFT mission [9]. Constraints on the EoS, inferred from GW events GW170817 [10] and GW190425 [11] that are associated with binary NS mergers, are important confirmations of the great potential that multi-messenger observations offer.

X-ray and $\gamma$-ray burst phenomena are clearly explosive events that perturb their sources [12] and have been generally associated with neutron stars (see, e.g. [13] and references therein). In a manner analogous to terrestrial seismology, observational methods and techniques in asteroseismology are using the frequency of seismic waves rippling throughout stars with the aim of probing their internal structure and thermodynamic properties. [14] The frequencies appearing in the spectrum of NS oscillations, which sensitively rely on accurate stellar models, are matched to the observed frequencies. The period of stellar oscillations for non-relativistic stars are in the range of minutes, whilst for NSs the periods are much shorter; they typically range from about 0.2 to 0.9 milliseconds (cf. eq. (88)). [15] These oscillations occur when a star is perturbed away from its dynamical equilibrium and a restoring force tries to return it back to that equilibrium state. Among the various types of oscillation modes, in this paper we focus on the three lowest-frequency modes (fundamental and the next two higher modes) of radial oscillations where the pressure provides the dominant restoring force that produces oscillations. In this mode the Lagrangian displacement of fluid elements in the star is purely radial and the spherical symmetry is preserved. [16, p. 55]

As the radial modes are the simplest oscillation modes of NSs, they have been comprehensively investigated over the past half-century. In his pioneering papers [17], Chandrasekhar introduced a variational method to impose a sufficient criterion for the dynamical stability of radial and non-radial stellar oscillations. The identification of stable modes has been in the focus of interest ever since and methods for obtaining spectra of oscillation modes have been thoroughly investigated by various authors (e.g [18, 19, 20, 21, 22] and references therein). The first extensive analysis of radial modes for stellar models with various zero-temperature EoSs was presented by Glass & Lindblom in Ref. [21]. Although Glass & Lindblom used appropriate equations, the numerical results for the eigenfrequencies were essentially flawed, as Väth & Chanmugam later pointed out in Ref. [22]. The correct eigenfrequencies were computed for six EoSs by Ref. [22] and verified by the two different numerical schemes of Ref. [23]. In addition to having the earlier studies re-examined, Kokkotas & Ruoff also included further six zero-temperature EoSs in Ref. [23].

Our present paper, together with its accompanying paper [24], is devoted to study the radial oscillations of NS models. In the earlier paper [24], we proved that the fundamental equation for radial pulsation expressed by a set of effective variables which involve dissipative terms can be recast in a self-adjoint form. In contrast to the common non-dissipative case, the associated Sturm–Liouville eigenvalue problem (SL-EVP) is generalized for a discrete set of eigenfunctions with complex eigenvalues which correspond to the squared frequencies of the oscillation modes and the imaginary part corresponds to the damped solution. The imaginary part directly determines the minimum period of observable pulsars, hence it is imperative to make an accurate estimation on the relative timescale of thermal conductivity. In addition to providing a formulation of the dynamical equations that govern the radial oscillations of neutron stars through a perturbation scheme, an order-of-magnitude estimation was given for the timescale of energy dissipation, rather than a precise one, by an analytical approximation method without relying on explicit numerical computations. In the present paper, the SL-EVP for radial-oscillation modes is converted to a system of finite difference equations where we implement a second-order accurate differencing scheme so the resulting system of finite-difference equations emerges as a tridiagonal matrix eigenvalue problem. In a manner similar to the approach of Kokkotas & Ruoff [23], we compute the four lowest-frequency radial-oscillation modes of neutron stars constructed from various representative EoS models considered by Özel & Freire [2]. The rapidly growing list of observational data combined with advances in EoS modeling has motivated us to investigate neutron stars in a wide range of total gravitational mass and for various types of representative EoS models, listed in Table 1, based on nucleonic and hybrid nucleon–hyperon–quark matter models. The algorithm yields zero-frequency modes at the maxima and minima of the mass curves, as emphasized by [22, 23], while the adiabatic index of the equilibrium state characterizes the stiffness of the EoS at a given density.

Although proto-NSs [25] and strange stars [22, 26] were also studied with finite-temperature EoS, we resort to the study of zero-temperature EoSs only in this paper. Some of our equations in this paper involve dissipative terms which will be applicable in a following research, where advanced theoretical data on the heat transfer of dense matter shall be implemented in order to accurately compute the rate at which viscosity and thermal conductivity drain the energy deposited in oscillation modes. Properties of the bulk and the shear viscosities, and of the thermal conductivity in NSs have been studied more in detail by a number of recent of works [27, 28, 29, 30]. The zero-temperature approximation is evidently no longer sufficient in case one investigates damped oscillations, due to internal friction or viscous effects, as a result of the finite temperature inside the star. The density and temperature dependence of transport coefficients in the crust and in the core have been described for different EoSs by [27]. Thermal conductivity and shear viscosity of nuclear matter arising from nucleon–nucleon interaction in non-superfluid neutron-star cores were considered by [28], whereas those arising from the collisions among phonons in superfluid neutron stars were considered by [30]. The conversion of kinetic energy into heat and effects of viscosity on stellar pulsations in general have been addressed first by [31, 32, 33]. And the effects of radial heat flow on the dynamical instability and on the thermal evolution of NSs have been extensively studied by Herrera and his colleagues in Refs. [34, 35, 36].

II Realistic equation of state and thermodynamic properties of matter

Internal structure and macroscopic properties of NSs are strongly correlated with the EoS of dense matter, even though the exact EoS remains exceedingly uncertain especially at high densities. Although, the latest discovery of high-mass NSs PSR J1614-2230 [5] and PSR J0348+0432 [6] has ruled out several realistic models of EoS, suggesting that the maximal mass for NSs has to be larger than $M\sim 2M_{\odot}$ for a given EoS, but the number of candidate models with maximal mass below this limit is still considerably large. In NS cores, the temperature of matter is far below the Fermi energy of its constituent particles and its particular thermodynamic state at $T\simeq 0$ is accurately described by the isentropic one-parameter equation of state

[TABLE]

relating the pressure $p$ and energy density $\epsilon$ to the rest-mass density $\rho$ which exceeds nuclear density [37]

[TABLE]

In fact, densities in the cores are expected to be as large as $\rho\sim 5-10\rho_{0}$, where the nuclear matter at saturation (i.e. at the minimum of the energy per nucleon) has the density $\rho_{0}\simeq 2.8\times 10^{14}\text{ g/cm}^{3}$ or $n_{0}\simeq 0.16\text{ fm}^{-3}$, where the baryon-number density is related to the baryon-mass density as $n_{B}=\rho_{B}/m_{u}$ and $m_{u}=931.494$ MeV is the atomic mass unit. Given that neutrons geometrically overlap at $\rho\sim 4\rho_{0}$ and with increasing overlap between nucleons, transitions to non-nucleonic phase are expected. [2] It is possible for ultra-dense matter to contain hyperon, pion or kaon condensates. [38] Some of the possibilities considered to date also include free quarks or color superconducting phases. [39]

II.1 Realistic nuclear-theory-based EoS models

With the intention of covering a wide range of potential types of representative EoS models and generation methods, here we consider four EoSs of cold nucleonic matter (pure $npe\mu$ matter, i.e. the hypothetical components composed of neutrons, protons, electrons, and muons) and we follow the widespread naming convention of Refs. [2, 3, 4]:

• •

APR4 was derived by a variational method with modern nuclear potentials [40];

• •

MPA1 was derived by a relativistic Brueckner–Hartree–Fock theory [41];

• •

MS1 was derived by a relativistic mean-field theory [42];

• •

SLy4 was derived by a potential method [43];

and we also include three EoSs of ‘exotic’ non-nucleonic matter (i.e. the hypothetical components composed of hybrid nucleon–hyperon–quark matter):

• •

SQM1 is a hybrid EoS which describes relativistic non-interacting gas mixed with strange quark matter [44];

• •

H4 is a stiffer hyperon-based EoS derived by a relativistic mean-field theory [38];

• •

ALF1 is a hybrid EoS which describes a APR4 nuclear matter for a low density and a color–flavor-locked quark matter for a high density with the transition density $3\rho_{0}$, where $\rho_{0}\simeq 2.8\times 10^{14}\text{ g/cm}^{3}$ [39].

Due to the fact that the realistic EoSs are given in tabulated forms only, one has to interpolate between the given values of pressure and energy density to obtain stellar models with continuous functions of radius $r$. [45, p. 24] We have simultaneously interpolated the pressure and its first and second derivatives as functions of energy density. We used different degrees of spline interpolants, including linear, quadratic and cubic spline curves, to find the smoothest possible interpolating function at the tabulated set of points. The following densities are determined at the boundary of two neighboring pieces for the zero-temperature EoSs that are describe in Sec II.2 by piecewise-polytropes with $n=3$ pieces: $\rho_{1}\simeq 10^{14}\text{ g/cm}^{3},\ \rho_{2}\simeq 5.012\times 10^{14}\text{ g/cm}^{3},\ \rho_{3}\simeq 10^{15}\text{ g/cm}^{3}.$ ALF1 has the lowest pressure among the above considered types of EoSs and thus, making it the softest one. APR4, MPA1 and Sly4 have also relatively small pressure as in the case of ALF1 for a low-density region $\rho_{1}\leq\rho\leq\rho_{3}$, but for $\rho_{2}\lesssim\rho\leq\rho_{3}$, the pressure is higher than that for ALF1. Thus, for $\rho<\rho_{3}$, which NSs of canonical mass $1.3\text{ -- }1.4M_{\odot}$ have, these EoSs are soft as far as the canonical NSs are concerned. It is worthy of note that for a relatively small value of $p_{2}$, the adiabatic index, as illustrated in Fig. 2, is as large as $\Gamma_{2}\sim 3$, owing to the fact that the maximal mass of NS has to be $M_{\text{max}}\lesssim 2M_{\odot}$ for a given EoS. Thus, an EoS that is soft at $\rho=\rho_{2}$ has to be in general stiff for $\rho\gtrsim\rho_{3}$. However, MPA1 has pressure that even exceeds that of H4 for a high-density region $\rho\gtrsim\rho_{3}$. By contrast, H4 and MS1 have pressure higher than the rest for $\rho\lesssim\rho_{3}$, although the EoS becomes softer for a high-density region $\rho\gtrsim\rho_{3}$. In particular, MS1 has extremely high pressure among many other types of EoSs for $\rho\lesssim\rho_{3}$, and on that account, it is the stiffest EoS. All the distinguishing features mentioned above are reflected in Fig. 1 which displays the pressure in NS as a function of the baryon-number density or of rest-mass density. Table 1 lists the representative EoSs, the constituent particles, the maximal mass and circumferential radius of NSs of total mass $1.4M_{\odot}$, respectively, in descending order of compactness.

II.2 One-piece and piecewise-polytropic EoSs for heating and cooling processes

Although the tabulated, nuclear-theory-based EoS models listed in Table (1) are more realistic, they cannot be given in analytic terms over the whole range of energy densities inside NSs. A very common closed-form EoS is the polytropic one,

[TABLE]

which describes a non-interacting, degenerate matter. In general, the ‘polytropic constant’ $K=K(s)$ depends on the entropy, however, the degenerated matter dynamics in zero-temperature approximation can be modeled as an adiabatic flow with a constant $K$. The internal energy is given by the first law of thermodinamics for adiabatic process $(\delta Q=0)$, which can be integrated to obtain

[TABLE]

where we have imposed that $\lim\limits_{\rho\to 0}\epsilon/\rho=1$. [46] With the adiabatic assumption, eqs. (3–4) represent a barotropic fluid where the pressure is just a function of $\rho$. The adiabatic index $\Gamma_{1}$, defined by

[TABLE]

is an important dimensionless parameter characterizing the stiffness of the EoS (3) at a given density. [47, p. 190] For instance, a non-relativistic degenerate Fermi gas is reasonably well described by a polytropic EoS that scales as $p\propto\rho^{5/3}$, and for highly relativistic degenerate Fermi gases, $p\propto\rho^{4/3}$. Generally, $\Gamma_{1}$ and the speed of sound $c_{s}$ depend on the dynamical regime determined by $\rho$ (or $\epsilon$) as presented in Fig. 2 for the representative EoSs listed in Table 1. The non-monotonic variation of $\Gamma_{1}$ with baryon-number density (or rest-mass density) is revealed by Figs. 2 and 7, respectively. As a sudden increase of $\Gamma_{1}$ comes about at $\rho\simeq 1.3\times 10^{13}\text{ g/cm}^{3}$ for the stiffer MPA1 and MS1, but for most EoSs, it only occurs at $\rho\simeq 10^{14}\text{ g/cm}^{3}$. An abrupt change occurs in the stiffness of the matter when the Fermi energy of free neutrons exceeds the rest-mass energy of the neutrons. It is quite apparent that the matter undergoes phase transitions at this point which is associated with the ‘neutron drip point’ and only slightly depends on the particular EoS model. [23] The decay rate of the adiabatic index after its maximum also varies from one EoS to another, cf. Ref. [48]. The following two conditions must be satisfied by these EoSs:

The thermodynamic stability requires the EoS be monotonic ($dp/d\rho\geq 0$ and $dp/d\epsilon\geq 0$), and therefore $\Gamma_{1}$ must be positive. 2. 2.

The causality requires the speed of sound $c_{s}$ be less than the speed of light which is expressed by

[TABLE]

The constraint of eqs. (5–6) on the pressure-averaged adiabatic index, defined by

[TABLE]

demand that

[TABLE]

within a dynamically stable star of radius $R$. [49] Expressly, $\bar{\Gamma}_{1}<4/3$ is a sufficient condition to ensure the dynamical stability of spherical stars in Newtonian gravity (as we have outlined in Appx. B). However, in the stronger gravity of General Relativity (GR), even models with the stiffest EoSs become unstable for some value $R/M>9/8$. The more rigorous constraint on $\bar{\Gamma}_{1}$ for a star is to be stable against radial perturbation

[TABLE]

where $K$ is positive and of order of unity. For dynamical oscillations of neutron stars, the adiabatic index $\Gamma_{1}$ does not coincide with the polytropic one $\Gamma$, and the critical point (see Sec. V.1) implies a secular instability whose growth time is long compared to the dynamical timescale of stellar oscillations (88). [49]

It has been demonstrated by Ref. [4] that a piecewise-polytropic EoS with three pieces $(n=3)$ above the nuclear density approximately reproduces most properties of the representative EoSs listed in Table 1. These models of nuclear-theory-based zero-temperature EoSs at high density are modeled with a small number of parameters and the expression for pressure (3) is written in a parameterized form as

[TABLE]

where $n$ is the number of the pieces used to parameterize a EoS at high density, $\rho_{i}$ is the rest-mass density at the boundary of two neighboring $(i-1)$-th and $i$-th pieces, $K_{i}$ is the polytropic constant for the $i$-th piece, and $\Gamma_{i}$ is the adiabatic index for the $i$-th piece. Here, $\rho_{0}=0$, $\rho_{1}$ denotes a nuclear density evaluated in eq. (2), and $\rho_{n+1}\to\infty$. Other parameters $(\rho_{i};\,K_{i};\,\Gamma_{i})$ are determined by fitting with an EoS. Requiring the continuity of the pressure at each $\rho_{i},\,2n$ free parameters, say $(K_{i};\,\Gamma_{i})$, determine the EoS completely. [50]

One of the most serious drawbacks of the polytropic EoS is that, although (3) is a good approximation for a ‘cold’ star, there are extremely energetic processes, like the merger of stars or accretion from a disk, which can increase enormously the temperature and a simple polytrope will not provide a physical description. A more realistic EoS in closed form can be obtained by a combination of the polytripic EoS to describe the cold part and an ideal EoS for the thermal one, allowing for fluid heating due to shocks. The hybrid EoS is given by

[TABLE]

with an adiabatic thermal index $\Gamma_{\text{th}}$ that can be different from the adiabatic cold index $\Gamma$. The internal energy can be split into a thermal and a cold part,

[TABLE]

The total internal energy density $\epsilon$ can be obtained from the evolution of the conserved quantities, whilst the cold part is described by (4), leading to the explicit expression

[TABLE]

It is possible to extend this approach by using a collection of continous piecewise-polytropes (10) in hybrid EoS models, which in turn allows an accurate match with any type of realistic nuclear-theory-based EoS at high density.

II.3 Stress–energy tensor for dissipative fluids

Neutrino emission processes are supposed to be the main sources of energy loss in the stellar core in the later stages of stellar evolution. For the reason, the equations of relativistic fluid dynamics to describe energy–momentum conservation are written as

[TABLE]

where $Q_{\nu}$ is the total neutrino emissivity of all processes outlined in Table 2, and the stress–energy tensor for dissipative fluids can be written as the sum of three components

[TABLE]

where

[TABLE]

are the perfect fluid, heat flux and viscosity stress–energy tensors, respectively. [52] Note that despite its causality and stability problems, the above description of stress–energy tensor has been widely used in Eckart’s theory of relativistic irreversible thermodynamics [53]. As in eq. (1), the variables $p$ and $\epsilon$ represent the isotropic pressure and energy density measured by a comoving observer with velocity $u^{\alpha}$, which satisfies $u^{\alpha}u_{\alpha}=1$ with $u^{0}>0$, and

[TABLE]

is the standard projection tensor onto 3-space normal to flow. The symmetric trace-free spatial shear tensor is defined as

[TABLE]

and expansion scalar (or dilatation rate)

[TABLE]

is associated with the convergence (or divergence) of the fluid world lines. The heat flux density in eqs. (14–16) written as

[TABLE]

is a spacelike vector, $Q^{\alpha}u_{\alpha}=0$, that describes the flow of thermal energy per unit of area along spatial cordinate axis $x^{\alpha}$ per unit of time. The first term in eq. (20) corresponds to the non-relativistic Fourier’s law of heat conduction, the second term takes into account the relativistic effect of isothermal heat flux due to the inertia of energy with $a_{\beta}=u^{\gamma}u_{\beta;\gamma}$ being the acceleration of fluid. The negative sign indicates that heat flows from higher to lower temperature regions.

In eqs. (16) and (20), $\eta,\,\zeta$, and $\kappa$ are collectively called transport coefficients (or dissipation coefficients). The bulk viscosity coefficient $\zeta$ defines the resistance of the medium to gradual uniform compression or expansion; and $\kappa$ is non-negative and accounts for the thermal conductivity, respectively. [54] The shear (also called as ‘common’ or ‘dynamic’) viscosity coefficient $\eta$ describes the fluid’s resistance to gradual shear deformation which is assumed to be equal to the electron shear viscosity $\eta_{e}$ in the stellar core; and we take it from [52] in this paper. Shternin and Yakovlev (2008) [27] concluded that the shear viscosity of neutrons and of protons can be neglected for the reason that it depends strongly on the nuclear interaction model and the many-body theory.

III Equilibrium stellar models

We consider a static spherically symmetric star, described by the metric

[TABLE]

where $t$ and $r$ are the time and radial coordinates, $d\Omega$ is a solid angle element in a spherical frame with the origin at the stellar center and

[TABLE]

are the metric functions. The later function is often replaced by the expression

[TABLE]

where the gravitational mass contained within the radius $r$ is given as

[TABLE]

Outside the star the metric reduces to the Schwarzschild form:

[TABLE]

The characteristic scale factor, defined by

[TABLE]

is referred to as Schwarzschild radius or ‘gravitational radius’ and the total gravitational mass of the star, defined by $M=m(R)$, which acts as the source of the gravitational field outside the star $(r>R)$. The crust corresponds to the layer $r_{cc}<r<R$, where $r_{cc}$ determines the crust-core interface. The proper depth below the stellar surface is defined by

[TABLE]

as the proper radial distance between the stellar surface $R$ and a given layer of radius $r$. The proper depth and the metric functions of NSs of total mass $M=1.4M_{\odot}$ with SLy4 EoS in their cores are shown in Figs. 3–4 as function of radius in units of km. (Cf. Figure 37 in Ref. [55] and Figure 4.1. in Ref [56].)

As in Ref. [24], the physical variables for energy density $\epsilon$ and for isotropic pressure $p$ are replaced by the corresponding effective variables

[TABLE]

that incorporate time-dependent dissipative contributions of the stress–energy tensor (16). From now on, we will only use these effective variables and for the sake of simplicity, we shall not put any bar over them. The structure of compact stars in hydrostatic equilibrium is then constrained by the following set of equations:

[TABLE]

where overdots denote partial differentiation with respect to $t$. [22] The mass-radius equation (29a) and the Tolman–Oppenheimer–Volkoff (TOV) equation (29b) constraint the stellar structure in static gravitational equilibrium, whereas eq. (29c) determines the metric function $\nu(r)$. Supplemented with a given EoS (1), the set of equations (29) can be integrated from the center of the star with initial condition $m(0)=0$ and an arbitrary value for the central density $\rho_{c}=\rho(0)$, until the pressure $p(r)$ will vanish at some radius $R$ which corresponds to the star’s surface.111 It is clear from eq. (29b) that $p(r)$ is strictly decreasing with increasing $r$. However, for some of the representative EoSs, $p$ may not drop to zero by the time $dp/dr$ becomes negligible. Therefore the numerical integration has to be stopped when the decrease of $dp/dr$ comes to a halt and the value of $p$ has fallen sufficiently low. Even though in Sec. II.1 we interpolate between the tabulated values of pressure and energy density such a way that produce the smoothest possible EoS, the stellar models obtained from the integration of eqs. (29) are not as smooth as those of polytropic ones. This hindrance is especially evident in the vicinity of the stellar surface, where the pressure is zero, but it gets immensely increased as one moves into the stellar interior. (See footnote 1 and Fig. 5.) The underlying reason is that $\rho(r)$ and other matter functions contain unwanted bumps and sharp bends as the matter undergoes phase transitions for increasing pressure. [45, p. 24] As stated in Sec. II.2, these sharp phase transitions from hadronic state into another (e.g. hyperons, pion or kaon condensates) change the state and the stiffness of matter quite abruptly, together with the adiabatic index $\Gamma_{1}$. The stiffening of EoS plays a crucial role in setting the maximal mass: the maximal mass increases as the EoS gets stiffer at high densities. [2] The extra pressure from the repulsive force raises the maximal mass. Considering that among the EoS models present in this paper, MS1 allows the largest maximal mass to pile up, Fig. 5 displays the radial profile of isotropic pressure for different central densities of NSs with MS1 EoS in their core.

For each particular EoS there is a unique family of equilibrium configurations which is parameterized by $\rho_{c}$. Figs. 6a–6b show the corresponding mass–radius and mass–central energy density relations, respectively, which can be interpreted as sequences of stellar models $M=M(R)$ and $M=M(\rho_{c})$. $M(R)$ curves demonstrate that low-mass configurations, having low central densities, mainly below $\rho_{0}$, depend dominantly on the low-density part of the EoS, whereas the maximal mass is hardly affected (see also [12]). On the other hand the Harrison–Zel’dovich–Novikov criterion [18, 19] for static stability of compact stars implies that $\partial M(\rho_{c})/\partial\rho_{c}\geq 0$. As mass piles up with increasing $\rho_{c}$, the $M$–$R$ relation probes the higher density part of the EoS, in the range $\rho_{0}\lesssim\rho_{c}\lesssim 3\rho_{0}$. Lattimer & Prakash [3] showed that across various EoS models, the scaling of radius is approximately $R\propto p^{1/4}$ at a central density $\rho_{c}\backsimeq 1.5\rho_{0}$, so that a radius measurement is probing the EoS, just above nuclear density (2). The $M(R)$ curves that have attained sufficient mass exhibit vertical segments with radii varying from $10$ to $16$ km. The vertical segments can be understood from the scaling of hadronic EoS models which is typically close to $p\propto\rho^{2}$ for higher densities (as demonstrated by $\Gamma_{1}\backsimeq 2$ in Fig. 2). In this case, the radius is nearly independent of mass. Those models of hadronic EoS with extreme softening (due to a kaon or pion condensate, high abundances of hyperons, or a low-density quark-hadron phase transition) do not have pronounced vertical segments. [59]

IV Adiabatic radial perturbations and speed of their propagation

To move beyond the state of equilibrium, we now introduce small, adiabatic perturbations of eqs. (29) in such a way that does not violate the spherical symmetry of the star.222Note that the inclusion of transverse perturbations results in the appearance of non-radial oscillation modes. Although non-radial modes of NSs (such as g-modes and r-modes, being in the sensitivity bands of aLIGO and aVirgo detectors) are of great importance for being accompanied by the emission of gravitational waves, in order to make the problem more easily tractable, we specialize to purely radial modes in this paper. The metric functions from the line element (21) can then be recast as

[TABLE]

where a small parameter $\delta\ll 1$ denotes the ratio between the scale of variation of the metric functions $(\nu_{0},\,\lambda_{0})$ corresponding to the unperturbed equilibrium stellar model and $(\delta\nu,\,\delta\lambda)$ are metric perturbations due to the radial pulsations. To obtain the equations that govern the radial pulsations, let us consider a fluid element displaced radially outward from an initial position at $r$ to $r+\delta r$. The fluid element expands (or, if displaced inward, contracts) with the radial velocity

[TABLE]

as its pressure adjusting with the group velocity (50) across the fluid to the external pressure, given by

[TABLE]

Note that the term involving $\delta T_{\text{visc}}$ and $\delta T_{\text{heat}}$ represents the radial component of linear momentum density and is present only if viscosity or heat transfer occurs across the fluid elements, cf. (74). We assume that the motion of the fluid element, moving with a four-velocity defined by

[TABLE]

is faster than the time for heat to flow out of the fluid element. The perturbation is thus adiabatic, and from eqs. (5–6) the change in the energy density of the fluid element due to the perturbation is

[TABLE]

where the adiabatic index of the perturbation $\Gamma$ is assumed to be equal to the adiabatic index of the equilibrium configuration $\Gamma_{1}$ in eq. (5).

Let us consider the radial displacement of the fluid element $\delta r$ which satisfies a linear wave equation being expanded in an asymptotic series

[TABLE]

where the angular frequency and the wavenumber (or spatial frequency), defined as

[TABLE]

are the coordinate components of the four-wavevector, respectively. [60] The higher-order terms $X_{n}$ in the series (35) become successively smaller in the small parameter $\delta\ll 1$ where $\partial X_{n}/\partial t=\partial X_{n}/\partial r=0$. The equation for $X_{0}^{2}=X_{0}X_{0}^{*}$ can be obtained from the Lagrangian

[TABLE]

where $G(\omega,k)=0$ is the desired dispersion relation. The Lagrangian for the pulsation equation (52) is

[TABLE]

where one can take the first-order approximation of the radial displacement (35):

[TABLE]

and regard the $\sin^{2}$ and $\cos^{2}$ terms as time averaged to 1/2. [60] As reported by Ref. [61, p. 392–399], the dispersion relation can be identified from the Lagrangian as

[TABLE]

Now, let us assume, instead of (35) and (39), that the radial displacement possesses a simple harmonic time-dependence of the form

[TABLE]

where $X(r)$ is the radially-dependent amplitude and $\omega$ is a complex characteristic eigenfrequency which consists of the so called ‘damped frequency’ and ‘relaxation time’,

[TABLE]

respectively. The aforementioned quantities are related to the natural frequency $\omega_{n}$ of the undamped system and the damping ratio $0\leq\zeta_{d}\leq 1$ which describes the rate at which the normal oscillation modes decay and critically determines the dynamical behavior through the energy-dissipation rate (44). [24] Being bilinear in the adiabatic perturbations, the total-energy density $\epsilon$ has a time dependence $\exp[-2\operatorname{Im}(\omega)t]$. [62] Subsequently, the time derivative of the total-energy density implies that

[TABLE]

which together with the energy-dissipation rate

[TABLE]

directly determines the characteristic damping time of perturbations as

[TABLE]

In order to discuss causality, eqs. (36) must be transformed to a local Lorentz frame. With respect to $(e^{-\nu_{0}/2}\partial_{t},\,e^{-\lambda_{0}/2}\partial_{r})$ the transformed quantities are

[TABLE]

The transformed dispersion relation is

[TABLE]

where

[TABLE]

and

[TABLE]

The group velocity with respect to local light cones is

[TABLE]

which is always greater than the local speed of sound and therefore $c_{s}<1$ is merely a weak condition (the very least required by causal wave propagation). [21] In line with Jeans’ criterion for general relativistic spherical stars, for all sufficiently long wavelengths such that

[TABLE]

a dynamical instability sets in. In the Newtonian limit, $\mu=2(Gm/r^{3})^{1/2}$ gives the spherical version of Jeans’ classic result for an infinite homogenous medium. [63] Note that there is a local, purely general relativistic instability at any point in the sphere where $2m(r)$ approaches $r$. This supports the well-known result that configurations at equilibrium must have $2m<r$ at all interior points.

V Linear pulsation equation and associated eigenvalue problem

Chandrasekhar [17, eq. (61)] showed that, with the assumption of harmonic time-dependence (41), the fundamental equation for radial pulsation can be cast in the form

[TABLE]

where the functions of equilibrium stellar configuration $\{p(r),\,\rho(r),\,m(r),\,\nu(r),\,\lambda(r)\}$ in eq. (29), the effective adiabatic index $\Gamma(r)$ and sound speed $c_{s}$ in eqs. (5–6) and

[TABLE]

are computed for specific stellar models. [23]

The fundamental equation for radial pulsation (52) together with its boundary conditions (56–57) constitutes a Sturm–Liouville eigenvalue problem (SL-EVP) for a discrete set of scalar-valued eigenfunctions of radial displacement $\{X_{0}(r),\,X_{1}(r),\,\ldots,\,X_{j}(r),\,\ldots\}$ with their respective eigenvalues $\{\omega^{2}_{0},\,\omega^{2}_{1},\,\ldots,\,\omega^{2}_{j},\,\ldots\}$. Successive eigenvalues correspond to squared complex-valued frequencies of the oscillation modes, arranged in ascending order, where the index $j=0,\,1,\,\ldots$ represents the number of nodes of each respective mode inside the star. The smallest eigenvalue $\omega_{0}^{2}$ is associated with the fundamental-mode frequency of radial oscillations which has no nodes between the center and the stellar surface, whereas the first excited mode $(j=1)$ has a node, the second one $(j=2)$ has two, and so forth. The normalized eigenfunctions functions of radial displacement are presented in Fig. 8. In terms of a superposition of sinusoids, the $\omega_{0}$ is the lowest-frequency sinusoidal in the sum (see dynamical stability in Sec. V.1). Owing to the fact that the general solution to eq. (52) is always a superposition of both growing and damped oscillation modes, the definitions (42) provide the means to separate the real and imaginary parts of the complex frequency:

[TABLE]

thus the imaginary part of inhomogenity of the eigenvalue problem appearing on the right-hand side of eq. (52) reduces to

[TABLE]

where the functions $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are expressed by eq. (75), in terms of transport coefficients $\eta,\,\zeta$, and $\kappa$. The derivation of eq. (52), which has been extended to involve transport coefficients in neutron-star matter through $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, has been outlined in detail in our preceding paper [24] and shall not be repeated here.

The homogeneous SL-EVP (52) is associated with harmonic oscillations and can be converted to a self-adjoint form which poses a boundary-value problem. [24] The eigenvalues and corresponding eigenfunctions of the given boundary-value problem are then uniquely determined in Sec. VI by the finiteness of $X(r)$ and $dX(r)/dr$ everywhere and by satisfying the following pair of boundary conditions of the third type333The Robin boundary condition or third-type boundary condition is a specification of a linear combination of the values of a function and its derivative on the boundary. It is commonly used in solving SL-EVPs.: [21]

The fluid at the center of the star is assumed to remain at rest; thus,

[TABLE] 2. 2.

The perturbed pressure is required to vanish on the perturbed boundary of the star; thus the Lagrangian change in the pressure $\delta p=e^{\nu_{0}/2}r^{-2}\Gamma p_{0}\frac{d}{dr}\left(e^{-\nu_{0}/2}r^{2}X\right)$ has to vanish at the surface; thus,

[TABLE]

Let us keep in mind that even the homogeneous SL-EVP cannot be solved analytically in general. However, a rough order-of-magnitude estimate can be made. Provided that $X_{0}(r)$ is nearly linear in $r$ (cf. blue line in Fig. 8), $\omega_{0}\propto(G\hat{\rho})^{-1/2}$ is a universal estimate for the f-mode, irrespective of the particular EoS, where $\hat{\rho}$ is the mean density of a star undergoing homologous contraction (see eq. (88) in Appendix B). Detweiler & Lindblom [15] has found that for most stellar models, the oscillation periods of the fundamental mode

[TABLE]

typically range from about 0.2 to 0.9 milliseconds, depending on the particular EoS and the central density. The ‘ordinary’ or ‘temporal’ frequency $\nu_{0}$ and angular frequency $\omega_{0}$ are typically $\sim 10^{4}$ $\text{s}^{-1}$ for neutron stars, but two quantities differ by a factor of $2\pi$. Accordingly, $\tau_{d}$ in eq. (88), which is in the range of 0.1 – 0.3 seconds, can also serve as an estimate of a hydrodynamical timescale for a given star. [64, p. 291] In order to compute the three lowest-frequency oscillation modes with high accuracy for each particular EoS in Table 1, a finite-difference method will be used in the subsequent Sec. VI.

V.1 Dynamical stability criteria based on the fundamental mode

In Sec. III the hydrostatic stability of equilibrium stellar configurations was shown to be firmly ensured by the TOV equation (29b). On the other hand, the dynamical stability of stellar configurations with respect to small radial perturbations is determined by the fundamental mode:

If $\omega_{0}^{2}>0$, then all the higher eigenfrequencies of the spectrum are real, which indicates that the equilibrium stellar model is dynamically stable with respect to small perturbations, and they form a lower-bounded infinite sequence $\omega_{1}\leq\omega_{2}\leq\omega_{3}\ldots$. In this case the small radial displacements (35) exhibit purely harmonic oscillations in the f-mode. [18] 2. 2.

If $\omega_{0}^{2}<0$, then $\omega_{0}$ becomes purely imaginary, which indicates that the equilibrium stellar model comes to be dynamically unstable (at least in the f-mode). In this case the small radial displacements (35) exponentially grow or damp with time (depending on the value of the damping ratio $\zeta_{d}$). As discussed in detail by our preceding study [24], the highest instability increment, associated with the most rapid exponential growth or dampening, is provided by the f-mode. Instability will indeed occur at a specific central energy density $\epsilon_{0}>\epsilon_{\text{crit}}$ in the crust of those neutron stars with total mass $M$ that exceeds the maximal-mass configuration of a family of non-rotating NS models for the particular EoS. Such NSs would eventually collapse into a black hole. [19] 3. 3.

Finally, at $\omega_{0}=0$ the equilibrium stellar model is at neutral stability; neither stable nor unstable. In this case $\epsilon_{0}=\epsilon_{\text{crit}}$ as marked in Fig. 6b by the symbols $\diamond$ for states at neutral stability.

The dynamical stability is also inherently related to the shape of the $M(R)$ curve. [64, p. 293] As we already stressed in Sec. III, $dM/d\epsilon_{c}>0$ is necessary but not sufficient condition for dynamical stability. At each critical point where $dR/d\epsilon_{c}<0$ holds an even number $n$ of radial oscillation mode changes its stability while for $dR/d\epsilon_{c}>0$ an odd-numbered mode changes its stability. [18] The dynamically stable and unstable branches of $M(R)$ curves are presented in Fig. 9 by solid lines and dashed lines, respectively. The left panel schematically illustrates the spiral structure of a general $M(R)$ curve for dynamically unstable states whereas the right panel depicts the unstable branches for three selected EoS (APR4, H4, SLy4) with particularly small radii ($R\simeq 8-9$ km) at their respective final critical points. The letters $C_{1},\,C_{2},\,C_{3}$ refer to critical points (turning points) which divide the $M(R)$ curve into four segments, with an increasing number of unstable oscillation modes. As we have seen in Fig. 6a, the central energy density $\epsilon_{c}$ is increasing from large radii and small masses until a certain critical central energy density $\epsilon_{c}^{\text{crit}}$ is reached at the maximal mass and the minimal radius. All radial oscillation modes are stable in the $\epsilon_{c}<\epsilon_{c}^{\text{crit}}$ regime, but once the central energy density is increased beyond the critical point $C_{1}$, it does not restore the stability. [65] At $C_{1}$ the fundamental mode becomes unstable. At $C_{2}$, on account of $(dR/d\epsilon_{c})_{C_{2}}>0$ the fundamental mode remains unstable and an additional $n=1$ mode becomes unstable. At $C_{3}$, all three oscillation modes become unstable.

VI Finite-difference method for computing eigenfrequencies

Various numerical methods have been investigated to accurately compute the eigenfrequencies of radial pulsation for particular nuclear-theory-based EoS models. The scholarly literature [21, 22, Kokkotas, 60] refers to two most widespread approach to solve the SL-EVP by Runge–Kutta integration as ‘shooting method’ and ‘finite-difference method’, respectively. We shall adopt the later method by converting the linear ODE (52) into a system of algebraic equations by replacing the derivatives with finite-difference approximations in the form of a tridiagonal matrix that can be solved by matrix algebra techniques, ideally suited to modern numerical analysis.

Let some grid functions $\{X_{0},\,X_{1},\,\dots,\,X_{N},\,X_{N+1}\}$ be defined on a grid with uniformly-spaced points $r_{j}=j\Delta r$ with a separation $\Delta r=R/(N+1)$ between adjacent grid points. An illustration of thin spherical shells of matter with thickness $\Delta r$ at successive grid points $r_{j-1},\,r_{j},\,r_{j+1}$ is shown in Fig. 10. We consider the finite-difference approximation $X_{j}$ of some eigenfunction $X(r_{j})$ for any positive integer $j=0,\,\ldots,\,N$. The discretization and evaluation of eq. (52) on the $j$th shell is formally expressed by

[TABLE]

where $\omega^{2}_{j}$ is the squared frequency of the $j$th oscillation mode, $[D_{r}X]_{j}$ and $[D^{2}_{r}X]_{j}$ are the first- and second-order discrete derivatives of the grid function $X(r)$ at a corresponding grid point $r=r_{j}$ for $j=1,\,\ldots,\,N$, and the coefficients $a_{j},\,b_{j},\,c_{j}$ are approximated by

[TABLE]

Let us consider approximation of the derivatives of grid functions in (59) by the following finite central differences:

[TABLE]

with the truncation error $\epsilon_{t,D}=\epsilon_{t,D^{2}}=\mathcal{O}(\Delta r^{2})$ is to be assessed in Sec. VI.1. By re-expressing the finite differences (61), eq. (59) reduces to a set of algebraic equations

[TABLE]

for $j=1,\,\ldots,\,N$ which represent the three successive non-zero entries of the $j$th row of a non-singular tridiagonal matrix $[\Sigma_{ij}]\in\mathbb{R}^{N\times N}$, defined by

[TABLE]

where the superdiagonal, diagonal and subdiagonal entries are as follows:

[TABLE]

Although $b_{j}$ prevents the matrix (63) from being symmeric, but due to its tridiagonal structure (inhereted from the central differencing scheme which involved only the closest pair of points), zero entries can be discarded and a significant amount of computational cost is saved. Consequently, a matrix equation of the form

[TABLE]

is constructed from the linear system of $N$ equations (62) for the central differencing scheme where $I$ denotes the $N$-by-$N$ identity matrix, and

[TABLE]

is the vector which consists of the unknown values $\{X_{1},\,\dots,\,X_{N}\}$ of the eigenfunction of a specific oscillation mode at corresponding grid points $\{r_{1},\,\dots,\,r_{N}\}$. The first and last entry of $\vec{X}$ involve the values $X_{0}$ and $X_{N+1}$, respectively, as a result of the pair of boundary conditions (56–57) which imply $X_{0}=X_{N+1}=0$.

In linear algebra, Cramer’s rule states that the homogeneous linear system of equations (65) has a non-trivial solution if and only if the determinant of the corresponding characteristic matrix is zero, expressly

[TABLE]

The characteristic polynomial of the tridiagonal matrix $\Sigma$, denoted by $p_{\Sigma}(\omega^{2})$ and defined by

[TABLE]

is a monic polynomial of degree $N$ which is invariant under similarity transformation and has the eigenvalues as roots. In other words, eq. (67) is referred to as the characteristic equation, considering that the roots of $p_{\Sigma}(\omega^{2})=0$ are exactly the eigenvalues of $\Sigma$. Harrison et al. [18] have argued that a sufficiently accurate numerical algorithm always yield a zero-frequency mode at that specific central density that corresponds to the maximal-mass configuration of a family of stellar models for the particular EoS. In fact, by computing $p_{\Sigma}(\omega^{2})$ for each particular EoS at various values of central density, all the eigenfrequencies and corresponding eigenmodes are simultaneously determined in a single run. Kokkotas & Ruoff [23] pointed out that although this method is more time consuming than searching for each eigenvalue separately by shooting method, but one may smoothly refine its accuracy by setting up grids with optimal size. In the next section the truncation error in the eigenvalues due to the finite-difference approximations (59) will be estimated. We demonstrate that it does not exceed ca. $10^{-4}$ Hz for setting up regular grid of size $N\simeq 3000$.

VI.1 An error estimate of the finite-difference approximation

Two potential sources of numerical errors occur in computing our finite-difference approximation to the solution of the SL-EPV. Errors due to inexactness in the floating-point representation of numbers and the arithmetic operations are called round-off errors and they are determined by the machine precision (roughly $10^{-15}$). Moreover, local truncation errors (LTE) incur due to the inaccuracy in the approximation (61) for derivatives of grid functions. The later one, however, can be controlled by varying the spacing of the radial grid on which equilibrium stellar models and the eigenfunctions of radial displacement are determined. [66]

As stated in the previous section, the discrete derivative $[D_{r}X]_{j}$ is not exactly equal to the derivative $X^{\prime}(r_{j})$. The truncation error in the first-order centered finite-difference approximation (61a), expressed by

[TABLE]

can be computed in the following manner by expanding $X_{j+1}$ and $X_{j-1}$ in a Taylor series around the point $r_{j}$:

[TABLE]

Now, by inserting the quadratic approximations (70) into the finite difference (61a) and collecting terms in $[D_{r}X]_{j}-X^{\prime}(r_{j})$, the result for the truncation error is

[TABLE]

an expression given in terms of a power-series in $\Delta r$ where the lowest power is proportional to $\Delta r^{2}$. [67, p. 419] The truncation error for the approximant of the second-order derivative (61b) is

[TABLE]

The truncation error, therefore, in both expressions is said to be of second order in $\Delta r$, as we stated after eq. (61). The accuracy of our computations was tested by varying the spacing of grid points, defined by $\Delta r\equiv R/(N+1)$. Our grids were set up over $N\simeq 3000$ points for each stellar models where the circumferential radius ranged from about $8$ to $14$ km, therefore the spacing varied between $2.666$ and $4.665$ metres (or equivalently between $8.892$ and $15.561$ nanoseconds of oscillation period $T$). Subsequently, our tests indicated that the truncation error (71) in the computation of eigenfrequencies did not exceed $\epsilon_{t}\sim 10^{-4}$ Hz.

VI.2 Numerical results and interpretation

This subsection discusses the numerical results obtained by the implementation of the finite-difference method for the homogeneous SL-EVP (52). We present the first three radial-mode frequencies as a function of central energy density in the form of tables and accompanying figures for each realistic EoS. In each panel of Table 3 we list the central energy density ($\epsilon_{c}$), the circumferential radius ($R$), the gravitational mass ($M$), the frequencies of the fundamental mode ($\nu_{0}$) and the first two lowest-frequency excited modes ($\nu_{1}$ and $\nu_{2}$) of the respective stellar model.444Let us note that the eigenfrequencies in Sec. VI were in fact angular frequencies $(\omega)$, whereas the use of ‘ordinary’ or ‘temporal’ frequencies $(\nu)$, which differ only by a factor of $2\pi$ in eq. (58), are more conventional in asteroseismology. In addition to the stable model, we also included a stellar configuration just beyond the limit of dynamical stability where, as stated above, the total mass $M$ as function of the central energy density $\epsilon_{c}$ exhibits a local maximum. The f-mode frequency of such configurations is marked by an asterisk. Considering that the f-mode oscillation period (58) relies extremely sensitively on the particular equilibrium stellar model (and thus inherently on the choice of interpolants used to fill the gaps in the associated discontinuous EoS, the smoothness of $\Gamma(\epsilon)$, etc.), slight deviances between each numerical predictions of the scholarly literature [20, 21, 22, 23] and also our results for the eigenfrequencies are anticipated. Ref. [22] pointed out that the frequencies obtained by different authors are not substantially different at low central densities, where effects of GR remain substantially negligible, but at higher central densities where the redshift of the star $z=(1-2GM/Rc^{2})^{-1/2}-1$ exceeds $z\simeq 0.1$, discrepancies occur which may be assigned to the slight deviations in the equilibrium stellar models.555As we started in Sec. III, equilibrium configurations are generally parameterized by $\rho_{c}$, $\epsilon_{c}$ or equivalently by $M$. However, the oscillation spectra of stellar models constructed from different EOSs are more directly comparable with a parametrization by the redshift parameter $z$ than by either of the earlier parameters (see [21, 1])..

On each panel of Fig. 11, the three lines represent the three lowest-frequency eigenmodes of radial oscillation as functions of central energy density. On each panel, the f-mode frequency drops towards zero as the particular stellar model approaches its dynamical stability limit which, indeed, is indicated by the presence of an eigenmode with zero frequency (cf. [45, p. 21]). The dynamical instability in stars with MS1 and APR4 EoSs is exposed by the presence of a very low frequency of the f-mode, which has dropped to less than $5\%$ of that of the first excited mode, at central energy densities associated with the maximal-mass stable configurations. The decay of the lowest-frequency eigenmodes appears to be a general feature of Fig. 11, irrespective of the particular EoS. The oscillation frequency of higher modes is always larger than that of a lower stable mode and for all modes it appears to decrease as the central energy density approches the smallest possible value $\epsilon_{\text{min}}$ of the particular stellar model. The explanation for our observations can thus be traced back to the fact that NSs at high densities can be approximated as being homogeneous [68] and the $\omega_{0}^{2}$ can be given by the relation

[TABLE]

Ref. [22, p. 253] confirms that when the central energy density of NSs is approaching $\epsilon_{\text{min}}$, such compact objects become approximately homogeneous and due to their small mass, non-relativistic [69]. The relation 73, which shall be explicitly explained by eq. 86, expresses that the eigenvalue $\omega_{0}^{2}$ for the fundamental mode is proportional to the squared average density, as well as that the value of $\omega_{0}^{2}$ scales with the adiabatic index $\Gamma$.

It is also apparent from Fig. 11 that stellar models of softer EoSs have higher frequencies in the f-mode than the stiffer ones for the same central density. As reported by Ref. [12], stellar models of softer EoSs are generally associated with more centrally condensed stars with larger average densities. Having the highest pressure among the considered types of EoSs at a high-density region $\epsilon\lesssim 0.57\text{ GeV fm}^{-3}$, MS1 and H4 are the stiffest EoSs and producing NSs with comparatively low-frequency oscillations. When the central energy density exceeds $\epsilon\lesssim 0.42\text{ GeV fm}^{-3}$, the frequency for stars with H4 EoS decreases by about $15\%$ in the fundamental mode compared to the frequency of the MS1 EoS. This is due to the fact that H4 being a stiffer mixture of nuclear matter at low densities and hyperons at high densities (cf. Sec. II.1), a phase transition to hyperon-admixed matter takes place in the core of NSs. As reflected in the decreasing frequency, the phase transition causes a drop of adiabatic index with respect to stars with MS1 EoS, which do not contain hyperons but are otherwise mostly similar.

Appendix A Components of the stress–energy tensor

The non-vanishing components of the viscosity and heat flux stress–energy tensors defined by (16) and quantities involving them shall be explicitly evaluated below. Here, we only repeat those quantities from our preceding paper [24] that appear in the present paper. The non-vanishing components of the unperturbed and perturbed states are written in the notation we used in this paper as

[TABLE]

respectively. The inhomogenity of the Sturm–Liouville equation (52) that determines the damped solution is expressed by

[TABLE]

where $\delta u_{r}$ is the radial velocity defined by eq. (31) and

[TABLE]

Appendix B Dynamical stability of non-rotating spherical stars in Newtonian gravity

This section provides a Lagrangian description of radial oscillations of stars in the framework of classical hydrodynamics in order for one to find a linear approximation for the dynamical times cale of these oscillations. Let us consider a mass element $dm$ at a radius $r$ inside the star at a given moment $t=t_{0}$ contained in a thin spherical shell with a thickness $dr$ (see Fig. 10). Considering that the spatial coordinate of the given mass element does not vary in time, it is more convenient to take a Lagrangian coordinate, which is connected to the mass element

[TABLE]

instead of the radius $r$. [70] The spatial behavior of the radial function $r(m,t)$ is then described by differential equation

[TABLE]

in accordance with the definition of gravitational mass (24). The Navier–Stokes equation of hydrodynamics with spherical symmetry provide the greatly simplified expression

[TABLE]

for the equation of motion of mass elements, where

[TABLE]

are the gravitational force gravity and pressure gradient, respectively. [71] In order to prevent the mass elements of the shell from being accelerated by the gravitational force towards the center of the star (as indicated by the minus sign), a net force of the same absolute value due to pressure acts on the shell in an outward direction (see Fig. 10). In the state of hydrostatic equilibrium $\partial^{2}r/\partial t^{2}=0$ and the sum of the forces arising from the gravity and pressure has to be zero as well, which yields the ODE of stellar structure

[TABLE]

where the pressure, defined by

[TABLE]

is consequently always greater at the lower boundary than at the outer boundary of the shell, expressively, $p_{j}>p_{j-1}$ (see Fig. 10). The pair of basic equations (78) and (81) in the Lagrangian description corresponds to the conservation of mass and of momentum, respectively.

Let us now recall the small, adiabatic perturbations of the equilibrium stellar configuration described in Sec. IV. We have seen that the change in the radius of every mass shell is $r=r_{0}(1+\delta)$, where the small parameter $\delta\ll 1$ may vary only with time (see eq. (41)). Mass conservation expressed by eq. (78) entails that the change in the rest-mass density is $\rho=\rho_{0}(1+3\delta)$. Considering that the perturbation is adiabatic, the definition (5) for the adiabatic index $\Gamma_{1}$ implies that

[TABLE]

so the change in the pressure is $p=p_{0}(1-3\delta)$, where the pressure and density perturbations are assumed to be very small: $|(p-p_{0})/p_{0}|\ll 1,\quad|(\rho-\rho_{0})/\rho_{0}|\ll 1$. [72]

The perturbed quantities no longer satisfy the equation of hydrostatic equilibrium (81), but rather the equation of motion (79). As the only time variable quantity is $\delta$ we may write the equation of motion in the form

[TABLE]

where any quantity with the subscript ‘0’ satisfies the equation of hydrostatic equilibrium (81). Therefore, combining equations (78) and (84) we obtain

[TABLE]

Now we make another crude approximation, $M/r_{0}^{3}\approx\bar{\rho}=\text{const.}$ where $M$ is the total mass of the star, and we obtain

[TABLE]

This equation has a solution

[TABLE]

If $\omega^{2}<0$, then the motion is oscillatory, i.e. the star is dynamically stable. If $\omega^{2}>0$, then there is a solution that increases exponentially, i.e. the star is dynamically unstable. The proper analysis leads to the requirement, which was already pointed out in Sec. II.2, that the average value of $\Gamma$ within the star has to be less than $4/3$ for the star to be dynamically unstable. The linear instability of an oscillating dynamical system occurs on a dynamical timescale

[TABLE]

where $\hat{\rho}$ is the mean density of a star. [64] $\tau_{d}$ is often referred as the hydrostatic timescale, considering that it represents the typical time after which a slightly perturbed star can return to equilibrium. The critical value for the adiabatic exponent, $\Gamma_{cr}=4/3$, is a consequence of the gravitational force varying as $r^{-2}$. If the gravitational acceleration was given as $-GM/r^{2+\alpha}$, then the critical value would be $\Gamma_{cr}=(4+\alpha)/3$. One of the effects of GR is to make gravity stronger than Newtonian near a black hole, but the effect is there even when the field is weak. For $r\gg r_{s}=2GM/c^{2}$ we have roughly $\alpha\approx r_{s}/r$. This means that dynamical instability may occur when $\Gamma$ is just very close to, but somewhat larger than $4/3$.

Acknowledgements.

This work was supported by the National Research, Development and Innovation Office of Hungary (NKFIH, by Hungarian abbreviation) via grants NKFIH K–124366 and K–124508. I also acknowledge the partial support of networking activities by the “PHAROS” EU COST Action CA16214. I wish to express my gratitude and appreciation to Prof. Kostas Kokkotas, director of Institute for Astronomy and Astrophysics at the University of Tübingen (IAAT), for his helpful discussions and expert guidance in dealing with numerical issues in computing accurate eigenfrequencies of the fundamental and higher-order radial modes of neutron-star oscillations. I would also like to extend my gratitude to my former doctoral advisor, Dr. Mátyás Vasúth, head of Gravitational Physics Research Group at the Wigner Research Centre for Physics, for his critical comments and for reviewing some of the most important calculations.