GRAVITATIONAL RADIATION FROM HYDRODYNAMIC TURBULENCE IN A DIFFERENTIALLY ROTATING NEUTRON STAR

and

Published 2009 December 29 © 2010. The American Astronomical Society. All rights reserved.
, , Citation A. Melatos and C. Peralta 2010 ApJ 709 77 DOI 10.1088/0004-637X/709/1/77

0004-637X/709/1/77

ABSTRACT

The mean-square current quadrupole moment associated with vorticity fluctuations in high-Reynolds-number turbulence in a differentially rotating neutron star is calculated analytically, as are the amplitude and decoherence time of the resulting, stochastic gravitational wave signal. The calculation resolves the subtle question of whether the signal is dominated by the smallest or largest turbulent eddies: for the Kolmogorov-like power spectrum observed in superfluid spherical Couette simulations, the wave strain is controlled by the largest eddies, and the decoherence time approximately equals the maximum eddy turnover time. For a neutron star with spin frequency νs and Rossby number Ro, at a distance d from Earth, the root mean square wave strain reaches hrms ≈ 3 × 10−24 Ro3s/30 Hz)3(d/1 kpc)−1. Ordinary rotation-powered pulsars (νs ≲ 30 Hz, Ro ≲ 10−4) are too dim to be detected by the current generation of long-baseline interferometers. Millisecond pulsars are brighter; for example, an object born recently in a Galactic supernova or accreting near the Eddington rate can have νs ∼ 1 kHz, Ro ≳ 0.2, and hence hrms ∼ 10−21. A cross-correlation search can detect such a source in principle, because the signal decoheres over the timescale τc ≈ 1 × 10−3 Ro−1s/30 Hz)−1 s, which is adequately sampled by existing long-baseline interferometers. Hence, hydrodynamic turbulence imposes a fundamental noise floor on gravitational wave observations of neutron stars, although its polluting effect may be muted by partial decoherence in the hectohertz band, where current continuous-wave searches are concentrated, for the highest frequency (and hence most powerful) sources. This outcome is contingent on the exact shape of the turbulent power spectrum, which is modified by buoyancy and anisotropic global structures, such as stratified boundary layers, in a way that is understood incompletely even in laboratory situations.

Export citation and abstract BibTeX RIS

1. INTRODUCTION

Shortly after the discovery of radio pulsars, speculation arose that the superfluid interior of a differentially rotating neutron star is turbulent (Greenstein 1970). Since then, the theme has resurfaced intermittently during the quest to understand pulsar rotational irregularities, such as glitches and timing noise (Anderson et al. 1978; Tsakadze & Tsakadze 1980). Neutron star turbulence can be hydrodynamic, taking the form of a Kolmogorov-like cascade of macroscopic eddies at high Reynolds numbers (Peralta et al. 2005, 2006b, 2008; Melatos & Peralta 2007). Complicated vorticity patterns of this sort are observed in terrestrial experiments on spherical Couette flow, which undergo transitions to non-axisymmetric flow states at high Reynolds numbers, e.g., spiral, shear herringbone, or Taylor–Görtler vortices (Belyaev et al. 1979; Bühler 1990; Junk & Egbers 2000; Sha & Nakabayashi 2001; Nakabayashi et al. 2002a, 2002b; Nakabayashi & Tsuchida 2005; Peralta et al. 2006a, 2008), re-laminarization (Nakabayashi & Tsuchida 2005), or Stewartson layer disruption (Hollerbach 2003; Hollerbach et al. 2004, 2006; Wei & Hollerbach 2008). Neutron star turbulence can also be quantum mechanical, comprising a self-sustaining tangle of quantized microscopic vortices (Glaberson et al. 1974; Jou & Mongiovi 2004; Jou & Mongiovì 2006), excited by bulk two-stream instabilities (Andersson et al. 2007), interfacial two-stream instabilities (Blaauwgeers et al. 2002; Mastrano & Melatos 2005), or meridional circulation (Peralta et al. 2005, 2006b; Melatos & Peralta 2007). In general, macroscopic and microscopic superfluid turbulence appear to trigger each other; it is an unsolved, chicken-or-egg question as to which comes first (Barenghi et al. 2001; Tsubota 2009).

Turbulence powered by differential rotation is axisymmetric when averaged over long times but non-axisymmetric instantaneously. Turbulent flows therefore emit stochastic gravitational waves. In an incompressible fluid, the waves arise mainly from current quadrupole (and higher multipole) source terms. In a compressible fluid, the mass multipoles also contribute; indeed, they can dominate, e.g., during post-glitch Ekman pumping in a neutron star (van Eysden & Melatos 2008). Recently, Peralta et al. (2006a) pointed out that there exists a fundamental theoretical uncertainty regarding the shape and strength of the gravitational wave signal emitted by hydrodynamic turbulence. The mechanical stress–energy in Kolmogorov-like turbulence is contained mostly in large eddies near the stirring scale. Naively, therefore, one might expect the gravitational wave signal to look like a "dirty sinusoid," which reflects circulation on the largest scales and decoheres in approximately one rotation period. However, the instantaneous wave strain is proportional to the second time derivative of the stress–energy tensor, and this quantity is greatest for small eddies near the dissipation scale, which turn over most quickly. If the latter effect dominates, one might expect the signal to resemble white noise. Of course, large eddies match better to low-order multipoles than small eddies, and low-order multipoles typically dominate the gravitational wave strain far from the source I. M. Wasserman (2009, private communication). A careful calculation is therefore required to select between the various possibilities and reliably estimate the detectability of the signal.

In this paper, we undertake such a calculation by combining the formalism of Kosowsky et al. (2002) and Gogoberidze et al. (2007), developed to calculate the gravitational radiation from a turbulent, first-order phase transition in the early universe, together with the formalism of I. M. Wasserman (2009, private communication), developed to calculate the gravitational radiation from non-axisymmetric vorticity fluctuations in neutron stars, e.g., due to clusters of quantized superfluid vortices. In Section 2, we analyze global hydrodynamic simulations of incompressible, shear-driven neutron star turbulence to extract the vorticity correlation function, which feeds into the statistics of stress–energy fluctuations in the source. We then calculate analytically the current multipole moments, root mean square (rms) wave strain, and decoherence time of the resulting, stochastic gravitational wave signal in Section 3. The results are applied in Section 4 to estimate the detectability of the signal, e.g., with long-baseline interferometers like the Laser Interferometer Gravitational-Wave Observatory (LIGO), and its polluting effect on continuous-wave searches currently under consideration. Hydrodynamic turbulence imposes a fundamental, quantifiable noise floor on gravitational wave observations of neutron stars. Astrophysical implications, including the rate of gravitational wave braking in different types of neutron stars, are briefly canvassed.

2. TURBULENT VORTICITY CORRELATIONS

Let ω(x, t) be a turbulent vorticity field which fluctuates stochastically with position x and retarded time t in the source. The gravitational wave strain generated by the (l, m)th current multipole is proportional to the lth time derivative of ω(x, t) integrated over the source volume, as discussed by Thorne (1980) and in Section 3.2. If the turbulence is stationary, the wave strain averages to zero over times longer than the maximum eddy turnover time. However, the rms wave strain is not zero. It is proportional to the autocorrelation function 〈ωi(x, t)ω*j(x', t')〉 integrated over the source volume. Here and elsewhere, angular brackets denote the usual ensemble average. Expressing ω(x, t) in terms of its spatial Fourier transform, we can write

Equation (1)

where v(x, t) is the turbulent velocity field satisfying ω(x, t) = curl v(x, t). For stationary turbulence, 〈ωi(x, t)ω*j(x', t')〉 depends on t and t' only through the combination τ = t' − t.

Global, three-dimensional, numerical simulations of neutron star turbulence driven by differential rotation indicate that the turbulence is approximately isotropic and stationary once the Reynolds number Re exceeds ∼104 (Peralta et al. 2005, 2006b; Melatos & Peralta 2007; Peralta et al. 2008). Figure 1 presents meridional streamlines of the viscous (left panel) and inviscid (right panel) components of a two-component, incompressible,4 Hall–Vinen–Bekarevich–Khalatnikov (HVBK) superfluid (Hills & Roberts 1977) in a differentially rotating shell, with dimensionless thickness δ = 0.3, Rossby number Ro = 0.1, and Reynolds number Re = 3 × 104, showing a snapshot of the flow at time t = 4.8 Ro−1Ω−1, where Ω denotes the angular velocity of the stellar surface. The simulation parameters are defined precisely in Peralta et al. (2008), where the numerical algorithm (pseudospectral collocation) is also laid out in detail.5 Although the Reynolds number in Figure 1 is ∼107 times less than in a realistic neutron star and only just above the threshold for turbulence, it is already possible to see that departures from isotropy are limited to the largest scales, i.e., ∼R*δ, where R* is the stellar radius. This is true even when the shear is stronger; Rossby numbers as high as 0.3 have been investigated. Additional pictorial examples appear in Peralta et al. (2008). Isotropy is expected to increase with Re, as in other turbulent systems, but simulations with Re ⩾ 105, which would test this claim, are not feasible computationally at present.

Figure 1.

Figure 1. Hydrodynamic turbulence in a two-component, incompressible, HVBK superfluid in a differentially rotating shell, with dimensionless thickness δ = 0.3, Rossby number Ro = ΔΩ/Ω = 0.1, and Reynolds number Re = 3 × 104. Parameters are defined in Peralta et al. (2008). The figure shows a snapshot of the meridional streamlines of the two components at time t = 48 Ω−1, taken from the numerical simulations in Peralta et al. (2008). The rotation axis points vertically. Left panel: viscous component. Right panel: inviscid component.

Standard image High-resolution image

The turbulence in Figure 1 is stationary to a good approximation. The streamline pattern reorganizes stochastically on the timescale Ω−1, and the velocity components at a fixed point alternate in sign, in such a way that the vorticity averages to the rigid body value $2{ \bOmega }$ over the long term. This behavior is summarized in Figure 2. The left panel shows the meridional streamlines of the viscous HVBK component at four instants in time, each separated by 2Ω−1. The eddies in the flow change noticeably in shape, size, and position throughout the simulation (and in the inviscid component, which is not plotted). In the right panel, we graph all three vector components of the vorticity at a mid-latitude point versus time, after subtracting the rigid body term $2{\bOmega }$. The mean of each component is plotted as a horizontal line for comparison. After initial transients die away, i.e., for t ≳ 20Ω−1, the turbulent vorticity fluctuates stochastically and without bias about its mean value and is independent of t. It has zero mean after adjusting for the residual differential vorticity $2(\Delta \Omega){\bf \hat{z}}$.

Figure 2.

Figure 2. Stationarity of the turbulence in Figure 1, for the same simulation parameters. Left panel: snapshots of the meridional streamlines of the viscous HVBK component at (a) t = 50 Ω−1, (b) t = 52 Ω−1, (c) t = 54 Ω−1, and (d) t = 56 Ω−1. The streamlines reorganize stochastically on the timescale Ω−1. Right panel: vorticity at a fixed point (r = 0.85R*, θ = π/4, ϕ = 0) as a function of time (150 ⩽ Ωt ⩽ 250). The spherical polar vorticity components are plotted after subtracting the rigid body rotation, viz., ωr − 2Ωr (solid curve), ωθ − 2Ωθ (dashed curve), and ωϕ (dotted curve). The horizontal lines represent the mean value of each component over the plotted range. The vector sum of the means equals $2(\Delta \Omega)\,{\bf \hat{z}}$, the residual differential vorticity.

Standard image High-resolution image

Isotropic, stationary turbulence has a velocity correlation function $\langle v_m({\bf k}) v_q^\ast ({\bf k}^{\prime }) \rangle = V (2\pi)^3 {\hat{P}}_{mq}(k) P(k) \delta ({\bf k} - {\bf k}^{\prime })$, where V is the total volume of the system (and drops out at the end of the calculation of any physical observable), ${\hat{P}}_{mq}(k)=\delta _{mq}-{\hat{k}}_m{\hat{k}}_q$ is a projection operator perpendicular to the wave vector k, and P(k) is the power spectrum of the turbulence, usually a power law. However, for gravitational wave problems, where we are interested in temporal fluctuations of the mean-square wave strain (and hence mean-square multipole moments), we are obliged to work with unequal time (tt') correlators of the kind appearing in Equation (1). As a working hypothesis, in this paper, we assume the standard Kraichnan form for high-Re turbulence in three dimensions (Kraichnan 1959; Kosowsky et al. 2002), viz.,

Equation (2)

with

Equation (3)

and

Equation (4)

In Equations (3) and (4), ε is the energy dissipation rate per unit enthalpy (units: m2 s-3), and η(k)−1 is the eddy turnover time at wavenumber k; that is, turbulent motions with wavelength 2π/k in any given local region decohere over time intervals longer than η(k)−1.

What is P(k) for neutron star turbulence? Alas, there is little one can say with confidence about this question, given the impossibility of direct measurements and the worrying experience with terrestrial systems, where idiosyncratic features are often imprinted on P(k) by boundary layers and other unavoidable global structures even in simple systems (see next paragraph). Nevertheless, in order to make progress, we assume that P(k) is a power law, P(k) ∝ kα, and that the power-law exponent is close to the Kolmogorov value for isotropic, high-Re turbulence, α = −11/3 (Kraichnan 1959; Gogoberidze et al. 2007). Following the standard normalization recipe, we can then write

Equation (5)

The power law stretches across an inertial range extending from the wavenumber corresponding to the stirring scale,

Equation (6)

up to the wavenumber corresponding to the viscous dissipation scale,

Equation (7)

where ν denotes the kinematic viscosity.

Direct numerical simulations provide reasonable support for the Kolmogorov scaling (Peralta et al. 2008). In Figure 3, we construct P(k) for the velocity field in Figure 1 from the simulation data by summing the squared pseudospectral coefficients |Cnlm|2 corresponding to each value of k for the top 103 modes (Peralta et al. 2005, 2008). That is, we plot P(k) as a function of kR*/2π = (n2 + l2 + m2)1/2, where n, l, and m denote radial, latitudinal, and toroidal mode indices, respectively. The viscous and inviscid HVBK components are analyzed in the left and right panels, respectively. We first note that the toroidal contribution |vϕ(k, t)|2 (open circles) dominates P(k), especially at small k, and adheres closely to the Kolmogorov scaling (solid line). The poloidal contributions |vr(k, t)|2 (open squares) and |vθ(k, t)|2 (open triangles) deviate from the Kolmogorov scaling at small k because isotropy breaks down for the largest eddies, as noted before; but, in any case, P(k) is dominated by |vϕ(k, t)|2 at small k. A least-squares fit to P(k) over the inertial range in Figure 3 yields α = −3.52 ±  0.35 for the viscous HVBK component (8 ≲ kR*/2π ≲ 38) and α = −3.55 ± 0.25 for the inviscid HVBK component (16 ≲ kR*/2π ≲ 53). These results agree surprisingly well with the Kolmogorov value α = −11/3, even though the turbulence in Figure 1 is not fully developed, the inertial range stretches over less than one decade in the simulations, and there are departures from isotropy at small k. We have verified that spectral filtering, which is implemented in the numerical solver to enhance its stability (Peralta et al. 2005, 2008), does not warp P(k) significantly for Re = 3 × 104.

Figure 3.

Figure 3. Turbulent power spectrum P(k) (in arbitrary units) for the simulation parameters in Figure 1, viz., δ = 0.3, Ro = 0.1, and Re = 3 × 104. The power spectrum is proportional to ∑n,l,m|Cnlm|2, where (n, l, m) are mode indices in the pseudospectral expansion; for a given wavenumber bin [k, k + dk], the sum runs over all indices satisfying k ⩽ (n2 + l2 + m2)1/2k + dk (among the top 103 modes). A snapshot of the spectrum at t = 50 Ω−1 is plotted as a function of normalized wavenumber kR*/2π, with logarithmic binning. It contains three contributions, |vr(k,  t)|2 (open squares), |vθ(k,  t)|2 (open triangles), and |vϕ(k, t)|2 (open circles), whose sum is proportional to P(k). The Kolmogorov scaling (solid line), with arbitrary normalization, is overplotted for comparison. Slightly different bins are used for the three terms. The steady-state differential rotation, contained in C010, lies off the scale. Left panel: viscous HVBK component. Right panel: inviscid HVBK component.

Standard image High-resolution image

It is important to reiterate that the scaling in Equation (5) applies to isotropic turbulence in the bulk, e.g., grid turbulence far from any walls. It is known from many laboratory experiments, e.g., in wind tunnels, that P(k) is modified by the presence of anisotropic global structures like boundary layers,6 to the point where it may not even be a power law (Saddoughi & Veeravalli 1994). The modifications are not merely of academic interest; we show in Section 3 that the amplitude of the gravitational wave signal from turbulence is sensitive to the form of P(k), e.g., through α. Unfortunately, calculating P(k) accurately is a formidable undertaking even in terrestrial situations, where conditions can be controlled, let alone in a neutron star. A voluminous literature exists on turbulent cascades in shear and viscous boundary layers; see, for example, the review by Robinson (1991). Experiments that use stereoscopic particle image velocimetry to measure the instantaneous velocity field and hence P(k) find that the boundary layer is populated by coherent structures, like hairpin vortices (Ganapathisubramani et al. 2003) or wall-wake flows (Perry & Marusic 1995; Marusic 2001), and anomalous Reynolds stresses (Wietrzak & Lueptow 1994; Ganapathisubramani et al. 2003), which are inconsistent with the Kolmogorov model (Saddoughi & Veeravalli 1994). The role of stratification (Fernando 1991), important in a neutron star, and the interplay between shear and buoyant convection (Moeng & Sullivan 1994), are also under active investigation. Resolving these matters lies far outside the scope of this paper, but it is important to recognize them and work toward a better understanding over time, e.g., by improving upon the pioneering superfluid spherical Couette simulations of Peralta et al. (2008).

Buoyancy suppresses radial motion in a neutron star, arguably reducing the turbulence to two dimensions. Kraichnan (1967) postulated that two-dimensional turbulence develops two inertial ranges: a −5/3 cascade (α = −11/3), which conserves kinetic energy and runs to low k, from the stirring scale up to the system scale; and a −3 cascade (α = −5), which conserves mean-square vorticity and runs to high k, from the stirring scale down to the dissipation scale. In a neutron star, the stirring and system scales are approximately the same, so the −3 cascade notionally covers a wider k range than the −5/3 cascade. However, laboratory experiments indicate that the situation is more complicated. For example, Iida et al. (2009) found P(k) ∝ k−4 for horizontal k and P(k) ∝ k−5 for vertical k. Vertical motions are suppressed, but the vertical cascade still plays a key role in routing energy through the system, by mediating the formation of sheared stacks of pancake-like structures. Sommeria (1986) showed experimentally that, even when the largest scales are pancake-like, smaller scales remain isotropic and Kolmogorov-like, perhaps due to the action of internal gravity waves (which are weakly damped at the Prandtl numbers Pr ≫ 1 expected inside a neutron star). As a rule, stratified turbulence fossilizes into two dimensions when the activity parameter I = ε/(νN2) drops below ≈7 (Iida et al. 2009), where N ≈ 500 rad s-1 is the Brunt–Väisälä frequency in a typical neutron star in beta equilibrium (van Eysden & Melatos 2008). Some of the candidate astrophysical sources considered in Section 4 satisfy this inequality, while others do not. We confirm, with the help of order-of-magnitude estimates in Section 3.3, that the gravitational wave results in Section 3 and Section 4 are probably insensitive to the dimensionality of the turbulence, at least for the values of α that one might reasonably expect in a neutron star. Nonetheless, we emphasize that the question is far from settled and needs to be studied more carefully with more sophisticated, compressible numerical simulations.

3. GRAVITATIONAL WAVE SIGNAL

3.1. Current Multipole Moments

The gravitational wave strain measured by an observer at a distance d from a source can be written in the transverse traceless gauge as a linear combination of gravitoelectric ("mass") and gravitomagnetic ("current") multipoles (see Equation (4.3) of Thorne 1980). The latter components take the form

Equation (8)

In Equation (8), Slm(t) denotes the (l, m)th current multipole moment, written as a function of the retarded time t, and TB2,lmjk denotes the associated gravitomagnetic tensor spherical harmonic, which describes the angular dependence (or beam pattern) of the radiation field. For a Newtonian source (slow internal motions, weak internal gravity), like a differentially rotating neutron star, the current multipole moments assume the form

Equation (9)
Equation (10)

where Yl−1,lm is a vector spherical harmonic of pure orbital type, Ylm is a scalar spherical harmonic, and ρ(x) is the fluid mass density. Equation (9) corresponds exactly to Equation (5.27b) in Thorne (1980). Equation (10) is derived from Equation (9) by expressing the vector harmonic in terms of gradients of scalar harmonics, viz.,

Equation (11)

and then integrating by parts I. M. Wasserman (2009, private communication). Physically, therefore, the current multipole arises from the magnetic component of the velocity field (Equation (9)) or, equivalently, from the radial component of the vorticity field (Equation (10)). We only consider incompressible turbulence in this paper (see Section 2), for which the mass multipoles vanish.

For the sake of simplicity, we take ρ to be uniform, i.e., ρ = 3M*/(4πR3*), where M* is the total mass of fluid in the stellar interior (cf. I. M. Wasserman 2009, private communication). A more realistic assumption is that ρ is incompressible (subsonic flow) but stratified gravitationally. However, we avoid stratification in this first pass at the problem because we wish to exploit the scalings for isotropic Kolmogorov turbulence in Section 2, which assume uniform ρ. As noted in Section 2, the stratified problem is much harder.

3.2. Root Mean Square Wave Strain

The mean wave strain at the observer is zero for stationary, isotropic turbulence, as discussed in Section 2. Therefore, to assess detectability, we compute the autocorrelation function

Equation (12)

which is positive definite for t = t', reduces to the rms wave strain hrms = 〈|hTTjk(t)|21/2for t = t', and is a function of t and t' only through the combination τ = t' − t for stationary turbulence. For the sake of simplicity, we compute C(τ) for a single multipole (l, m). By doing so, we circumvent the following complication: an observer at a particular position, doing a realistic detection experiment, sees cross terms in C(τ) which mix multipoles together (e.g., S20S21*). The cross terms only vanish when averaged over all possible observer positions (by the orthonormality of TJS,lmjk; see Equation (2.36) of Thorne 1980). In this sense, our detectability estimates are conservative; the signal from one multipole is obviously a lower bound on the total signal.

The wave strain autocorrelation function is the ensemble average of a product of time derivatives of Slm. In the special case of stationary turbulence, where it is always possible to (notionally) fix t (t') at some instant during the average and exchange d/dt' (d/dt) with d/dτ (−d/dτ), it is possible to simplify the ensemble average using the following formula from turbulence theory (e.g., Equation (12) of Gogoberidze et al. 2007):

Equation (13)

Upon combining Equations (8), (10), (12), and (13), and noting that |TB2,lmjk|2 ⩽ 1, we obtain the maximum autocorrelation an optimally situated observer can detect:

Equation (14)

Upon substituting Equations (1)–(3) into Equation (14), to deal with the ensemble average, and performing the k' integral over the delta function, we arrive at the expression

Equation (15)

Equation (15) can be simplified analytically in a number of ways. Here, we elect to expand the two plane-wave factors in terms of scalar spherical harmonics, viz.,

Equation (16)

and exploit the orthogonality properties of the spherical harmonics to make progress. In Equation (16), jL denotes a spherical Bessel function of the first kind. The six angular integrals (over ${ \hat{\bf x}}$, ${ \hat{\bf x}}^{\prime }$, and ${ \hat{\bf k}}$) now factorize easily, and we arrive at the final expression for C(τ),

Equation (17)

with

Equation (18)
Equation (19)

One can evaluate Ki and Nij by writing ${\hat{x}}_i$ and ${\hat{k}}_i$ in terms of Y10 and Y1,±1and using Clebsch–Gordan coefficients to evaluate the triple products. For example, Ki is nonzero only if L = l ± 1 and |M + m| ⩽ 1. In what follows, however, we do the integrals directly with the help of a symbolic algebra package.

3.3. Quadrupole

Let us begin by specializing to the case l = m = 2, where the sums in Equation (17) are nonzero for L = 1, 3 and L' = 1, 3 only. (The same is true for (l, m) = (2, 1).) Doing the r and r' integrals and τ derivatives, we find

Equation (20)

The integral defined by

Equation (21)

when combined with P(k) through the k integral in Equation (20), expresses mathematically the fact that the wave strain is an incoherent sum of randomly phased eddy motions, a subset of which are wavenumber matched to the multipole moment under consideration. It can be shown that |C(τ)|1/2 is roughly proportional to the square root of the mean number of eddies at the relevant scale I. M. Wasserman (2009, private communication).

The maximum rms wave strain is obtained for τ = 0, when there is no turbulent dephasing. The k integral in Equation (20) runs from the stirring to the dissipation scale, i.e., from ks in Equation (6) to kd in Equation (7). The ψ function is oscillatory, but its envelope grows ∝k3 for a = 5. Hence, if η(k) and P(k) are given by Equations (4) and (5) respectively, the integrand scales ∝k−3 overall in the limit τ → 0. The rms wave strain is therefore dominated by motions at the stirring scale. Under these circumstances, with kskd, the result is

Equation (22)

If the turbulence is powered by differential rotation, the specific (per unit mass) energy input per unit time is ε = R2*(ΔΩ)3 (Landau & Lifshitz 1959),7 from which we obtain

Equation (23)

The Rossby number is defined as Ro = ΔΩ/Ω, as in Section 2.

The above results confirm that we can reliably use the quadrupole moment to estimate detectability. As hrms is dominated by ks, and the stirring scale is well matched to l = 2, we are not missing a lot of power emitted by eddies at large k and showing up at proportionally higher l. On the debit side of the ledger, the assumption of isotropy is least valid at the stirring scale, as discussed in Section 2. Equation (20) makes it clear why this may ultimately turn out to be a serious flaw: hrms depends quite sensitively on the shape of P(k). For example, if the exponent of the power spectrum satisfies α> − 5/3, the dissipation scale governs hrms, not the stirring scale, and Equation (23) would exhibit different scalings with Ω and ΔΩ, plus an explicit dependence on the kinematic viscosity of the stellar interior.

Buoyancy arguably reduces the turbulence to two dimensions by suppressing radial motion, as discussed in Section 2. However, order-of-magnitude estimates suggest that the gravitational wave signal is insensitive to the dimensionality of the turbulence. Evaluating Equations (17) and (20) for α = −5, appropriate for the vorticity-conserving −3 cascade postulated by Kraichnan (1967), we find that hrms changes by less than 20% for kskd. In reality, because hrms is dominated by motions near ks, where the forward and reverse cascades cross, α effectively lies closer to −5/3 than to −3, even in the two-dimensional model of Kraichnan (1967), and the change is even smaller. In all cases, the power-law dependences on ε and the other quantities in Equation (22) are universal, as long as we have α < −5/3; the shape of the turbulent spectrum affects just the numerical prefactor in Equation (22), and even then only weakly, unless we have α> − 5/3, whereupon the dissipation scale dominates the signal. Theory and experiment are united in deeming spectra shallower than P(k) ∝ k−11/3 to be extremely rare in nature; certainly, buoyancy acts in the contrary direction to steepen P(k). As noted in Section 2, although neutron stars are strongly stratified, the activity parameter I = ε/(νN2) ≈ 40(ΔΩ/1 rad s-1)3(ν/10 m2 s-1)−1(N/500 rad s-1)−2 falls below the fossilization threshold I ≈ 7 (Iida et al. 2009) for ΔΩ < 0.6(ν/10 m2 s-1)1/3(N/500 rad s-1)2/3   rad s-1. Some of the candidate sources discussed below in Section 4 do not satisfy this inequality (see Table 1); that is, hydrodynamic turbulence in these sources may be effectively three dimensional despite strong stratification, lending extra support to the results in this section.

Table 1. Candidate Source Categories

Type d (kpc) Ω (rad s-1) Ro hrms τc (s) $|\dot{\Omega }_{\rm GW}|$ (rad s-2)
Protoneutron Star
(a) Core (30 km) 10 2 × 103 1 9 × 10−21 1 × 10−4 3 × 100
(b) Envelope (300 km) 10 1 × 103 0.1 5 × 10−22 3 × 10−3 1 × 10−6
Glitching Pulsar
(a) Vela-like 1 1 × 102 10−2 5 × 10−31 3 × 10−1 3 × 10−27
(b) Crab-like 1 2 × 102 10−4 4 × 10−36 1 × 101 4 × 10−41
Accretor
(a) Millisecond pulsar 1 3 × 103 10−2 1 × 10−26 9 × 10−3 7 × 10−17
(b) White dwarf 0.1 0.2 0.1 2 × 10−28 1 × 101 3 × 10−29
Nearby Pulsar
(a) Fast rotator 0.01 3 × 103 0.1 1 × 10−21 9 × 10−4 7 × 10−9
(b) Slow rotator 0.01 2 × 102 0.1 4 × 10−25 1 × 10−2 4 × 10−17

Download table as:  ASCIITypeset image

3.4. Higher Multipoles

As a general rule, the current multipole moment depends increasingly strongly on kd as l increases, because higher multipoles match better to small-scale eddies. This effect is communicated through the ψ function combined with P(k), e.g., as in Equation (20). Whenever l increases by 1, the τ derivatives in Equation (17) bring down two extra powers of η(k), the r and r' integrals contribute an extra factor of k2, and the ψ functions contribute new k factors too. For l = 3, we find

Equation (24)
Equation (25)

with kskd. The integrand in Equation (24) scales ∝k−5/3, so hrms for l = 3 is dominated by motions at the stirring scale, just like for l = 2. The ratio of the l = 3 and l = 2 wave strains is 0.83(R*Ω/c)Ro, implying that the l = 3 radiation is significantly weaker for any realistic neutron star rotating slower than centrifugal breakup.

For l ⩾ 4, one finds that hrms is dominated by motions at the dissipation scale, not the stirring scale. Consequently, the dependence of hrms on Ω and ΔΩ also changes, and a new, explicit dependence on the kinematic viscosity appears. For example, for l = 4, one obtains hrms ∝ ρR5*ε5/3k1/3dd−1 ∝ ν−1/4, and the ratio of the l = 4 and l = 2 wave strains is ≈(R*Ω/c)2Ro2(kdR*)1/3. In view of the Kolmogorov scaling kd/ks = Re3/4 (see Equation (17) of Kosowsky et al. 2002), we conclude that the l = 4 radiation is weaker than the l = 2 radiation for Re ≲ (R*Ω/c)−8Ro−8, which is always satisfied except possibly in a strongly sheared millisecond pulsar, e.g., one born recently in a supernova.

3.5. Decoherence Time

The decoherence time τc is the time that must elapse before the instantaneous wave strains hTTjk(t) and hTTjk(t') become statistically uncorrelated. We define it to be the value of τ = t' − t at which C(τ), defined as the ensemble-averaged correlator in Equation (12), decreases to some fixed fraction (say, one quarter) of its maximum (at τ = 0).

The results in Section 3.3 and Section 3.4 demonstrate that hrms is generated predominantly by the stress–energy in motions at the stirring scale. It is tempting, therefore, to predict that the signal decoheres on the eddy turnover time at the stirring scale, i.e., the rotation period of the star multiplied by Ro (Kosowsky et al. 2002). In fact, the situation is potentially more complicated. Looking at Equation (20), which describes how C(τ) decreases with τ, we see that the integrand contains three terms, which scale ∝k−3, k−5/3, and k−1/3. The first, positive term is the only one which contributes to the maximum of hrms, i.e., Equation (22); its integral is dominated by ks, as discussed above. The second, negative term is the only one which can reduce C(τ) and cause decoherence. Its integral is also dominated by ks. The integral of the third, positive term involves both ks and kd but is dominated by the former (see next paragraph). Hence the final result for τc contains information about both ks and kd in general. But, for the l = 2 and l = 3 signals, it turns out that the dependence on kd is very weak, unlike for higher l. This result is crucial for the question of detectability, as explained in Section 4.

By graphing C(τ) numerically in Figure 4, we find that it falls to one quarter of its maximum over a time comparable to the eddy turnover time at the stirring scale, η(ks)−1. This graphical task is tricky, because the square of the ψab functions oscillates rapidly for kdks. However, we can obtain an excellent analytic approximation by averaging over many cycles of the fast oscillation. Upon writing [ψ51(x) +  ψ53(x)]2 ≈ 25x6/2 plus fast oscillations ∝cos(2x) in the regime x ≫ 1, we arrive at the formula

Equation (26)

where Erf(x) symbolizes the error function. Equation (26) is accurate to ≈12% across the full range of τfor the parameter range under consideration. The dependence on kd/ks ≫ 1 in the third term is weak, although this stops being true for the rare case of non-Kolmogorov spectra with α> − 5/3, where kd dominates (see Section 3.3). From Equations (4), (6), (7), and (26), with ε = R2*(ΔΩ)3, we arrive at the following expression for the decoherence time corresponding to the half-strain point:8

Equation (27)
Equation (28)

The shear viscosity in a neutron star is given by the neutron–neutron scattering formula derived by Cutler & Lindblom (1987), viz., ν = 10(ρ/6 × 1017 kg m-3)5/4(T/108 K)−2 m2 s-1, where T is the temperature of the stellar interior. From Equation (7), we obtain kdks for fiducial values of ρ and T. In a protoneutron star, discussed in Section 4, the bulk viscosity (∝T8) exceeds the shear viscosity (Sawyer 1989). Under these conditions, small compressions enhance the turbulent dissipation rate, effectively reducing kd. Numerical simulations of a compressible HVBK superfluid (outside the scope of this paper) are needed to quantify this effect properly. However, the foregoing gravitational wave calculations continue to hold to a good approximation, provided that the stirring and dissipation scales remain moderately well separated (e.g., kd/ks ≳ 5 in Equation (26)).

Figure 4.

Figure 4. Wave strain autocorrelation function C(τ), defined by Equation (26), normalized by the mean-square wave strain h2rms, as a function of the time lag τ, normalized by the maximum eddy turnover time η(ks)−1 at the stirring scale, for kd/ks = 103.

Standard image High-resolution image

4. DETECTABILITY AND ASTROPHYSICAL IMPLICATIONS

In this paper, we calculate analytically the current multipole moments generated by vorticity fluctuations in high-Re turbulence in a differentially rotating neutron star. We derive an analytic expression (Equation (20)) for the wave strain autocorrelation function C(τ) of the resulting, stochastic gravitational wave signal to leading (quadrupole) order, in terms of the turbulent power spectrum and eddy turnover time spectrum, and show that C(τ) is governed by motions at the stirring scale. From the same equation, we also compute the root mean square wave strain hrms and the decoherence time τc of the signal and show that τc approximately equals the eddy turnover time at the stirring scale (usually much longer than the rotation period). Convenient formulae for hrms and τc, scaled in terms of astrophysical parameters, are presented in Equations (23) and (28), respectively.

We can use these formulae to assess qualitatively whether the stochastic signal can be detected by existing and planned long-baseline gravitational wave interferometers. Clearly, neutron star turbulence imposes a fundamental, unavoidable, astrophysical noise floor on continuous-wave searches for neutron stars, e.g., Abbott et al. (2009). But how seriously does it pollute such searches? To answer this question, we note two things. First, most neutron stars are either rotating too slowly or with too little shear to cause trouble, according to Equation (23), at least at the sensitivities anticipated for Enhanced and Advanced LIGO. Second, for the small subset of neutron stars that are potentially powerful stochastic emitters, the signal decoherence time is relatively short; indeed, for ΔΩ ∼ Ω and the fastest rotators, τc approaches (without ever quite reaching) the sampling time of LIGO-like interferometers. Equations (23) and (28) imply that τc decreases as hrms increases.

Table 1 illustrates these conclusions by listing hrms and τc for several realistic categories of neutron star sources.

  • 1.  
    Protoneutron stars. The stellar angular velocity profile is a key output of radiation (magneto)hydrodynamics simulations of core-collapse supernovae (Ott et al. 2006; Burrows et al. 2007). All progenitor models tested so far lead to strong differential rotation. Looking, for example, at the fastest rotators in Figures 8, 9, and 15 of Ott et al. (2006) or Figure 14 of Burrows et al. (2007), captured 0.2 s after core bounce in a two-dimensional, unmagnetized explosion, we see that Ω decreases gradually from ≈3 × 103 rad s-1 at r = 3 km to ≈1 × 103 rad s-1 at r = 30 km (enclosed mass ≈1.2 M), before dropping steeply to ≈20 rad s-1 at r = 300 km (enclosed mass ≈1.7M).9 Conservatively, for a hypothetical supernova in the Milky Way, the above figures imply hrms = 9 × 10−21and τc = 0.1 ms  (the first line of Table 1). LIGO is capable of detecting such stochastic emission, because τc is greater than the sampling time of the interferometer (≈60 μs), even though the signal partially decoheres. The prospects are brighter if the protoneutron star rotates slower and with less shear, and the gravitational radiation emanates mainly from the extended envelope (second line of Table 1), so that lower hrms (5 × 10−22) is traded for higher τc (3 ms). Detection by Advanced LIGO is then possible in principle, e.g., via a cross-correlation search (Dhurandhar et al. 2008). The signal persists for ≳102 s ≫ τc before the differential rotation dissipates and/or the protoneutron star spins down, e.g., due to r-modes, fallback, or a magnetized wind (Lai et al. 2001; Ott et al. 2006).
  • 2.  
    Glitchers. Discontinuous spin-up events ("glitches") observed in some rotation-powered pulsars are adduced as evidence that neutron stars rotate differentially; the nuclear lattice spins down electromagnetically, lagging the superfluid due to vortex pinning (Anderson & Itoh 1975). The fractional jump in angular velocity ranges from 10−11 to 10−4 across the pulsar population. The surprising absence of a "reservoir effect" (i.e., glitch size ∝ time elapsed since the preceding glitch) in most objects implies that the observed spin-ups are a small percentage of the underlying shear (Melatos et al. 2008; Warszawski & Melatos 2008). It is therefore safe to expect Ro ⩾ 10−4in some pulsars, although the proportion is hard to quantify. Two examples are given in Table 1: an adolescent, Vela-like pulsar (age ∼104 yr), which spins relatively slowly (∼10 Hz) but undergoes relatively large glitches (fractional jump ∼10−6); and a young, Crab-like pulsar (age ∼103 yr), which spins faster (∼30 Hz) but undergoes smaller gliches (fractional jump ∼10−8). Assuming that the typical glitch resets ∼0.01% of the underlying shear, we find that the decoherence time (0.3 s to 10 s) matches well to the LIGO pass band, but the signal is too weak (hrms ≲ 10−30) to be detected by interferometers under development, or to pollute continuous-wave searches targeted at glitching pulsars (van Eysden & Melatos 2008).10
  • 3.  
    Accretors. The hydromagnetic accretion torque acting on a compact star in a mass-transfer binary is often comparable to, or greater than, the electromagnetic spin-down torque acting on an isolated pulsar (Hartman et al. 2008) and fluctuates by several percent daily, causing X-ray variability. Accreting objects are therefore likely to rotate differentially, with Rossby numbers comparable to, or greater than, those of isolated glitching pulsars. Two examples are given in Table 1: a standard, accreting millisecond pulsar (fifth line), and an accreting white dwarf with an ONeMg core (sixth line) in a system on the verge of accretion-induced collapse (Blackman & Yi 1998; Dessart et al. 2007; Metzger et al. 2008). The angular velocity profile of the white dwarf before (and after) collapse is plotted in Figure 12 of Dessart et al. (2006). In both cases, hrms ≲ 10−26 falls below the threshold for detection by Advanced LIGO but may approach the sensitivity of the next generation of interferometers.
  • 4.  
    Nearby neutron stars. The nearest radio pulsar discovered to date is PSR J0108 − 1431, with d = 85 pc (Tauris et al. 1994). The nearest millisecond pulsar is PSR J0437 − 4715, with d = 157 pc (Verbiest et al. 2008). However, many radio-quiet neutron stars with d < 85 pc should reside in the solar neighborhood; the evidence is both observational (e.g., radio-quiet X-ray point sources like the "Magnificent Seven," some with parallaxes; see Popov et al. 2003 and references therein) and theoretical (e.g., population synthesis models predict ∼109compact objects in the Milky Way; see Kiel et al. 2008 and references therein). If the nearest objects have d ∼ 10 pc and rotate reasonably fast, they represent promising LIGO candidates.11 For example, Table 1 quotes hrms for two nearby isolated pulsars with Ro = 10−1. The millisecond pulsar in particular is a bright source, although it suffers from the usual drawback of fast rotators: τc is short. Although it is sometimes assumed that millisecond pulsars do not rotate differentially, because they are not seen to glitch, it is equally possible that they glitch strongly but infrequently, because the electromagnetic spin-down torque is weak (Melatos et al. 2008).

Even before detecting gravitational waves, we can use existing radio timing observations of neutron stars to test and constrain the model in this paper. The radiation emitted by hydrodynamic turbulence exerts a reaction torque, which fluctuates in sign instantaneously but drains rotational kinetic energy from the fluid over time. The angular velocity of the star decreases at the average rate (Thorne 1980; I. M. Wasserman 2009, private communication)

Equation (29)
Equation (30)
Equation (31)

Equation (31) follows from Equation (29) via the sequence of steps in Sections 3.13.3. In the last column of Table 1, we compute $|\dot{\Omega }_{\rm GW}|$ for the source categories discussed above. In several instances, the fiducial value of $|\dot{\Omega }_{\rm GW}|$ is close to, or even greater than, the value of $|\dot{\Omega }|$ measured in radio timing experiments. Already, this places constraints on the source parameters, as the inequality $|\dot{\Omega }_{\rm GW}| \le |\dot{\Omega }|$ must always hold for any specific object. For example, the oldest radio millisecond pulsars are measured to have $|\dot{\Omega }| \sim 10^{-15}\,{\rm rad\,s^{-2}}$ (see, for example, Table 10.1 in Lyne & Graham-Smith 1998). Accreting millisecond pulsars in low-mass X-ray binaries have $|\dot{\Omega }| \lesssim 10^{-13}\,{\rm rad\,s^{-2}}$during X-ray outbursts and $|\dot{\Omega }| \sim 10^{-16}\,{\rm rad\,s^{-2}}$ between X-ray outbursts (Hartman et al. 2008). If these data are representative of the entire millisecond pulsar population, then we can start to rule out the existence of large shears such as those quoted in the fifth and seventh entries of Table 1. Indeed, Equation (31) and the data combine to yield a direct, observational upper limit on the rotational shear in any neutron star, viz.

Equation (32)

Equation (32) is an important result. It allows the theory in this paper to be falsified, if a glitching pulsar is discovered whose fractional angular velocity jumps exceed the right-hand side of Equation (32). We will investigate this application more thoroughly in a forthcoming article. Figure 5 gives a taste of what is possible. In the left panel of Figure 5, we plot Romax (i.e., the right-hand side of Equation (32)) versus spin-down age $\Omega /(2|\dot{\Omega }|)$ for all objects with Romax ⩽ 1 in the Australia Telescope National Facility (ATNF) Pulsar Catalog (Manchester et al. 2005).12 Clearly, millisecond pulsars with ages ≳108 yr place strict limits on Ro, with Romax ∼ 10−2 in some cases. Furthermore, the pulsars marked with diamonds have been observed to glitch at least once. If the theory in this paper is correct, then the largest angular velocity jump observed, (ΔΩ/Ω)g,max, cannot exceed the underlying shear, ΔΩ/Ω, in each of the glitching pulsars; in fact, it is expected to be much smaller, given the absence of a reservoir effect. As a test, we plot (ΔΩ/Ω)g,max/Romax versus spin-down age in the right panel of Figure 5. In all cases, the plotted ratio is much smaller than unity, consistent with theory.13 Future observational campaigns to simultaneously monitor large groups of pulsars, e.g., with multibeam radio telescopes such as the Square Kilometer Array, will challenge the theory more keenly.

Figure 5.

Figure 5. Upper limit on the angular shear in rotation-powered pulsars, assuming gravitational wave spin-down due to hydrodynamic turbulence. Left panel: maximum Rossby number, Romax, given by the right-hand side of Equation (32), vs. characteristic age, $\Omega /(2| \dot{\Omega } |)$ (in yr), for objects with Romax ⩽ 1 in the ATNF Pulsar Catalog. The points marked with diamonds indicate pulsars with a history of glitch activity. Right panel: maximum observed glitch size, (ΔΩ/Ω)g,max, expressed as a fraction of Romax, plotted vs. characteristic age, for the points marked with diamonds in the left panel. The ratio should not exceed unity for any pulsar.

Standard image High-resolution image

In summary, hydrodynamic turbulence imposes a fundamental noise floor on gravitational wave observations of neutron stars. Its polluting effect is muted by partial decoherence in the hectohertz band, where current continuous-wave searches are concentrated, but only for the fastest rotators with the strongest shear. In addition, the mechanism sets a fundamental lower limit on the spin-down rate $|\dot{\Omega }|$ and hence an observational upper limit on the Rossby number Ro, when combined with pulsar timing data. These conclusions hold subject to one major caveat, which is discussed at length in Section 2, to wit, that hrms and τc depend somewhat on the exact shape of the turbulent power spectrum, P(k). It is known from laboratory experiments that P(k) is modified away from its isotropic, Kolmogorov form by the presence of anisotropic global structures like turbulent boundary layers, which are profoundly difficult to model in terrestrial contexts, let alone in a neutron star. Likewise, buoyancy modifies P(k) by creating unequal turbulent cascades parallel and perpendicular to the stratification direction, in a fashion that is still not understood completely even in controlled laboratory experiments. The order-of-magnitude estimates in Section 3.3 provide some comfort that the gravitational wave results are insensitive to the above considerations, but we emphasize that much work (e.g., compressible HVBK simulations) still needs to be done to clarify the issue.

More complete detectability estimates referring to specific search pipelines lie outside the scope of this paper. Likewise, we defer calculating the detailed frequency spectrum of the stochastic signal, a substantial task.

The authors thank Ira Wasserman for privileged early access to an illuminating preprint (I. M. Wasserman 2009, private communication), Sterl Phinney for directing us to the work of Kosowsky and collaborators, Scott Hughes for a valuable discussion on multipole expansions, and the anonymous referee for spotting an error in our original calculation of the decoherence time and for generally improving the manuscript. We also acknowledge the substantial computing resources and support provided by the Victorian Partnership for Advanced Computing and the Damiana Cluster at the Albert-Einstein-Institut, which made the simulations in this paper possible. A.M. thanks the Astronomy Department at Cornell University and the LIGO Data Analysis Group at the California Institute of Technology for their hospitality during the period when this work began, and the LIGO Visitor Program for financial support. C.P. acknowledges the support of the Max-Planck Society (Albert-Einstein Institut).

Footnotes

  • The incompressible approximation is acceptable when the turbulent motions are subsonic, as in a neutron star (Peralta et al. 2005).

  • Such simulations typically adopt no-slip and no-penetration boundary conditions for the viscous component, perfect slip or no slip for the inviscid component, a Stokes flow initially, and a mutual friction force of the Gorter-Mellink form appropriate for a quantized vortex tangle (Gorter & Mellink 1949; Glaberson et al. 1974; Peralta et al. 2005; Andersson et al. 2007).

  • Boundary layers in a neutron star are thin, e.g., Re−1/2 (surface Ekman layer), Re−1/3 (Stewartson layer tangent to the core), or Re−2/5 (equatorial Ekman layer) in units of R* (Peralta & Melatos 2009). Nevertheless, they influence a large volume of fluid by partitioning the flow globally into cells, thereby shaping P(k) at low k. Laboratory experiments also reveal transient, ribbon-like streamers, which resemble boundary layers, throughout the body of otherwise isotropic Kolmogorov turbulence. The streamers generate anomalous Reynolds stresses and significant instantaneous departures from isotopy at low k (Marusic 2001; Ganapathisubramani et al. 2003).

  • This choice of ε is conservative. Other possible choices, e.g., ε = R2*Ω2ΔΩ, imply a higher gravitational wave strain.

  • It should be noted that η(k) does not always follow Equation (6) near the stirring scale (Gogoberidze et al. 2007), as discussed in Section 2. We aim to refine Equation (27) in future work.

  • Strongly magnetized models rotate ∼10 times slower (Heger et al. 2005; Burrows et al. 2007); cf. millisecond protomagnetar engine for long-duration gamma-ray bursts (Bucciantini et al. 2007).

  • 10 

    The persistent, stochastic, gravitational wave emission from shear-driven hydrodynamic turbulence is unrelated to the putative burst emission from the glitches themselves.

  • 11 

    Our inability to measure an ephemeris from radio timing observations does not affect cross-correlation searches for the stochastic radiation from neutron star turbulence, although it is a major obstacle to coherent continuous-wave searches.

  • 12 

    Objects with Romax>1 do not place a useful limit on the shear.

  • 13 

    The limits derived from Figure 5 are conservative. In the standard vortex unpinning paradigm, the fractional angular velocity jump observed in a large glitch can be related to the internal shear according to ∼(Is/Ic)(Δr/R*)(ΔΩ/Ω), where Is/Ic ∼ 102 is the ratio of the superfluid and crustal moments of inertia, and Δr/R* ∼ 10−6 is the normalized radial distance moved by the unpinned vortices (Alpar et al. 1986). Arguably, therefore, the theory is challenged if (ΔΩ/Ω)g,max/Romax exceeds ∼10−4 rather than unity, a stronger test.

Please wait… references are loading.
10.1088/0004-637X/709/1/77