Is there a spectral turnover in the spin noise of millisecond pulsars? Free

ABSTRACT

1 INTRODUCTION

It has long been proposed that pulsars can be used to detect gravitational waves in the nHz band (Sazhin 1978; Detweiler 1979; Hellings & Downs 1983). Millisecond pulsars, first discovered in Backer et al. (1982), provide promising prospects for gravitational wave detection thanks to their exceptional rotational stability. The concept of a pulsar timing array (PTA), long-term monitoring of pulse arrival times from a spatial array of millisecond pulsars, was conceived three decades ago (Romani 1989; Foster & Backer 1990). Currently, several collaborations are conducting PTA observations, including the Parkes Pulsar Timing Array (Manchester et al. 2013), the European Pulsar Timing Array (Kramer & Champion 2013), and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (McLaughlin 2013). A consortium of these collaborations is called the International Pulsar Timing Array (IPTA, Hobbs et al. 2010; Perera et al. 2019).

The first gravitational-wave signal detected with PTAs is likely to be a stochastic gravitational-wave background, formed by a cosmic population of supermassive binary black holes (Rosado, Sesana & Gair 2015). Apart from the detection of gravitational waves, PTAs also offer the opportunity to establish a pulsar-based time standard (Hobbs et al. 2012), to study the Solar system (Caballero et al. 2018), the interstellar medium (Coles et al. 2015) and the Solar wind (Madison et al. 2019), and to constrain ultralight dark matter candidates (Porayko et al. 2018).

The science output of PTA data relies on how well we model noise. Incorrect noise models can also lead to false detection in gravitational-wave searches (Arzoumanian et al. 2018b; Hazboun et al. 2020). At low frequencies, where we are most sensitive to the stochastic gravitational-wave background, some millisecond pulsars, primarily studied in the PTA context, have measureable levels of red noise (Coles et al. 2011; Reardon et al. 2015; Lentati et al. 2016; Caballero et al. 2016; Arzoumanian et al. 2015, 2018a). The red noise power spectrum is modelled by either a power law, or the broken power law, which introduces a corner frequency below which the noise power spectrum plateaus. Additional opportunities also include the free spectral model (see e.g. Lentati et al. 2013) and the power-law model with deviations at each frequency bin (Caballero et al. 2016). One particular source of red noise is the spin noise, which might be associated with pulsar rotational irregularities (see e.g. Shannon & Cordes 2010). While some young pulsars show hints of a spectral turnover at low frequencies (Parthasarathy et al. 2019), it has not yet been found for millisecond pulsars. If the typical time-scale of a spectral turnover for millisecond pulsars is on the order of years or shorter, it reduces the red noise in the most sensitive frequency band of PTAs, yielding a faster detection of a stochastic gravitational-wave background. Implications of how a spectral turnover will affect times to detection of a stochastic background were discussed in Lasky et al. (2015). One of the conclusions of Lasky et al. (2015) is that the gravitational wave power spectrum will only surpass the steeper timing noise spectrum if the latter flattens below some frequency.

Moreover, pulsar timing red noise provides interesting prospects for studying neutron star physics. A range of mechanisms have been proposed to explain pulsar red noise, including switching between two different spin-down rates (Lyne et al. 2010), recovery from a glitch – a sudden increase in the rotational frequency (Johnston & Galloway 1999), a cumulative effect of frequent microglitches (Cordes & Downs 1985; D’Alessandro et al. 1995; Melatos, Peralta & Wyithe 2008), variable coupling between the crust and liquid interior (Alpar, Nandkumar & Pines 1986; Jones 1990), influence of planets (Cordes 1993) and asteroids (Shannon et al. 2013). Nevertheless, there are not many models that link power spectral density model parameters to physical features. One such model by Melatos & Link (2013), which we explore in this paper, predicts a superfluid turbulence in neutron star interiors as the origin of red noise. The turbulent process exerts a torque on the star’s crust, where the external magnetic field of the star is produced. The model features a spectral turnover.

In this work, we employ Bayesian inference to search for evidence of spectral turnover in pulsar spin noise in the first data release (DR1) of the IPTA (Verbiest et al. 2016). We discuss our data analysis methods in Section 2. Our simulation study is presented in Section 3. We describe the noise processes of the first IPTA data release in Section 4. We present the results in Section 5, and discuss our conclusions in Section 6.

2 METHOD

2.1 Bayesian methodology in pulsar timing

Some timing model processes are covariant with red noise. In particular, in analyses by Coles et al. (2011) and Reardon et al. (2015), the least-squares timing model fit absorbs some red noise. This absorption of power causes an apparent visible turnover in the measured spectra of red post-fit residuals, which is why the model with the broken power law was used for these analyses. In Caballero et al. (2016), the regular power-law was used, as the effects of timing model fitting were taken into account. In our analysis, we employ analytical marginalization over the uncertainty of timing model parameters in equation (4), which is equivalent to the simultaneous fitting of the timing model parameters and the red noise parameters, under the assumption that non-linear dependencies of the likelihood on the timing model parameters are negligible. This avoids the problem of detecting a spectral turnover that is actually due to the timing model fit, and makes it possible to target the spectral turnover in the spin noise itself. During marginalization, one loses sensitivity at low frequencies, especially at frequencies ≤1/T_obs, due to taking the uncertainty of the timing model into account.

2.2 Modelling stochastic processes

2.3 Red noise models

In Fig. 1, we plot examples of models of spin noise power spectral density. Note, at high frequencies, Model M with two parameters asymptotically approaches Model PL with fixed γ = 2 and only one free parameter (amplitude), so parameters η(R⁻¹) and λ of Model M become degenerate. In order to break this degeneracy, and to distinguish Models PL and M, one must observe a spectral turnover. This conclusion will be important later when we find pulsars that prefer Model M over Model PL, but realize that at the current stage of observations the performance of Model M is largely determined by the consistency of Model PL’s estimate of γ with 2.

Figure 1.

Models for pulsar red noise power spectral density. The blue solid line represents Model PL (equation 13) and the orange dashed line represents Model BPL (equation 14). For both of them, we chose A = 2 × 10⁻¹³ and γ = 2. For the orange dashed line |$f_{\text{c}} = {{0.5}{\, \mathrm{yr}^{-1}}}$|⁠. The green dotted line represents Model 15 (equation 15) with |$\eta (R^{-1}) = {{0.5}{\, \mathrm{yr}^{-1}}}$|⁠, λ = 0.5, assuming pulsar spin period of |${1}\, {\mathrm{ms}}$|⁠.

Open in new tabDownload slide

In our analysis, we model N_F = 30 Fourier components of red noise processes. For power-law P(f), the fraction of the signal power above 1/T_obs that is fit with N_F components is equal to |$1 - N_{\text{F}}^{1-\gamma }$| when γ > 1. As an example, for a typical γ = 3, with 30 Fourier components we take into account |$99.9 {{\ \rm per\ cent}}$| of the red noise power above 1/T_obs. Below γ = 1.5, where 30 Fourier components take into account |$81.7 {{\ \rm per\ cent}}$| of the red noise power above 1/T_obs, it is better to use more Fourier components. In reality, after we calculate this fraction for the power up to the sampling frequency, this fraction will be greater. Nevertheless, for PSR J2145−0750, where in Lentati et al. (2016) it has been estimated that γ = 0.6 ± 0.2, we use 100 Fourier components (107 components were used in Lentati et al. 2016). We model remaining pulsars with 30 Fourier components, which is a reasonable and computationally cheap approximation. More comments on the consequences of this choice are provided in Section 5.

2.4 Software

We estimate the design matrix using the designmatrix plugin in tempo2 (Hobbs, Edwards & Manchester 2006). We simulate data and access tempo2 using libstempo (Vallisneri 2013). We construct our models and likelihood, and do parameter estimation using enterprise (Ellis et al. 2019). We perform likelihood sampling using the ptmcmcampler (Ellis & van Haasteren 2017) for IPTA DR1 data. For simulations, we use a nested sampler dynesty (Speagle & Barbary 2018), and we use bilby (Ashton et al. 2019) to access the Dynesty sampler.

3 SIMULATION STUDY

We perform a simulation study to demonstrate our ability to do Bayesian model selection. We also demonstrate some potential subtleties in recovering a low-frequency turnover. We simulate ToAs, ToA errors, and timing residuals for the pulsar PSR J0711−6830, using ephemerides from the ATNF Pulsar Catalogue (Manchester et al. 2005). We simulate ToAs evenly sampled once every 30 d between MJD 53000 and 56650, which is roughly consistent with the average cadence of a typical IPTA observatory (see Verbiest et al. 2016, table 1). We assume ToA errors to be 0.5 μs, which is within the range of ToA errors as found in the DR1 of the IPTA. These parameters are applied to all simulations described in this section of the paper. In our noise simulations we only assume one observing system, one observed radio frequency, and only red and white noise. The red noise parameters chosen for simulations are described in the following subsections. We choose them, so that they are approximately consistent with noise parameters of the real data (see e.g. Lentati et al. 2016, table 6). The parameter values recovered from simulations in this section have been confirmed to be consistent with injected values.

3.1 Red noise in an ensemble of pulsars

We simulate 50 mock pulsars with different random realizations of Model PL red noise and white noise. Then, we perform model selection between Models PL and BPL. The simulated white noise parameters throughout the subsection are EFAC = 1 and EQUAD  = 0.1 μs. According to section 3.3 of Verbiest et al. (2016), these are the typical EFAC and EQUAD values found in IPTA DR1. The simulated red noise amplitude is different for the three cases we describe in this subsection, while the priors for red noise power-law index and corner frequency are |$\pi (\gamma) = \mathcal {U}(2,5)$| and |$\pi (f_{\text{c}}) = \text{log}_{10}~\mathcal {U}(10^{-10},10^{-6})$|⁠. Here |$\mathcal {U}$| stands for a uniform distribution, and |$\text{log}_{10}~\mathcal {U}$| stands for a uniform in log₁₀ distribution. We use the same red noise priors for A and γ for Models PL and BPL, for both injection and recovery.

First, we simulate Model PL with a prior |$\pi (A) = \log _{10}~\mathcal {U}(10^{-14},10^{-11})$|⁠. The prior range for noise amplitude is chosen such that red noise is overall stronger than white noise. As a result, with all simulated pulsars, we obtain |$\log \mathcal {B}^{\text{BPL}}_{\text{PL}} = -30.8$|⁠. Hence, Model PL is correctly preferred over Model BPL.

Second, we demonstrate that we do not prefer the wrong model if the red noise is overall much weaker than white noise. The prior for simulation and recovery of red noise amplitude is reduced to |$\pi (A) = \log _{10}~\mathcal {U}(10^{-17},10^{-14})$|⁠. Now, |$\log \mathcal {B}^{\text{BPL}}_{\text{PL}} = 1.0$|⁠. Therefore, if the red noise is too weak, we cannot distinguish between two models, as expected.

Finally, we demonstrate that, when the data from multiple pulsars are injected with Model BPL, our algorithm prefers Model BPL over Model PL. To do this, we use the following prior on red noise amplitude |$\pi (A) = \log _{10}~\mathcal {U}(10^{-14},10^{-11})$|⁠. Now we obtain |$\log \mathcal {B}^{\text{BPL}}_{\text{PL}} = 96$| favouring the correct model. Our results for this subsection are summarized in Table 1. All injected signals were successfully recovered.

3.2 Prior mismatch in simulations

Most of the IPTA pulsars from DR1 are dominated by white noise Lentati et al. (2016). In this subsection, we perform simulations that demonstrate that model selection for red noise in data, dominated by white noise, can lead to the false detection of a spectral turnover, if we do not carefully choose our prior. We perform simulations of only white noise with EFAC = 1 and EQUAD  = |${{0.1}\, {\mu \mathrm{s}}}$|⁠. We perform model selection between Models BPL and PL. We observe that evidence for the absence of red noise (Model |$\varnothing$|⁠) is always the strongest, while either Model PL or BPL may be preferred, depending on our prior on f_c parameter. As we allow our prior on f_c to include only low values less than around 1/T_obs, we cannot distinguish Models PL and BPL. As we allow our prior on f_c to include only frequencies higher than our sampling frequency, we cannot distinguish between Models BPL and |$\varnothing$|⁠, and model selection between PL and BPL prefers BPL. This is not surprising, as white noise and Model PL are limiting cases of Model BPL. Therefore, for the case of the DR1 analysis, when the true distribution of spin noise parameters is unknown, we propose to account for this effect by including in equation (8) only pulsars having |$\text{log}~\mathcal {B}_{\varnothing ,i}^{\text{PL}} \ge 5$| or |$\text{log}~\mathcal {B}_{\varnothing ,i}^{\text{BPL}} \ge 5$|⁠. This way we exclude pulsars with no evidence of any spin noise and do not obtain false positives in favour of either a spectral turnover or its absence. Another solution to this problem is to fit the priors using the hierarchical inference (MacKay 2003), which we defer to a future work.

3.3 The effect of sample variance in recovery of high amplitude red noise

In this subsection, we find that with a PTA observation time of 10 yr, we are unlikely to resolve a turnover in the red noise process of any particular pulsar, assuming a fiducial f_c = 10 nHz. The is because factors a_i and b_i in equation (9) become a source of noise themselves, and we do not have a data span long enough to effectively probe residuals spectra at frequencies around the turnover.

To demonstrate this, we simulate 1000 pulsars with red noise Model PL amplitude |$\pi (A) = \log _{10}\mathcal {U}(10^{-15};10^{-11})$| and γ = 3, and simulate additional 1000 pulsars with red noise Model BPL with the same parameters and a corner frequency f_c = 10 nHz. As the amplitude of the red noise in the set of simulated pulsars increases, the average |$\text{log}\mathcal {B}_i$| in favour of the correct model plateaus. This is demonstrated in Fig. 2. We can see that, at some point, increasing |$\log \mathcal {B}(f)$| starts slightly favouring the correct model, but then saturates, so that increasing the amplitude of the red noise does not help to resolve a low-frequency turnover. In this medium-to-strong red noise regime, some realizations of Model PL may favour the Model BPL hypothesis, and vice versa. However, the mean |$\log \mathcal {B}^{\text{BPL}}_{\text{PL},i}$| (red line in Fig. 2) favours the correct model.

Figure 2.

The demonstration of the effect of sample variance on the recovery of a spectral turnover. Each point represents |$\log \mathcal {B}^{\text{BPL}}_{\text{PL},i}$|⁠. The top plot with blue points is for different realizations of a power law, Model PL (equation 13), while the bottom plot with orange points is for different realizations of a broken power law, Model BPL (equation 14). The injection parameters, except red noise injection amplitude A (horizontal axes), are the same for both plots. As the amplitude of the red noise is increased, the evidence in favour (bottom plot) and against (top plot) the spectral turnover plateaus. Red lines are mean values for every 200 simulations.

Open in new tabDownload slide

4 SOURCES OF NOISE IN THE FIRST IPTA DATA RELEASE

In this section, we describe sources of noise in the IPTA DR1 data set. We use Lentati et al. (2016) as a guide for choosing what noise terms to include in our model. In Table 2, we list the prior distributions for parameters used in our models. Then, we perform Bayesian inference of these parameters and model selection for millisecond pulsar spin noise.

4.1 White noise

IPTA pulsars are often monitored by several radio observatories. The raw voltages from each telescope are processed by different hardware. Each observing system has different measurement errors, contributing to measured white noise. Noise parameter EFAC, introduced in equation (12), accounts for ToA uncertainty, associated with errors during the process of cross-correlation of pulse profile templates with observed pulse profiles. Parameter EQUAD is introduced to account for stochastic variations in both phase and amplitude of radio pulse profiles. These variations are called ‘pulse jitter’ (Osłowski et al. 2011; Shannon et al. 2014). Parameters EFAC and EQUAD are introduced for each backend system that processes raw telescope data, in accordance with equation (12). In NANOGrav data, one epoch of observations with wide-band receivers is split into multiple ToAs, corresponding to different radio frequencies, or sub-bands. Thus, for NANOGrav data, ECORR parameters are introduced to account for correlations between sub-banded ToAs at each epoch (Arzoumanian et al. 2018b).

4.2 DM noise

Dispersion measure (DM) is the electron column density, integrated along the line of sight to a pulsar. Stochastic variations in dispersion measure result in DM noise. We model DM noise as a power law with A_DM and γ_DM, where |$\kappa _j = K^2 \nu _j^{-2}$| in equation (9). So, both κ_j and F_{i, j} depend on the radio frequency ν_j (Hz) of the j’th ToA. A constant K = 1400 MHz can be thought of as a reference radio frequency. We account for DM variations for every pulsar in IPTA analysis.

4.3 Band noise and system noise

Lentati et al. (2016) found that specific IPTA pulsars show evidence of band noise and system noise, which introduces additional red noise in some observing systems and radio frequency bands. In order to separate band noise and system noise from spin noise, we add a separate power law with A_BS and γ_BS on specific radio frequency bands and systems for specific pulsars where band and system noise for IPTA data release 1 has been found (see table 4 in Lentati et al. 2016, for details).

4.4 Spin noise

We model spin noise as a common red noise process between all observing systems and radio frequencies. Model PL depends on parameters A_SN and γ_SN, and Model BPL depends on an additional parameter f_c. We refer to a hypothesis that no spin noise is present in the data, as to Model |$\varnothing$|⁠. In this work, we are mostly interested in resolving a spectral turnover in spin noise, characterized by the parameter f_c in Model BPL. We are also interested in Model M with parameters M_SN and t_c.² When carrying out model selection between Models M and PL, we chose our prior on Model PL amplitude A to match the range of spin noise amplitudes that is allowed by our priors for η(R⁻¹) and λ in Model M. Otherwise, the model with a wider prior range on spin noise amplitude would be incorrectly penalized when calculating a Bayes factor.

4.5 Transient noise events

5 RESULTS

We perform parameter estimation and model selection for pulsars from the IPTA DR1. A summary of our analysis for individual pulsars is given in Table 3. The first column contains pulsar names and the second column contains observation spans. The next two columns, log₁₀A_SN and γ_SN, represent parameter estimates for Model PL with errors, based on 16 per cent and 84 per cent levels of marginalized posteriors. The last two columns contain the results of spin noise model selection. From the seventh column, we see that specific pulsars do not show support in favour of a spectral turnover because |$| \log \mathcal {B}^{\text{BPL}}_{\text{PL},i} | \lt 2$| for all pulsars.

Next, we employ equation (8), in order to use all available data for model selection. We perform our analysis with five different Solar system ephemeris models, as it has been found that errors in Solar system ephemerides contribute to pulsar red noise (Caballero et al. 2018; Arzoumanian et al. 2018b; Guo et al. 2019). We find that data favour neither Model PL, nor Model BPL. This result is summarized in Table 4.

Note, Tables 3 and 4 contain only results from seven pulsars where |$\log \mathcal {B}^{\text{PL}}_{\varnothing }\gt 5$| or |$\log \mathcal {B}^{\text{BPL}}_{\varnothing }\gt 5$|⁠. In table 6 in Lentati et al. (2016), authors present eight pulsars that show evidence for spin noise in their analysis. Seven of them can be found in our Table 4: PSRs J0613−0200, J0621+1002, J1713+0747, J1824−2452A, J1939+2134, J2145−0750, and J1024–0719. In the remaining PSR J1012+5307, we did find some evidence of spin noise, |$\text{log}~\mathcal {B}_{\varnothing ,i}^{\text{PL}} = 4.3$|⁠, assuming the default Solar system ephemeris DE421. However, PSR J1012+5307 did not satisfy our formal criteria to be included in Table   4. It is worth noting that in Lentati et al. (2016) pulsar J2145−0750 is found to have the most shallow power-law index γ_SN = 0.6 ± 0.2. For the reasons discussed in Section 2.3, PSR J1012+5307 only showed evidence of spin noise in our analysis after we changed a number of Fourier components N_F from 30 to 100 for this pulsar.

The last column in Table 3, |$\log {\mathcal {B}^{\text{M}}_{\text{PL},i}}$|⁠, presents log Bayes factors in favour of Model M over Model PL. We find that no pulsars show a strong support for Model M. However, PSRs J1939+2134, J1024–0719, and J1713+0747 disfavour Model M with |$\log {\mathcal {B}^{\text{M}}_{\text{PL},i}} \lt -4$|⁠.

Figure 3.

Hyperposteriors |$\mathcal {P}(\xi)$| for DR1 pulsars. Orange lines are posteriors |$\mathcal {P}^{\text{BPL}}_{\text{PL}}(\xi)$| for the fraction of pulsars that are described by Model BPL (equation 14), assuming other pulsars are described by Model PL (equation 13). Green lines are posteriors |$\mathcal {P}^{\text{M}}_{\text{PL}}(\xi)$| for a fraction of pulsars that are described by Model M (equation 15), assuming other pulsars are described by Model PL. For solid lines, we assume that spin noise in PSR J1024−0719 is intrinsic to the pulsar. For dashed lines, we assume that the apparent spin noise in J1024−0719 is caused by the second spin frequency derivative of the pulsar induced by gravitational interaction of PSR J1024−0719 with a binary companion star (Kaplan et al. 2016).

Open in new tabDownload slide

6 CONCLUSIONS

We perform Bayesian model selection to search for a spectral turnover in pulsar spin noise using the DR1 of the IPTA. We find support, with a log Bayes factor above 4, for spin noise in eight pulsars, which is consistent with Lentati et al. (2016). However, we find no evidence for a spectral turnover either in individual pulsar data or by combining different pulsars. We also fit the data to the superfluid turbulence model proposed by Melatos & Link (2013). Our results show that whereas this model is indistinguishable from the power-law model for most pulsars, it is strongly disfavoured by three pulsars, especially PSR J1939+2134 with a log Bayes factor of 110.

Based on a range of simulations, we find that one is unlikely to resolve a spectral turnover with a fiducial corner frequency of 10 nHz in any pulsar with ≈10 yr of observations. Longer data spans are required to increase the detection confidence of a spectral turnover in individual pulsars, while a larger number of pulsars with red noise can help to resolve the presence of a spectral turnover in a population of pulsars. A follow-up study using longer data sets and a larger sample of pulsars, for example the IPTA second data release (Perera et al. 2019), will prove useful in not only understanding the nature of red noise in millisecond pulsars but also in evaluating the realistic prospect of gravitational-wave detection. A more detailed simulation study is required to explore PTA configurations that would resolve spectral turnover in the individual pulsars. Whereas our simulation study assumed a pulsar with observation span of 10 yr, two pulsars from the DR1 of the IPTA have observations spans above 25 yr. At the same time, next-generation PTAs based on MeerKat, the Five-hundred-meter Aperture Spherical Telescope (FAST), and the Square Kilometre Array (SKA), will have a reduced radiometer noise. Both greater observation spans and reduced white noise levels will increase the sensitivity of a PTA to the spectral turnover, and the future study could help to estimate by how much. Simulations that attempt to provide a precise answer to these questions for the real data ought to include all other noise sources (i.e. DM noise and band noise), multiple observing backends with realistic observation cadences. Another interesting future simulation study would determine whether the broken power-law model would be favored over the power-law model when the superfluid turbulence model is simulated.

ACKNOWLEDGEMENTS

We thank the anonymous referee for valuable feedback. We thank Nataliya Porayko, Yuri Levin, Daniel Reardon, and Paul Lasky for useful discussions. BG, XZ, and ET are supported by the Australian Research Council (ARC) Centre of Excellence CE170100004. ET is additionally supported by the ARCFuture Fellowship FT150100281.

DATA AVAILABILITY

The DR1 of the IPTA underlying this article is available at gitlab.com/IPTA/DR1. The simulated data underlying this article will be shared on a reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: EXPLICIT FORM OF MODEL M POWER SPECTRAL DENSITY

APPENDIX B: POSTERIOR PROBABILITY DISTRIBUTION EXAMPLES

In Fig. B1, we demonstrate posterior distributions for spin noise parameters of two pulsars, J0621+1002 and J1939+2134, where the highest and the lowest |$\log \mathcal {B}_{\text{PL}}^{\text{BPL},i}$| are found (see Table 3 for details).

Figure B1.

This figure represents posterior distributions for spin noise parameters for J0621+1002 (left, B1a) and J1939+2134 (right, B1b). Vertical dashed lines represent 1/T_obs. For J1939+2134, with least evidence for the spectral turnover in Table 3, measurement of f_c does not affect the measurement of the amplitude and the slope of spin noise. However, for J0621+1002, with the highest evidence for the spectral turnover in Table 3, measurement of f_c does affect measurement of the power-law index.

Open in new tabDownload slide