Correcting for the solar wind in pulsar timing observations: the role of simultaneous a nd l ow-frequency observations

TL;DR

This paper investigates how simultaneous low-frequency and wide-bandwidth observations can improve correction of solar wind effects in pulsar timing, which is crucial for gravitational wave detection.

Contribution

It demonstrates that a single wide-bandwidth receiver can effectively model and remove solar wind effects, highlighting limitations for pulsars near the Sun.

Findings

Wide-bandwidth receivers can model solar wind effects.

Single observations may be insufficient near the Sun.

Two-telescope observations are often inadequate for close-passing pulsars.

Abstract

The primary goal of the pulsar timing array projects is to detect ultra-low-frequency gravitational waves. The pulsar data sets are affected by numerous noise processes including varying dispersive delays in the interstellar medium and from the solar wind. The solar wind can lead to rapidly changing variations that, with existing telescopes, can be hard to measure and then remove. In this paper we study the possibility of using a low frequency telescope to aid in such correction for the Parkes Pulsar Timing Array (PPTA) and also discuss whether the ultra-wide-bandwidth receiver for the FAST telescope is sufficient to model the solar wind variations. Our key result is that a single wide-bandwidth receiver can be used to model and remove the effect of the solar wind. However, for pulsars that pass close to the Sun such as PSR J1022+1022, the solar wind is so variable that observations at two telescopes separated by a day are insufficient to correct the solar wind effect.

Tables

Table 1. Table 1: A list of pulsars that FAST is expected to time with 100 ns timing precision. The second column is their ecliptic latitude. The third and the fourth columns is the peak amplitude of the solar wind signal in our simulation for the low frequency (270  MHz) data and high frequency (1470  MHz), respectively.

Pulsar	Ecliptic latitude(deg)	Peak amplitude for 270 MHz	Peak amplitude for 1470 MHz
J1713$+$0747	30.70	6.75($\mu$s)	0.29($\mu$s)
J1744$-$1134	11.81	22.5($\mu$s)	0.9($\mu$s)
J1911$-$1114	11.09	31.8($\mu$s)	1.2($\mu$s)
J1939$+$2134	42.30	7.7($\mu$s)	0.32($\mu$s)
J2229$+$2643	33.29	10.4($\mu$s)	0.51($\mu$s)

Equations

$$t_{\mathrm{DM}}={\lambda }^2\left(\frac{e^2}{2m_e\pi c^3}{\int }_0^d{n}_{e}dl\right)$$

$$h=\left(\frac{A^2}{12{\pi }^2}\right){\left(\frac{f}{f_y}\right)}^{2\alpha -3}$$

Full text

1 Introduction

Pulsar timing arrays (PTAs) projects, such as the Parkes Pulsar Timing Array111The Parkes Pulsar Timing Array (PPTA) in Australia is an implementation of the PTA concept based on observations with the Parkes 64-m radio telescope. 24 millisecond pulsars is being observed at three radio-frequency bands, 50  cm ($\sim$ 700  MHz), 20  cm ($\sim$ 1400  MHz) and 10  cm ($\sim$ 3100  MHz).(Manchester et al. 2013), aim to detect ultra-low-frequency gravitational waves by observing a large number of millisecond pulsars (MSPs) (Shannon et al. 2015). To do this, these pulsars are observed every few weeks over decades. Differences between the measured and predicted pulse times of arrival (ToAs) are referred to as “timing residuals” (Hobbs et al. 2006). Gravitational waves are usually not included in the timing model and therefore the existence of such waves will mean that the ToAs cannot be perfectly predicted with the model and hence gravitational waves will induce timing residuals. Other phenomena can also lead to discrepancies between the measurements and the predictions. For instance, the timing model may be incorrect, the clock used for determining the ToAs may have an error and the pulsars themselves may not be perfectly stable (Hobbs et al. 2004). In many cases the dominant noise process are variations in a pulsar’s dispersion measure (DM) caused by the fluctuations in the interstellar medium or in the solar wind (You et al. 2007a). The challenge for the PTA experiments is therefore to identify the signature caused by the gravitational wave signal in the presence of these other noise processes. This can be done by looking for correlations between the timing residuals for different pulsars (Shannon et al. 2015). However, Tiburzi et al. (2016) showed that, for a realistic array of pulsars, phenomena such as errors in the solar system ephemeris, instrumental errors or uncorrected solar wind effects potentially can lead to false detections of the gravitational wave signal.

The dispersive group delay for a pulsar signal, $t_{\rm DM}$, is related to the integral of the electron density, $n_{e}$, from the pulsar to the Earth:

[TABLE]

where ${\rm DM}=\int_{0}^{d}n_{e}\,dl$ is the ‘dispersion measure (DM)’, $\lambda$ is the observing wavelength, $c$ is the speed of light, $e$ is the electron charge, $m_{e}$ is the mass of an electron and $d$ is the distance between the Earth and the pulsar. The DM at a given time can be determined if pulse arrival times have been measured at multiple observing frequencies. The precision with which the DM can be measured is determined by the arrival time uncertainties and the frequency coverage of the observations. If the DM variations can be measured with sufficient precision then such variations can be accounted for in the pulsar timing model. Methods to do this for the PPTA sample of pulsars were described by You et al. (2007b), Keith et al. (2013) and Reardon et al. (2016)222PTA experiments in Europe and in North America have also studied DM variations; see Desvignes et al. (2016) and Lam et al. (2016).. These papers showed that the DM variations can often, but not always, be modelled by assuming Kolmogorov turbulence in the interstellar medium. However, Keith et al. (2013) showed that, without care, the method used to account for DM variations can remove an underlying gravitational wave signal. By simultaneously modelling the frequency-dependent signal (the DM variations) along with a frequency-independent signal (known as the “common mode” and would include the gravitational wave signal), it was shown that the gravitational wave signal would not be affected. This method works well when the DM variations are relatively smooth and the changes occur over the timescale of many observations (i.e., weeks to months).

You et al. (2007a) and You et al. (2012) studied observations of PPTA pulsars whose line-of-sight pass close to the Sun. They demonstrated that, at an observing frequency of 1400 MHz, the solar wind can contribute ToA errors of 100 ns for sources 60∘ from the Sun and more than 1$\mu s$ within 7∘. In contrast to the variations caused by the turbulence in the interstellar medium, the DM variations caused by the solar wind depend strongly on the line-of-sight angle to the Sun and change on various time-scales including day-to-day variability and multi-year variability because of the solar cycle.

The standard pulsar timing software packages such as tempo and tempo2 (Hobbs et al. 2006) include simple models to account for the solar wind. The default model is a time-independent, spherically symmetric model completely defined by the angle of the line-of-sight to the Sun and a single electron density parameter (NE1AU: the electron density at 1 A.U.). You et al. (2007b) demonstrated that this parameterisation was not suitable for modelling the actual observations333As part of our work, we updated this model to allow changes in the electron density parameter as a function of time, but we determined that even this model was too simplistic to model real changes in the solar wind.. You et al. (2007b) attempted to model the measured DM changes using observations from the Wilcox Solar Observatory444The Wilcox Solar Observatory (WSO), located in California in US has been collecting solar polar field data since 1975. The solar field data can be used to predict the magnitude of the next cycle and the peak of the current cycle.. The WSO models do predict variations in the DM changes caused by the solar wind that are similar to what we observe, but as shown by You et al. (2007b), they are not sufficient for correcting a given pulsar data set.

To date, the PPTA observations have been carried out either with a relatively small bandwidth centred in the 20 cm observing band, or using a dual-band (10/50  cm) receiver. In a typical observing session, a given pulsar is first observed using one of these receivers and then (perhaps a few days later) observed using the other receiver. Although this is not regularly done, we can also obtain estimates of a pulsar’s DM at a given time using low frequency telescopes, such as the Murchison Widefield Array (MWA) (Bhat et al. 2014)555The Murchison Widefield Array (MWA) is a Square Kilometre Array (SKA; Dewdney et al. 2009) precursor telescope at low radio frequencies in western Australia.. This provides very low frequency (80-300 MHz) observations, but not necessarily at the same time as the Parkes observations. We are also currently commissioning an ultra-wide-bandwidth receiver for Parkes that should enable us to simultaneously observe from 700 MHz to 3GHz. Many other observatories are also building ultra-wide-bandwidth receivers. For instance, the Five-hundred-meter Aperture Spherical Telescope666The Five-hundred-meter Aperture Spherical Telescope (FAST), situated is GuiZhou province in China, is the world’s largest single dish radio telescope and opened on 25th September 2016. will have a receiver providing a band from 270 MHz to 1620 MHz (Zhang et al. 2016).

In this paper, we:

• •

demonstrate how realistic solar wind effects can be simulated within the tempo2 software package

• •

simulate (1) realistic observations for the ultra-wide-band receiver with the FAST telescope and (2) observations in the low-frequency band with MWA and in the 20 cm band with Parkes telescope.

• •

trial different methods for measuring and removing the DM variations and test whether any underlying gravitational wave signal would be affected.

• •

discuss the implications of these results for FAST and the PPTA.

2 Simulation

We used the ptaSimulate software package to simulate pulsar parameters and ToAs. This software package was developed in order to make realistic simulations of pulsar timing array data sets. For our work, we choose only to make a simple simulation of the radiometer noise (leading to the errors on the ToAs), but realistic simulations of an isotropic, stochastic gravitational waves and of the solar wind. For the radiometer noise we simply assume that the ToA uncertainties are 100 ns which is typical of a “good” pulsar in the PPTA. For the gravitational waves we use the routines described by Hobbs et al. (2009) to simulate a background defined by:

[TABLE]

where $\alpha=-2/3$, $f_{y}=1/{\rm 1yr}$ and we choose $A=10^{-14}$. This amplitude is much larger than our current upper bounds (Shannon et al. 2015), but is used to demonstrate the possibility that the gravitational wave signal can be affected by the correction method used to model the solar wind. The solar wind is modelled based on the tempo2 implementation described in You et al. (2007a) which uses WSO data as input. We have simulated various pulsars in the PPTA sample, but for this work, we concentrate on PSR J1022$+$1001 which has the lowest ecliptic latitude ($-0.06^{\circ}$) of the pulsars in the sample. We use ptaSimulate to provide three noise models: (N1) includes the radiometer noise and gravitational waves, (N2) the radiometer noise and the solar wind and (N3) all three noise processes.

For each of the noise processes we make two data sets:

• •

(DS1): we simulate the FAST telescope and the ultra-wide-band receiver covering a range from 270 to 1620  MHz and we assume that we obtain a single ToA for every 300 MHz interval.

• •

(DS2): we simulate the Parkes telescope observing in the 20  cm and the MWA telescope observing in lower band.

As shown in Figure 1, for DS1, we simulate 270, 570, 870, 1170 and 1470  MHz as suggested by Table 2 in Zhang et al. (2016). We assume that we obtain a single ToA for every 300 MHz interval777Pulsar timing has not yet been carried out with such a wide instantaneous bandwidth and so it is not clear exactly how the band will be divided. We note that the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) project (Arzoumanian et al. 2016) obtain a ToA for each 4MHz subband. However, processing such a large number of ToAs is impractical with the computing resources that we have available and therefore we simulated ToAs every 300 MHz.. For DS2, we simulate 150 and 180 MHz for MWA according to the Table 1 in Tingay et al. (2013) and Table1 in Loi et al. (2015) and 1300,1400 and 1500 MHz for Parkes according to the Table 1 in Manchester et al. (2013).

Even though the Parkes and MWA observations provide a frequency span from 80 to 3100  MHz, it is not common for both telescopes to observe simultaneously. Therefore, for DS2 we have also tried different time offsets between the MWA and the Parkes observations: (TO1) the observations were simultaneous; (TO2) the MWA observation occurred one day after the Parkes observations. All of our simulated data sets have a data span of ten years and an observing cadence of 14 d. We simulated from 2005 to 2015 as we use real WSO data sets for the Solar wind modelling - clearly the actual FAST and Parkes/MWA data do not cover this time range, but the results will not depend upon the exact dates chosen.

In order to illustrate the effect of the solar wind we have also simulated the highest and lowest (270 MHz 1470 MHz) observing bands for FAST, but with a 1 day sampling. In Figure 2a we show a typical realisation of the gravitational wave background signal over the 10 year data span. In Figure 2b we show the simulated signal with the gravitational wave background and the solar wind. The solar wind clearly produces timing residuals many orders of magnitude larger than the gravitational wave signal. As expected the signal has an annual periodicity representing the times when the line-of-sight to the pulsar goes close to the Sun. The exact signal depends upon the solar cycle and the exact data sampling. In the inset (b1) and (b2) to Figure 2b we show a zoom-in around year 2008 for 1470  MHz and 270  MHz. In this panel we have not attempted to remove the effect of the solar wind using any method. In Figure 2c we show the same data, but, this time, use the standard tempo/tempo2 spherically symmetric solar wind model in which we assume that the electron density at 1AU is $10$ e cm*-3*. Clearly, significant (both positive and negative) residuals are still seen and so this simple model both over-corrects and under-corrects the data at different epochs and does not provide a sufficient means to remove the solar wind variations.

3 Method and Results

In Figure 3, we show the power spectral density of the pre-fit residuals for DS2. The green line shows the spectral density for the timing residuals affected by the solar wind and radiometer noise. The purple line is the same, but for the simulation containing only the gravitational wave and radiometer noise. The dashed line is the analytic representation of the power spectrum of the simulated gravitational waves. Producing a single spectrum using an identical method is challenging in this case as the gravitational wave signal has a steep power-law signal which requires that a pre-whitening method to mitigate the effects of spectral leakage. However, the solar wind signal does not follow a steep power-law spectrum and therefore we should not use pre-whitening methods. Our analysis here is to remove the effect of the solar wind without removing the gravitational wave signal. We therefore present these power spectra in a way that is optimal for the steep gravitational wave power law spectrum and note that the spectra of the solar wind will not be perfect (particularly at low frequencies). Throughout this paper we therefore input a set of pulsar timing residuals and form a power spectrum assuming a gravitational wave background power law spectrum using the “analyticChol” and “cholSpectra” tempo2 plugins888These two plugins use Cholesky fitting routines to obtain a power spectral of the timing residuals. Note that this same procedure was followed by Tiburzi et al. (2016). Clearly, without correction the solar wind effect dominates the entire spectrum and makes the gravitational waves undetectable. Our aim in this section is to determine whether existing methods can be used to remove the effect of the solar wind, without also removing the gravitational wave signal.

In order to mitigate the effect of DM variations caused by the interstellar medium, Keith et al. (2013) used the simultaneous fitting of a frequency-dependent and frequency-independent function to the residuals. The function is sampled at specified time intervals. The user is able to choose the time sampling, but they reported optimal time intervals to correct DM variations caused by the interstellar medium. These optimal time intervals are long and, for the PPTA pulsars, are listed in their Table 2. We have trialled 30, 100, 200 and 300 days sampling. However, as seen from the inset to Figure 1b, the solar wind induces timing residuals on much shorter time scales. As we expected (and shown in Figure 4), the correction is poor999Here we only show the results corresponding to the 200 days sampling, but the other sampling intervals give similar results.. It is clear that the solar wind features in the power spectrum have not been removed and the power spectrum remains orders-of-magnitude higher than the power spectrum of the simulated gravitational wave signal.

The Keith et al. (2013) method does not require a common gridding between the frequency-dependent and the frequency-independent functions. We expect that the frequency-dependent function (representing the solar wind) will vary quickly, but the frequency-independent function (such as the gravitational wave signal) can be sampled on a much sparser grid. We therefore fitted using a grid spacing of 10 days for the DM (frequency-dependent) function combined with a 300 days sampling for the common mode (frequency-independent) function.

We trialled this method on both DS1 and DS2 and present the resulting power spectra in Figure 5. We demonstrate that using a sampling (14 days) for the frequency-dependent function that implies that we can model any DM variations occurring during a single observing epoch. We simultaneously fit a 300 days frequency-independent function. In Figure 5a we present the pre- and post-fit spectra for our simulation of FAST (DS1)101010Note that the fitting did not converge with a single iteration and as the gravitational wave signal induced residuals significantly larger than the 100 ns error bars and it was necessary to use the Cholesky fitting routines as described by Coles et al. (2011) when carrying out the fit.. In Figure 5b we show the same, but for the Parkes/MWA data set in which both telescopes are assumed to observe at the same time. Both panels (a) and (b) indicate that this process works well; we can remove the signature of the solar wind without removing the gravitational wave signal. However, in Figure 5c we show the resulting power spectra if the Parkes and MWA observations are offset in time by 1 day. In this case it is clear that we are unable to completely remove the solar wind variations (because the solar wind varies on a time scale faster than a day; see Figure 5 in (You et al. 2007b) for a demonstration of this with actual observations) although we note that the lowest frequency channels in the power spectrum are now at the expected level of the gravitational wave signal, also highlight that our simulated gravitational wave signal is at least one-order-of-magnitude higher than the currently expected level.

4 Discussion and conclusion

Tiburzi et al. (2016) highlighted that when the line-of-sight to pulsars at low ecliptic latitudes passes close to the Sun then the effects of the solar wind can become large. Typically the PTA experiments simply exclude observations within a certain angle from the Sun, but with sufficiently high precision (as will be achievable with telescopes like FAST and the SKA) simply excluding observations will not always be practical. Tiburzi et al. (2016) attempted to model the solar wind using a 100 days gridding, but noted that this was ineffective as the variations caused by the solar wind are narrow, cuspy and change from year to year. Other correction/modelling methods, such as the Bayesian analysis methods (Lentati et al. 2014) also assume simplistic models for the solar wind.

Our results indicate that the Keith et al. (2013) procedure can be used to model and remove the solar wind effect assuming that the frequency-dependent function has the same gridding as the observing cadence and a wide-bandwidth receiver is used. This is promising for FAST assuming that the wide-bandwidth receiver is available for high-precision pulsar timing. As the primary goal for early science operations for FAST is to carry out a survey with a multibeam, 20 cm receiver, it is possible that the wide-bandwidth receiver will not be commonly available. In this case it will become much harder to correct for the effect of the solar wind and as shown in the Parkes-MWA simulations it will not be trivial to model the solar wind sufficiently using other telescope observations (unless they are obtained simultaneously with the FAST data).

The sampling required for the solar wind correction is much faster than the optimal sampling given by Keith et al. (2013) for modelling and removing DM variations caused by the interstellar medium. When processing actual PTA data sets it will therefore be necessary to measure the DM variations on an uneven grid spacing. At times when the line of sight to the pulsar goes close to the Sun then a fast sampling will be required. At other times the Keith et al. (2013) sampling will need to be used. As the solar wind is variable and also the interstellar medium variations can have abrupt changes (for instance, Coles et al. (2015); reports observations of two extreme scattering events) it is likely that the grid spacing will need to be chosen carefully and manually.

We note that the pulsar we chose in this paper (PSR J1022$+$1001) is the PTA pulsar that has the smallest ecliptic latitude and therefore the solar wind will have the largest effect for this pulsar. Zhang et al. (in preparation) suggests that five IPTA pulsars (PSRs J1713$+$0747, J1744$-$1134, J1911$-$1114, J1939$+$2134 and J2229$+$2643) are likely to have ToA uncertainties smaller than 100 ns during a typical FAST timing observation. We have re-simulated our daily-sampled FAST data set for these pulsars. Clearly the signal strength amplitude is related to the ecliptic latitude. However, as shown in Table 1, for all pulsars the solar wind effect is detectable in the simulated data sets and for some pulsars (such as PSRs J1911$-$1114 and J1744$-$1134) the signal is significantly higher than the 100 ns level.

We also note that the amplitude of gravitational waves we set is much larger than our current upper bounds (Shannon et al. 2015). To make sure our method is effective for detecting gravitational wave, we re-simulate data for Figure 5b with a lower value of the background amplitude, $A=10^{-15}$. The results which is shown in Figure 6 and we note that the algorithm presented here works well even with the lower gravitational wave background amplitude.

For this paper we have studied the use of the Keith et al. (2013) method for modelling and removing the effect of the solar wind from pulsar timing array data sets. Other methods are available. For instance, the “DMX” procedure (e.g., Demorest et al. 2013) can be used in a similar manner, but the common-mode signal must be accounted for in the fitting process to ensure that the gravitational wave signal is not removed. Similarly, Bayesian methods exist (e.g., Lentati et al. 2014), but producing a model for the solar wind is non-trivial and, to our knowledge, has not been attempted within the Bayesian frameworks.

In conclusion, the solar wind will induce timing residuals that will need to be modelled and/or removed when analysing PTA data. Our work suggests that this can be done with a wide-enough, simultaneous observing frequency coverage. We also note that the solar wind is intrinsically interesting. Instead of simply removing and discarding the solar wind signature, it can be used (as shown by You et al. 2012) to improve the accuracy of magnetic field estimates of the solar wind near the sun.

5 Acknowledgements

This work is supported by West Light Foundation of CAS(No.XBBS201322) and National Natural Science Foundation of China(Nos.11403086 and U1431107), the Strategic Priority Research Programme(B) of the Chinese Academy of Sciences(No.XDB23010200). The FAST FELLOWSHIP is supported by Special Funding for Advanced Users, budgeted and administrated by Center for Astronomical Mega-Science, Chinese Academy of Sciences (CAMS). We thank Prof. Zong-hong Zhu for facilitating this project and allowing Ze-xi Niu to spend time in Australia working with the Parkes Pulsar Timing Array project team. We specifically thank Dr R. Shannon for useful comments on this project. We also thank Dr. Xing-Jiang Zhu for his great help.