Gravitational waves from nonradial oscillations of stochastically accreting neutron stars

Oscillating neutron stars are sources of continuous gravitational waves. We study analytically the excitation of stellar oscillations by the mechanical impact on the stellar surface of ''clumps'' of stochastically accreted matter. We calculate the waveform and spectrum of the gravitational wave signal emitted by the accretion-driven pulsations. Results are generated for an idealised model of a nonrotating, unmagnetised, one-component star with uniform polytropic index $n_{\rm poly}$ assuming Newtonian gravity and the Cowling approximation. We find that the excited mode amplitudes grow with increasing $n_{\rm poly}$ and mode order $n$. The gravitational wave signal forms a sequence of amplitude-modulated packets for $n_{\rm poly}=1$, lasting $\sim 10^{-3}$s after each impact. The gravitational wave strain increases with increasing $n_{\rm poly}$, but decreases with increasing $n$ and increasing multipole order $l$ for $n_{\rm poly}=1$. In the observing band of current long-baseline interferometers, $g$-modes emit higher, narrower peaks in the amplitude spectral density than $f$- and $p$-modes, with the highest peaks reaching $\sim 10^{-26}$Hz$^{-1/2}$ for modes with damping time $\tau_{nl} \sim 10^{8}$yr. The root-mean-square strain $h_{\text{rms}}$, calculated by summing over modes with $2\leq l\leq4$ and $\tau_{nl} \leq 10^{8}$yr, spans the range $10^{-33} \leq h_{\text{rms}} \leq 10^{-32}$ for $1\leq n_{\text{poly}}\leq 2$.

Previous studies of neutron star oscillation modes have analysed ★ E-mail: wddong@student.unimelb.edu.au† E-mail: amelatos@unimelb.edu.auvarious excitation mechanisms, including thermal excitation by type I X-ray bursts (Heyl 2004;Watts 2012;Chambers et al. 2018), excitation during magnetar flares (Levin 2006;Levin & van Hoven 2011), resonant tidal excitation of -modes during binary inspirals (Lai 1994;Reisenegger & Goldreich 1994;Lai & Wu 2006;Sullivan et al. 2023), opacity excitation (Dziembowski & Cassisi 1999), conversion from radial to nonradial oscillations (Van Hoolst et al. 1998;Dziembowski & Cassisi 1999), and excitation of transient superfluid modes by a rotational glitch (Andersson & Comer 2001a;Sidery et al. 2010;Yim & Jones 2020).Less attention has been paid to excitation by the mechanical impact of infalling matter, including stochastically accreted matter.Numerical simulations of radially accreting quadrupolar mass shells indicate that  -modes dominate the generation of gravitational radiation (Nagar et al. 2004).Relatedly, de Araujo et al. (2006) argued that  -modes in low-mass X-ray binaries are continuously excited by thermonuclear burning of accreted matter in addition to mechanical impact.Otherwise, to the authors' knowledge, detailed calculations of neutron star oscillation modes excited by mechanical impact have not been attempted in the literature.Irrespective of the excitation mechanism, accretion-driven oscillations contribute to the phenomenon of spin stalling, when the gravitational radiation reaction torque balances the accretion torque (Wagoner 1984;Bildsten 1998;Wagoner et al. 2001;Riles 2023).
In this paper, we analyse the mechanical excitation of neutron star oscillations by stochastic accretion and calculate the GW signal thereby generated.The first goal is to predict the characteristic gravitational wave strain ℎ 0 as a function of the statistics of the accretion process, in order to estimate the signal's detectability.The second goal is to predict the waveform in the time and frequency domains, as a guide to the development of search algorithms for this signal class.The paper is structured as follows.In Section 2, we describe the neutron star model and its equilibrium state.In Section 3, we briefly review the linear theory of free stellar oscillations.In Section 4, we calculate the forced oscillatory response to the mechanical impact of accreting matter using a Green's function framework established in Appendix A. We present Molleweide surface maps for the modes that are predicted to dominate the generation of gravitational radiation.We then compute the waveform, temporal autocorrelation function, characteristic wave strain, and amplitude spectral density of the GW signal in Sections 5 and 6.In doing so, we calculate the ensemble statistics of the GW signal as a function of the statistics of the accretion flow, e.g. the size and impact frequency of accreting "clumps".Astrophysical implications are canvassed briefly in Section 7.

NEUTRON STAR MODEL
The neutron star model analysed in this paper is kept simple deliberately, in order to focus attention on the novel elements of the calculation, namely mode excitation by mechanical impact and the implication for GWs.We assume that the star is made up of a onecomponent fluid obeying a polytropic equation of state (EOS).We also assume Newtonian gravity and zero rotation.Free, nonradial oscillations of this idealised stellar model have been analysed in detail previously by many authors (Cowling 1941;Ledoux & Walraven 1958;Dziembowski 1971;Cox 1980).We neglect the corrections introduced by a rigid crust with nonzero shear modulus (McDermott et al. 1985;McDermott et al. 1988;Passamonti & Andersson 2012;Passamonti et al. 2021), multiple fluid components or superfluidity (McDermott et al. 1985;McDermott et al. 1988;Lindblom & Mendell 1994;Andersson & Comer 2001b;Prix & Rieutord 2002;Passamonti & Andersson 2012;Gualtieri et al. 2014), nonpolytropic EOS (Lindblom & Detweiler 1983;McDermott et al. 1988;Mc-Dermott 1990;Strohmayer 1993), rotation (Papaloizou & Pringle 1978;Bildsten et al. 1996;Friedman & Stergioulas 2013;Andersson & Gittins 2023), relativistic gravity (Thorne & Campolattaro 1967;Thorne 1969a;McDermott et al. 1983;Finn 1988), and magnetic fields (Carroll et al. 1986;Lee 2007;Sotani et al. 2007).These and other corrections can be accommodated within the theoretical framework set out in Sections 3 and 4, once GW detections of such sources are made.They will be important when analysing actual GW data but are premature for the order-of-magnitude detectability estimates in Sections 5 and 6.
The free and driven oscillations analysed in Sections 3 and 4 respectively are linear perturbations of a nonrotating steady state which is calculated by standard methods.We take the EOS to be polytropic (Lindblom et al. 1998), with where  is the pressure,  is the density,  is a proportionality constant, and  poly is the polytropic index.The exponent Γ =  log / log  = 1 + 1/ poly differs in general from the adiabatic index Γ 1 as discussed in equation ( 12); the latter refers to an isentropic motion without any heat exchange.The difference in Γ and Γ 1 leads to -modes among other phenomena as discussed in Section 3.2.Mass continuity implies where  () is the mass enclosed within radius .

SMALL-AMPLITUDE STELLAR OSCILLATIONS
Stellar oscillations in the Newtonian regime are well-explored (Dziembowski 1971;Cox 1980;McDermott et al. 1988).In this section, we summarise the main results from the theory of stellar oscillations as they apply to perturbations of the equilibrium state defined in Section 2. The oscillations are assumed to be linear, nonradial, and adiabatic.The equations of motion, adiabatic condition, eigenmode decomposition, and boundary conditions are presented in Sections 3.1-3.4respectively.Examples of eigenvalues and eigenfunctions for the  -, -, and -modes relevant to the application in this paper are presented in Section 3.5.

Equations of motion
The mass and momentum conservation laws and Poisson's equation for the gravitational field are given by v where v, , and Φ represent the flow velocity, viscous part of the stress tensor, and gravitational potential respectively.Let  denote Eulerian perturbed quantities, e.g.v = v 0 + v and v = / + (v 0 • ∇), where  (x, ) is the Lagrangian displacement of a fluid element from its equilibrium position.Setting v 0 = 0, we write The linearised versions of ( 4)-( 6) read v We simplify ( 8)-( 10) in two important ways in this paper.First, we neglect dissipation and set  = 0. Second, we make the Cowling approximation (Cowling 1941) and neglect ∇Φ with respect to the other terms in (9).With the two simplifications above, the equation of motion for  becomes

Adiabatic condition
The Lagrangian perturbations for pressure and density are assumed to be related by the adiabatic condition where Γ 1 is the adiabatic index (Chandrasekhar 1957).We adopt the convention that Δ denotes a Lagrangian perturbation, with Δ = +  • ∇ for a quantity .The adiabatic index Γ 1 does not necessarily equal Γ.We restrict attention to Γ 1 ≥ Γ to avoid destabilising the -modes, as discussed below.
The adiabatic, Eulerian perturbation for the pressure of an ideal fluid is Upon substituting ( 8) and ( 13) into (11), we arrive at the pulsation equation where r is the radial unit vector.Equation ( 14) agrees with equation (6) in McDermott et al. (1988), when the perturbed gravitational potential and shear modulus vanish.The adiabatic assumption is justified within the region  ≳ 10 4 g cm −3 , where the thermal conduction timescale over a pressure scale height /(|∇Φ 0 |) at a temperature of ∼ 10 8 K is longer than the period of the surface  = 1 -modes (∼ 10 −2 s) (Bildsten & Cutler 1995).For a typical polytropic EOS ( poly = 1.5), the nonadiabatic surface layer is thin enough (∼ 10 −2 cm), that its effect on the oscillation spectrum is negligible.

Eigenmode decomposition
The normal modes of oscillation of the linear system (14) exhibit a harmonic time dependence  ∝    , because the coefficients of  and its derivatives on the right-hand side of ( 14) do not depend on .In general, the eigenvalue  =  +  is a complex number, whose imaginary part, , describes mode damping or growth.Under adiabatic conditions (see Section 3.2), we have  = 0.
Two classes of oscillation modes exist, namely spheroidal and toroidal modes (Aizenman & Smeyers 1977).They have different angular structures.An inviscid, nonrotating, unmagnetised neutron star does not support toroidal motions in the Cowling approximation due to the absence of shear stress.Henceforth, we concentrate on spheroidal modes.These take the form where   (, ) denotes a spherical harmonic (Jackson 1998).
The pulsation equation ( 14) resolves into two independent, ordinary differential equations when evaluated for the ansatz (15), viz.
The oscillation spectrum of ( 16) and ( 17) consists of three modes: acoustic -modes, fundamental  -modes, and buoyancy-driven modes.We follow convention and designate modes by the spherical harmonic index  and the overtone number , which counts the number of radial nodes.We denote -modes and -modes by    and    respectively.One  -mode exists for each , denoted by   (Cowling 1941). -modes can be thought of as -modes with no radial nodes.The azimuthal harmonic  does not appear, because the modes are degenerate with respect to , when rotation and magnetic fields are absent (Perdang 1968).

Boundary conditions
We impose central and surface boundary conditions on equations ( 16) and ( 17).The central boundary condition follows from requiring that the Taylor expansions of  1 and  3 are finite at  = 0.This regularity condition reduces to The surface boundary condition requires the Lagrangian perturbation Δ to vanish at the surface.This condition reduces to We normalise the amplitude of the free oscillations to the oscillatory moment of inertia  =  *  2 * according to where the volume integral is over the entire undisturbed star,  = (, , ) labels the modes, and the superscript * denotes the complex conjugate.The orthogonality implied by equation ( 21) was first derived by Chandrasekhar (1963).The energy contained in the free oscillation labelled by  under this normalisation, scales with  2  .The amplitudes of the driven modes are calculated in terms of a Green's function for excitation by accretion-driven mechanical impact in Section 4, while the normalisation in equation ( 21) is reserved for the undriven eigenfunctions.The energy in a mode driven by the mechanical impact of accreting gas is much smaller than  2  /2 for any astrophysically plausible model of the impact process, as shown in Sections 4-6.

Numerical results for eigenmodes
We integrate equations ( 16) and ( 17) using a fourth-order Runge-Kutta scheme starting from the centre and matching (20) at the surface Table 1.Eigenmode spectrum for three EOS with 1 ≤  poly ≤ 2 and 5/3 ≤ Γ 1 ≤ 2. Columns from left to right per panel: mode designations, squared eigenfrequencies  2 in units of   * / 3 * , eigenfrequencies /2  in units of Hz, eigenfunction free amplitudes (dimensionless) at the surface   / * and  ⊥ / * , dimensionless overlap integral   , and the gravitational radiation damping time-scale   in seconds for the first few  = 2 eigenmodes.We focus on  = 2 modes as they emit quadrupole GWs, which dominate the total GW signal.The free amplitudes are normalised according to the normalisation condition (21) and are much greater than the driven amplitudes resulting from an astrophysically plausible accretion process (see Sections 4-6).Top left:  poly = 1, Γ 1 = 2; -modes are not present in this model.Top right:  poly = 1.5, Γ 1 = 5/3; -modes are not present in this model.Bottom left and right:  poly = 2, Γ 1 = 5/3; -modes are present.Fixed stellar parameters:  * = 1.4 ⊙ ,  * = 10 4 m.using a Newton-Raphson shooting scheme.Eigenvalues in the shortwavelength limit (see below) serve as initial guesses.The approach is standard except for one subtlety: the 1/ r terms in equations ( 16) and ( 17) and the vanishing pressure and density at the surface cause  1 and  3 to diverge at r = 0 and r = 1.The divergences are resolved by shifting the boundaries slightly.Specifically, we apply equations ( 19) and ( 20) at r = r1 = 10 −10 and  =   = 10 7 g cm −3 respectively, noting that other physics (e.g. the magnetic stress) becomes important at  <   (Bildsten & Cutler 1995) and is neglected in equations ( 16) and ( 17).We verify numerically that the eigenvalues computed for   = 10 7 g cm −3 differ by one part in 10 6 from those computed for   = 10 4 g cm −3 (say) as a comparator, confirming that the results are insensitive to the choice of   ≤ 10 7 g cm −3 .Usually, about ten iterations are enough to reach the targeted accuracy  = |Δ|/|| ≤ 10 −11 for  -and -modes, and  = |Δ σ| ≤ 10 −11 for -modes, where Δ (Δ σ) is the change in  ( σ) between successive iterations.We choose  unusually small in order to handle the cancellation arising from the oscillatory integrand in the overlap integral introduced in Section 5 [see equation ( 43)] and Appendix B. Table 1 lists the eigenfrequencies calculated for -,  -, and modes with different EOS, in order to validate the above procedure against existing results in the literature (e.g.Robe 1968, Table 2; note the unit difference). -modes oscillate at ≳ kHz depending on the EOS.Most -modes oscillate at ≳ 10 kHz, except for the first overtone.Less centrally condensed stars (with lower  poly ) have lower frequency  -modes and higher-frequency overtones.The frequencies of the first ten -modes run from ∼ 0.4 kHz to 1 kHz.Note that the frequency of a  -mode calculated under the Cowling approximation can be 50 per cent higher than the exact value, whereas the frequencies of -modes are nearly unchanged (Robe 1968).Therefore, all the -modes listed in the table, some  -modes, and the first overtone of some -modes fall within the LVK observing band (Finn & Chernoff 1993;Buikema et al. 2020).The free mode amplitudes   and  ⊥ in Table 1 are normalised according to (21) and are much greater than the driven amplitudes resulting from an astrophysically realistic accretion process, as shown in Sections 4-6.
In the short-wavelength limit, we can obtain an analytic expression for the dispersion relation at high orders  or  (Finn 1988) through a local analysis, assuming  ∼   and   * ≫ 1 (McDermott et al. 1983(McDermott et al. , 1988)).The -and -modes satisfy respectively The sign of  2 determines the stability of -modes via (24).For a polytropic star, one finds  2 ∝ Γ 1 − Γ < 0 for Γ 1 < Γ, leading to instability.The values of  2  ∝ Γ 1 − Γ in the bottom right panel of Table 1 are relatively low, because one has Γ 1 ≈ Γ typically in a neutron star.The slowly varying wavenumber satisfies a stationarity condition (Lai 1994;Aerts et al. 2010), where  is a constant of order unity arising from the phase change at the boundaries, whereupon ( 23) and ( 24) take the form Real neutron stars are viscous.Viscosity modifies the dispersion relation significantly for modes above a certain order.This happens for ∇ 2 ∼  / , i.e.  2 /(4 2 ) ∼ 1 where  is the shear viscosity.For a star with temperature  * = 10 9 K, viscosity modifies -modes with order  ≳ 10 6 (/10 25 g cm −1 s −1 ) and -modes with order  ≳ 10 4 (/10 25 g cm −1 s −1 ).One gets similar estimates for the spherical harmonic order  (Cutler & Lindblom 1987;Friedman & Stergioulas 2013).We show in Section 5 that modes of such high order emit negligible power in GWs.

EXCITATION BY ACCRETION
In this section, we calculate the oscillatory response of the star to the mechanical impact of gas accreting onto the stellar surface.The response is a linear combination of the eigenmodes in Section 3, whose coefficients are determined by the near-surface condition at the point of impact of the accreting gas.We formulate the inhomogeneous boundary value problem in terms of a Green's function in Section 4.1.We then solve for the response to the impact of a single clump of gas (Section 4.2) and a stochastic sequence of equally sized clumps with random arrival times (Section 4.3).

Green's function
We address the inhomogeneous excitation problem where L is the spatial operator in the pulsation equation ( 14) (see also Appendix A) and F  (x, ) is the -th Cartesian component of the mechanical force density from accretion.Equation ( 28) has the general solution where the subscript  denotes the -th Cartesian component, the Einstein summation convention on repeated indices is adopted, and   is the Green's function derived in Appendix A. In (30),  (•) denotes the Heaviside step function, and  sums over all modes.Equation ( 30) is equivalent to the general solution derived by Schenk et al. (2001), viz.equation (2.24)-(2.26) in the latter reference.The solution in Schenk et al. (2001) follows from the mode decomposition formalism and does not introduce a Green's function explicitly.The formalism applies to both rotating and nonrotating stars.
In this paper, we analyse the situation, where the accretion occurs stochastically through a sequence of discrete "clumps".Hence we model the force density from the mechanical impact of a clump at the surface phenomenologically by where   is -th Cartesian component of the momentum of the clump, x f is the impact point on the boundary of the star (taken to be the same for every clump for simplicity), and T () specifies the temporal history of the impact, including its duration.Strictly speaking, the idealised neutron star model in this paper has a vanishing mass density at the surface; that is, there is nothing at the surface that feels the force of the impact.This artificial feature is resolved self-consistently by moving the outer boundary condition inward slightly, as when applying the zero-pressure boundary condition to the eigenfunctions in Section 3.5, so that the star feels the force from the clump, when the clump reaches  =  B (see Section 3.5).

Single clump
We start by considering the impact of a single clump.If the impact occurs, such that the force density remains constant throughout the interaction between the clump and the stellar surface, we can write where  is the duration of a collision starting at  s .Upon substituting (32) into ( 29) and ( 30), we obtain for  s ≤  <  s + , and for  >  s + .In ( 33) and ( 34), we follow standard practice (Echeverria 1989;Kokkotas et al. 2001;Ho et al. 2020) and incorporate phenomenological damping factors ∝ exp[−( −  s )/  ] to capture the back reaction from gravitational radiation.A self-consistent treatment of gravitational radiation reaction lies outside the scope of this paper (Wheeler 1966;Thorne 1969a,b).We assume that the gravitational radiation reaction operates slowly over many mode cycles, viz.    ≫ 1, in order to simplify (33) and (34).Typical gravitational damping timescales (see Table 1 and Section 5.4) for  -modes and2  1 -modes are ∼ 10 −2 s; for -modes, we have ≳ 10 2 s.
The factor p •  ( ) (x f ) in ( 33) and ( 34) implies that the amplitude of an excited mode depends on the angle of infall. -and -modes are more strongly excited when the clump strikes the surface radially, while -modes are strongly excited, when p has a sizable tangential component.For a self-consistent nonrotating picture, we take p ∝ r, i.e. p is always radial, unless otherwise specified.
We see from equation (34) that the oscillation resembles that excited by a delta-function impulse ( −  s ) for    ≪ 1 and is suppressed for    ≫ 1. Physically, this occurs because the top-hat impact (32) averages over high-frequency modes.The height of the accretion column  is assumed to be small ( ≲ 1 km) for a neutron star with  * = 1.4 M ⊙ and  = 10 4 m (Mushtukov et al. 2018).Accretion columns with  ∼  * are possible (Poutanen et al. 2013) but are not considered in this paper.Taking  ≳ / (see Appendix C), modes with frequencies   ≳ 10 5 rad s −1 are suppressed, i.e. (  ) −1 ≲ 1.However, with top-hat averaging, the system prefers to excite high-order modes, because p ) increases with ascending order.This is because higher- modes require less energy to reach the same amplitude than lower- modes, as higher- modes involve larger motion in the surface layer, which is weakly bound gravitationally.Preferential excitation of high-order modes is also present in slowly pulsating B-type stars (Dupret et al. 2008;Aerts et al. 2010), as well as in Sun-like stars that are excited by turbulent convection (Dziembowski et al. 2001;Chaplin et al. 2005). 2 A more realistic collision model, e.g., the turbulent plume interaction studied by Colgate & Petschek (1981), would involve extended radial and angular profiles of the force density, smooth (x − x f ) spatially in (31), and average spatially the shorter-wavelength modes through (29).To capture this consequence qualitatively without complicating the calculation unnecessarily, we cut off the  sum in ( 33) and (34) at || ≤  c ≈ 5.The cutoff at || ≈ 5 is motivated by the computational cost required to compute the mode sum in equations ( 33) and (34) to the desired accuracy.The overlap integral [equation ( 44)] decreases monotonically with .At the same time, the grid resolution required to calculate equation ( 44) increases with .Therefore, above a certain  value, additional orders do not contribute to the desired accuracy of the mode sum.In this paper, the desired accuracy is 0.1%, which implies a cutoff satisfying 3 ≲  ≲ 7, depending on the physical parameters in the problem.
Molleweide plots of the evolution of  ( =  * , , , ) at a sequence of four snapshots in time are presented in Figure 1.The snapshots are timed to illustrate the spreading of the oscillation.The oscillation amplitude is relatively small at the beginning as it is not fully excited before the impact ends.The surface is pushed inwards in response to the impact and bounces back within a few multiples of the collision duration.The striped response pattern, for 2 ≤  ≤ 4, is always centred on the impact point.The (dimensionless) radial surface displacement, | surf |/ * =  (, x,  s , )/ * for |x| =  * , during the spreading motion is bounded by where  is the mass of the clump, and  is its infall speed.The single-clump response is to be compared with the typical ellipticity ≲ 10 −8 predicted theoretically to arise from steady-state elastic and hydromagnetic mountains on the stellar surface (Ushomirsky et al. 2000;Melatos & Payne 2005;Haskell et al. 2006;Priymak et al. 2011).

Stochastic sequence of clumps
In reality, accretion onto neutron stars is a stochastic process.Complicated and possibly chaotic dynamics are imprinted on the mass accretion rate by continuously excited, three-dimensional, hydromagnetic instabilities at the disk-magnetosphere boundary (Romanova et al. 2003;Romanova & Owocki 2015).In this subsection, we extend the discussion about a single clump in Section 4.2 to a random sequence of clumps.In the idealised phenomenological model described by ( 31) and ( 32), there are several quantities which may be converted into random variables: the clump infall velocity v, the clump mass , the position x f where the accreted matter falls, and the duration of the collision .The randomness in v could come from turbulent or thermal motion inside the clump, three-dimensional twisting of the accretion column [e.g.magnetic Rayleigh-Taylor fingers; see Romanova et al. (2003)], or "bouncing" of the stellar surface in response to previous clumps.However, for an accretion rate  ≈ 1 × 10 −8  ⊙ yr −1 with mean clump impact frequency matching the LVK band,  acc ∼ kHz, the speed of the oscillating surface,  surf ∼ |v| | | 2 /(2  acc  ) ≈ 10 −15 , and the thermal velocity,  drift ≈ 10 −3 ( * /10 7 K), are negligible compared to the clump's terminal velocity, which satisfies |v| ∼ .Moreover, if the stellar magnetic field guides the clumps onto the magnetic polar cap, the impact zone is comparable to or smaller than the wavelength for most of the modes in Table 1 (e.g. with  ≲ 5), and the factor (x − x f ) with fixed x f in (31) serves as a fair approximation. 3 The main motivation for choosing  acc ∼ 1 kHz is to match the observing band of ground-based detectors.This value is broadly consistent with what is known about clumpy accretion streams.We expect the mean difference in the arrival times of consecutive clumps,  −1 acc , to be approximately the fragmentation period of the accretion stream.In turn, the fragmentation period is comparable to the Kepler period at the disk-magnetosphere boundary (at radius   ), under a range of plausible physical conditions (Romanova & Owocki 2015).For example, one obtains  −1 acc ≈ 1ms, for   ≈ 2 * .Simulations of matter fragmenting as it flows through the magnetosphere with 10 * ≤   ≲ 100 * at different accretion rates show variability on time-scales between ∼ 10 ms and ∼ 100 ms (Romanova & Owocki 2015;Mushtukov et al. 2024, e.g., Fig. 30 in the first reference).
In this paper, for the sake of simplicity, we assume that every clump has the same mass, and every impact has the same duration.We do so, because X-ray observations of accreting neutron stars are unable to resolve the impacts of individual clumps and hence measure  and  ∼  −1 acc , where the latter scaling assumes a quasi-continuous flow.Fortunately, many of the order-of-magnitude estimates in Sections 5 and 6 depend on  and  acc chiefly through the time-averaged accretion rate  =   acc rather than  and  acc separately.It is therefore important to choose  and  acc consistently with astrophysically realistic values of .With the above approximations, we allow the time of impact  s to fluctuate randomly, i.e.
where { s } is a sequence of impact epochs following a Poisson counting process, with mean frequency  acc , and  () is the total (and random) number of impacts up to time .
For a canonical polytropic neutron star with  poly = 1.5 and Γ 1 = 5/3 accreting at  acc ∼ kHz, the upper bound of  surf estimated using the Cauchy-Schwarz inequality is ( surf / * ) max =  clumps  surf / * , where  clumps =  acc  is the mean number of clumps striking the surface during a mode-weighted average damping time .The typical value of the amplitude of  surf / * excited randomly by multiple clumps is expected to scale with  clumps 1/2 and can be written as Equations ( 37) should be compared to equation ( 35) for a single clump, with  = /  acc .
Figure 2 shows one realisation of the evolution of  ( * , , , ) in a sequence of four snapshots (up to 8.3 ms) in response to eight stochastic impacts obeying (36).In Figure 2, modes excited shortly before the second snapshot are out of phase with modes excited before the first snapshot.We see that the amplitude is lower during the impact and the oscillation spreads out.The third snapshot at the bottom left highlights the radially inward displacement around the impact point, shortly after the latest excitation.The last snapshot displays evidence of destructive interference of modes excited by randomly timed clumps.The upper bound of  surf / * in Figure 2 agrees roughly with equation (37) for  clumps = 8.

GRAVITATIONAL WAVEFORM
In this section, we calculate the gravitational waveforms emitted by stellar oscillations excited by single and multiple clumps as in Section 4. In Section 5.1, we calculate the mass quadrupole GW strain emitted by an arbitrary stellar mode (labelled by , as above) in terms of an overlap integral, whose integrand is proportional to  ( ) multiplied by   .In Section 5.2, we calculate the waveform generated by four types of stochastic impact: single delta function, single top hat, multiple periodic top hats, and multiple Poisson top hats.Energy spectra of GW signals generated by multiple Poisson top-hat impacts with various durations are presented in Section 5.3, in order to identify the energetically dominant modes.The GW damping timescale for each mode is calculated a posteriori in Section 5.4.

Strain
The metric perturbation ℎ TT   in the transverse traceless gauge measured by an observer at a distance  generated by the time-varying mass multipoles of a Newtonian source is given by (Thorne 1980) with where  is the retarded time,  is the Eulerian density perturbation from equation ( 8), and  E2,   is the gravitoelectric tensor spherical harmonic which describes the angular radiation pattern.Only the mass multipoles emit GWs in a nonrotating and inviscid neutron star, because the current multipoles ∝ ∇ ×  0  ( ) have no radial component for spheroidal modes.
The time derivative in (38), after integrating (39) by parts, can be written as plus the time derivative of a surface integral which vanishes for ( * ) = 0.In the above notation, we have which equals the mode amplitude excited by a clump striking at time  s .We also have while   is the dimensionless overlap integral (Press & Teukolsky 1977;Lai 1994;Reisenegger & Goldreich 1994;Kokkotas & Schäfer 1995;Sullivan et al. 2023) If the accreting gas is guided magnetically onto the magnetic polar cap, then the impact footprint is approximately axisymmetric, and modes with  = 0 dominate. 4In what follows, we present results for the special case of even multipoles (where  ≥ 2 is an even integer) to illustrate by way of example some of the diverse waveforms that are possible.With     ≫ 1 and writing ℎ TT   () = ℎ 0 () E2,0   , we obtain Interestingly, for  = 2,   appears explicitly in the phase but not the amplitude of ℎ 0 (), for a top-hat impact.However, the amplitude does depend on  2 implicitly.One can see this by replacing the sum over  with its continuous analog (i.e. an integral), whereupon one obtains ℎ 0 ∝   −1 2 × (terms depending on ) approximately.

Types of impact
In this subsection, we evaluate ℎ 0 () for an illustrative selection of the impact models in Section 4: a single delta-function in time, a single top hat, a sequence of top hats which is periodic, and a sequence of top hats which is Poisson (uniformly) distributed in time, following (36).
Figure 3 shows the dominant  = 0 quadrupolar gravitational waveforms from accretion-induced oscillations.Each quadrant (a)-(d) of the figure contains four subpanels, which display ℎ 0 () versus  for the four kinds of impact defined in the previous paragraph, ordered from top to bottom.The top two quadrants (a) and (b) correspond The waveforms produced by a top-hat impact and a deltafunction impact are similar in shape and amplitude.There is a burst of activity lasting ≲ 10 −3 s after the impact, during which an amplitudemodulated "packet" of shorter-period modes is excited.The packet evolves into a longer-period oscillation for  ≳ 10 −3 s.The similarity between delta-function and top-hat impacts arises, because both are shorter than the periods of modes that dominate the GW emission.
The strain from periodic impacts is only ∼ 2 times greater than the strain produced by one impact, because all modes have frequencies different from  acc and therefore are not excited resonantly.Periodic impacts form amplitude-modulated packets, whose envelopes oscillate with period ∼ 10 −2 s, whereas the waveform produced by Poisson impacts is shaped irregularly.The strain from Poisson impacts is 2-4 times greater than from periodic impacts due to the randomness in  s , which disrupts the cancellation involved in non-resonant periodic excitation.
The GW signal from -modes, plotted separately for clarity in quadrant (d), is considerably weaker (by about four orders of magnitude) than from  -and -modes for three reasons.First, -modes have smaller   because of more complete cancellation in the oscillatory integrand in (44).Second, the accreting material is assumed to strike the surface radially, and -modes have smaller r •  * ( ) (x f ) than modes.Third, the GW signal is proportional to the time derivatives of   , and -modes oscillate slower than  -and -modes.
Polytropes with stiffer EOS lead to smaller mode amplitudes for the same impact.The overlap integral   is also smaller (except for  -modes), because the wave strain contribution from expanding and contracting mass shells cancel more completely in a less centrally condensed star.Hence, polytropes with smaller  poly emit weaker GW, as confirmed by Figure 3.We find ℎ 0 ≈ 1 × 10 −33 , 5 × 10 −33 , and 2 × 10 −32 in quadrants (a), (b), and (c) respectively.Multipole order plays a significant role in setting the shape and amplitude of the emitted waveform.The discussion in the previous paragraphs applies principally to the  = 2 multipole, which is drawn as a blue curve throughout Figure 3.We overlay the  = 3 and  = 4 signals as orange and green curves respectively.The "packets" produced by periodic impacts have similar periods for  = 3 and  = 4.The single-impact waveforms show that  = 2 dominates usually.However, the quadrupole GWs emitted by  -and -modes are damped until they reach similar magnitude to higher orders within ≈ 3  for 2  -modes (∼ 0.02 ms).The amplitude of the strain emitted by  -and -modes from  = 3 and  = 4 multipoles is an order of magnitude smaller than from  = 2 for  poly = 1.0 and of the same order for  poly = 1.5 and 2. The quadrupole GWs emitted by -modes always dominate over higher-order multipoles.

Dominant modes
In this subsection, we explore the energy spectra of GW signals generated by multiple Poisson top-hat impacts with different , in order to compare with the simulated GW signals from a single, infalling, quadrupolar mass shell in Nagar et al. (2004).
Figure 4 2004) is ∼ 10ms.Figure 4a shows that  1 contains most of the energy during the initial 10 ms for  = 10 −5 s, whereas the  -mode dominates in the simulations by Nagar et al. (2004).The difference arises because -modes with  ≥ 2 are suppressed for  = 10 −5 s through top-hat averaging, except for the  1 mode.In Figure 4b, the  -mode emerges as the highest peak, because the  1 mode is suppressed for  = 10 −4 s, which is of the same order of magnitude as  in Nagar et al. ( 2004) (see the text above equation ( 43) in the latter reference).
Another essential distinction between this paper and Nagar et al. ( 2004) is that this paper focuses on long-lasting accretion, while Nagar et al. (2004) applies to transient core-collapse scenarios.In Figure 4c, we plot  (2,0) / for a longer spectral window of 100 ms. 1 re-emerges with the highest peak.Figure 4b and 4c imply that more energy accumulates in higher- modes, as the spectral window lengthens (see Section 6.1).

Damping timescale
A variety of damping mechanisms exist in realistic neutron stars, including but not limited to gravitational radiation reaction, neutrino emission, nonadiabaticity, internal friction, and electromagnetic radiation.Nonadiabaticity and viscosity are insignificant as explained in Section 3. Electromagnetic damping and neutrino damping can be important for certain modes (McDermott et al. 1988).However, gravitational radiation reaction dominates the damping for the  -,-, and -modes we study in this paper.
We estimate the damping timescale a posteriori via   = 2  / gw , where  gw is the power of the gravitational radiation, a standard approach (Thorne 1969b;Balbinski & Schutz 1982;Mc-Dermott et al. 1988;Owen et al. 1998).In our notation, the -folding timescale reads The damping modification of   is of order O  −2   −2  and is therefore insignificant for     ≫ 1, as typically occurs.Numerical results for the damping timescale are given in Table 1.  ∝  −2  grows exponentially with  because   decreases exponentially with  as demonstrated in Appendix B. Therefore, gravitational radiation reaction does not damp high- modes efficiently.

DETECTABILITY
In this section, we calculate the temporal autocorrelation function of the stochastic GW strain calculated in Section 5 with the aid of the analytic formulation in Appendix D. The autocorrelation function is important for two reasons in this paper: (i) when evaluated at zero time lag, it returns the mean square amplitude of the GW signal, a key quantity when assessing detectability; and (ii) when Fourier transformed, it returns the power spectral density of the GW signal, which contains important information about the source physics, e.g. the modes excited by accretion, and the statistics of the accreted clumps.In Section 6.1, we calculate the temporal autocorrelation function of ℎ 0 ().We focus on the autocorrelation of modes with themselves and de-emphasise the cross-correlation between different modes, which dephase due to frequency mismatch over time lags that are relevant to GW observations.The detectability of the GW signal is investigated in both time and frequency domains, via the rootmean-square wave strain and the amplitude spectral density (ASD) in Sections 6.2 and 6.3 respectively.

Temporal autocorrelation function
We define the autocorrelation function as where angular brackets denote an ensemble average.For the sake of simplicity, we focus on the dominant multipole (, ) = (2, 0).For a stationary process, we have  (,  ′ ) =  (| |), where  =  ′ −  is the time lag between measurements.We find with r() =    () ( ) .The indices (, ) = (2, 0) are suppressed for convenience in   ,   , and   for the remainder of this subsection.The statistical properties of   ()-like random variables, including their autocorrelation, are derived in Appendix D. From Appendix D, we have The double integral in the first line of ( 49) is zero as it involves the derivatives of constant mean amplitudes.The second line in ( 49) is negligible, except when two modes have the same frequencies.Integrating the second line in (49) with  =  ′ , we isolate the dominant term Equation ( 50) suggests that the amplitude autocorrelation increases with , because   increases with ; that is, higher-order modes accumulate more energy from multiple impacts before being damped.
Upon substituting ( 46) and ( 50) into ( 48) and neglecting the cross terms,  () simplifies to with   = 1/ √︁ 1 − (|v|/) 2 .Interestingly,  () does not depend on   , because   cancels   and we have  2 / 2 ∝  2 .The cancellation occurs when GW damping dominates; a different scaling applies for other damping mechanisms, to be studied in future work.We redefine the mode cutoff  c to incorporate this   -independence, ensuring   ≲  ⊙ /  (e.g. c ≈ 15 for  poly = 1.0).Modes with  >   generate negligible ℎ 0 because the accumulation of energy from multiple impacts cannot compensate for the exponentially decaying   .To understand this concept mathematically, consider the rate of change of the energy   in mode , viz.
where  in  is the input power.Consider a higher-order mode A with  A >  0 and a lower-order mode B with  B <  0 , where  0 is the lifetime of accretion episode.The ratio of the emitted GW power at Since the damping time-scale increases exponentially with increasing mode order , and  in does not increase as fast, we have ( in A / in B )( 0 / A ) ≪ 1.When other damping mechanisms (e.g., neutrino emission or viscosity) are considered, the power of the signal is scaled by ∼  total / GW , with 1/ total = 1/ GW + 1/ other , where  GW and  other are the damping time-scales via gravitational radiation and another mechanism, respectively.Therefore, modes with  other ≪  GW <  0 also contribute negligibly to the long-term signal.Unlike (34), the factor | () |/( * σ2  ) in (52) increases with  for a relatively soft EOS (e.g. poly = 1.5 or 2.0), but decreases with increasing  for a relatively stiff EOS (e.g  poly = 1.0).

Root-mean-square wave strain
As a first pass, we assess the detectability of the GW signal in the time domain.We calculate the root-mean-square characteristic wave strain ℎ rms in terms of the temporal autocorrelation function at zero lag, viz.ℎ rms =  ( = 0) 1/2 and hence Equation ( 53) follows immediately from ( 52), except that the sum in the last line is over all  ≥ 2, not just  = 2.Note that equation ( 53) implies ℎ rms ∝  −1 → ∞ if one naively takes the limit  → 0. The divergence occurs, because one has T ∝  −1 nominally in Equation ( 36) for a top-hat impact.The divergence disappears, if one directly computes the ordered limit lim  →0 lim →0  (, ) in ( 51).Note also that  * ,  * and σ are related through the EOS.For example, one must simultaneously scale  * and  * according to the EOS, in order to leave σ unchanged.Equation ( 53) implies that the GW signal emitted by accretiondriven stellar oscillations is too weak to be detected by the current generation of LVK interferometers (Watts et al. 2008;Haskell et al. 2015).The value of ℎ rms in ( 53) is lower than the upper bounds ℎ 95% 0 ∼ 10 −26 (at 95% confidence) inferred from recent searches of continuous waves from accreting systems (Middleton et al. 2020;Zhang et al. 2021;Covas et al. 2022;LIGO Scientific Collaboration et al. 2022a,b).
radiation torque in the mechanical excitation scenario is negligible, compared to the accretion torque.

ASD
A complementary approach to assessing the detectability of the GW signal is to compute its ASD,  ℎ (  ) 1/2 , in the frequency domain and compare it with the ASD of the LVK interferometers.This is important, as some of the signal power ∝ ℎ 2 rms falls outside the LVK observing band and cannot be detected.Conversely, predicting the ASD theoretically becomes important, once a signal is detected, because the locations of the spectral peaks contain important information about the stiffness of the EOS.For example, high-order modes are prominent, if the EOS is relatively soft.
The one-sided ASD is related to the Fourier transform of  () by where  > 0 is the Fourier frequency in Hz.Examples of the ASD are plotted in Figure 5 for three EOS with  poly = 1.0, 1.5, and 2.0 (top to bottom subpanels), and modes with  = 2 and  ≤   in the upper LVK observing band.The peaks reach  ℎ (  ) 1/2 ∼ 10 −34 Hz −1/2 for  -and -modes, and 10 −32 Hz −1/2 ≲  ℎ (  ) 1/2 ≲ 10 −25 Hz −1/2 for -modes.For  poly = 1.0 and 1.5, only the  -mode falls within the calibrated observing band6 (10 Hz ≲  ≲ 5 kHz), and there are no -modes (Γ 1 = Γ).For low-order modes ( ≲ 5), the predicted ASD drops below the sensitivity curve of the current generation of long-baseline interferometers, e.g. the LVK detectors (Cahillane & Mansell 2022); see Figure 4 in the latter reference.However, the ASD increases with mode order in general, mainly because   increases with mode order.Moreover, -modes generally attain larger  ℎ (  ) 1/2 than  -and -modes because they have higher   and lower   .The peaks in Figure 5 grow narrower, as  increases.
Figure 5 demonstrates how various approximations, made for computational convenience, affect the ASD.Each panel displays three curves for quadrupolar GWs: the exact ASD, that includes all terms and cross-correlations in  () (dashed black curve); the ASD that excludes cross-correlations in  () (dash-dot purple curve); and the ASD that includes only the dominant terms of  (), i.e. equation (52) (solid blue curves).We see that the solid blue curves align well with the dash-dot purple curves, except for  ≲ 1 kHz.The cross-correlation terms also do not modify  ℎ (  ) 1/2 much near the peaks (  ≈   ), which implies that ( 52) is a good approximation.They do modify  ℎ (  ) 1/2 between the -mode peaks, e.g. for 5 × 10 3 Hz ≲  ≲ 10 4 Hz in the middle panel of Figure 5.This occurs because we have ⟨  ()   ′ ( ′ )⟩ > ⟨  ()   ( ′ )⟩ for  ′ ≠  and | ′ − | ≫   (see Section 6.1).
Figure 6 presents a magnified comparison between the sensitivity curve of a typical LVK detector during the third observing run and the ASD of higher-order ( ≥ 6) -modes for  = 2, 3, and 4, and  poly = 2.0.The peaks of  ℎ (  ) 1/2 approach the sensitivity curve, as  ≤   increases, but are bounded by  ℎ (  ) 1/2 ≲ 10 −26 Hz −1/2 .However, the peaks are narrow, so their signal-to-noise ratio (SNR) ∝ ∫   ℎ (  ) 1/2 is low, and they remain difficult to detect, in keeping with Section 6.2.A rigorous calculation of the SNR requires information about the specific detection algorithm.It goes beyond the scope of this paper and will be tackled in future work.

CONCLUSION
In this paper, we study the gravitational radiation emitted by nonradial oscillations of a neutron star excited by the mechanical impact of stochastically accreting clumps of matter.We calculate analytically the star's oscillatory response to a single impulsive, delta-function clump in terms of a Green's function and hence calculate the gravitational wave signal emitted by multiple delta-function or top-hat clumps, which recur periodically or with Poisson statistics.The mode calculations are done for a nonrotating, unmagnetised, onecomponent star with a polytropic EOS assuming Newtonian gravity and the Cowling approximation and reproduce standard results in the literature for the eigenfrequencies and eigenfunctions.The outputs of the calculations include gravitational waveforms as functions of the EOS and type and duration of impact (Figure 3), which will help to guide future work on detection algorithms, as well as predictions of the signal's amplitude and spectral properties, such as the temporal autocorrelation function [Equation ( 52)], root-mean-square wave strain [Equation ( 53)], and ASD (Figure 5 and 6), which are important for assessing detectability.
The main quantitative results are summarised as follows.(i) The gravitational waveforms take the form of sequences of amplitudemodulated packets, which last ∼ 10 −3 s for  poly = 1 after each excitation.Poisson impacts lead to more irregular waveforms and slightly increase the signal strength compared to periodic impacts.(ii) The signal amplitude ℎ rms increases with  poly and decreases with  for  poly = 1.Summing over modes with damping time   ≲ 10 8 yr and multipole order 2 ≤  ≤ 4, ℎ rms ranges from 10 −33 to 10 −32 for 1 ≤  poly ≤ 2. (iii) The  -mode dominates the energy spectral density  (2,0) / in simulations by Nagar et al. (2004) of transient, infalling, quadrupolar mass shells during supernova core collapse.In this paper, where the modes are excited continuously and stochastically by accretion, higher- modes dominate, as the spectral window lengthens.(iv) -modes produce higher and narrower peaks than and -modes in the ASD in the LVK observing band.The highest peaks ∼ 10 −26 Hz −1/2 approach the sensitivity of LVK detectors for modes with   ∼ 10 8 yr.(v) Overall, the signal is about four orders of magnitude weaker than LVK detectors can detect at present.
The calculations in this paper are idealised.They ignore rotation, magnetic fields, general relativistic corrections, and multiple stellar components (e.g., crust, superfluid).To further guide searches involving real data, it is worth extending the theory to treat a realistic EOS, model the size and location of the impact region realistically, and consider toroidal oscillations such as -modes (Lindblom et al. 1998;Owen et al. 1998).

APPENDIX B: NUMERICAL EVALUATION OF THE OVERLAP INTEGRAL
The overlap integral   in (43) decreases quickly as  and  increase.The integrand in (43) oscillates rapidly for large  and , making it difficult to compute   accurately.Moreover, as noted by Reisenegger (1994), small mixing of  -modes in the numerical computation of eigenfunctions can introduce errors in   , even if the eigenfrequencies are calculated accurately.If the numerical eigenfunction  () acquires a numerical error of the form with | | ≪ 1 and | 0 | ≪ 1, the leading-order error in computing   of ξ () is ∼  0 , because  -modes mostly dominate with   ≈ 0.5 (except in the most centrally condensed stars).Therefore, one needs  0 /|  | ≪ 1, in order to evaluate   accurately (Reisenegger 1994).One way to estimate  0 is to check the orthogonality conditions (Lai 1994).In this work, we find  0 ∼ 10 −9 , which is of the same order of magnitude as the precision in solving the Lane-Emden equation.
The three panels in Figure B1 are plotted on log-linear axes.We find that   decreases approximately exponentially with  for 1 ≤  poly ≤ 2, in line with the trend for modes with  ≤ 4 in Passamonti et al. (2021) for example.

APPENDIX C: INTUITIVE SPRING-BASED MODEL OF IMPACT
In this appendix, we develop intuition about the energy and momentum transfer that occur, when a clump of matter strikes the star's surface, by analysing a toy model of an idealised, double-spring system.Consider a heavy spring (the neutron star) with variable length  and fixed spring constant , attached to a light spring (the clump) with initial length  and spring constant .For simplicity, we ignore the length dependence in , which describes the deformation of the clump during an impact.We model the system as in Figure C1, where the springs are massless but are attached to blocks of mass  and  ≪ , corresponding to the star and clump respectively.The selfgravity of  is neglected.The clump initially has a speed  and experiences an acceleration of magnitude  towards .
Let  1 denote the variable length of the light spring and  2 denote the displacement of  from it equilibrium location.The equations of motion for this system are with initial conditions  1 (0) = ,  1 (0) = −,  2 (0) = 0, and  2 (0) = 0.An overdot denotes a derivative with respect to time.

Figure 1 .
Figure 1.Molleweide plot for the (dimensionless) radial displacement of the stellar surface at different times in response to a hypothetical clump with mass  = 6 × 10 11 kg striking with speed |v| = 0.4 at  = 0.The collision duration  is set to be  = 10 −5 s.The red and blue regions indicate outwards and inwards radial motions respectively.Dark colours indicate greater displacement (absolute value), as per the colour bar.The impact point, marked with a cross, is at   = /3,   = 0.  -and -modes are summed for 2 ≤  ≤ 4. Damping is not included.Top left: Snapshot at 2 × 10 −6 s before the collision ends.Top right: Snapshot at 10 −5 s exactly when the collision ends.Bottom left: Snapshot at 3 × 10 −5 s, three times the collision duration.Bottom left: Snapshot at 5 × 10 −5 s, five times the collision duration.Stellar parameters:  * = 1.4  ⊙ ,  * = 10 km,  poly = 1.5, and Γ 1 = 5/3.
displays the energy spectral density  (,) / of the GW signals emitted by oscillations excited by a sequence of Poisson top-hat impacts for (, ) = (2, 0).The time span used calculate  (,) / by Fourier transforming ℎ 0 () is termed the spectral window. 5In panels (a)-(c), we present, from top to bottom,  (2,0) / for (a)  = 10 −5 s with a spectral window of 10ms, (b)  = 10 −4 s with a spectral window of 10ms, and (c)  = 10 −4 s with a spectral window of 100ms respectively.The spectral window used by Nagar et al. (

Figure B1 .
Figure B1.Overlap integral |  | for modes with  poly = 1.0 (top left panel), 1.5 (top right panel), and 2.0 (bottom panel).The blue, orange, and green dots are for multipole orders  = 2, 3, and 4 respectively.The trends in every panel are exponential (log-linear axes).Every second -and -mode in the bottom panel is not labelled for readability.

Figure C1 .
Figure C1.Schematic illustration of the impact of an accreted clump on the stellar surface, modelled as an idealised mechanical system comprising two springs (spring constants  and ) connecting two masses ( and  ≪ ).