- Split View
-
Views
-
Cite
Cite
Andrew Melatos, Bennett Link, Pulsar timing noise from superfluid turbulence, Monthly Notices of the Royal Astronomical Society, Volume 437, Issue 1, 01 January 2014, Pages 21–31, https://doi.org/10.1093/mnras/stt1828
- Share Icon Share
Abstract
Shear-driven turbulence in the superfluid interior of a neutron star exerts a fluctuating torque on the rigid crust, causing the rotational phase to walk randomly. The phase fluctuation spectrum is calculated analytically for incompressible Kolmogorov turbulence and is found to be red; the half-power point is set by the observed spin-down rate, the crust-superfluid lag and the dynamical response time of the superfluid. Preliminary limits are placed on the latter quantities using selected time- and frequency-domain data. It is found that measurements of the normalization and slope of the power spectrum are reproduced for reasonable choices of the turbulence parameters. The results point preferentially to the neutron star interior containing a turbulent superfluid rather than a turbulent Navier–Stokes fluid. The implications for gravitational wave detection by pulsar timing arrays are discussed briefly.
INTRODUCTION
Timing noise is a type of rotational irregularity observed in all isolated radio pulsars, in which pulse arrival times wander stochastically about the fitted ephemeris (Boynton et al. 1972; Cordes & Helfand 1980; Cordes & Downs 1985; D'Alessandro et al. 1995; Baykal et al. 1999; Hobbs, Lyne & Kramer 2010; Shannon & Cordes 2010). It is characterized as a random walk in the rotational phase, angular velocity or torque; the Fourier spectrum is always red, implying a process autocorrelated on a time-scale of hours to years (Cordes & Helfand 1980). Timing noise has been attributed to various mechanisms: microjumps akin to small glitches (Cordes & Downs 1985; Janssen & Stappers 2006; Melatos, Peralta & Wyithe 2008), recovery from unseen glitches (Johnston & Galloway 1999), quasi-periodic changes in magnetospheric structure (Lyne et al. 2010), variable coupling between the crust and liquid interior (Alpar, Nandkumar & Pines 1986; Jones 1990), stochastic variations in the star's shape (Cordes 1993) and fluctuations in the spin-down torque (Cheng 1987a,b; Urama, Link & Weisberg 2006). The most recent comprehensive survey of timing irregularities in 366 pulsars over time-scales longer than a decade found that a low-frequency noise process cannot explain all the observations on its own (Hobbs et al. 2010), suggesting that more than one physical mechanism contributes. In some pulsars, the pulse times-of-arrival are correlated over weeks with the distinctive signature of a relaxation process, e.g. damping of internal differential rotation (Price et al. 2012).
Recent research is lending growing support to the hypothesis, propounded originally by Greenstein (1970), that the superfluid interior of a neutron star is turbulent. The turbulence takes one of two forms: macroscopic, Kolmogorov-like eddies driven by crust–core shear at high Reynolds numbers, possibly involving unstable structures like Stewartson layers (Peralta et al. 2005, 2006a,b; Melatos & Peralta 2007; Peralta & Melatos 2009; Melatos 2012); and microscopic tangles of self-regenerating, reconnecting, quantized vorticity, driven by Kelvin-wave instabilities like those seen in terrestrial superfluids (Peralta et al. 2006a; Andersson, Sidery & Comer 2007) or dissipative instabilities arising from perfect or imperfect pinning in the inner crust or outer core (Link 2012a,b). It has been suggested that these forms of turbulence develop in a sustained manner and contribute to the stochastic spin variations observed in radio pulsars (Greenstein 1970; Peralta 2007; Link 2012b).
In this paper, we calculate from first principles the torque statistics and phase wandering produced by superfluid turbulence. The effect sets a timing noise floor, on top of which other processes like magnetospheric state changes add their contributions. In Section 2, we calculate the autocorrelation function for angular momentum fluctuations in the context of an idealized neutron star model. In Section 3, we calculate the power spectral density of the phase residuals and present convenient analytic formulas for the roll-over-frequency and zero-frequency normalization. In Section 4, we compare the predicted spectrum with timing data from a few representative objects and show how to place limits on quantities of fundamental physical importance, like the dynamical response time of the neutron superfluid and the moment of inertia of the stellar crust, mindful that a low-frequency noise process cannot explain all the irregularities observed (Hobbs et al. 2010) and that a comprehensive comparison with more data must still be done. The calculations are closely related to techniques developed to compute the stochastic gravitational radiation emitted by a turbulent neutron star (Melatos & Peralta 2010) and by phase transitions in the early Universe (Kosowsky, Mack & Kahniashvili 2002; Gogoberidze, Kahniashvili & Kosowsky 2007).
TURBULENT TORQUE
Idealized neutron star model
We start by considering an idealized model of a neutron star as two coupled subsystems. The first subsystem is the rigid crust and charged electron–proton fluid, which we assume are locked together by magnetic stresses and hence corotate (Alpar, Langer & Sauls 1984; Ruderman, Zhu & Chen 1998).1 The second subsystem is the inviscid neutron condensate. Our interest is in the scenario where the condensate is turbulent, driven by one or more of the processes referenced in Section 1. For simplicity, we assume that the mass density ρ of the star is uniform.
The two subsystems couple through friction. The exact nature of the friction is unimportant for this paper, but we now describe some relevant processes to give physical context. In the outer core, the dominant contribution to friction appears to be the scattering of electrons off vortices in the neutron condensate, whose cores are magnetized by entrainment of the neutron and proton mass currents (Alpar et al. 1984). This interaction could be modified significantly by vortex clustering (Sedrakian & Sedrakian 1995; Sedrakian et al. 1995). The protons are expected to form a type II superconductor, with the magnetic flux confined in flux tubes, which are frozen to the highly conducting, charged plasma. The vortices of the neutron superfluid pin to the flux tubes, primarily through a magnetic interaction, with pinning energies as high as ∼102 MeV per vortex–flux-tube junction (Srinivasan et al. 1990; Jones 1991; Mendell 1991; Chau, Cheng & Ding 1992; Ruderman et al. 1998; Glampedakis, Andersson & Samuelsson 2011; Link 2012a). Pinning partly decouples the neutron and charged fluids, increasing the coupling time-scale (Link 2012c) and sustaining an angular velocity difference, as the charged component of the star is spun down by the external, electromagnetic torque. While pinning energies remain rather uncertain, the conclusion that vortices pin to outer-core flux tubes appears to be increasingly likely (provided the outer core is indeed a type II superconductor); this point is discussed further in Section 2.2.
For simplicity, we assume that angular momentum transport through the neutron condensate occurs instantaneously. This approximation is good, even if angular momentum is stored temporarily in a third subsystem, e.g. Kelvin waves propagating along superfluid vortices pinned to the inner crust or hydromagnetic-inertial and cyclotron-vortex waves propagating through the charged fluid (Easson 1979a; Mendell 1998; Melatos 2012); the associated wave-crossing time-scales are still fast (≲102 s) in a typical star.
We suggest as a useful mechanical analogy a boiling pot of water on a frictionless stove. As the water (the turbulent condensate) boils, the total angular momentum of the water fluctuates, and a stochastic torque is applied to the pot (the charged component). This paper is concerned with analysis of the stochastic torque in the context of a neutron star.
In what follows, we regard the turbulent condensate as driving the crust, not vice versa; that is, the angular random walk executed by the crust does not feed back to modify the turbulence, at least on observational time-scales of decades. We justify this approximation quantitatively a posteriori in Section 4.1.
Turbulent condensate
The physics of the turbulent condensate has been examined previously by many authors (Peralta et al. 2005, 2006b; Andersson et al. 2007; Melatos & Peralta 2010; Link 2012a,b). Uncertainties remain. Large-scale simulations in the non-linear regime have been performed for some driving mechanisms, but even so the limited dynamic range means that important physics is not always captured. Here, we review the main possibilities briefly, emphasizing those aspects that motivate the idealized model developed in Sections 2 and 3. In essence, the model postulates the existence of turbulence with Kolmogorov-like statistics in some fluid component that couples to the crust. There are many ways to realize this scenario, and we now describe some of them.
Turbulence driven by vortex instabilities can arise in the core and/or the inner crust. In the core, where the vortices may be pinned to flux tubes, imperfect pinning destabilizes the vortex lattice; the source of free energy is the differential rotation between core neutrons and the proton–electron fluid, which locks magnetically to the crust (Link 2012a). In the absence of turbulence, the angular velocity lag is ∼0.1 rad s−1; in its presence, the lag and hence the steady-state injected power per unit enthalpy (labelled ε in Section 2.3) remain unknown, because the non-linear saturation time-scale has not yet been calculated. Pinning-driven vortex instabilities can also occur in the inner crust (Link 2012b), with pinning by flux tubes replaced by pinning at nuclear lattice sites. The two scenarios are essentially identical with regard to the calculations in this paper, the main difference being the inertia carried by the turbulent condensate (e.g. core neutrons versus inner crust neutrons) and associated entrained components.
Turbulence driven by meridional circulation (e.g. unstable Ekman pumping at high Reynolds number) can also occur in the core and/or inner crust. Again, the angular velocity difference and steady-state power have not been calculated self-consistently in the literature; the angular velocity of the outer crust is specified by fiat in simulations, without adjusting for the back-reaction torque from the viscous proton–electron component (Peralta et al. 2005, 2006b). Nevertheless, for realistic neutron star parameters, the differential velocity projected along the rotation axis arising from spin-down-powered Ekman circulation greatly exceeds the Donnelly–Glaberson instability threshold (Peralta et al. 2005; Andersson et al. 2007) in large parts of the core and inner crust, generating islands of tangled vorticity and patchy mutual friction (Peralta et al. 2006b).2 Liquid helium experiments show that turbulent velocity spectra in a superfluid are Kolmogorov-like (over two decades in wavenumber) in various grid, wake and ‘chunk’ flows (Salort et al. 2010). The chief theoretical input, ε, is left as a free parameter, to be constrained by pulsar timing data, although it can also be estimated robustly from the spin-down rate and angular velocity difference by an energy balance argument.
In this paper, we assume that the decelerating crust comes into dynamical equilibrium with the Kolmogorov cascade in the condensate. The existence of a continuously driven, statistically steady, turbulent state (driven here by electromagnetic spin-down) has been confirmed experimentally in Navier–Stokes and superfluid turbulence (Sagaut & Cambon 2008; Salort et al. 2010). Stratification can quench the turbulence intermittently under certain conditions (Chung & Matheou 2012; Lasky, Bennett & Melatos 2013), but such quenching is incompletely understood even in terrestrial contexts and falls outside the scope of this paper.
Angular momentum fluctuations
Two-point velocity fluctuations decorrelate exponentially with τ according to equation (5). There is some latitude inherent in the functional form: laboratory data and numerical simulations variously point to an exponential or Gaussian cut-off in Navier–Stokes turbulence (Comte-Bellot & Corrsin 1971; Dong & Sagaut 2008), while analogous measurements in a superfluid have never been done (Salort et al. 2010) and the role played by long-duration intermittency is still unclear (Mercado et al. 2012; Zrake & MacFadyen 2013). In the absence of definitive experiments, we adopt the exponential form here, anticipating the empirical finding that the spectrum of pulsar timing noise is observed to be red (Hobbs et al. 2010) with a power-law high-frequency tail.3 The decorrelation time-scale η(k)−1 and kinetic energy per unit wavenumber k2P(k) are proportional to γ−1 and γ2, respectively, as discussed in Section 2.2; in general, γ itself may also be a function of k. In the limit γ → 0, the turbulence is quenched. As superfluid turbulence remains poorly understood (Salort et al. 2010), especially in rotating systems, it is hard to compute γ reliably from first principles. Instead, as foreshadowed in Section 2.2, we keep it as a model parameter (always appearing in the combination ε1/3γ) and explain in Section 4 how pulsar timing noise measurements constrain it.
PHASE RESIDUALS
COMPARISON WITH OBSERVATIONS
We now undertake some preliminary comparisons between theory and data to lay the groundwork for more comprehensive population studies in the future. In Section 4.1, we verify that the theory predicts roughly the correct normalization and shape of Φ(f) for two representative objects with well-measured spectra, given sensible choices of the underlying physical variables. This is just a rudimentary consistency check; the constraints thereby derived on λ and ε1/3R−2/3γ are indicative only; the ultimate goal is to place unified constraints on these quantities across the pulsar population. In Section 4.2, we begin the latter task by examining time-domain, root-mean-square measures of timing noise in a sample of 366 objects, most of which do not yet have Φ(f) measured. We find that the theory predicts Ω and |$\dot{\Omega }$| scalings in accord with the data. The residual scatter may contain clues about how γ (and hence the physics of pinning) varies across the pulsar population. It deserves further study. We caution that coefficients like γ are governed by non-equilibrium transport processes, so the existence of a simple, one-parameter family of models (indexed by stellar mass or temperature, for example) is not guaranteed.
Power spectral density
Existing data already permit consequential tests of the theory. As an example, Fig. 1 displays the timing noise spectra of two representative millisecond pulsars, one quiet (PSR J1909−3744; lower, purple curve) and one noisy (PSR J1939+2134; upper, blue curve). The power spectral density Φ(f) (vertical axis) is plotted in units of year; we convert from the units of yr3 favoured elsewhere (Hobbs et al. 2010) by dividing by the spin period squared, so that multiple objects can be compared meaningfully on the same plot. The flat portions of the two zig-zag curves correspond to white Gaussian noise arising from measurement errors and the ephemeris fitting process, as well as possibly a component intrinsic to the pulsar. The data are post-processed by jointly whitening the low-pass-filtered phase residuals and timing model by applying a Cholesky transformation to the covariance matrix to compensate for correlated noise (Coles et al. 2011). The whitened correlations arise chiefly from inadequate calibration of the raw observations and imperfect correction for variations in the interstellar dispersion, i.e. they are predominantly extrinsic (Coles et al. 2011). The left-hand portion of the blue curve is genuine timing noise, with a red spectrum below f ≲ 3 yr−1. Overplotted are theoretical curves for λ = 3 × 10−2 and four values of ε1/3R−2/3γ specified in the caption.
A striking feature of Fig. 1 is that the predicted phase noise amplitude is high; superfluid turbulence can perturb the rotational phase of the crust at an observable level. The top (green) theoretical curve, plotted for η(R−1) = 1 yr−1, lies well above the data. Generally, at a particular observation frequency f0, the theoretical spectral power peaks for η(R−1) ≈ 1.6f0, i.e. when the spectrum rolls over near f0. The theoretical peak amplitude, |$\Phi (f_0) \approx 2 \times 10^{-7} \lambda ^{-2} f_0^{-1}$|, typically exceeds the observed spectral power by a wide margin for traditional values of the crust's moment of inertia, viz. 10−2 ≲ λ ≲ 0.5 (Lyne, Shemar & Smith 2000; Lattimer & Prakash 2007; van Eysden & Melatos 2010; Haskell, Pizzochero & Sidery 2012).
In order to pull the theoretical curve below the observations, the decorrelation frequency η(R−1) must fall well below or well above the observation band. For slow decorrelation, i.e. η(R−1) ≪ f0, we obtain a red spectrum of the form Φ(f) ∝ f−2 within the observation band. This scenario corresponds to the two middle, diagonal curves, whose parameters are chosen to match the red noise signal measured in PSR J1939+2134 [blue curve; η(R−1) = 2 × 10−4 yr−1] and to lie underneath the white noise background measured in PSR J1909−3744 to give an upper bound [purple curve; η(R−1) ≤ 7 × 10−8 yr−1]. The agreement with PSR J1939+2134 is excellent given the simplicity of the model, and the inferred limits on ε1/3R−2/3γ are reasonable for both objects, as discussed below. For fast decorrelation, i.e. η(R−1) ≫ f0, we obtain a white spectrum Φ(f) ≈ constant within the observation band. This scenario cannot explain the red noise in PSR J1939+2134, which would then arise from a different physical process, but it still constrains the turbulence parameters usefully: the bottom (brown, horizontal) theoretical curve yields an approximate lower bound η(R−1) ≥ 5 × 106 yr−1 for both objects. It is straightforward to compute the above bounds as functions of the moment-of-inertia ratio λ.
The above analysis can be extended fruitfully to other objects. High-quality power spectral density curves like those displayed in Fig. 1 are challenging to generate. We plan to undertake a systematic analysis of more objects in the near future, as more data flow out of pulsar timing array projects searching for gravitational waves (Hobbs et al. 2010; Yardley et al. 2011; Manchester et al. 2013) and timing noise experiments targeting young pulsars (Zhang et al. 2012).7
Amplitude versus spin-down rate
In the time domain, the theory predicts how far the phase wanders stochastically from the underlying, deterministic ephemeris over the observation time Tobs. The wandering is quantified (i) cumulatively, in terms of the cubic Taylor series term |$(12\pi )^{-1} \ddot{\Omega } T_{\rm obs}^3$| left over after subtracting Ω and |$\dot{\Omega }$| from ϕ(t) or (ii) progressively, in terms of the root-mean-square phase residuals σ = 〈δϕ(t)2〉1/2 for 0 ≤ t ≤ Tobs. Time-domain tests are more ambiguous than a straight measurement of the power spectral density, because there are many competing ways to subtract polynomial and/or harmonic terms from the time series, each introducing a degree of whitening that is difficult to quantify. Still, despite the risk of ambiguity, time-domain tests have certain advantages: they are quick, they can be attempted on many objects with existing data, and they are independent of frequency-domain tests, in the sense that they address the non-stationary component of the timing noise, which is explicitly subtracted to get Φ(f) in Section 3.
Fig. 2 displays σz(10 yr) = σz(Tobs = 10 yr) as a function of |$\dot{\Omega }$| for the pulsar sample analysed by Hobbs et al. (2010). In the top panel, we plot the raw data. In the bottom panel, we plot the normalized quantity |$\sigma _z(10\,{\rm yr}) (\Omega /2\pi \,{\rm s^{-1}}) (\dot{\Omega } / 10^{-13}\,{\rm s^{-2}})^{-1}$|. According to equation (24), the normalized σz should be independent of Ω, |$\dot{\Omega }$| and Tobs; its scatter should reflect the scatter in the nuclear-related quantities λ, ΔΩ and γ across the pulsar population. Do the data support this? On balance, yes. The raw σz values span more than seven decades and display a clear trend with |$\dot{\Omega }$| (Pearson correlation coefficient 0.83).9 By contrast, the normalized data span four decades and do not display a statistically significant trend with |$\dot{\Omega }$| (Pearson correlation coefficient −0.0052), as predicted by the theory. This is encouraging, given how little is known about the precise form of the temporal correlations in Kolmogorov turbulence even in terrestrial experiments, let alone a neutron star superfluid. The spread in the normalized data (355 out of 366 points in the bottom panel of Fig. 2 lie between 10−11 and 10−8 on the vertical axis) is consistent with γ varying moderately by a factor of ∼10 across the population. The rough proportionality between σz and |$\dot{\Omega }$| is also in accord with many previous studies (Cordes & Helfand 1980; Arzoumanian et al. 1994; Matsakis et al. 1997; Hobbs et al. 2010; Shannon & Cordes 2010); for example, Hobbs et al. (2010) found |$\sigma _z \propto \Omega ^{-0.40} \dot{\Omega }^{0.75}$| with Pearson correlation coefficient 0.77. Similar conclusions hold for the stability parameter Δ8 (Arzoumanian et al. 1994), which expresses |$\langle \ddot{\Omega }^2 \rangle ^{1/2}$| for Tobs = 108 s in terms of the logarithm of a dimensional quantity (essentially σzTobs).
Fig. 3 displays the whitened root-mean-square residuals as a function of |$\dot{\Omega }$| for the pulsar sample analysed by Hobbs et al. (2010). In the top panel, we plot the raw data in units of ms, i.e. the quantity labelled σ3 in Hobbs et al. (2010). In the bottom panel, we plot the normalized residuals |$\sigma _3 (\Omega /2\pi \,{\rm s^{-1}}) (\dot{\Omega }/10^{-13}\,{\rm s^{-2}})^{-1/6} (T_{\rm obs}/ 10\,{\rm yr})^{-1/2}$|, converting σ3 into a dimensionless quantity, which can be compared directly with σ in equation (25). It is difficult to disentangle the instrumental and intrinsic components of σ without further investigation (e.g. altering the instrumental configuration). Neither the raw nor the normalized data exhibit a trend with |$\dot{\Omega }$| in Fig. 3 (Pearson correlation coefficients −0.047 and −0.033, respectively), in keeping with the prediction of equation (25) but also with what one expects if the noise is instrumental. The normalized data span three decades, making it unlikely that intrinsic noise dominates instrumental noise in every object in Fig. 3; otherwise, equation (25) would imply that λ−1(ΔΩ)1/6γ1/2 spans three decades too, which is conceivable but unlikely in the light of independent empirical studies of glitch recovery time-scales (van Eysden & Melatos 2010; Espinoza et al. 2011; Yu et al. 2013), the time-averaged spin-up rate due to glitches (Lyne et al. 2000; Espinoza et al. 2011) and nuclear physics calculations (Lattimer & Prakash 2007).10 On the other hand, there is not enough evidence to support the opposite conclusion, namely that instrumental noise dominates intrinsic noise in every object.
We check equation (25) for consistency by inferring limits on λ from the data and then asking whether they are sensible on theoretical grounds. At one somewhat unlikely extreme, if every observed object is dominated by intrinsic noise, the data imply σ ≤ 8 × 10−2 (upper envelope of the points in the bottom panel of Fig. 3) and hence (λ/0.03)−1(ΔΩ/1 s−1)1/6(γ/10−4)1/2 ≤ 1.4 from equation (25). This lower bound agrees well with independent empirical and theoretical studies (Lyne et al. 2000; Lattimer & Prakash 2007; van Eysden & Melatos 2010; Espinoza et al. 2011) and is already astrophysically interesting. Reducing the instrumental component of σ3 will tighten the bound. At the other extreme, if all the observed noise is instrumental, the data imply σ ≤ 1 × 10−4 and hence (λ/0.03)−1(ΔΩ/1 s−1)1/6(γ/10−4)1/2 ≤ 2 × 10−3, requiring γ ≲ 10−10. It is intriguing to speculate whether future observations will reduce the instrumental component of σ3 or whether we are starting to see an intrinsic white noise floor. Reducing the instrumental noise ultimately creates an opportunity to falsify the turbulence model, at least in its present idealized form, if it proves possible to push the measured σ and hence the inferred γ well below a physically reasonable value, after allowing for the ambiguities inherent in the whitening process. Further, theoretical work is required to determine from first principles what the lower limit on γ should be.
CONCLUSION
In this paper, we calculate analytically the statistics of the rotational phase fluctuations produced by superfluid turbulence in a neutron star in terms of two fundamental parameters: the non-condensate fraction of the moment of inertia, λ = Ic/I0, and the decorrelation time-scale, η(R−1)−1, which depends on the steady-state angular velocity shear and the dynamical response time of the superfluid. The calculation is idealized, in the sense that the turbulence is assumed to obey the isotropic Kolmogorov law, without allowing for the undoubtedly important but poorly understood effects of buoyant stratification, fast rotation, two-component superfluidity and turbulent hydromagnetic stresses (Melatos & Peralta 2010; Salort et al. 2010; Melatos 2012). Simple formulas are given for the autocorrelation function of the phase residuals in the time (equation 14) and frequency (equation 16) domains. It is shown that the spectrum is red, consistent with radio pulsar timing data. Simple recipes are also presented for how to extract λ and η(R−1) from the half-power point of the spectrum (equations 17 and 18) or place a limit on their product from the f−2 tail (equation 19). Steeper tails can be accommodated within the theory by modifying slightly the exponential temporal decorrelation function in equation (5), a generalization that will be considered in future work.
The theory is applied to data from a representative group of ordinary and millisecond pulsars to illustrate in a preliminary fashion how the theory can be tested; a full comparison will be undertaken in a future paper. For the objects studied, the decorrelation frequency is bounded by η(R−1) ≲ 10−3 yr−1 or η(R−1) ≳ 106 yr−1, and the pinning turbulence parameter satisfies γ ≲ 10−4 or (for a limited subclass of objects) γ ∼ 1, consistent with other work (Link 2012a,b). Superfluidity enters the theory purely through γ; one has γ = 1 for Navier–Stokes turbulence and γ ≪ 1 for pinned superfluid turbulence. Hence, the preference for γ ≪ 1 implied by the data amounts to indirect yet independent evidence for superfluidity in neutron stars and warrants further study. Good agreement is obtained with popular measures of the root-mean-square phase residuals like the Allan variance σz, both with respect to the overall normalization and the spin-down trend. We show that whitened phase residuals can be used to place astrophysically interesting bounds on λ−1γ1/2. The results may find practical application to experiments currently underway to detect gravitational radiation with pulsar timing arrays (van Haasteren et al. 2011; Yardley et al. 2011; Manchester et al. 2013), chiefly by clarifying the relative strength of the reducible and irreducible components in timing noise.
Additional observational tests are needed, starting with extending the preliminary tests in this paper to more objects. Direct measurements of Φ(f) are the cleanest signature of the stationary component of the red noise, but they also require the greatest effort. Root-mean-square residuals carry time-integrated information about the low-frequency, high-power component, which cannot be resolved spectrally with existing, multidecade data sets. The challenge is to construct a stable root-mean-square statistic, which does not depend on how the ephemeris is subtracted, as many authors have noted previously (Cordes 1980; Arzoumanian et al. 1994; Matsakis et al. 1997; Hobbs et al. 2010; Shannon & Cordes 2010; Zhang et al. 2012). Equation (14), which gives the phase autocorrelation function, lends insight into what additional tests are likely to be profitable. One approach is to study the angular velocity residuals instead of the phase residuals, since the former, unlike the latter, constitute a stationary time series, with 〈δΩi(t)δΩj(t′)〉 ∝ Gij(τ) via equation (8). Baykal et al. (1999) constructed spectra for angular velocity residuals by removing quadratic and cubic trends simultaneously from pulse-frequency and time-of-arrival data. They found scalings of the form f−q, with 0.4 ≲ q ≲ 2.4, in the tail of the spectra of four pulsars with anomalous braking indices, PSR B0823+26, PSR B1706−16, PSR B1749−28 and PSR B2021+51, but with error bars on q of between ±0.5 and ±1.4, i.e. consistent with Gij(f) ∝ f−2 but inconclusive.11 It is worth testing, perhaps via Monte Carlo simulations, whether the advantage of stationarity enjoyed by 〈δΩi(t)δΩj(t′)〉 outweighs the disadvantage of differentiating numerically the ϕ(t) time series generated by the timing software. Finally, whatever the technique, the theory can be tested by observing for longer and extending the spectrum to lower frequencies, where it is predicted to rise to Φ(0) ≈ 1.2 × 10− 2(λ/0.03)− 2|$(1-\lambda )^{-1/3}(\dot{\Omega } / 10^{-13}\,{\rm s^{-2}})^{-1/3}(\Delta \Omega / 1\,{\rm s^{-1}})^{-1/3} (\gamma /10^{-4})^{-1} \, {\rm yr}$|.
Time-integrated braking indices also contain information about timing noise (Johnston & Galloway 1999; Hobbs et al. 2010). They have anomalous absolute values as large as ∼104, which manifestly do not describe magnetic dipole braking. Evidence exists that pulsars younger than ∼105 yr have predominantly positive braking indices dominated by glitch recoveries, whereas the braking indices of pulsars older than ∼105 yr are positive or negative with roughly equal likelihood and reflect some non-glitch, non-magnetic process, possibly superfluid turbulence (Urama et al. 2006; Hobbs et al. 2010). More work is needed to determine how to extract from equation (14) a time-integrated braking index, which is directly comparable to the available data.
The theory presented in this paper can be extended in several ways. First, an improved description of superfluid turbulence is required: on a local level, to calculate γ from simulations that account for the pinning microphysics, and on a global level, to account for stratification, hydromagnetic stresses and multiple superfluid components, which influence the Kolmogorov physics as well as γ. Work is underway along these directions, but the problem is formidable even under terrestrial conditions and is unlikely to be solved soon (Salort et al. 2010). Turbulence itself alters transport coefficients like the viscosity, both macroscopically through mixing length physics and microscopically through scattering in a vortex tangle. Secondly, off-axis torque fluctuations cause the rotation axis to precess (cf. Chandler wobble), with the angular displacement set by the dissipation physics (cf. Section 3). The rigid crust and corotating charged fluid are asymmetric in general under the action of elastic and hydromagnetic stresses (Melatos 2000; Mastrano et al. 2011), the pinned superfluid vorticity induces gyroscopic precession on the time-scale 2πλ− 1Ω− 1 (Shaham 1977) and secular and/or stochastic torques do not necessarily average to zero over many precession cycles (Melatos 2000). It is interesting to speculate whether the absence or presence of precession explains the different types of timing noise observed in individual pulsars, characterized as phase, frequency and torque noise in the literature (Cordes & Helfand 1980; Cordes & Downs 1985), and whether there is any correlation with pulse profile/polarization variations (Osłowski et al. 2011). Thirdly, it is worth asking whether the theory in this paper can help relate the physics of timing noise and glitches, in pulsars where both phenomena are present. For example, if the observed post-glitch recovery reflects the dynamics of the core superfluid, temporarily decoupled from the crust by strong flux tube pinning (Link 2012c), then η(R−1)−1 can be identified approximately with the recovery time-scale (van Eysden & Melatos 2010), and the normalization of Φ(f) is proportional to the recovery time-scale through equation (16), a testable prediction.
Finally, we emphasize that there is compelling evidence that timing noise is dominated by magnetospheric state switching in certain pulsars (Lyne et al. 2010) and that this phenomenon is not incorporated in the theory presented here.
The authors thank V. Ravi, G. Hobbs, W. Coles, R. Shannon and R. Manchester for stimulating discussions and for preparing the observational data plotted in the figures. The power spectral density curves for PSR J1909−3744 and PSR J1939+2134 in Fig. 1 are provided by the Parkes Pulsar Timing Array Project and have been published in fig. 11 of Manchester et al. (2013). This research was supported by the Australian Research Council and the USA National Science Foundation (Award AST-1211391), and by NASA (Award NNX12AF88G).
Crust–core corotation is not guaranteed. The hydromagnetic coupling is subtle (Melatos 2012); it is weakened by buoyancy (Mendell 1998), type I superconductivity (Haskell, Pizzochero & Sidery 2012), thermally activated vortex creep (Link 2012c), and in non-trivial magnetic geometries (Easson 1979b).
The mutual friction force is ‘patchy’ in the sense that it takes different forms locally. Specifically, it is isotropic where the vortices are tangled, anisotropic in a rectilinear vortex array and ∼103 times weaker in the former configuration than in the latter for typical neutron star parameters (Peralta et al. 2006b).
The two secular terms combine with the other three terms to give zero drift as t, t′ → 0, as required.
The observed statistics depend on the time origin of the measurements, because the random walk prior to t = t′ = 0 adulterates the random walk at t, t′ > 0 by randomizing δϕ(t) and its derivatives at t = 0, as proved in section 3c of Cordes (1980).
Landau & Lifshitz (1959) proposed ε ∝ (ΔΩ)3 for Navier–Stokes turbulence driven by a constant shear ΔΩ, which produces a more energetic flow (and hence stronger timing noise) than |$\varepsilon \propto \dot{\Omega } \Delta \Omega$| for typical pulsar parameters. Here, we stick with the latter alternative to be conservative, noting only that there is legitimate debate around what form of ε suits the boundary conditions best; see also Melatos & Peralta (2010).
Shannon, private communication.
Hobbs, private communication.
Errors in |$\dot{\Omega }$| leak into |$\ddot{\Omega }$| to leading order, when |$\dot{\Omega }$| is subtracted from the ephemeris ϕ(t), and hence may explain part of the trend in σz versus |$\dot{\Omega }$|.
The range covered by the normalized σ3 does not change significantly, when we exclude the 25 ms pulsars with periods shorter than 10 ms, whose residuals are systematically lower (σ3 ≤ 0.09 ms), and the pulsars which are known to glitch (Espinoza et al. 2011), whose quasi-exponential recoveries may pollute σ3.
Magnetar torque spectra may be analysed too, e.g. figs 10 and 11 in Woods et al. (2002). Magnetic stresses change the character of the turbulence, e.g. its effective dimensionality, a topic for future work.