Gravitational waves from r-mode oscillations of stochastically accreting neutron stars

ABSTRACT

1 INTRODUCTION

r-modes are inertial modes in rotating neutron stars, whose restoring force is the Coriolis force. Gravitational waves (GWs) generated by r-mode oscillations have attracted significant attention in continuous wave searches (Riles 2023; Wette 2023). This is because r-modes grow unstably through the emission of gravitational radiation (Andersson 1998; Friedman & Morsink 1998; Lindblom, Owen & Morsink 1998; Owen et al. 1998). Recent targeted and all-sky searches for GWs from r-modes report no detections and infer upper bounds on the characteristic wave strain satisfying |$h_{0}^{95 {{\ \rm per\ cent}}} \lesssim 10^{-25}$| at |$95 {{\ \rm per\ cent}}$| confidence, using data from the first three observing runs of the Laser Interferometer Gravitational-wave Observatory (LIGO), Virgo, and Kamioka Gravitational Wave Detector (KAGRA) (Fesik & Papa 2020; Middleton et al. 2020; Rajbhandari et al. 2021; Abbott et al. 2021a, b; Covas et al. 2022; LIGO Scientific Collaboration and KAGRA Collaboration 2022; LIGO Scientific Collaboration and Virgo Collaboration 2022; Vargas & Melatos 2023). Indirect upper limits on |$h_0^{95 {{\ \rm per\ cent}}}$| from r-modes can also be inferred from electromagnetic observations of neutron stars, assuming that the gravitational radiation reaction torque causes the star to spin down (Owen 2010; Glampedakis & Gualtieri 2018; Caride et al. 2019; Middleton et al. 2020) or balances the accretion torque (Wagoner 1984; Bildsten 1998; Glampedakis & Gualtieri 2018; Middleton et al. 2020), or that photon cooling is in thermal equilibrium with the viscous damping of r-modes (Schwenzer et al. 2017).

Unstable r-modes grow exponentially from small, random, initial fluctuations, until non-linear mode coupling saturates their growth (Schenk et al. 2001; Arras et al. 2003; Bondarescu, Teukolsky & Wasserman 2007). The initial fluctuations are sometimes attributed to non-specific origins (e.g. thermal) and sometimes modelled in detail. For example, buoyant r-modes may be excited in the surface ocean during thermonuclear type I X-ray bursts (Strohmayer & Lee 1996; Heyl 2004; Chambers & Watts 2020) and can lead to coherent X-ray oscillations during a superburst, as claimed in XTE J1751–305 (Strohmayer & Mahmoodifar 2014). The excited surface r-mode is not related directly to the r-modes in the bulk of the star, which themselves may be amplified near the surface (Lee 2014). In binary systems accreting at present, oscillation modes can also be excited by repeated mechanical impacts of clumps of stochastically accreted matter (Nagar et al. 2004; Dong & Melatos 2024). Previous work along these lines has focused on non-radial f-, p-, and g-modes in a non-rotating star with a polytropic equation of state (EOS). The root-mean-square GW strain increases with increasing polytropic index, attaining |$10^{-33} \le h_{\rm rms} \le 10^{-32}$| under astrophysically plausible conditions (Dong & Melatos 2024).

In this paper, we extend previous calculations to include rotation and investigate non-radial r-mode oscillations excited mechanically by stochastic accretion. We calculate the amplitude of the excited r-modes and the root-mean-square amplitude and Fourier spectrum of the resulting gravitational radiation, assuming that the star rotates slowly enough not to trigger the Chandrasekhar–Friedman–Schutz (CFS) instability. The paper is organized as follows. In Section 2, we introduce the stellar model and briefly review the linear perturbation theory of free r-mode oscillations. In Section 3, we calculate the impulsive r-mode response to the mechanical impact of accreted matter, first from a single clump and then from a stochastic sequence of clumps, as implied by modern numerical simulations (Romanova & Owocki 2015). We compute the GW signal from the excited r-modes as a function of the ensemble statistics of the accretion clumps in Section 4. We review the conditions for exciting the CFS (in)stability and discuss the astrophysical implications of the cyclically triggered CFS instability in Section 5. Section 6 is the conclusion.

2 r-MODE OSCILLATIONS IN A SLOWLY ROTATING STAR

The neutron star model in this paper is kept as simple as possible to focus on the new feature of the analysis: stochastic excitation of r-modes by the mechanical impacts of clumps of accreting matter. We model the star as a one-component ideal fluid rotating uniformly with angular velocity |$\boldsymbol {\Omega }$| in Newtonian gravity. We also assume that the EOS is barotropic. Equilibria and r-mode oscillations of this idealized model have been examined in detail in the literature (Papaloizou & Pringle 1978; Tassoul 1978; Saio 1982), as well as their generalized counterparts in the broader inertial modes family (Lindblom & Ipser 1999; Lockitch & Friedman 1999). We neglect the corrections from rapid rotation (Lindblom, Mendell & Owen 1999), shear viscosity at the crust–core boundary (Bildsten & Ushomirsky 1999; Levin & Ushomirsky 2001), stratification (Yoshida & Lee 2000; Passamonti et al. 2009; Andersson & Gittins 2023; Gittins & Andersson 2023), non-linear mode coupling (Schenk et al. 2001; Arras et al. 2003; Bondarescu et al. 2007), superfluidity (Lindblom & Mendell 2000; Andersson & Comer 2001), and magnetic fields (Rezzolla, Lamb & Shapiro 2000; Morsink & Rezania 2002). Although pure r-modes do not exist in general relativistic barotropes, relativistic corrections to the Newtonian r-modes have also been calculated using axial-led inertial modes in relativistic barotropic stars (Lockitch, Andersson & Friedman 2000; Lockitch, Friedman & Andersson 2003; Idrisy, Owen & Jones 2015; Ghosh, Pathak & Chatterjee 2023). For relativistic non-barotropic stars, however, the eigenvalue problem of pure r-modes is singular (Kojima 1998). These and other realistic corrections can be incorporated into the theoretical framework Sections 2 and 3, when the need arises. They are essential for a detailed analysis of GW observations, once detections are made.

In this section, we review the linear perturbation theory of r-mode oscillations in a slowly rotating star. In Sections 2.1 and 2.2, we set out the equilibrium configuration and the linearized equations of motion, respectively. The r-mode oscillations are assumed to be adiabatic. We introduce the slow-rotation expansion and derive the eigenfrequencies and eigenfunctions for r-modes in Section 2.3. Orthogonality of the eigenmodes is discussed in Section 2.4.

2.1 Equilibrium

The mass continuity, Euler, and Poisson equations in the corotating frame are

where |$\boldsymbol {u} = \boldsymbol {v}-\boldsymbol {\Omega } \times \boldsymbol {x}$| is the flow velocity in the corotating frame, |$\boldsymbol {v}$| is the flow velocity in the inertial frame, |$\boldsymbol {x}$| is the position vector, |$\rho$| is the density, P is the pressure, and |$\Phi$| is the gravitational potential.

At equilibrium, the star obeys the unperturbed stationary Euler equation with |$\boldsymbol {u}_0 = 0$|⁠, viz.

where |$\Phi _{\rm rot} = - (\boldsymbol {\Omega }\times \boldsymbol {x})^2/2$| is the centrifugal potential. The subscript zero denotes an equilibrium quantity. For slowly rotating stars satisfying |$|\boldsymbol {\Omega }| \ll (\pi G \rho _0)^{1/2}$|⁠, the centrifugal potential is negligible compared to the gravitational potential, and the background star is approximately spherically symmetric.

2.2 Linearized equations of motion

We describe the linear perturbations of the star by Eulerian perturbed quantities, which are denoted by the operator prefix |$\delta$|⁠, e.g. the pressure perturbation is |$\delta P$|⁠. Lagrangian perturbations, denoted by the prefix |$\Delta$|⁠, are related to Eulerian perturbations via |$\boldsymbol {\xi }(\boldsymbol {x}, t)$|⁠, the Lagrangian displacement of a fluid element from its equilibrium position. Mathematically, we write

With |$\boldsymbol {u}_0 = 0$| and |$\delta \boldsymbol {u} = \Delta \boldsymbol {u}$|⁠, we obtain

The linearized forms of equations (1)–(3), expressed in terms of Eulerian perturbations, are then given by

where |$\bf{B}$| has Cartesian components

|$\epsilon _{ijk}$| is the Levi-Civita symbol, the Einstein summation convention applies to repeated indices, and |$\bf{C}$| satisfies (Lynden-Bell & Ostriker 1967; Schenk et al. 2001)

Note that the centrifugal term in equation (2) does not appear in equation (8), because |$\Delta [\boldsymbol {\Omega } \times (\boldsymbol {\Omega } \times \boldsymbol {x})]$| cancels with |$(\boldsymbol {\xi } \cdot \nabla)(-\nabla P_0/\rho _0-\nabla \Phi _0)$| on the right-hand side of equation (2).

We assume that the perturbations are adiabatic, satisfying

where |$\Gamma _1$| is the adiabatic index. In general, |$\Gamma _1 = \Gamma _1(\boldsymbol {x})$| is a function of the background structure of the star. Equation (8) together with equation (12) reduces to equation (14) in Dong & Melatos (2024) for |$\boldsymbol {\Omega }=0$|⁠. With equation (5), equation (12) can also be written as

under the assumption that the Schwarzschild discriminant |$\boldsymbol {A}$| satisfies

equation (15) holds for barotropic stars, where P is a function of |$\rho$| only [|$P = P(\rho)$|]. In other words, a barotrope has no buoyancy, because the squared Brunt-Väisälä frequency |$N^2 \propto \boldsymbol {A}\cdot \nabla P_0$| vanishes. The barotropic approximation is valid, when the Coriolis force dominates buoyancy, with |$N \lesssim |\boldsymbol {\Omega }|$|⁠.

The Brunt-Väisälä frequency can be estimated by |$N \sim (x/2)^{1/2} g c_{\text{s}}^{-1}$|⁠, where |$x = 6 \times 10^{-3} \rho /\rho _{\text{nuc}}$| is the local ratio of protons to neutrons in chemical equilibrium, |$\rho _{\text{nuc}} = 2.7 \times 10^{14} \, \text{g cm}^{-3}$| is the nuclear saturation density, g is the local gravitational acceleration, and |$c_{\text{s}}$| is the adiabatic sound speed (Reisenegger & Goldreich 1992; Lockitch et al. 2003). In a recent study, Altiparmak, Ecker & Rezzolla (2022) reported |$c_{\text{s}}^2/c^2 \gt 1/3$| at |$\rho \gtrsim 2\rho _{\text{nuc}}$| for most of the EOSs they consider, which parametrize |$c_{\text{s}}$| as a function of the chemical potential in the range |$1.1 \rho _{\text{nuc}} \lt \rho \lesssim 40 \rho _{\text{nuc}}$|⁠. This implies |$150\, \text{Hz} \lesssim N \, (\rho /\rho _{\text{nuc}})^{-1/2} (g/10^{14} \text{cm s}^{-2})^{-1} \lesssim 260$| Hz. Therefore, the barotropic approximation |$N \lesssim |\boldsymbol {\Omega }|$| is valid at the margin for accreting millisecond X-ray pulsars with |$|\boldsymbol {\Omega }|/2\pi \gtrsim 200$| Hz (Patruno & Watts 2021). This is consistent with the findings of Passamonti et al. (2009), who demonstrated numerically that the difference between |$l=|m|$|r-mode frequencies in barotropic and non-barotropic stars is negligible. In this paper, we assume that the barotropic condition is fulfilled for r-modes for simplicity. A more realistic study of the general problem in non-barotropic stars is left for future work.

2.3 Slow-rotation expansion of eigenmodes

The normal modes of oscillation |$\boldsymbol {\xi }_{\alpha }(\boldsymbol {x}, t)$| of equation (8), summing to give the total displacement vector |$\boldsymbol {\xi } = \sum _{\alpha } \boldsymbol {\xi }_{\alpha }$|⁠, exhibit a harmonic time dependence,

because operators |$\bf{B}$| and |$\bf{C}$| are time-independent. The subscript |$\alpha$| labels the eigenmodes. The eigenvalue |$\sigma _{\alpha } = \omega _{\alpha } + i/\tau _{\alpha }$| is complex in general, and its imaginary part |$1/\tau _{\alpha }$| describes the damping or growth rate of the mode. As the system described by equations (1)–(3) involves no dissipation (e.g. ideal fluid, adiabatic oscillations), the eigenvalues |$\sigma _{\alpha } = \omega _{\alpha }$| are real.

We adopt the slow-rotation approximation and expand the eigenfrequencies and mode functions in orders of |$\Omega =|\boldsymbol {\Omega }|$| in the corotating frame (Schenk et al. 2001), viz.

where the bracketed superscript denotes the order in |$\Omega$|⁠. The zeroth-order terms correspond to eigenfrequencies and mode functions in the non-rotating limit. The operator |$\bf{B}$| by definition [equation (10)] is of the order of |$\Omega$|⁠, i.e. |$\bf{B} = \bf{B}^{(1)}$|⁠. The operator |$\bf{C}$| is expanded as |$\bf{C} = \bf{C}^{(0)} + \bf{C}^{(2)} + \mathcal {O}(\Omega ^4)$|⁠, as implied by equation (4).

The frequencies of r-modes vanish in the non-rotating limit, i.e. |$\omega _{\alpha }^{(0)} = \bf{C}^{(0)} \cdot \boldsymbol {\xi }^{(0)}_{\alpha } = 0$|⁠. The non-trivial equation at order |$\Omega ^2$| reads

The component of equation (19) projected into the subspace spanned by zero-frequency modes |$\boldsymbol {\xi }_{\alpha }^{(0)}$| in a non-rotating barotrope gives (Schenk et al. 2001)

which is sufficient to determine |$\omega _{\alpha }^{(1)}$| and |$\boldsymbol {\xi }_{\alpha }^{(0)}$|⁠. Equation (20) can also be derived equivalently from vorticity conservation (Friedman & Stergioulas 2013).

In slowly rotating barotropic stars, r-modes exist only for |$l=|m|$| in Newtonian gravity (Friedman & Stergioulas 2013), where |$(l,m)$| are the orders of the spherical harmonic |$Y_{lm}$|⁠, and we write |$\alpha = (l, m)$|⁠. For purely axial r-modes, one obtains

where |$U_{lm}$| is a radial eigenfunction, and |$\hat{\boldsymbol {r}}$| is the radial unit vector. The rotational distortion of the stellar surface into an ellipse affects |$U_{lm}$| at order |$\Omega ^2$|⁠. Upon substituting equation (21) into equation (20), we arrive at

and

for |$l =|m|$|⁠, and |$U_{lm} = 0$| otherwise. The derivation is detailed in Appendix  A. Following the convention in Owen (2010), we normalize the Eulerian velocity perturbation of r-modes, |$\delta \boldsymbol {u}_{\!\alpha }(\boldsymbol {x}) = i \omega _{\alpha } \boldsymbol {\xi }_{\alpha }(\boldsymbol {x})$|⁠, according to

where |$R_*$| is the radius of the star.

2.4 Orthogonality

Eigenfunctions in rotating stars are not orthogonal in the conventional sense familiar from non-rotating stars. We define a symplectic product W of two vector functions (Friedman & Schutz 1978a)

where the angular brackets denote the inner product defined by

The symplectic product |$W(\boldsymbol {\xi }_{\alpha }, \boldsymbol {\xi }_{\beta })$| is time-independent (⁠|$\text{d}W/\text{d}t = 0$|⁠) for |$\boldsymbol {\xi }_{\alpha }$| and |$\boldsymbol {\xi }_{\beta }$| satisfying equation (8). When one has |$\omega _{\alpha } \ne \omega _{\beta }$| for |${\alpha }\ne {\beta }$|⁠, |$\text{d}W(\boldsymbol {\xi }_{\alpha }, \boldsymbol {\xi }_{\beta })/\text{d}t=0$| implies the modified orthogonality condition (Schenk et al. 2001; Pnigouras et al. 2024)

where |$\delta _{\!\alpha \beta }$| is the Kronecker delta, |$\mathcal {N}_{\alpha }$| is a normalization constant, and the Einstein summation convention is suppressed temporarily on the right-hand side of equation (27). Similarly, we have |$W(\boldsymbol {\xi }_{\alpha }^{*}, \boldsymbol {\xi }_{\alpha }) = 0$|⁠. Following Friedman et al. (2017), we write |$\mathcal {N}_{\alpha }$| as

with

For r-modes, the leading slow-rotation order of |$\kappa _{\alpha }$| is |$\kappa _{\alpha }^{(0)} = 1/2$| (see Appendix  A), and the leading order of |$\mathcal {N}_{\alpha }$| is

where |$\langle \cdot \, ,\, \cdot \rangle ^{(0)}$| denotes the inner product in the non-rotating limit, as introduced in Appendix  A, and the Einstein summation convention is suppressed temporarily again in equation (30).

3 EXCITATION BY ACCRETION

In this section, we calculate the response of r-modes to the mechanical impact of accreted matter. We assume that the accretion occurs stochastically as a random sequence of discrete clumps of matter, as implied by modern numerical simulations (Romanova & Owocki 2015). In Section 3.1, we formulate the associated inhomogeneous boundary value problem and solve it in terms of a Green’s function. The force densities corresponding to a single clump and a stochastic sequence of clumps with random arrival times are defined in Sections 3.2 and 3.3, respectively. We also calculate the temporal autocorrelation function of the stochastic response in Section 3.4, to lay the foundation for calculating the root-mean-square wave strain and power spectral density of the emitted gravitational radiation in Section 4.

3.1 General solution to the inhomogeneous problem

The inhomogeneous excitation problem can be formulated as

where |$\boldsymbol {\mathcal {F}}(\boldsymbol {x}, t)$| is the mechanical force density from accretion impacts. The standard mode decomposition of |$\boldsymbol {\xi }(\boldsymbol {x}, t)$| (Schenk et al. 2001; Dong & Melatos 2024; Pnigouras et al. 2024) couples together the coefficients in equation (31) for distinct modes |$\alpha$| and |$\alpha ^{\prime } \ne \alpha$| in a rotating star. Therefore, it is convenient to recast equation (31) as (Dyson & Schutz 1979)

Hence the orthogonality condition (27) for the column matrix |$\boldsymbol {\chi }_{\alpha } = [\boldsymbol {\xi }_{\alpha }, \boldsymbol {\dot{\xi }}_{\alpha }]^{\text{T}}$| becomes

where the superscript |$^{\text{T}}$| denotes the transpose, the dagger denotes the conjugate transpose, and the symplectic matrix |$\bf{W}$| is defined by (Dyson & Schutz 1979)

We expand |$\boldsymbol {\chi }$| according to the following ansatz (Schenk et al. 2001; Pnigouras et al. 2024),

where c.c. denotes the complex conjugate. Upon substituting equation (35) into equation (32) and applying equation (33), we arrive at the following equations of motion for |$c_{\alpha }(t)$| and |$c_{\alpha }^{*}(t)$|⁠,

equation (36) is diagonalized, in the sense that there is no coupling between the equations governing distinct modes |$\alpha$| and |$\alpha ^{\prime } \ne \alpha$|⁠. The solution to equation (36) is given by

assuming that the initial condition is |$c_{\alpha }(t) = 0$| for |$t\lt t_0$|⁠. We set |$t_0 = 0$| without loss of generality. The solution |$c_{\alpha }^{*}(t)$| is obtained by taking the complex conjugate of equation (37). Equation (37) gives an explicit expression for the mode amplitude |$c_{\alpha }(t)$| excited by the force density |$\boldsymbol {\mathcal {F}}(\boldsymbol {x}, t)$| and is a key input into the GW calculations in Section 4.

Real neutron stars have viscosity. We follow standard practice (Owen et al. 1998; Kokkotas, Apostolatos & Andersson 2001; Ho et al. 2020) and incorporate phenomenological damping or growth factors |$\propto \exp [-(t-t^{\prime })/\tau _{\alpha }]$| in equation (37) to account for viscous damping or growth due to the back reaction from gravitational radiation; see Section 4.

3.2 Force density from a single clump

We denote the force density from the impact of a single clump starting at |$t_{\rm s}$| as |$\boldsymbol {\mathcal {F}}(\boldsymbol {x}, t; t_{\rm s})$|⁠. For simplicity, we model |$\boldsymbol {\mathcal {F}}(\boldsymbol {x},t; t_{\rm s})$| to be a Dirac delta function in position (i.e. the idealized clump strikes the surface at a single point) and constant in time during its interaction with the stellar surface; see section 4 in Dong & Melatos (2024) and references therein for a physical justification. We then write

where |$\boldsymbol {p}_{\rm s}$| is the momentum of the clump in the corotating frame, |$T_{\rm s}$| is the duration of the impact, |$\boldsymbol {x}_{\rm s}$| is the impact point on the surface of the star, and |$H(\cdot)$| denotes the Heaviside step function.

The force density model in equation (38) is highly idealized. For example, non-polar-cap accretion flows and inhomogeneous hotspot geometries are expected to produce more complicated spatial profiles of the force density (Romanova et al. 2004; Kulkarni & Romanova 2008; Romanova & Owocki 2015). Improvements to the model could include extended radial and angular profiles corresponding to the aforementioned geometries, as well as oscillatory temporal profiles. The extended radial and angular profiles are likely to reduce the r-mode amplitude due to spatial averaging. The reduction is predicted to be smaller for r-modes, whose eigenfunctions are nodeless and monotonic in r, than for high-radial-order p-modes, whose eigenfunctions oscillate in r. Temporal oscillations affect the r-mode amplitude negligibly, as long as |$T_{\rm s} \lesssim \omega _{\alpha }^{-1}$| is satisfied. This is predicted to be the case even near resonance, where the oscillation frequency of the force density is close to |$\omega _{\alpha }$|⁠, because the impact duration is too short for the oscillatory force density to go through many full cycles. Equation (38) with |$T_{\rm s} \sim 10^{-4}\,$| s replicates the f-mode dominance in the GW energy spectrum emitted by an infalling, quadrupolar mass shell simulated by Nagar et al. (2004) (see section 5.3 in Dong & Melatos 2024). Nevertheless, it is straightforward to generalize the analysis to a more distributed (and hence more realistic) |$\boldsymbol {\mathcal {F}}(\boldsymbol {x},t; t_{\rm s})$| through equation (37), which expresses the mode response in terms of the Green function response to an impulsive excitation convolved with |$\boldsymbol {\mathcal {F}}(\boldsymbol {x},t; t_{\rm s})$|⁠.

3.3 Stochastic sequence of clumps

In real astrophysical systems, accretion is stochastic. The aperiodic X-ray flux is observed to fluctuate at frequencies |$\sim \text{kHz}$| in accreting neutron stars (Sunyaev & Revnivtsev 2000), largely due to mass accretion rate fluctuations arising from clumps in the flow. Clumps arise naturally in inhomogeneous accretion flows through self-healing magnetic and radiation-driven Rayleigh–Taylor instabilities (Morfill et al. 1984; Wang 1987; Demmel, Morfill & Atmanspacher 1990; Spruit & Taam 1993; D’Angelo & Spruit 2010). Quasi-periodic oscillations in X-ray light curves may be associated with beat frequencies between the angular velocities of the neutron star and the clumps (Lamb et al. 1985). Correlations between data in the near-IR and optical bands provide evidence of accreted plasma clumps in PSR J1023+0038 (Shahbaz et al. 2018). Typical X-ray pulsars have root-mean-square fractional variabilities of |$\sim 20 {{\ \rm per\ cent}}$| in their X-ray fluxes (Belloni & Hasinger 1990). The observed power spectral density of the X-ray flux in many objects is white below the turnover frequency and red above the turnover frequency, with power-law indices |$\approx 1$|–1.5 (Hoshino & Takeshima 1993; Lazzati & Stella 1997). The turnover frequency approaches |$\Omega$| for systems with magnetospheric radius close to the corotation radius (Revnivtsev et al. 2009).

State-of-the-art, 3D, magnetohydrodynamic simulations show that the accretion rate fluctuates on a time-scale of milliseconds in response to self-healing Rayleigh–Taylor and Kelvin–Helmholtz instabilities at the disc–magnetosphere boundary (Romanova, Kulkarni & Lovelace 2008; Romanova & Owocki 2015; Mushtukov et al. 2024). The millisecond variability of the accretion rate is consistent with the expected time-scale of the fragmentation of the clumpy accretion stream at the Keplerian period |$\Omega _{\rm K}^{-1}$| near the disc–magnetosphere boundary |$r = r_{\rm m}$|⁠, with |$\Omega _{\rm K}(r_{\rm m})^{-1} = (GM/r_{\rm m}^3)^{-1/2} = 2.7 (M/\text{M}_{\odot })^{-1/2} (r_{\rm m}/100\, \text{km})^{3/2}$| ms. It is also broadly consistent with the propagating fluctuation model for a geometrically thick accretion disc (⁠|$H/r_{\rm m} \sim 1$|⁠), in which clumps propagate inwards on the viscous time-scale in the inner emission region (Lyubarskii 1997). The viscous time-scale at |$r = r_{\rm m}$| is estimated as |$t_{\rm visc}(r_{\rm m}) \propto \alpha _{\rm disk}^{-1} (H/r_{\rm m})^{-2} \Omega _{\rm K}(r_{\rm m})^{-1}$|⁠, where |$\alpha _{\rm disk} \le 1$| is the viscosity parameter, and |$H/r_{\rm m}$| is the disc scale height to radius ratio (Lyubarskii 1997). The stochasticity enters equation (38) through random arrival times |$t_{\rm s}$|⁠, random impact points |$\boldsymbol {x}_{\rm s}$|⁠, the random impact duration |$T_{\rm s}$|⁠, and the random clump momentum |$\boldsymbol {p}_{\rm s}$| (Dong & Melatos 2024). However, it is difficult to resolve individual clumps observationally. In this paper, for the sake of simplicity, we assume that every clump has the same momentum and impact location, and that every impact has the same duration. We model the stochasticity in |$t_{\rm s}$| according to

where |$\lbrace t_{\rm s}\rbrace$| is a sequence of impact epochs drawn from a homogeneous Poisson point process with mean clump impact rate |$f_{\rm acc}$|⁠, and |$N(t)$| counts the total (and random) number of impacts up to time t. Equation (39) is consistent with the standard shot noise model (Terrell 1972; Lamb et al. 1985; Lochner, Swank & Szymkowiak 1991), with shot noise profiles in accretion discs replaced by |$\boldsymbol {\mathcal {F}}(\boldsymbol {x}, t; t_{\rm s})$|⁠. The shot noise model reproduces the observed red-noise power spectrum in some (although not all) sources (Lazzati & Stella 1997) and is representative enough to be used as a first pass to study the stochastic excitation of oscillations by accretion. It is also worth noting that the shot noise model struggles to reproduce the observed proportionality between rms variability and X-ray flux on short time-scales in some sources (Uttley & McHardy 2001). Generalizing the analysis to more sophisticated models, such as the propagating fluctuation model (Lyubarskii 1997), is left for future work.

3.4 Temporal autocorrelation function

The autocorrelation function of |$c_{\alpha }(t)$| excited by the force density in equation (39) can be calculated from

where angular brackets with the subscript ‘s’ denote the ensemble average over all random realizations of the stochastic impact process. Equation (40) is derived in appendix D in Dong & Melatos (2024). For |$T_{\rm s}^{-1} \gg \omega _{\alpha }, \omega _{\beta }$|⁠, the first term on the right-hand side of equation (40) reduces to

Notably, |$T_{\rm s}$| does not appear in equation (41), as an impact with |$T_{\rm s}^{-1} \gg \omega _{\alpha }, \omega _{\beta }$| approximately behaves, as if the mode is excited by a delta-function impulse. Equation (41) tends to a constant for |$\tau _{\alpha } \gt 0$| and |$\tau _{\beta } \gt 0$|⁠, as t and |$t^{\prime }$| go to infinity. The sign of |$\tau _{\alpha }$| depends on whether or not the CFS instability is triggered; see Section 5. The second term on the right-hand side of equation (40) is negligible, unless the two modes have the same frequency. With |$T_{\rm s}^{-1} \gg \omega _{\alpha }$| and |$\alpha = \beta$|⁠, the second term on the right-hand side of equation (40) reduces to

where we denote the time lag between t and |$t^{\prime }$| by |$\zeta =|t-t^{\prime }|$|⁠. Equation (42) is stationary, i.e. it depends only on |$\zeta$|⁠, in the limit |$\min (t, t^{\prime })\rightarrow \infty$|⁠. The root-mean-square mode amplitude |$c_{\alpha ,\text{rms}}$| can be evaluated by setting |$\alpha =\beta$| and |$\zeta = 0$| in equation (40).

We note for completeness that the following result is also useful for calculating the root-mean-square GW strain in Section 4:

equation (43) is related closely to equation (42), except that we do not take the complex conjugate of |$c_{\alpha }(t^{\prime }; t_{\rm s})$|⁠.

4 GRAVITATIONAL RADIATION

In this section, we explore the stochastic GW signal emitted by the accretion-excited r-modes calculated in Section 3. We calculate the wave strain generated by modes excited by a single clump in Section 4.1. In Sections 4.2 and 4.3, we calculate the root-mean-square wave strain and the power spectral density of the emitted gravitational radiation, respectively, in order to investigate its detectability.

4.1 Instantaneous wave strain

The metric perturbation |$h_{jk}^{\rm TT}$| in the transverse traceless gauge from Newtonian r-modes, measured by an observer at a distance d, is given by a linear combination of the time-varying current multipole moments (Thorne 1980),

with (Melatos & Peralta 2009)

where t is the retarded time, and |$T^{\text{B}2,lm}_{jk}$| is the gravitomagnetic tensor describing the angular beam pattern associated with the |$(l, m)$| spherical harmonic. Equation (46) follows from equation (45), because one has |$\delta \boldsymbol {u} = \delta \boldsymbol {v}$| and |$\delta \rho = 0$| for purely axial r-modes.

We write equation (44) in terms of |$h_0(t; t_{\rm s}, l)$|⁠, the strain generated in the l-th multipole moment by a mode |$c_{\alpha }(t; t_{\rm s})$| excited by a single clump striking at time |$t_{\rm s}$|⁠, as

with

Interestingly, |$h_0(t; t_{\rm s}, l)$| does not depend on the EOS, because |$\delta S_{lm}(\delta \boldsymbol {u}_{\!\alpha })$| cancels out |$1/\mathcal {N}_{\alpha }^{(1)}$| in |$c_{\alpha }(t; t_{\rm s})$|⁠, and |$\rho$| does not enter the integrand in equation (37). The cancellation occurs physically, because purely axial r-modes in a barotropic star have power-law radial profiles |$U_{lm}(r)$|⁠. Evaluating equation (48) assuming equation (38), as well as |$T_{\rm s}^{-1} \gg \omega _{\alpha }$| and |$\omega _{\alpha } \tau _{\alpha } \gg 1$|⁠, the strain |$h_0(t; t_{\rm s}, 2)$| from the dominant |$l=2$| multipole simplifies to

where we write |$\gamma _{v} = [1-(|\boldsymbol {v}|/c)^2]^{-1/2}$|⁠, |$\nu _{\rm s}$| is the spin frequency, |$\dot{M}$| is the accretion rate, |$\boldsymbol {v}$| is the clump velocity in the corotating frame, and |$\omega _{\alpha , \text{i}}$| represents the angular velocity in the inertial frame. In deriving equation (49), we assume that the accretion rate is constant over the accretion episode and write |$\boldsymbol {p}_{\rm s} = \gamma _{v} m \boldsymbol {v}$| in the force density [equation (38)], where |$m = \dot{M}/f_{\rm acc}$| is the typical mass of a clump.

The strain in equation (49) scales |$\propto \nu _{\rm s}^{2}$|⁠. This differs from the |$\nu _{\rm s}^{3}$| scaling in Riles (2023), because the impact-excited r-mode amplitude scales |$\propto \nu _{\rm s}^{-1}$|⁠. Physically, this occurs because the |$t^{\prime }$| integral in equation (37) is independent of |$\nu _{\rm s}$| for |$T \ll \omega _{\alpha }^{-1}$|⁠, and we have |$c_{\alpha } \propto \mathcal {N}_{\alpha }^{-1} \propto \omega _{\alpha }^{-1} \propto \nu _{\rm s}^{-1}$|⁠, i.e. the mode energy |$E_{\alpha } \propto |c_{\alpha }|^{2} \nu _{\rm s}^{2}$| is independent of |$\nu _{\rm s}$|⁠. The impact location |$\boldsymbol {x}_{\rm s}$| not only affects the amplitude, but also the phase through the |$e^{2i\phi }$| term in |$Y_{22}$|⁠. When |$\boldsymbol {x}_{\rm s}$| is fixed, the phase is set to zero without loss of generality.

4.2 Root-mean-square wave strain

In Sections 4.2 and 4.3, we assume that the viscous damping time-scale |$\tau ^{\text{visc}}_{\alpha }$| and gravitational radiation reaction time-scale |$\tau ^{\text{gw}}_{\alpha }$| are related by |$\tau ^{\text{visc}}_{\alpha } \lt \left|\tau ^{\text{gw}}_{\alpha }\right|$|⁠, so that the CFS instability is not triggered. We define the temporal autocorrelation function of the strain as |$C_{\alpha }(t, t^{\prime }) = \left\langle h_{0}(t; \alpha) h_{0}^{*}(t^{\prime }; \alpha) \right\rangle _{\rm s},$| with |$h_0(t; \alpha) = \sum _{t_\text{s} \le t} h_0(t; t_{\rm s}, \alpha) .$| We calculate the root-mean-square characteristic wave strain, |$h_{\text{rms}}$|⁠, by evaluating |$C_{\alpha }(t, t^{\prime })^{1/2}$| at zero lag (⁠|$\zeta = 0$|⁠), in order to assess the detectability of the stochastic signal.

In the regime |$\omega _{\alpha } \tau _{\alpha } \gg 1$|⁠, the dominant term of |$C_{\alpha }(t, t^{\prime })$| reads

with

Physically, |$K_{\alpha }(\zeta =0; t)$| equals the mean number of clumps striking the stellar surface up to time t for the mode |$\alpha$|⁠. In the limit |$t \rightarrow \infty$|⁠, |$K_{\alpha }(\zeta =0; t)$| tends to |$f_{\rm acc} \tau _{\alpha }$|⁠, the saturated number of clumps during the damping time-scale |$\tau _{\alpha }$|⁠. In the regime |$\min (t, t^{\prime }) \ll \tau _{\alpha }$|⁠, |$K_{\alpha }(\zeta =0; t)$| reduces to |$2 f_{\rm acc} t$|⁠, and |$h_{\rm rms}$| depends on the duration of any particular uninterrupted accretion episode, with |$\Delta t_{\rm acc} \ll \tau _{\alpha }$|⁠. Hence, we have

One can easily obtain the saturated |$h_{\rm rms}$|⁠, by replacing |$\Delta t_{\rm acc}$| with |$\tau _{\alpha }/2$| in equation (52).

The impact velocity direction |$\hat{\boldsymbol {v}}$| and the impact location |$\boldsymbol {x}_{\rm s}$| play important roles in determining |$h_{\rm rms}$|⁠. For example, when taking |$\hat{\boldsymbol {v}} = \hat{\boldsymbol {\theta }}$| in equation (52), we find |$\left| \hat{\boldsymbol {v}} \cdot \left[(\boldsymbol {x} \times \nabla)Y_{22}(\hat{\boldsymbol {x}})\right]_{\boldsymbol {x} = \boldsymbol {x}_{\rm s}} \right|$| attains its maximum value of 0.77 at the equator. As |$v_{\phi }$| increases, |$h_{\rm rms}$| decreases. The impact location |$\boldsymbol {x}_{\rm s}$| that maximizes |$h_{\rm rms}$| shifts towards the pole for |$\hat{v}_{\phi } \gt \hat{v}_{\theta }$|⁠. It maximizes |$h_{\rm rms}$| at |$\theta = \pi /4$|⁠, with |$\left| \hat{\boldsymbol {v}} \cdot \left[(\boldsymbol {x} \times \nabla)Y_{22}(\hat{\boldsymbol {x}})\right]_{\boldsymbol {x} = \boldsymbol {x}_{\rm s}} \right| = 0.39$|⁠, when |$\hat{\boldsymbol {v}} = \hat{\boldsymbol {\phi }}$|⁠. Equation (52) implies that GWs emitted by impact-excited r-mode oscillations in neutron stars are too weak to be detected by the current generation of LIGO detectors. In an optimistic scenario, with |$v_{\rm s}=0.5$| kHz and |$\Delta t_{\rm acc}=10^{6}$| yr, one finds |$h_{\rm rms} \lesssim 1 \times 10^{-28}$|⁠, which corresponds to |$|c_{\alpha }| \lesssim 10^{-8}$| (Riles 2023). These values of |$h_{\rm rms}$| and |$|c_{\alpha }|$| are below the upper bounds |$h_{0}^{95 {{\ \rm per\ cent}}} \sim 10^{-25}$| (at 95  per cent confidence level) inferred from recent searches of continuous waves from accreting systems (Middleton et al. 2020; Covas et al. 2022; LIGO Scientific Collaboration 2022) and the theoretically predicted amplitude set by non-linear mode coupling, viz. |$10^{-7} \lesssim |c_{\alpha }| \lesssim 10^{-4}$| (Bondarescu et al. 2007; Gusakov, Chugunov & Kantor 2014).

The impact-excited f- and p-modes in a non-rotating, barotropic star emit GWs with |$h_{\rm rms} \lesssim 10^{-33}$| (Dong & Melatos 2024). The f- and p-modes are strongly averaged during the top-hat impact. On the other hand, the g-modes, like r-modes, are less affected by the top-hat averaging. The GW strain emitted by g-modes varies from |$10^{-34}$| to |$10^{-32}$|⁠, depending on |$\Delta t_{\rm acc}$| and |$\hat{\boldsymbol {v}}$| (Dong & Melatos 2024). There are also differences between the EOS dependence of f- and p-modes and that of r-modes in a barotropic star. For instance, the GWs from f- and p-modes are weaker for a stiffer EOS (larger |$\partial \log {P}/\partial \log {\rho }$|⁠) (Dong & Melatos 2024). In contrast, the GWs from impact-excited r-modes are approximately independent of the EOS, as explained in Section 4.1, except that the EOS enters |$h_{\rm rms}$| implicitly through |$\tau _{\alpha }$|⁠.

4.3 Power spectral density

A complementary way to assess the detectability of the GW signal is to calculate its power spectral density, |$S(f)$|⁠, where f is the Fourier frequency, and compare it with the sensitivity curves of the LIGO detectors. The power spectral density measures how radiated power is distributed in the frequency domain and helps to guide narrowband searches. In the stationary regime, i.e. |$\tau _{\alpha } \ll \min (t, t^{\prime })$| and |$C_{\alpha }(t, t^{\prime }) = C_{\alpha }(\zeta)$|⁠, we compute the one-sided |$S(f)$| by Fourier transforming equation (50) with respect to |$\zeta$|⁠, viz.

Examples of |$S(f)^{1/2}$| in the LIGO observing band are plotted in Fig. 1, for three spin frequencies, viz. |$\nu _{\rm s}=10$|⁠, 100, and 716 Hz¹, and r-modes with |$l=m=2$|⁠, 3, and 4. Each peak is actually 601 peaks, with |$\tau _{\alpha }$| spanning |$10^{-6} \le \tau _{\alpha }/(1\, \text{yr}) \le 10^{6}$|⁠, which almost overlap. The colour gradient along each peak, defined by the colour bar at the right of the figure, describes how |$S(f)^{1/2}$| near |$f\approx \omega _{\alpha }/(2\pi)$| varies with |$\tau _{\alpha }$| for multiple peaks. The difference in |$S(f)^{1/2}$| away from the peak for different |$\tau _{\alpha }$| is not visually apparent, due to the broad logarithmic scale. The damping time-scale |$\tau _{\alpha }$| serves as a proxy for the EOS, as discussed in Section 4.2. We see that the peaks grow narrower and decrease in height as |$l=m$| increases. The highest peak reaches |$\max _{f} S(f)^{1/2} \approx 1.3 \times 10^{-21} \,\text{Hz}^{-1/2}, 1.4 \times 10^{-22}\, \text{Hz}^{-1/2}$|⁠, and |$1.7 \times 10^{-23}\, \text{Hz}^{-1/2}$| for |$l=m=2$|⁠, 3, and 4, respectively, at |$\nu _{\rm s} = 716$| Hz and |$\tau _{\alpha } = 10^{6}$| yr. The peak height scales |$\propto \nu _{\rm s}^{2} \tau _{\alpha }$|⁠, for |$\omega _{\alpha } \tau _{\alpha } \gg 1$|⁠.

Figure 1.

Amplitude spectral density |$S(f)^{1/2}$| versus Fourier frequency f within the LIGO observing band, for |$l=m=2$|⁠, 3, and 4 and |$\nu _{\rm s} = 10 \,$|⁠, 100, and 716 Hz. Modes are labelled by {|$(l,m), \nu _{\rm s}$|}. The peaks of the amplitude spectral density are colour-coded by the damping time-scale |$\tau _{\alpha } \in [10^{-6}, 10^{6}]$|  yr. The colour gradient along each curve is an illusion; every curve has a single colour, but multiple curves overlap indistinguishably except near their peaks. The sensitivity curve achieved during the third LIGO observing run (blue curve) (LIGO Scientific Collaboration 2021) and projected sensitivity curves for the Einstein Telescope (baseline 10 km, red curve) (Hild et al. 2011) and Cosmic Explorer (baseline 40 km, orange curve) (Evans et al. 2021) are also shown for comparison. The parameter |$\left| \hat{\boldsymbol {v}} \cdot \left[(\boldsymbol {x} \times \nabla)Y_{lm}(\hat{\boldsymbol {x}})\right]_{\boldsymbol {x} = \boldsymbol {x}_{\rm s}} \right|$| is set to unity as an upper bound. Other parameters: |$R_{*} = 10\,$| km, |$d = 1\,$| kpc, |$\dot{M} = 10^{-8}\, \text{M}_{\odot } \, \text{yr}^{-1}$|⁠, |$f_{\rm acc} = 1\,$| kHz, and |$|\boldsymbol {v}| = 0.4c$|⁠.

Open in new tabDownload slide

For modes with |$l=m=2$| or 3, the peaks of |$S(f)^{1/2}$| exceed the LIGO sensitivity (blue curve) in the O3 observing run (LIGO Scientific Collaboration 2021) for |$\tau _{\alpha } \gtrsim 10^{5}$| yr. However, this does not imply automatically, that the GW signal is detectable. One needs an integration time |$T_{\rm obs} \gtrsim \tau _{\alpha }$| to detect the peak associated with the mode |$\alpha$|⁠. As one example, for the fastest spinning neutron star at |$\nu _{\rm s} = 716$| Hz and |$d\approx 5.5$| kpc (Hessels et al. 2006; Ortolani et al. 2007) with |$\tau _{\alpha } = 10^{6}$| yr, one needs |$T_{\rm obs} \gtrsim 10^{3}$| yr to reach the LIGO sensitivity curve. Therefore, the GW signal from impact-excited r-modes is not detectable by the current generation of LIGO interferometers, consistent with Section 4.2. On the other hand, the observing time requirements are less demanding for the Einstein Telescope (Hild et al. 2011) or Cosmic Explorer (Evans et al. 2021), with |$T_{\text{obs}} \sim 1$|  yr, for the parameters above. Signals from low-mass X-ray binaries with |$\tau _{\alpha } \lesssim 10^{5}$|  yr remain challenging to detect even with the next generation of interferometers.

5 CFS INSTABILITY

In the preceding sections, we assume that the r-modes are CFS-stable. The assumption is motivated observationally by the non-detection of GWs from r-modes, which implies that the dimensionless r-mode amplitude satisfies |$\alpha _{r} \lesssim 10^{-9}$| for millisecond pulsars at frequencies inside the nominal instability window (Ho, Heinke & Chugunov 2019; Boztepe et al. 2020). The CFS instability is triggered, when |$\Omega$| exceeds a critical value |$\Omega _{\rm CFS}$|⁠, leading to |$\tau _{\alpha } \lt 0$| for |$\Omega \gt \Omega _{\rm CFS}$|⁠. The nominal instability window is computed, assuming that |$\tau _{\alpha }$| depends only on the gravitational radiation reaction, bulk viscosity from modified URCA reactions, and shear viscosity from neutron–neutron and electron–electron scattering (Friedman & Stergioulas 2013; Glampedakis & Gualtieri 2018; Ho et al. 2019). Generalizations of the minimal damping model, which lie outside the scope of this paper, include Ekman layer damping at the crust–core interface (Bildsten & Ushomirsky 1999; Levin & Ushomirsky 2001; Bondarescu et al. 2007), mutual friction (Haskell, Andersson & Passamonti 2009), and hyperons or strange quarks (Nayyar & Owen 2006; Alford, Mahmoodifar & Schwenzer 2012). The reader is directed to section 9.6.2 of Friedman & Stergioulas (2013) and sections V.E– V.G of Glampedakis & Gualtieri (2018) for further details and additional references. It is instructive to analyse CFS-stable r-modes on an equal footing with the stably excited f-, p-, and g-modes, extending the study by Dong & Melatos (2024). The resulting, stochastic GW signal represents quiescent-like emission, which complements the quasi-monochromatic emission from CFS-unstable r-modes analysed by other authors (Andersson 1998; Friedman & Morsink 1998; Owen et al. 1998).

Once the instability switches on, the r-mode amplitude grows exponentially. The saturation of the r-mode amplitude at a value |$|c_{\alpha , \text{sat}}|$| through non-linear mode couplings has been analysed in detail in the literature (Schenk et al. 2001; Arras et al. 2003; Bondarescu et al. 2007). The instability is quenched, as the gravitational radiation carries away angular momentum, lowering |$\Omega$| below |$\Omega _{\rm CFS}$|⁠.

In this section, we briefly discuss how the temporal autocorrelation function |$C_{\alpha }(t, t^{\prime })$| is modified, when the CFS instability switches on and off. We consider a single on-off episode, noting that an ongoing cycle of on-off episodes is likely to occur in reality (Levin 1999; Heyl 2002; Arras et al. 2003; Bondarescu et al. 2007; Glampedakis & Gualtieri 2018). For a single on-off episode, we write

where |$t_{\rm on}$| and |$t_{\rm off} \gt t_{\rm on}$| are the on-off trigger epochs, respectively, and |$\tau _{\alpha , \text{on}}$| and |$\tau _{\alpha , \text{off}}$| are the damping time-scales when the CFS instability is on and off, respectively. Equation (54) is a simplified model, which relies on two assumptions: (i) |$\tau _{\alpha , \text{off}}$| and |$\tau _{\alpha , \text{on}}$| are constant, with an instantaneous transition between the CFS-stable and the CFS-unstable regime; and (ii) |$\tau _{\alpha , \text{off}}$| is assumed to be the same for |$t \lt t_{\rm on}$| and |$t \gt t_{\rm off}$|⁠. The former assumption is justified by treating |$\tau _{\alpha , \text{off}}$| and |$\tau _{\alpha , \text{on}}$| as averaged, effective damping/growth time-scales in the CFS-stable and CFS-unstable regimes, respectively. The latter assumption is likely to be inaccurate immediately after the onset of the CFS-stable regime. In future work, the calculation will be improved by replacing |$\tau _{\alpha }$| with |$\tau _{\alpha }(t)$| calculated numerically to account for spin and temperature evolution (Levin 1999; Bondarescu et al. 2007).

The dominant contribution to |$C_{\alpha }(t, t^{\prime }) \propto K_{\alpha }[\zeta , \min (t, t^{\prime })]$| is given by the second term in equation (40), and reads

for |$t^{\prime } \gt t \gt t_{\rm off}$|⁠, where a tilde above a temporal quantity indicates, that the quantity is normalized by |$\tau _{\alpha , \text{off}}$|⁠. The duration of the CFS instability is denoted by |$T_{\rm CFS} = t_{\rm off}-t_{\rm on}$|⁠. Equation (55) reduces to equation (51) for |$\tau _{\alpha , \text{on}} = \tau _{\alpha , \text{off}}$| or |$T_{\rm CFS} = 0$|⁠, i.e. when the CFS instability is not triggered.

Interestingly, equation (55) can increase or decrease with t depending on the values of |$\tilde{t}_{\rm off}, \tilde{T}_{\rm CFS}$|⁠, and |$\tilde{\tau }_{\alpha , \text{on}}$|⁠, whereas equation (51) always increases with t. From equation (55), we find

In deriving equation (56), we assume |$\left|\text{d} \ln {\nu _{\rm s}}/\text{d} t\right| \ll \left|\partial \ln {K_{\alpha }(\zeta ; t)}/\partial t\right|$|⁠, which is valid throughout the observation span, |$\left|\text{d} \ln (\dot{M}|\boldsymbol {v}|)/\text{d} t\right| \ll \left|\partial \ln {K_{\alpha }(\zeta ; t)}/\partial t\right|$|⁠, and |$\left|\text{d} \ln {f_{\rm acc}}/\text{d} t\right| \ll \left|\partial \ln {K_{\alpha }(\zeta ; t)}/\partial t\right|$|⁠. For |$\tilde{T}_{\rm CFS} \gt \tilde{T}_{\rm CFS, crit} =|\tilde{\tau }_{\alpha , \text{on}}| \ln [(1 +|\tilde{\tau }_{\alpha , \text{on}}|) /|\tilde{\tau }_{\alpha , \text{on}}|]/2$|⁠, one always has |$\partial C_{\alpha }(\zeta ; t)/\partial t \lt 0$|⁠, for all |$\tilde{t}_{\rm on} \ge 0$|⁠. Physically, this is because the CFS instability, coexisting with impact-excited r-modes, lasts long enough, so that |$C_{\alpha }(\zeta ; t \gt \tilde{t}_{\rm off})$| exceeds the limit |$C_{\alpha , \text{sat}}(\zeta)$| set by the saturated number of clumps during |$\tau _{\alpha }$|⁠, regardless of the initial amplitude. For |$\tilde{T}_{\rm CFS} \lt \tilde{T}_{\rm CFS, crit}$| on the other hand, one has |$\partial C_{\alpha }(\zeta ; t)/\partial t \gt 0$| for |$\tilde{t}_{\rm on} \lt \tilde{t}_{\rm on, crit}$|⁠, and vice versa. These regimes, where |$\tilde{T}_{\rm CFS} \lt \tilde{T}_{\rm CFS, crit}$|⁠, correspond to |$C_{\alpha }(\zeta ; t \gt \tilde{t}_{\rm off}) \lt C_{\alpha , \text{sat}}(\zeta)$|⁠, if the initial amplitude is small enough, or |$C_{\alpha }(\zeta ; t \gt \tilde{t}_{\rm off}) \gt C_{\alpha , \text{sat}}(\zeta)$|⁠, if the initial amplitude is high enough. The critical value of |$\tilde{t}_{\rm on, crit}$| is given by |$\tilde{t}_{\rm on, crit} = - \ln \lbrace (1 +|\tilde{\tau }_{\alpha , \text{on}}|) [1-\exp (-2 \tilde{T}_{\rm CFS} /|\tilde{\tau }_{\alpha , \text{on}}|)]\rbrace /2$|⁠.

What are the characteristic astrophysical numbers associated with the above regimes? We estimate |$T_{\rm CFS}$| by the spin-down time-scale |$T_{\rm sd} = \Omega /|\dot{\Omega }| \approx 5|c_{\alpha , \text{sat}}|^{-2}|\tau _{\alpha , \text{on}}|$| (Levin 1999; Glampedakis & Gualtieri 2018). For a system that does not reach |$|\tau ^{\text{gw}}_{\alpha }| = \tau ^{\text{visc}}_{\alpha }$| for |$t \gt t_{\rm off}$| [i.e. a cycle with thermogravitational runaway, followed by saturation], viscous damping is dominated by the shear viscosity, with |$\tau _{22, \text{off}} \approx 1\, (T/10^{9}\, \text{K})^{2} \,$|  yr from electron–electron scattering for temperatures |$\lesssim 10^{9}\, \text{K}$| (Friedman & Stergioulas 2013). One then requires |$|c_{\alpha , \text{sat}}|$| of the order of unity in order to achieve |$\tilde{T}_{\rm CFS} \lt \tilde{T}_{\rm CFS, crit}$|⁠, for an optimistic value |$|\tau _{\alpha , \text{on}}| =|\tau ^{\text{gw}}_{\alpha }| \approx 7 \times 10^{-7} (\nu _{\rm s}/1 \, {\rm kHz})^{-6}\, \text{yr}$|⁠.² Therefore, a neutron star whose CFS instability evolves cyclically almost always has |$\partial C_{\alpha }(\zeta ; t)/\partial t \lt 0$|⁠. On the other hand, if a neutron star evolves towards equilibrium on the r-mode stability curve (Bondarescu et al. 2007), we have |$T_{\rm CFS, crit} \rightarrow \infty$| and |$\partial C_{\alpha }(\zeta ; t)/\partial t \propto 2 f_{\rm acc} \gt 0$|⁠.

Table 1 lists the sign of |$\partial C_{\alpha }(\zeta ; t\gt t_{\rm off})/\partial t$| in three regimes and the possible physical scenarios leading up to the observation span in each regime. The second order time derivative |$\partial ^2 C_{\alpha }(\zeta ; t)/\partial t^2 = -(2/\tau _{\alpha }) \partial C_{\alpha }(\zeta ; t)/\partial t$| can be used to discriminate whether the CFS instability is on or off during the observation span, with |$\tau _{\alpha } = \tau _{\alpha , \text{on}} \lt 0$| and |$\tau _{\alpha } = \tau _{\alpha , \text{off}} \gt 0$|⁠, respectively.

In summary, the slope and concavity of |$C_{\alpha }[\zeta , \min (t, t^{\prime })]$| can be used to infer observationally the coexistence or otherwise of the CFS instability with impact-excited the r-mode oscillations, at least in principle. That is, by autocorrelating the wave strain |$h(t)$| detected by an instrument like LIGO, we can infer whether the CFS instability is on or off before or during the observation span, assuming the change in |$\nu _{\rm s}$| is negligible. We can also infer upper/lower bounds on the duration of the CFS instability or the duration of the uninterrupted accretion episode, compared to the unknown damping time-scale |$\tau _{\alpha , \text{off}}$|⁠, even if these durations are longer than the observation span itself.

6 CONCLUSIONS

The gravitational radiation emitted by r-modes that are mechanically excited by the stochastic impact of clumps of matter on an accreting neutron star is studied. The standard results for the r-mode eigenfrequencies and eigenfunctions, in an unmagnetized, slowly rotating, one-component barotrope are reproduced, assuming Newtonian gravity. We then calculate analytically, using a Green’s function, the response of the r-modes to impulsive, top-hat impacts which obey a Poisson point process. The stochastic radiated signal is quantified in terms of the instantaneous and root-mean-square wave strain [equations (49) and (52), respectively], as well as the power spectral density (Fig. 1).

The instantaneous wave strain excited by a single clump, |$h_0(t; t_{\rm s})$|⁠, does not depend strongly on the EOS. The root-mean-square wave strain satisfies |$h_{\rm rms} = N_{\rm clumps}^{1/2} h_0(t; t_{\rm s})$|⁠, where |$N_{\rm clumps}$| is the characteristic number of clumps striking the stellar surface, while the r-modes oscillate, which depends implicitly on the EOS through the damping time-scale |$\tau _{\alpha }$|⁠. For an uninterrupted accretion episode with duration |$\Delta t_{\rm acc} \ll \tau _{\alpha }$|⁠, one has |$h_{\rm rms} \lesssim 10^{-35} (\nu _{\rm s}/ 10\, {\rm Hz})^{2} (\Delta t_{\rm acc}/1\, {\rm yr})^{1/2}$|⁠, with the exact value depending on the impact velocity direction |$\hat{\boldsymbol {v}}$| and the impact location |$\boldsymbol {x}_{\rm s}$|⁠. This is comparable to the strain from impact-excited f-, p-, and g-modes in a non-rotating star (Dong & Melatos 2024), and is too weak to be detected by the current generation of LIGO detectors. The amplitude spectral density peaks at |$\max _f S^{1/2}(f) \sim 10^{-22} \text{Hz}^{-1/2}$| and thereby exceeds nominally the LIGO sensitivity curve in the O3 observing run for |$l=m=2$|r-modes with |$\nu _{\rm s} = 100\,$| Hz and |$\tau _{\alpha } \sim 10^{6}$| yr. The peak amplitude spectral density scales approximately as |$\propto \nu _{\rm s}^{2} \tau _{\alpha }$|⁠. However, the signal is not detectable, as the required integration time, |$T_{\rm obs} \gtrsim \tau _{\alpha } \sim 10^{6}$| yr, is much longer than a typical LIGO observation, consistent with the conclusion drawn from |$h_{\rm rms}$|⁠.

The slope and concavity of the temporal autocorrelation function offer an observational route to test for the coexistence of the CFS instability with impact-excited r-modes. Specifically, coexistence strictly before the GW observation implies |$\partial C_{\alpha }(\zeta ; t)/\partial t \lt 0$|⁠, and coexistence strictly during the GW observation implies |$\partial C_{\alpha }(\zeta ; t)/\partial t \gt 0$| and |$\partial ^2 C_{\alpha }(\zeta ; t)/\partial t^2 \gt 0$|⁠.

ACKNOWLEDGEMENTS

This research is supported by the Australian Research Council (ARC) through the Centre of Excellence for Gravitational Wave Discovery (OzGrav) (grant number CE230100016). We thank the anonymous referee for their constructive feedback.

DATA AVAILABILITY

There are no data generated in this paper.

Footnotes

REFERENCES

APPENDIX A: DERIVATION OF r-MODE PROPERTIES

In this appendix, we derive in abridged form some key properties of purely axial r-modes, which serve as inputs into the GW calculations in Sections 4 and 5. Eigenfunctions and eigenfrequencies are discussed in Appendix A1. Normalization is discussed in Appendix A2. The gravitational wave growth time-scale due to the CFS instability is discussed in Appendix A3.

A1 Frequencies and eigenfunctions

We work in a basis of unit vectors (⁠|$\hat{\boldsymbol {r}}, \hat{\boldsymbol {\theta }}, \hat{\boldsymbol {\phi }}$|⁠) in spherical polar coordinates. Equation (21) expressed in component form is given by

We write the angular velocity as |$\boldsymbol {\Omega } = \Omega \hat{\bf z} = \Omega \cos {\theta } \hat{\boldsymbol {r}}-\Omega \sin {\theta } \hat{\boldsymbol {\theta }}$|⁠, and hence obtain

Upon substituting (A1) and (A2) into (20), the |$\hat{\boldsymbol {r}}$| component reduces to

After applying the orthogonality condition of the spherical harmonics, equation (A3) implies

for an individual |$l \ge |m|$|⁠, and |$U_{l} = 0$| otherwise. The |$\hat{\boldsymbol {\theta }}$| component reduces to

where the prime denotes a derivative with respect to r, and we exploit the recursion relations

with

equation (A5) agrees with equation (35) in Lockitch & Friedman (1999) for |$U_{l}(r) \ne 0$|⁠. Upon applying (A4) and the orthogonality condition of the spherical harmonics to (A5), we arrive at two equations for the r-mode eigenfunctions (Gittins & Andersson 2023),

In order to simultaneously satisfy equations (A9) and (A10), we must have |$l =|m|$|⁠, and

Therefore, we have shown that one has |$U_{l} \propto r^{l + 1}$| for |$l =|m|$|⁠, and |$U_{l} = 0$| otherwise (Friedman & Stergioulas 2013).

A2 Normalization

In this subsection, we derive the leading-order normalization constants |$\mathcal {N}_{\alpha }$|⁠, in order to quantify the gravitational response of the r-modes to the mechanical impact of accreted matter in Section 4.

We introduce the slow-rotation expansion of the equilibrium density |$\rho _0$| in the inner product (26) (Schenk et al. 2001) and express it as

The angular frequencies of r-modes is |$\sim \mathcal {O}(\Omega)$|⁠. Hence, to leading order, we obtain

with

Upon substituting (A1) and (A2) into (26), we obtain

To pass from (A17) to (A18), we integrate by parts. We then use (A4) and

to convert (A19) into (A20). Simplifying (A15) using (A20), we obtain |$\kappa _{\alpha }^{(0)} = 1/2$|⁠.

A3 Growth time-scale of the CFS instability

In this subsection, we briefly review the theory of the CFS instability of r-modes to support Section 5, where the effect of the CFS instability on GW signals from clump-excited r-modes is discussed.

The energy of the mode in the rotating frame is given by

where angular brackets denote an inner product defined in Section 2.4. The power radiated by a normalized r-mode due to gravitational radiation reaction in the inertial frame is (Thorne 1980)

The rotating frame power |$\dot{E}_{\alpha , \text{r}}$| is related to |$\dot{E}_{\alpha , \text{i}}$| via the gravitational radiation reaction torque |$\dot{J}_{\alpha }$|⁠, viz.

where we use |$\dot{J}_{\alpha } = - (m/\omega _{\alpha , \text{i}}) \dot{E}_{\alpha , \text{i}}$| and |$\omega _{\alpha , \text{i}} = \omega _{\alpha }-m\Omega$|⁠. We estimate the growth time-scale due to gravitational radiation reaction a posteriori by (Owen et al. 1998; Friedman & Stergioulas 2013)

When |$\omega _{\alpha , \text{i}}$| and |$\omega _{\alpha }$| have opposite signs, one obtains |$\dot{E}_{\alpha , \text{r}} \gt 0$| and |$\tau ^{\text{gw}}_{\alpha } \lt 0$|⁠, indicating that such modes are unstable. This is known as the CFS instability (Friedman & Schutz 1978b).

Upon combining equations (22)–(24), (46), and (A22)–(A26), we arrive at

For r-modes with |$l=m=2$| in a uniform density neutron star, equation (A27) reduces to (Friedman & Stergioulas 2013)