Effect of viscosity and thermal conductivity on the radial oscillation and relaxation of relativistic stars

TL;DR

This paper develops a linearized model for small radial oscillations in relativistic stars, incorporating viscosity and thermal conductivity effects, to analyze eigenfrequencies and relaxation times relevant for neutron star pulsations.

Contribution

It introduces a generic formulation of the dynamical equations including viscosity and thermal conductivity effects, and applies a variational principle to determine eigenfunctions and eigenvalues for neutron star oscillations.

Findings

Eigenvalues give squared frequencies and relaxation times.

Framework allows integration with various nuclear matter models.

Qualitative estimates of energy dissipation rates are provided.

Abstract

In this paper we present a generic formulation of the linearized dynamical equations governing small adiabatic radial oscillations of relativistic stars. The dynamical equations are derived by taking into consideration those effects of viscosity and thermal conductivity of neutron-star matter which directly determine the minimum period of observable pulsars. A variational principle is applied to determine a discrete set of eigenfunctions with complex eigenvalues. The real and imaginary parts of eigenvalues represent the squared natural frequencies and relaxation time of radial oscillations of non-rotating neutron stars, respectively. We provide a suitable framework which may be supplemented with various potential species of cold-nuclear-matter models to compute the spectra of the normalized eigenfrequencies with a certain numerical precision. In the last section, we provide a qualitative estimation of the rate at which viscosity and thermal conductivity drain the kinetic energy of radial oscillation mode in reasonably uniform neutron stars, without relying on explicit numerical computations.

Equations

$$d{s}^2=e^{\nu }d{t}^2-e^{\lambda }d{r}^2-r^2(d{\vartheta }^2+{\sin }^2\vartheta d{\phi }^2)$$

$${\nabla }_{\mu }{J}^{\mu }=0,$$

$${\nabla }_{\mu }{T}^{\mu \nu }=0,$$

$${\nabla }_{\mu }{S}^{\mu }\ge 0$$

$$J^{\mu }=J_{\text{PF}}^{\mu }+J_{\text{NPF}}^{\mu },T^{\mu \nu }=T_{\text{PF}}^{\mu \nu }+T_{\text{NPF}}^{\mu \nu },$$

$$J_{\mu }^{\text{PF}}=\rho u_{\mu },T_{\mu \nu }^{\text{PF}}=(ϵ+p)u_{\mu }{u}_{\nu }-p{g}_{\mu \nu },$$

$$u_{\mu }=(e^{\nu /2},0,0,0).$$

$$ϵ=\rho (1+\epsilon )$$

$$h=\frac{p+ϵ}{\rho }=1+\epsilon +\frac{p}{\rho }.$$

$$h=1+\epsilon +\frac{p}{\rho }\to 1.$$

$$S^{\mu }=s\rho u^{\mu }+q^{\mu }/T,$$

$$T=\frac{1}{\rho }{\left(\frac{\partial ϵ}{\partial s}\right)}_{\rho }$$

$$q^{\mu }=(0,q^i)$$

$$T_{\text{NPF}}^{\mu \nu }={\mathcal{S}}^{\mu \nu }+T_{\text{flux}}^{\mu \nu },$$

$${\mathcal{S}}^{\mu \nu }={\pi }^{\mu \nu }+\Pi h^{\mu \nu },$$

$$h_{\mu \nu }=g_{\mu \nu }-u_{\mu }{u}_{\nu },$$

$$T_{\text{flux}}^{\mu \nu }=u^{(\mu }{q}^{\nu )}$$

$$T^{\mu \nu }=(ϵ+p)u^{\mu }{u}^{\nu }+(p+\Pi )h^{\mu \nu }+{\pi }^{\mu \nu }+u^{(\mu }{q}^{\nu )}.$$

$$\begin{array}[]{l}\pi_{\mu\nu}=-2\eta\sigma_{\mu\nu},\\[10.0pt] \Pi=-\zeta\Theta,\\[10.0pt] q_{\mu}=-\kappa\left(h^{\nu}_{\;\mu}T_{,\nu}-Tu_{\mu;\nu}u^{\nu}\right),\end{array}$$

$${\sigma }_{\mu \nu }=h_{\mu }^{\alpha }{h}_{\nu }^{\beta }{u}_{(\alpha ;\beta )}-\frac{2}{3}\Theta h_{\mu \nu }$$

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

$$T^{\mu \nu }=(ϵ+p)u^{\mu }{u}^{\nu }+p{g}^{\mu \nu }-2\eta {\sigma }^{\mu \nu }-\zeta \Theta h^{\mu \nu }+u^{(\mu }{q}^{\nu )},$$

$$\begin{array}[]{l}\displaystyle\bar{\epsilon}=\epsilon+T_{\text{NPF}}{}^{0}_{\;0}\\[10.0pt] \displaystyle\bar{p}=p-T_{\text{NPF}}{}^{1}_{\;1}\end{array}$$

$$G_{\mu }^{\nu }=8\pi T_{\mu }^{\nu },$$

$$8\pi \overline{ϵ}=e^{-\lambda }\left(\frac{{\lambda }^{\prime }}{r}-\frac{1}{r^2}\right)+\frac{1}{r^2}$$

$$8\pi \bar{p}=e^{-\lambda }\left(\frac{{\nu }^{\prime }}{r}+\frac{1}{r^2}\right)-\frac{1}{r^2}$$

$$8\pi \tilde{p}=e^{-\lambda }\left(\frac{{\nu }^{\prime \prime }}{2}-\frac{{\lambda }^{\prime }{\nu }^{\prime }}{4}+\frac{{\nu }^{\prime 2}}{4}+\frac{{\nu }^{\prime }-{\lambda }^{\prime }}{2r}\right)+\frac{e^{-\nu }}{4}\left(2\ddot{\lambda }+\dot{\lambda }(\dot{\lambda }-\dot{\nu })\right)$$

$$8\pi T_0^1=e^{-\lambda }\dot{\lambda }/r,$$

$$m^{\prime }=4\pi r^2\rho ,$$

$$\nu ={\nu }_0+{\int }_0^R\frac{2(4\pi r^{3}p+m)}{r(r-2m)}dr.$$

$$m(R)=4\pi {\int }_0^R\rho (r)r^{2}dr$$

$$d{V}_0=\sqrt{\det ({\gamma }_{ij})}d{x}^3=4\pi e^{\lambda /2}{r}^{2}dr,$$

Full text

• August 1, 2019

1 Introduction

1.1 Astrophysical motivation and relation to earlier works

Neutron star oscillations may impact on a range of observations, involving in particular radio and X-ray timing and gravitational waves. The possibility that radial oscillations of neutron stars give rise to oscillations observed within radio subpulses of pulsars was proposed by V. Boriakoff [1] in 1976. X-ray and $\gamma$-ray burst events have been generally associated with neutron stars by many authors e.g. [2]. Periodicities have been observed in X-ray bursts, which has raised considerable interest in radial [3] and other types [4] of neutron-star oscillations since the early 1980s. Bursters exhibit periodic and rapid increases in luminosity (typically a factor of 10 or greater) when enormous amount of energy deposited in oscillation modes is released in a short period of time through heat outflow via neutrino emission. [5] Physically, some mass from the stellar interior is drawn toward the surface where the hydrogen fuses to helium which accumulates until it fuses in a burst, producing X-rays. After the subsequent emission of thermal photon from the surface, the crust thermally relaxes toward equilibrium with the core. X-ray observations from the recently lunched NICER mission [6] and from the upcoming LOFT mission [7] will yield the mass and radius of a few stars up to $\sim 5\%$ precision. The observations of binary NS and BHNS coalescence events by gravitational-wave detectors such as the transient signal GW170817 will also dramatically improve our understanding of ultra-dense matter in neutron stars. Further interest in the study has been stimulated by [8], where $\gamma$-ray bursts were assumed to originate as results of collisions between strange stars at cosmological distances. Even though radial oscillations of strange stars are expected to be damped rapidly [9], such strange stars during their short time-scales are likely to be promising targets for multi-messenger observations. For most stellar models, the periods (typically ranging from 0.2 to about 0.9 milliseconds) depend on the stellar model and its central density [10], while the relaxation time is in the range of $0.1-0.3$ seconds. Therefore the study of oscillation spectra and relaxation times of neutron stars could be very valuable as their dynamical behaviour may eventually be deduced from observations.

Theoretical interest in the dynamical stability of relativistic stars has arisen since 1964 from the seminal works of S. Chandrasekhar, R.F. Tooper and J.M. Bardeen [11, 12, 13, 14, 15] and a general stability criterion was formalized in the 1970s by J. L. Friedman and B. F. Schutz [16]. The stability of spherically symmetric stars under radial adiabatic perturbations has been extensively studied and reported in the literature (e.g [3, 4, 17, 18, 19]). Induced by the small radial perturbations, dynamical instability will intervene by radial oscillation before the star contracts. Several techniques for obtaining spectra of oscillation modes have been developed for various stellar equilibrium models, mostly with zero-temperature equation-of-state (EOS). We have applied Chandrasekhar’s linear varational method [11, 12] to formulate the variational principle which forms the basis for determining the characteristic eigenfrequencies and relaxation times of radial oscillations.

Concurrently, the conversion of kinetic energy into heat and effects of viscosity on stellar pulsations in general has been addressed first by [20, 21, 22]. Properties of transport coefficients (bulk viscosity, shear viscosity, thermal conductivity) in neutron stars have been studied more in detail by a number of recent of works [23, 24, 25, 26]. The density and temperature dependence of shear viscosity and of bulk viscosity in the crust and in the core, respectively, have been described for different EOS of neutron stars by [23]. Thermal conductivity and shear viscosity of nuclear matter arising from nucleon-nucleon interaction in non-superfluid neutron-star cores were considered by [24], whereas those arising do from the collisions among phonons in superfluid neutron stars were considered by [26]. An extension of [24] for different nucleon-nucleon potentials and different three-body forces in [25] found that the nucleon contribution dominates the thermal conductivity, but the shear viscosity is dominated by leptons. The most up-to-date review of the transport properties and the underlying reaction rates of dense hadronic and quark matter in the crust and the core of neutron stars is found in [27].

1.2 Objectives and content

The present paper, together with its accompanying paper [28], is part of a more extensive body of work investigating the radial oscillations of neutron stars affected by the viscosity and thermal conductivity of cold nuclear matter that brings about damping of stellar oscillations and directly determine a minimum period of observable pulsars. The separation of the analytical and numerical aspects of study into two distinct papers was felt desirable, partly in order to not to divert attention from the details of each respective aspect, partly for the inherently different approaches involved.

Here, in addition to providing a comprehensive review of the basic theory of relativistic stellar pulsations, we present an analytical formulation of the dynamical equations that governs the radial mode of linear adiabatic stellar oscillations through a perturbation scheme. We prove that, similarly to the non-dissipative case, the pulsation equation 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 (SLEVP) 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. However, the main novelty of this approach is the ability to directly relate the damping ratio of oscillations to the expressions $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, which stem from the viscous and heat-conductive contributions to the stress–energy tensor, without relying on explicit numerical computations. The scale of relaxation time, directly related to the damping ratio, is identical to approximate solution for the time-scale of energy dissipation given by [29] in an alternative way. The usefulness of our analytical approximation method is evidently restricted to providing qualitative and “order-of-magnitude” information about the dissipative time-scales in (92) rather than a precise one.

Conversely, the numerical solution of the eigenvalue problem is considered in the accompanying paper [28]. The SLEVP for the radial oscillation modes of stars 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 Kokkotas and Ruoff (2001) [30] approach, we compute the four lowest-frequency radial-oscillation modes of neutron stars constructed from various potential EOS of cold-nuclear-matter considered by Özel and Freire (2016) [31]. The algorithm yields zero-frequency modes at the maxima and minima of the mass curves while the equilibrium adiabatic index characterizes the stiffness of the EOS at a given density. Finally, we evaluate the rate at which viscosity and thermal conductivity drain energy from the radial oscillation mode.

The present paper is organized as follows: Mostly based on [32, p. 285–289], sec. 2 gives a brief relativistic description of perfect fluids and of non-perfect fluids (viscous and heat-conductive) in Eckart frame (or ’particle frame’). The Einstein field equations and the generalized Tolman–Oppenheimer–Volkoff equation are expressed in sec. 3 through effective variables that incorporate time-dependent dissipative contributions of the neutron-star matter. In sec. 4, we present a variational method to formulate a lowest-order asymptotic approximation of infinitesimal adiabatic radial oscillations. We adopt a linear perturbation theory from [11] for the Lagrangian displacement of fluid elements from the equilibrium state where the adiabatic perturbations are considered to be purely radial. Stellar equilibrium models are described in sec. 5 through a local adiabatic index which can be regarded as constant near the center, but in general, varies, depending on the dynamical regime. Sec. 6 discusses the second-order linear ordinary differential equation of radial pulsations (78), derived from the perturbation equations. Sec. 6.1 deals with this pulsation equation that imposes a regular SLEVP for natural frequencies of oscillation with separated boundary conditions (83) induced by evaluating and setting $\chi$ to zero at the boundary. The characteristic time-scale for the relaxation of radial oscillation mode due to viscous dissipation is discussed in sec. 6.2, based on [29]. Lastly, the principal results established in this paper are summarized in sec. 7.

2 Stress–energy tensor of perfect and of non-perfect fluids

We may first consider a finite spherical distribution of fluids bounded by a sperical surface $\Sigma$ and held together by its own gravitational attraction. The line element for such a system is given in Schwarzschild-like coordinates by

[TABLE]

with the metric potentials $\nu=\nu(t,r)$ and $\lambda=\lambda(t,r)$ being dependent on the temporal and radial coordinates $(t,\,r)$. The general-relativistic hydrodynamic equations for a generic fluid involve the equations of motion that are given by the conservation of rest mass and by the conservation of energy and momentum:

[TABLE]

[TABLE]

respectively, and the second law of thermodynamics

[TABLE]

that are relativistically consistent. Although the perfect-fluid approximation disregards scenarios when dissipation and energy fluxes are present, it works well for most fluids under generic conditions. However, it loses its validity when thermodynamic (i.e microscopic) time-scales are comparable to the dynamic (i.e. macroscopic) ones and thus when assumption of local thermodynamic equilibrium breaks down. The requisite extension of perfect-fluid description that accounts for dissipative terms and energy fluxes is non-perfect fluid. In general, one can assume the rest-mass density current and stress–energy tensor as the linear combination of two contributions:

[TABLE]

where the indices ”PF” and ”NPF” refer to the perfect and non-perfect fluid contributions, respectively.

2.1 Perfect fluids

For a system which consists of perfect fluid with total-energy density $\epsilon$, isotropic pressure $p$ and covariant metric elements $g_{\mu\nu}$ corresponding to the antecedent line element (1), one shall have

[TABLE]

where the spatial components of fluid four-velocity $u_{\mu}$ are zeros. Normalized to $u^{\mu}u_{\mu}=1$, it becomes

[TABLE]

By construction, the quantity $\epsilon$ introduced above in eq. (6) represents the total-energy density of the fluid, given by

[TABLE]

which consists of both the rest-mass density of the fluid $\rho$ and the specific internal-energy density $\varepsilon$, internal-energy density per unit rest mass or which in this case represents the thermal motion of the constituent fluid particles. [32, p. 98] Finally, we define the specific enthalpy $h$ as

[TABLE]

Now, recognizing that in a non-relativistic regime $\varepsilon\ll c^{2}$ (i.e., the energy density of the fluid is essentially given by the rest-mass density) and $p/\rho\ll c^{2}$ (i.e., the pressure contribution to the energy density is negligible), the Newtonian limit of the specific enthalpy is given by

[TABLE]

Note that there are two natural ways to define four-velocity $u^{\mu}$. One option, given by Eckart, uses a unit timelike vector $\textbf{{u}}_{\text{N}}$ parallel to J, whereas the other, suggested by Landau, defines a unit timelike vector $\textbf{{u}}_{\text{E}}$ parallel to $\textbf{{T}}\cdot\textbf{{u}}_{\text{E}}$. However, these two vectors are identical for a perfect fluid and parallel to entropy current S. The (maximum) entropy principle (4) implies a strict equality for perfect fluids, whose entropy current is then given simply as $S^{\mu}=s\rho u^{\mu}$. However, for relation (4) to be strictly non-zero, the entropy current must have an additional contribution from dissipative parts (first-order theories) with non-zero divergence such that

[TABLE]

where the temperature $T$, deduced from the first law of thermodynamics is given by

[TABLE]

and

[TABLE]

are the components of the heat-flux four-vector that describes the rate of energy flow per unit area along each spatial coordinate axis within the Eckart frame. We will adopt the Eckart frame in which the continuity equation retains the same equation (3) as for a perfect fluid and there is no dissipative contribution to the the rest-mass density current and to the energy density.

2.2 Non-perfect fluids

The dissipative contributions $T^{\mu\nu}_{\text{NPF}}$ to the stress–energy tensor in eq. (5) can be further decomposed into the form

[TABLE]

where $\mathcal{S}^{\mu\nu}$ is referred to as viscous stress tensor. An appropriate definition that avoids the possibility of superluminal propagation of speed is given by

[TABLE]

where $\pi^{\mu\nu}$ is the anisotropic stress tensor, $\Pi=p-p_{\text{eq}}$ is the viscous bulk pressure that measures the deviation of pressure from its equilibrium value, and

[TABLE]

is the standard projection tensor onto 3-space normal to flow; all tensors are symmetric. The heat-flow tensor

[TABLE]

accounts for the generation of energy fluxes. Altogether, $\mathcal{S}^{\mu\nu},\,\pi^{\mu\nu}$ and $\Pi$ are known as the thermodynamic fluxes and they account for the deviations of the fluid from a perfect fluid. Finally, the full stress–energy tensor will take the form

[TABLE]

All the different equations for the thermodynamical fluxes $(\pi_{\mu\nu},\,\Pi,\,q_{\mu})$ reduce under the condition [29] that

the dissipation coefficients are sufficiently small and 2. 2.

the time derivates of perturbed quantities remain sufficiently small

or in other words: provided that $T\nabla_{\mu}S>0$, i.e. the entropy four-divergence satisfy the condition of being non-negative [32], the relations of thermodynamical fluxes to thermodynamic forces $(\Theta,\,h^{\nu}_{\;\mu}\ln T_{,\nu}+u_{\mu;\nu}u^{\nu},\,\sigma_{\mu\nu})$ become linear:

[TABLE]

where the symmetric, trace-free, spatial shear tensor $\sigma_{\mu\nu}$ is defined by

[TABLE]

and expansion scalar (or dilatation rate)

[TABLE]

is associated with the divergence or convergence of the fluid world lines. [33] $\eta$ is the shear (also ’common’ or ’dynamic’) viscosity coefficient that describes the fluid’s resistance to gradual shear deformation; $\zeta$ is the bulk (or ’second’) viscosity coefficient that defines the resistance of the medium to gradual uniform compression or expansion; and $\kappa$ is non-negative and accounts for the thermal conductivity, respectively. Collectively, $\kappa,\,\zeta$, and $\eta$ are called transport coefficients or dissipation coefficients, while eqs. (19) are referred to as the constitutive equations of classical irreversible thermodynamics, also known as Eckart’s theory of relativistic irreversible thermodynamics [34]. Eqs. (19–21) give rise to an alternative form of the full stress–energy tensor,

[TABLE]

that is expressed through thermodynamic forces and transport coefficients. The introduction of the transport coefficients provides us an opportunity to shed light on the different definitions of fluids that are possible in combination with the different values assumed by $\eta,\,\zeta,\,\kappa$. [32]

3 Field equations, effective variables and the post-quasistatic approximation

It is convenient to replace the energy density $\epsilon$ and isotropic pressure $p$ with effective variables

[TABLE]

so that Einstein’s field equations

[TABLE]

upon satisfying the metric (1), provide five partial differential equations for each of the non-vanishing mixed-variance components of the stress–energy tensor, of which the four distinct ones are

[TABLE]

where $\tilde{p}=\bar{p}-T_{\text{NPF}}{}^{1}_{\;1}+T_{\text{NPF}}{}^{2}_{\;2}$ and overdots and primes denote partial differentiation with respect to $t$ and $r$, respectively. Another equation is proportional to (25c), thus it is needless to consider it separately. The effective variables $(\bar{\rho},\,\bar{p})$ satisfy the same Einstein’s field equations in the quasi-static regime as the corresponding physical variables $(\rho,\,p)$ (take account of the contribution of $T^{\mu\nu}_{\text{NPF}}$ to $p$ and $\epsilon$ in eq. (8)). Therefore, the effective and physical variables share the same radial dependence. [35] Owing to the fact that (25a) involves only $\lambda$ and $\rho$, it becomes

[TABLE]

once the new radial-dependent variable $m(r)=rt(1-e^{-\lambda})/2$ had been introduced. Similarly,

[TABLE]

Suppose that the radius extends to $R$, from (26) it is evident that the integral of the effective rest-mass density over the stellar interior

[TABLE]

can be interpreted as the ’gravitational mass’ of the system which includes all contributions to the relativistic mass (rest mass, internal energy, and the negative gravitational binding energy). However, integrating the total energy-density over the proper volume

[TABLE]

where the curvature of 3-space has been taken account of through $\lambda$, one obtains the mass

[TABLE]

that represents the sum of rest mass and internal energy. The difference between the two arises as a result of the mutual attraction of the fluid elements, called the binding energy, which is given by

[TABLE]

and exhibits the amount of energy required to disassemble a whole system into separate elements.

It is convenient not to use eq. (25c) directly, but instead substracting it from (25b). It can be replaced by the first-order ODE

[TABLE]

which stands for the condition for hydrostatic equilibrium, provided that the right-hand side of the equation is zero (cf. (46)). The first term on the right-hand side recovered from eq. (18) is simply the result of

[TABLE]

recovered from eq. (18). After the elimination of $\nu^{\prime}$, one obtains the generalized Tolman–Oppenheimer–Volkoff (TOV) equation

[TABLE]

for non-perfect fluids. Let us note that the non-zero terms on the left-hand side of eq. (34) are not contradictory to conservation of energy and momentum in (2–3). It will be shown that a set of two equations (46) takes over its role; one describing the hydrostatic equilibrium and another describing the perturbation-induced departure from the equilibrium configuration. In order to close the system of the field equations (25) and the TOV eq. (34), it must be supplemented with an EOS relating some fundamental thermodynamical quantities. In general, the EOS takes the form

[TABLE]

where $\bar{\epsilon}$ and $\bar{p}$ are the total-energy density and isotropic pressure, given by eq. (23).

4 Infinitesimal radial oscillations

Let us suppose an equilibrium configuration of non-perfect fluid governed by the eqs. (25) is subject to a small linear perturbation that does not violate its spherical symmetry. Let $\delta$ hereafter denote a small ratio between the scale of variation of the perturbed variables and the correspondig ones in unperturbed configuration. Any quantity associated with the unperturbed equilibrium state is denoted by the subscript ”0” whereas those that represent perturbations are equipped with the subscript ”1”. In respect to such a perturbation, motions in the radial directions arise. While formulating the equations governing the perturbed state, we shall ignore all quantities which are of the second or higher orders in motions. Consequently, the four-velocity of a fluid element will be expressed as

[TABLE]

where

[TABLE]

is the radial velocity with respect to the time coordinate $t$.

4.1 Perturbation equations for stellar oscillations

One way to describe perturbations is the ”microscopic” point of view where the observer follows individual fluid particles as they move through space and time. [36] It is convenient to introduce a displacement field $\xi_{1}$ (in our case, a scalar field) in the Lagrangian representation defined by

[TABLE]

which connects fluid elements in the equilibrium with corresponding ones in the perturbed configration. There is yet another, ”macroscopic” way of looking at the perturbations. [37] In the Eulerian representation of fluid motion, we simply consider changes in the variables $(\lambda,\,\nu,\,\epsilon,\,p)$ at a fluid element fixed in space and time. [38] This means that

[TABLE]

where the linear perturbations $(\delta\lambda_{1},\,\delta\nu_{1},\delta\epsilon_{1},\,\delta p_{1})$ are Euler changes. In this manner, retaining only those terms of the full stress–energy tensor (22) that do not involve second or higher orders in $\delta$. The PF and NPF parts of the tensor assume the form of

[TABLE]

respectively. Notice that in eq. (40a), the off-diagonal elements are the only elements where the leading-order terms are absent, whereas in eq. (40b) the leading-order terms are present only in these very same elements (see Appendix). In view of the above, field equations (25a) and (25b) hold true if the pair of static metric potentials $(\lambda_{0},\,\nu_{0})$ is replaced by $(\lambda,\,\nu)$. The corresponding linearized equations governing the perturbations are, thus,

[TABLE]

and

[TABLE]

The appropriately linearized form of eqs. (25d) and (3), respectively,

[TABLE]

and

[TABLE]

shall serve in place of the remaining two field equations. At this point we may recall the definitions (23) for the effective variables $(\bar{\epsilon},\,\bar{p})$ and introduce their respective linearized forms $(\bar{\epsilon}_{0}+\delta\bar{\epsilon}_{1},\,\bar{p}_{0}+\delta\bar{p}_{1})$. Eqs. (39–40) enable us to identify these new variables easily:

[TABLE]

Consequently, expressed in these more suitable variables, eq. (44) decomposes into a set of two equations

[TABLE]

where the first equation assumes the role of the constraint (34) for hydrostatic equilibrium while the second one governs the dynamics of perturbation-induced departure from the equilibrium configuration. As will be shown later on in sec. 6, fluctuations of the stellar radius exhibit oscillatory behaviour: expansions and contractions in the outer layers as a star pursues to maintain equilibrium. The remaining NPF contributions are gathered on the right-hand side to form the source term

[TABLE]

which stems from frictional forces in the fluid itself. As such, it is responsible for exponential growth or damping, depending on the friction coefficient.

Now, in the sense of the definition (38), the integration of eq. (43) yields

[TABLE]

or in consideration of the first eq. of (46),

[TABLE]

Eqs. (41) and (48) provide

[TABLE]

Substituting for $d\bar{p}_{0}/dr$ in the last equation from the first eq. of (46), we can also write

[TABLE]

Considering next equation (42) and substituting for $\delta\lambda_{1}$ in accordance with eq. (48), we obtain

[TABLE]

or in view of the first eq. of (46),

[TABLE]

4.2 Damping of stellar oscillations

Suppose that for normal modes of the fluid perturbations $(\delta\lambda_{1},\,\delta\nu_{1},\,\delta\bar{\epsilon}_{1},\,\delta\bar{p}_{1})$ possess a harmonic time-dependence of the form $\exp(i\Omega t)$ where

[TABLE]

is a complex characteristic frequency to be determined in sec. 6. Being subject to the damping effect of the dissipative forces (19), the fluid perturbations oscillate with a damped angular frequency (sometimes called pseudo-frequency)

[TABLE]

related to $\omega_{n}$ which is the natural frequency (or resonant frequency) of the undamped system. The rate at which the normal modes of radial oscillations are damped is characterized by the relaxation time (or damping time)

[TABLE]

and can be determined from energy-dissipation eq. (90). [39] Expressed by the damping ratio $\zeta$, a dimensionless non-negative parameter, the complex characteristic frequency is

[TABLE]

The value of $\zeta$ prescribes the frequency response and critically determines the dynamical behaviour. The oscillation is undamped $(\zeta=0)$, if the star oscillates with its natural angular frequency $\omega_{n}$. It is underdamped $(\zeta<1)$, if the star oscillates with a damped frequency $\omega_{d}$ and with the amplitude gradually decreasing with the rate of decay $1/\tau$. It is critically damped $(\zeta=1)$, if the star returns to steady state as quickly as possible without any oscillation.

With the harmonic time-dependence, we can rewrite the second equation of (46) in the form

[TABLE]

where we have substituted for $(\bar{p}_{0}+\bar{\epsilon}_{0})\nu^{\prime}_{1}$ in accordance with eq. (53). Also, it should be recalled that $\bar{\epsilon}_{1}$ is expressed in terms of $\xi_{1}$ and the perturbed variables by eq. (50). The additional term $\mathcal{S}_{1}$ is given by the expression (A.99).

5 The conservation of baryon number

The continuity equation (2) of the rest-mass density current given as in (6) involves the conservation of the baryon number, represented by the condition

[TABLE]

provided that $N$ is the number of constituent baryons per unit volume. Let the covariant derivative be written differently, in the form

[TABLE]

with the expression $\partial(\log\sqrt{-g})/\partial x^{\mu}$ being put in the place of the Christoffel symbols corresponding to the metric given in eq. (1) and

[TABLE]

being the determinant of the metric tensor. In the framework of the present linearized theory, it is quite reasonable to define the baryon number by

[TABLE]

so that eq. (60) endowed with the non-vanishing components of the four-velocity (36) emerges as

[TABLE]

With $v_{1}$ replaced by the Lagrangian displacement $\xi_{1}$ defined in eq. (38), the last equation brings about fluctuations in the baryon number:

[TABLE]

where the second term on the right-hand side vanishes on account of eq. (49) and

[TABLE]

represents the effect of dissipative terms on the baryon-number perturbation. Subsequently, eq. (64) reduced significantly and one obtains

[TABLE]

Provided that $N\equiv N(\epsilon,\,p)$ is an EOS that corresponds to (35), small linear perturbations in the energy density or pressure treated as variation are expected to induce baryon-number perturbations given by

[TABLE]

which in turn yields

[TABLE]

under the assumption that variables in the perturbed state relate to each other roughly the same way as the corresponding variables in equilibrium, that is $\partial\bar{N}_{1}/\partial\bar{\epsilon}_{1}\approx\partial\bar{N}_{0}/\partial\bar{\epsilon}_{0}$ and $\partial\bar{N}_{1}/\partial\bar{p}_{1}\approx\partial\bar{N}_{0}/\partial\bar{p}_{0}$. With $\delta\bar{\epsilon}_{1}$ and $\delta\bar{N}_{1}$ given by eqs. (50) and (66), respectively, eq. (68) comes to be

[TABLE]

where $\Gamma$ is identified with the adiabatic index given by

[TABLE]

Within a star that is dynamically stable against infinitesimal radial adiabatic perturbations, the average value for the adiabatic index is greater than or equal to $4/3$. Expressly,

[TABLE]

where $\epsilon$ denotes the total-energy density, $p$ the radial pressure, and $c_{s}$ the radial sound speed, respectively [40].

6 The pulsation equation and eigenvalue problem

With elimination of $\bar{\epsilon}_{1}$ and $\bar{p}_{1}$ from eq. (58), through eqs. (50) and (69), it emerges as

[TABLE]

where the third and forth terms on the right-hand side are the result of the transformation $\frac{\mathrm{d}f}{\mathrm{d}r}+f\frac{\mathrm{d}g}{\mathrm{d}r}=\exp(-g)\frac{\mathrm{d}}{\mathrm{d}r}\left[f\exp(g)\right]$, which applies to any two functions $(f,\,g)$ of the variable $r$ and $\mathcal{S}_{2}$ is written out in eq. (A.101) in appendix.

Substituting for $\bar{p}^{\prime}_{0}$ from the first eq. of (46) and applying the same transformation for the first two terms on the right-hand side of eq. (72) allow us to merge them with the next two terms in the expression

[TABLE]

which (cf. (56) in [11]), in turn, compared with the field equation (25c) restricted by equilibrium conditions,

[TABLE]

is reduced to three terms. Making use of the first eq. of (46) once again, the relation (72) appreciably reduces to give the pulsation equation

[TABLE]

associated to the class of second-order linear ODEs. Besides this, the definitions (54–57) make possible the separation of real and imaginary parts of the complex frequency squared:

[TABLE]

thus the last term of the right-hand side:

[TABLE]

For positive $\Omega^{2}$, the characteristic frequency $\Omega$ is real and thus, the solution is purely oscillatory. However, for $\Omega^{2}<0$, $\Omega$ contains an imaginary part, which corresponds to a damped solution. Since the general solution is a superposition of damped modes, the occurrence of a negative value of $\Omega^{2}$ corresponds to a secular instability whose growth time is long compared to the dynamical time of radial oscillations. The Harrison–Zel’dovich–Novikov criterion [17] for static stability of compact stars states that the total mass of such stars increases with the central density $\rho_{0}$ which implies that $dM(\rho_{0})/d\rho_{0}\geq 0$ for the stable region where $M(\rho_{0})$ is the function of total mass in terms of the central density. For neutron stars this will, indeed, happen for $\rho_{0}$ larger than the critical central density $\rho_{\text{crit}}$ at which the stellar mass $M(\rho_{0})$ as a function of $\rho_{0}$ has its maximum. In this case the star will ultimately collapse to a black hole. For $\rho_{0}=\rho_{\text{crit}}$ there must be a neutral mode with the corresponding eigenvalue $\omega^{2}=0$. [40]

6.1 Eigenfrequencies of radial pulsation

The real part of left-hand side of eq. (75) can be equated to the real part of right-hand side, comparably to the imaginary parts. Intrinsically, we can recast it in the so-called homogeneous Sturm–Liouville form

[TABLE]

with a free parameter

[TABLE]

denoting the eigenvalues and a set of coefficient functions

[TABLE]

hailing from eq. (75), are specified at the outset. The function $\mathcal{R}(r)$ is referred to as weighting function. The normalized Lagrangian displacement defined by

[TABLE]

is a scalar-valued function of the variables $(t,\,r)$. Provided that $\chi$ satisfies eq. (78), it is called a solution. Solutions of (78) are subject to the boundary conditions

[TABLE]

Together with the boundary condition (82), the pulsation equation (78) imposes a Sturm–Liouville eigenvalue problem (SLEVP), which seeks non-trivial solutions only for a countable set of real eigenvalues $\{\Lambda_{1},\Lambda_{2},\ldots,\Lambda_{n}\}$. The SLEVP is said to be regular if $\mathcal{P}>0$ and $\mathcal{R}>0$ for any $r\in[0,\,R]$, the functions $(\mathcal{P},\,\mathcal{P}^{\prime},\,\mathcal{Q},\,\mathcal{R},\,\mathcal{S})$ are continuous over the finite interval $[0,\,R]$, and the problem has separated boundary conditions of the form

[TABLE]

Sturm–Liouville theory states that the eigenvalues of the regular SLEVP are real and can be arranged in ascending order such that

[TABLE]

Corresponding to each eigenvalue $\Lambda_{n}$ is a unique (up to a normalization constant) eigenfunction $\chi_{n}(r)$ which has exactly $n-1$ zeros in $(0,R)$. Moreover, the normalized eigenfunctions form an orthonormal basis

[TABLE]

where $\delta_{mn}$ is the Kronecker delta and to each $\Lambda_{n}$ is associated with a single eigenfunction $\chi_{n}$.

6.2 Characteristic relaxation time of radial pulsation

In the Newtonian limit, the kinetic energy contained in these oscillations is given by

[TABLE]

an integralover an element of proper volume of fluid $dV_{0}$ (see eq. (29)), where $\delta v_{1}^{*}$ is the complex conjugate of the velocity perturbation $\delta v_{1}$. Associated with the radial displacement $\delta\xi_{1}=\delta r/r$, the later is given for radial oscillations by eq. (38) as

[TABLE]

where, in accordance with second expression of (A.95), we have substituted for $(\delta T^{\text{NPF}}_{1})^{1}_{\;0}$. The total energy in an oscillating star consists of kinetic and potential energy which are supposed to equally contribute to the total energy of harmonic oscillations, thus given by $E=2E_{k}$. On account of the density, $\rho$, is reasonably uniform in neutron stars [29], the total energy contained in the oscillation is given as

[TABLE]

by evaluating the integral (86) explicitly for an avaraged density $\bar{\rho}=3M/4\pi R^{3}$. Being bilinear in the fluid perturbations, $E$ has a time-dependence $\exp[-2\operatorname{Im}(\Omega)t]$. [39] Subsequently, its time derivative implies that

[TABLE]

which together with the energy-dissipation rate for the stress–energy tensor (22) as

[TABLE]

directly determines the dissipative time-scale of small perturbations of the fluid away from the equilibrium state as

[TABLE]

An approximate formula for each dissipative time-scales can be given by evaluating the corresponding dissipation integral (90). Cutler & Lindblom [29] have found that the following propotionalities hold for the dissipative time-scales:

[TABLE]

which reveals that the time-scale of shear viscosity is much shorter than either that of bulk viscosity or thermal conductivity. The imaginary part of the pulsation equation (75), together with eqs. (76–77) yields

[TABLE]

which implies that a unique damping ratio $\zeta_{n}$ corresponds to each eigenfrequency $\omega_{n}$ computed for undamped oscillations. From (93), it is evident that higher frequency components die out first. These damping ratios are identical with those associated with the dissipative time-scales (92) which, according with (56), are obtained from $\zeta_{n}=\tau_{n}\omega_{n}$.

7 Concluding remarks

A generic formulation of the dynamical equations governing small adiabatic radial oscillations of pulsating relativistic stars has been proposed in this paper through a perturbation scheme that, combined with the equations of viscous thermally-conductive fluids, constitutes an extension of radially pulsating perfect-fluid stellar models. We have proved that, similarly to the regular perfect-fluid case, the stellar pulsation equation (75) expressed by a set of effective variables (23) which involve dissipative terms, can be recast in a homogenous Sturm–Liouville form (78) with separated boundary conditions (83). In contrast to the regular perfect-fluid case, the associated eigenvalue problem is generalized for a discrete set of eigenfunctions with complex eigenvalues where the real and imaginary parts of the eigenvalues represent the squared natural frequency and relaxation time (or decay rate) of the oscillation, respectively. In the absence of dissipation, the discrete spectrum consists of real eigenvalues that form a complete set.

The main novelty of this approach is the ability to directly relate the damping ratio of oscillations to the expressions $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ in (A.100–A.101), which stem from the viscous and heat-conductive contributions to the stress–energy tensor, without relying on explicit numerical computations. An illustrative example set by [29] for neutron star with uniform density, allowed us in this paper to estimate the rate at which the viscosity and thermal conductivity of the nuclear matter drains energy from the oscillations. In accordance with the literature, the time-scale of shear viscosity is much shorter than either that of bulk viscosity or thermal conductivity and the imaginary part of the pulsation equation indicate that higher components vanish first from the frequency spectrum. The usefulness of our analytical approximation method is evidently restricted to providing qualitative and “order-of-magnitude” information about the dissipative time-scales in (92) rather than a precise one.

Appendix

We are to enumerate the non-vanishing components of the NPF stress–energy tensor $(T_{\text{NPF}})^{\mu}_{\;\nu}$ referred in eq. (40b). To explicitly evaluate them, the arithmetic operations in eqs. (20–22) have to performed. The non-vanishing components of the zeroth-order are

[TABLE]

those of the first-order are

[TABLE]

with an only radial-dependent coefficient

[TABLE]

In view of (A.95), it is evident that to abide by the nature of harmonic time-dependence, a relation of the form

[TABLE]

is implied. According to the definition (38), this implication requires

[TABLE]

to hold. In regard to to eq. (A.98), the expression $e^{\lambda_{0}-\nu_{0}}(\dot{T}^{\text{NPF}}_{1})^{1}_{\;0}$ is combined with $\mathcal{T}$ given by eq. (47) to compose

[TABLE]

which first appears in eq. (58). Substituted with the set of eqs. (A.94–A.98), the expression unfolds as

[TABLE]

Next to the $\mathcal{S}_{1}$, another quantity, expressed by

[TABLE]

appears in different forms the pulsation equation, namely in eqs. (75) and (72). Together, they constitute the expression $i(\Omega\mathcal{S}_{1}-\Omega^{-1}\mathcal{S}_{2})$ which represents the source of inhomogenity of the SL equation.

Acknowledgments

I would like to gratefully acknowledge the financial support from ”PHAROS”, COST Action CA16214. I am also grateful for K. Kokkotas for valuable discussions and for his generous hospitality at the University of Tübingen in April 2018.