Bayesian Solar Wind Modeling with Pulsar Timing Arrays

(The NANOGrav Collaboration)

Corresponding author: Jeffrey S. Hazboun

[email protected]

Introduction

Radio observations of distant astrophysical sources have long been used to study the content and characteristics of the solar wind. For example, Hewish & Dennison (1967) used observations of scintillation from quasars as a way to probe the structure of the solar wind. The dispersion of pulsed radio emission from pulsars, due to the diffuse ionized medium along the observational line of sight (LOS), has been known from their first radio observations [31]. The utility of these radio pulses for investigations of electron density was used in Counselman et al. [19], only two years after the discovery of pulsars, to measure the solar electron density ($10\ \mathrm{cm}^{-3}$ at 1 au) based on a spherically symmetric model for the solar wind.

Since these early investigations, pulsar astronomers have often included a model for the solar electron density as a part of pulsar ephemerides—see, for example, the three main pulsar timing software packages Tempo, Tempo2, and Pint [55, 33, 50]. In Lommen et al. [49] and Splaver et al. [60], the solar wind signal was shown to be highly covariant with astrometric components of the timing models for individual pulsars, in particular the parallax and sky position, because these parameters contain strong Fourier components at 1/yr and higher harmonics. Beyond the importance of a solar wind model for accuracy in pulsar astronomy, You et al. [77] and You et al. [76] showed the sensitivity of pulsar data sets to more complex (than the spherically symmetric $1/r^{2}$ wind) features in the solar electron density². More recent work [51, 66, 67] has shown the potential of pulsar timing arrays (PTAs) as independent probes of the solar wind and its behavior as a function of time.

In this paper we introduce new methods, expanding upon those used in Madison et al. [51] and Tiburzi et al. [67], to obtain information about the solar electron density from PTA data. Using a fully Bayesian framework, we show that much more information about the solar wind can be obtained from the same set of pulsar data used in Madison et al. [51], the NANOGrav 11-year Data Set Arzoumanian et al. [8]. As in these recent studies, it is important that we are using an array of pulsars. This allows us to separate the variations in electron density of the ionized interstellar medium (ISM) along the kiloparsec distances to pulsars from the local fluctuations in electron density due to variations in the solar wind.

The NANOGrav 11-year Data Set consists of high-precision time-of-arrival (TOA) measurements 45 millisecond pulsars spanning up to 11.4 years. Each pulsar was observed at approximately a monthly cadence, over widely separated radio frequencies in order to measure both TOAs and dispersion of the pulsar signal at each observing epoch. Observations were made using the Arecibo Observatory and the Green Bank Observatory. Typical TOA measurement precision was in the range 0.1 to 1 $\mu\mathrm{s}$. Further details of the data set are in NG11.

Monitoring of Solar Electron Density

Since the initial use of astrophysical radio sources to measure solar wind density via scintillation [32], a program for monitoring solar wind densities, especially at higher solar latitudes, has continued [16, 52, 68]. In addition to ground-based monitoring, a considerable number of resources have been used to study the solar wind from space. Many of these spacecraft have included instruments for electron density measurements, including Ulysses’s Solar Wind Observations Over the Poles of the Sun (SWOOPS; Bame et al. [13]), the Orbiter Retarding Potential Analyzer (ORPA) instrument mounted on the NASA Pioneer Venus Orbiter spacecraft (Knudsen et al. [39]) and the Advanced Composition Explorer (ACE) with its Solar Wind Electron Proton Alpha Monitor (McComas et al. [53]). Similar to pulsar timing measurements, the Viking mission used dual frequency delays to make early measurements of the latitudinal dependence of the solar electron density (Muhleman & Anderson [54]). The Parker Solar Probe (Bale et al. [11]) with its Solar Wind Electrons Alphas and Protons (SWEAP) instrument (Kasper et al. [38]) is revolutionizing our understanding of the solar wind in the inner solar system (Bale et al. [12]).

Apart from the Viking data, these space missions take in situ measurements of the solar electron density. These allow for one variety of long-term monitoring of the solar wind and Ulysses data have been used to study its structure and behavior at high to mid latitudes (Issautier et al. [36]). Obviously, this monitoring requires the use of costly spacecraft, and while the data from pulsar timing does not provide the same type of fine-grained local spatial information, the many lines of sight to pulsars allow for omnidirectional integral monitoring of the solar electron density from ground-based facilities and the scrutinization of various models for its structure and evolution. PTAs provide a probe of the solar wind that is distinct from and complementary to space missions.

Noise Mitigation in PTAs

The main goal of PTA experiments is to use precise long-term measurements of millisecond pulsars to observe gravitational-waves (GWs) in the nanohertz regime. Many pulsars are used since the unique noise properties of individual pulsars necessitates the corroboration of common³ astrophysical signals across multiple sources and because lower-amplitude common signals can be drawn out of the noise of more pulsars. Additionally, PTAs allow us to search for other common signals, for instance the motion of the solar system barycenter (Vallisneri et al. [69]) or errors in terrestrial time standards (Hobbs et al. [34]).

One motivation for this work is to develop a pan-PTA solar wind model as one component of the next generation of PTA noise models. Dispersion measure (DM) variations are one of the largest noise sources in PTA data, ([42, 18, 17, 37], and modeling these variations is an important part of PTA noise mitigation strategies—it is, in fact, a major driver of PTA observing strategies ([59, 44]). DM values in NG11 vary from $\sim3 - 300\ \mathrm{pc}/\mathrm{cm}^{3}$, while the variations are usually on the order of $10^{-3}\ \mathrm{pc}/\mathrm{cm}^{3}$. While a piecewise binned model of DM variations (DMX, see Arzoumanian et al. [7]) has served NANOGrav well, it removes up to a third of the power through the timing model transmission function ([29, 7, 37]), a problem endemic to any DM variation model that is part of the timing model fit. Additionally it has been shown that asynchronous multiband observations can lead to misestimation of DM variations (Lam et al. [41]) when binned together, either in DMX or an interpolation basis. Niu et al. [56] demonstrated that observations only a day apart can be insufficient to remove sharp solar wind effects, or “cusps”. In Hazboun et al. [30] and Lam et al. [47] the effect of sharp unmodeled cusps in DM variations was demonstrated. Mismodeling of solar wind cusps could then also present as broadband white noise in pulsar datasets. In particular, this could adversely affect GW searches for single sources which are especially dependent on the high frequency noise floor of PTAs ([40, 43] (Lam et al. 2018b)).

Perhaps most importantly for PTAs, in addition to the noise introduced by mismodeling DM variations, it was shown in Tiburzi et al. [65] that the solar wind can actually manifest spatial correlations amongst the pulsars, a potential source of confusion when searching for GWs, especially a GW background which presents in the data as a low-frequency spatially-correlated process. In this paper, we demonstrate how a fully Bayesian framework, built upon the extensive analysis infrastructure developed by PTA collaborations to search for GWs, can isolate the signal from the solar wind in PTA data. This allows for precise noise mitigation as well as monitoring of the overall solar electron density both as a function of time and observational LOS.

The paper is organized as follows. In § 2 the DM parameter will be derived along with various analytical models for a spherically symmetric and static solar electron density. In § 3 we present the full Bayesian model for pulsar timing data and how we incorporate our solar wind models. In § 4 three models for the behavior and structure of the solar electron density will be presented along with the results of using these models to analyze NG11. In § 5 we discuss the incorporation of these models into PTA data analyses as well as how these models might be used in the future to bolster efforts to study solar physics.

Solar Electron Density Models

The simplest and most common model for the solar electron density used in pulsar timing software packages like TEMPO, TEMPO2 and PINT ([55, 33, 50]) is the spherically-symmetric, time-independent expression $n_e(\vec{r}) = n_E (1\,\mathrm{au}/r)^2$, where $n_E$ is the electron density at 1 au, herein measured using units of $1/\mathrm{cm}^{3}$ to match pulsar timing software. The DM parameter used in pulsar timing is the column density, i.e. the integral of $n_e(\vec{r})$ along the line-of-sight to the radio source

where the first line is the generic expression, and the last equality uses the relationship $\frac{z_{\oplus}}{b}=\cot(\theta_i)=\tan\left(\frac{\pi}{2}-\theta_i\right)$ and $b=r_{\oplus}\sin\theta_i$. See Figure 1 for the definitions of these variables.⁴ The time delay is then dependent on the radio frequency and the DM by

DM is usually measured in units of $\mathrm{pc}/\mathrm{cm}^{3}$, and the constants of nature in front combine to $4.15\times10^{3}\ \mathrm{MHz}^{2}\,\mathrm{pc}^{-1}\,\mathrm{cm}^{3}\,\mathrm{s}$.

In You et al. [77] the authors use a two phase model assuming that the LOS to the pulsar crosses the fast and slow solar wind. Higher order terms, in addition to the simple $1/r^2$ model are used to ameliorate the effect of realistic solar wind time delays. The slow solar wind model⁵,

was developed in You et al. [77] by combining the near and far distance models in Allen [5] and Muhleman & Anderson [54]. A separate model was used for the fast solar wind

developed in Guhathakurta & Fisher [26, 27]. These models were then combined with measurements of the portion of the line-of-sight in the fast and slow wind to better estimate the cusps to the pulsar timing data. In You et al. [76] the fraction parameter for the two phases of the model is fit for as part of the analysis. Additionally a time dependent component is investigated. The model in You et al. [77, 76] has potential as a useful tool for removing the solar wind signal from pulsar timing data, especially for single pulsars, where disentangling the DM variations from the ISM is difficult.

Recently, however, Tiburzi et al. [66, 67] showed that the simple $1/r^{2}$ model for $n_{\mathrm{E}}$ commonly used in pulsar timing software packages and the You et al. [76] model are both insufficient for accurately modeling dispersion delays from the solar wind. Even with allowances for time variability, these models still fall short of the mitigation needed for the low frequency radio data taken by the LOw-Frequency ARray (LOFAR) [61, 71] used in the analyses.

In Aksim et al. [2] the authors compare their measurements of DM using Very Long Baseline Interferometry (VLBI) observations of quasars to both a spherically symmetric model and a numerical integral of the Alfvén Wave Solar Model (AWSoM) [70]. As in other studies the AWSoM matches well with the astrophysical measurements, however this model is currently too cumbersome to be used effectively by PTAs. Aksim et al. [2] also fit for a variable power spherically symmetric model and find reasonable agreement with their time delay measurements.

General Spherically Symmetric Solar Wind Density

Here we derive generic expressions for the dispersion measure for any spherically symmetric electron density fall off. While these expressions will assume the source of the electrons is the sun, they are general enough to be used for any spherically symmetric source of streaming electrons along the LOS to a pulsar, a stellar wind from a star orbiting a pulsar, for instance. Relaxing the assumption that the electron density around the sun drops off as $1/r^{2}$, we model the drop off with a more generic power law dependence. These calculations are similar to those in Aksim et al. [2], where the analytical expression for the time delay was calculated directly.

Assuming that $p$ is the exponent in the power law dependence we go through a similar calculation as in Equation (4),

where ${}_2F_1$ is a hypergeometric function. This expression can be further simplified by assuming the pulsar is sufficiently distant that we can take the limit $z_p \to\infty$. This gives

where $\Gamma(x)$ is the Gamma function. Compare this expression to the time delay equation given in Aksim et al. [2]. Equation (11) simplifies to Equation (4) when $p = 2$. See Appendix § 8 for expressions with higher integer values of $p$. The relationships below Equation (4) can be used to write this in terms of $\theta_i$. One can choose to build a solar electron density model using multiple summed terms, as in You et al. [77, 76]; Allen [5]; Muhleman & Anderson [54]; Guhathakurta & Fisher [26, 27], or allow the index to vary, as in Aksim et al. [2].

Bayesian Methods

The solar wind modeling presented here relies on the analysis infrastructure developed by PTAs for GW searches ([72, 48, 64, 74, 73]). While a gravitational wave background only causes residuals of tens of nanoseconds, the solar wind signal is much stronger ($\sim1.7\,\mu\mathrm{s}$ at 1 GHz with $\theta_i = 10^\circ$). This allows us to undertake an in depth study of the solar wind by modeling the solar electron density across all of the pulsars simultaneously.

To decouple the variations within the ISM along the LOS to each of the pulsars from the solar electron density that is varying locally, we use a deterministic solar wind model based on the expressions in § 2.1. The solar wind electron density is fit as a global parameter across all pulsars, while the geometry of the solar wind model is dependent on the details of the individual Sun-Earth-Pulsar angles, $\theta_i$, shown in Figure 1. In addition to this model each pulsar is also fit with a quadratic polynomial in time over the data set for the DM, encoded with the pulsar timing parameters, DM1 and DM2, as well as an additional Gaussian process model that emulates the variations of the ISM, which will be discussed in § 3.2.

The PTA analysis framework includes a timing model marginalization, based on a linearized timing model (van Haasteren & Levin [72]) and a full noise treatment that includes a generalized covariance matrix, parametrized by TOA errors, three additional white noise parameters and a power law red noise model (Lentati et al. [48]; Taylor et al. [64]). The PTA likelihood model is implemented using the ENTERPRISE software suite (Ellis et al. [23]) and the various models and extensions compiled in enterprise_extensions (Taylor et al. [63]). We use the Parallel Tempering Markov-Chain Monte Carlo Sampler, PTMCMCSampler, for numerically integrating our likelihood (Ellis & van Haasteren [22]).

It should be noted that while this formalism is focused on high precision millisecond pulsars, the code infrastructure can be used on data from any type of pulsar and its timing data. Observations of the more numerous canonical pulsars may be useful for these types of solar investigations in the future.

ACE SWEPAM Prior

As mentioned in § 1, a number of past and ongoing space missions have collected solar electron density data over the past few decades. In order to take advantage of this large base of knowledge about the solar electron density we used the solar electron density data from the ACE Solar Wind Electron Proton Alpha Monitor (McComas et al. [53]) to build an informative prior for $n_E$, the solar electron density at 1 au. We binned the data from the same time span as NG11, corrected the density from the L1 Lagrange point using a $1/r^2$ density model⁶ and used an empirical distribution to build a random variable in SciPy. The distribution and median are shown in Figure 2. The ACE SWEPAM in situ data are a useful prior for limiting the extent of parameter space at which we expect to be working, roughly $0.01$–$50\ \mathrm{cm}^{-3}$. However, the electron flow in the ecliptic plane is known to be different than the fast flow at other solar latitudes [54, 13, 26]. We do not expect the ACE measurements to match our results perfectly, but it is a good comparison in these analyses as a sanity check. In practice the prior was only informative when compared to the data in some individual pulsar runs — usually pulsars without strong solar wind signals in their TOAs. In tests, the full PTA runs retrieved the same posteriors from uninformative uniform priors, but we use the ACE prior throughout the work presented here to decrease convergence times in our Markov Chain Monte Carlo runs.

Uniform priors were used for the spectral index in all higher order searches and log-uniform priors were used for $n_{\mathrm{E}}^{(p)}$ priors when $p \ne2$.

Dispersion Measure Variation Model

Gaussian processes (GPs) have been used extensively to model time-correlated noise in pulsar datasets [10, 9, 25]. Commonly, they are used to model the GW background, achromatic red noise, such as that caused by intrinsic spin noise, and chromatic noise due to the bulk movement of the ISM [25, 48].

The standard prior function of GPs in these datasets is a red (negative spectral-indexed) power law model with a Fourier-basis. The solar wind has strong spectral characteristics at Fourier frequencies of $1/\mathrm{yr}$ and higher modes [51]. In early testing these modes showed themselves to be highly covariant with some frequencies in a Fourier-basis GP, so we have instead adopted for our DM variations model another common ansatz for GPs: a square exponential kernel,

where $\sigma$ is an overall variance, $\ell$ is a time scale of variation, and a constant variance is also added as a regularizing term for stable inversions and to ensure a minimal value of the variance. This kernel represents the realization average of covariance matrices described by this type of variation. They are similar to the structure functions used in the pulsar timing literature for characterizing noise due to variations in the ISM—for instance, see Equation (7.1) in Foster & Cordes [24] or Section 5 of Lam et al. [42]. Such structure functions were also used as a check in Tiburzi et al. [66] to ensure the iterative process used therein full disentangled the solar wind signal form the DM variations in the ISM. Our Bayesian framework allows for simultaneous fitting across multiple pulsars and frees us from the iterative process used in this latter publication. Tiburzi et al. [67] also uses a Bayesian analysis to separate these signals, but only within single pulsar datasets in order to keep track of solar wind dependencies on ecliptic latitude. The GP model above (along with our solar wind model) replaces the use of a short-time-span DMX model normally used in NANOGrav datasets [8, 3 b], since DMX cannot distinguish between the ISM and solar wind contributions to the measured DM value at each epoch. Here a standard 15 day linear interpolation basis in the time domain is used for constructing realizations of the DM variations.

Noise Analyses

Since we are implementing entirely new DM variation models from the original NG11 analysis, we redid single pulsar noise analyses using our solar wind and DM GP models. This allowed us to recover new values for the standard white noise parameters, EFAC, EQUAD and ECORR, (see e.g., [7, 46]) and power law achromatic red noise parameters for each pulsar. These single pulsar noise runs then consist of a set of three white noise parameters for each backend/receiver combination used for observations, two DM GP parameters, $\sigma$ and $\ell$, a solar electron density, $n_{\mathrm{E}}$, and two power law red noise parameters, $A_{\mathrm{RN}}$ and $\gamma$.

Using only one pulsar’s dataset it is impossible to fully disentangle the solar wind electron density from the ionized ISM. Any given observation’s DM variation can be attributed to either changes in the interplanetary or interstellar media. The key to these analyses is a global fit where the DM GP and red noise parameters are allowed to vary individually for each pulsar, while the solar electron density is set as a global parameter across all of the pulsars. Since pulsars have upwards of 20+ white noise parameters it is cumbersome to vary all of these parameters across the entire PTA, so we adopt the usual practice in PTA analyses [10] of setting their values constant in the full PTA analysis following single-pulsar noise modeling.

Bayesian Model Implementation and Analysis

We use the expressions in § 2.1 to build models for a spherically symmetric solar electron density and analyze NG11. We explore three different models for the behavior and structure of the solar electron density: higher order spherically symmetric solar wind models, time-binned solar wind density models, and continuous time-dependent solar wind perturbations. Each model demonstrates how our full Bayesian framework provides access to the unique probes of the solar wind available in PTA datasets.

Higher Order Models

Any number of higher order spherically symmetric models can be constructed using the tools developed here. We present only a few to demonstrate the feasibility of detecting these higher order signals using PTA datasets. In general we follow the methods of Aksim et al. [2] and fit for a varying power dependence, $\sim1/r^{p}$, of the solar electron density, though we still investigate models of the type in You et al. [77, 76] with multiple summed fixed-power terms. The fit for the proportion of the two-phases (slow and fast solar wind) in You et al. [76] is dependent on the individual lines-of-sight for the individual pulsars which is out of the scope of the current tests for these modeling tools.

In Table 1 we summarize results from a number of similar analyses. The aim here was not to exhaustively test the various powers from Equation (6) and Equation (7), as the real model is expected to be a mixture of the two [77], but rather to test the sensitivity of NG11 to some of these higher order terms. Table 1 shows a progression from the overly simplistic $1/r^2$ model, where we find $n_{\mathrm{E}}^{(2)} = 6.9^{+0.13}_{-0.13}\ \mathrm{cm}^{-3}$, to models with more terms, and more free parameters. Model 2, where the exponent in the density relation is allowed to vary, shows that there is support for a higher order model for the solar electron density, $p = 2.29^{+0.01}_{-0.02}$. This power is in agreement with the best fit power from Aksim et al. [2] of $p = 2.3$. Their value of $n_{\mathrm{E}}^{2.3} = 2.59 \pm0.13\ \mathrm{cm}^{-3}$ is smaller than our recovered value, but the Aksim et al. [2] value is from only one observation, while our value is an 11.4 year average. We will see in § 4.2 that time varying values can differ by up to a factor of 4 or more.

Model	First Term $p$	First Term $n_{\mathrm{E}}^{(p)}$ [cm$^{-3}$]	Second Term $p$	Second Term $\log_{10} n_{\mathrm{E}}^{(p)}$	Third Term $p$	Third Term $\log_{10} n_{\mathrm{E}}^{(p)}$
1	2	$6.9^{+0.13}_{-0.13}$	—	—	—	—
2	$2.29^{+0.01}_{-0.02}$	$4.7^{+0.20}_{-0.20}$	—	—	—	—
3	2	$2.5^{+0.56}_{-0.63}$	$2.4^{+0.04}_{-0.05}$	$0.38^{+0.10}_{-0.07}$	—	—
4	$2.30^{+0.02}_{-0.01}$	$4.5^{+0.18}_{-0.22}$	$4.3^{+0.58}_{-0.73}$	$-4.1^{95\%}$	—	—
5	2	$6.8^{+0.18}_{-0.19}$	4.39	$-2.87^{+0.03}_{-0.02}$	16.25	$-24.8^{95\%}$

Models 3 and 4 show that there is broad support for a second term in the density model when the power of the first term is kept constant, but the analysis only retrieves upper limits for a third term when the first two are kept constant. These analyses broadly support what has been previously shown in You et al. [76, 77] that pulsar timing data, especially for an array as we are treating here, is sensitive to a much more complicated model of the solar electron density than is routinely assumed. Model 5 shows the limits of the data set’s sensitivity, as we do not seem to detect a third term with the large index predicted in You et al. [77].

Time-Binned Solar Wind Density

As in Madison et al. [51] and Tiburzi et al. [67], we search for separate values of $n_{\mathrm{E}}$ with $p = 2$ during different periods of time by fitting values of $n_{\mathrm{E}}$ for discrete time bins of the dataset. In the Madison et al. [51] analysis the bins are set to be a year long and a fit is done to the DMX times series of the pulsars. DMX is a piece-wise time-binned analysis of the DM variations, (see, e.g. Arzoumanian et al. 2015). We present our results and those of Madison et al. [51] in Figure 3, along with binned ACE SWEPAM data for the same time period. The ACE data are in situ measurements of the local electron density, and so can be considerably noisier than the large effective average taken by the column density measurements to which the pulsar timing data are sensitive. In addition, since ACE orbits in the ecliptic plane, the SWEPAM instrument is really only sampling the slow solar wind, and would not track changes in the higher altitude fast wind. This is supported by scintillation studies that track high altitude solar wind density, see for instance Porowski et al. [57] for a model based on scintillation data. Their model for solar electron density shows strong positive correlations on the high latitude wind density with the solar cycle. Lastly, the ACE data are taken by the same instrument across the 11+ year time span shown here, while the NANOGrav PTA continually increased the number of pulsar in the array (hence LOSs) and sensitivity, including a “high cadence” observing campaign in the latter portion of the dataset Arzoumanian et al. [8].

The largest differences between our values and those of Madison et al. [51] occur in the first part of the dataset where the observation cadence was significantly less. The fit in Madison et al. [51] was done using a $\chi^2$-minimization, and did not include any priors for the solar electron density. Hence the unphysical negative value of the second bin is a result of the $\chi^2$-minimization. The large value in the fifth bin occurs during a time when there is a gap in the PTA data. This will be explored more in the finer binning used next. The Madison et al. [51] analysis is not as robust as the analysis here, since rather than fitting a secondary data product (DMX) we are fitting the solar wind as part of a full PTA-wide analysis, which includes a noise model fit and timing model marginalization.

In Figure 4 we show Bayesian binned results for bin sizes of three months again with $p = 2$, showing that NG11 has enough information to fit $n_E$ with finer resolution than is done in Madison et al. [51]. The 68% credible intervals for the binned values are small compared to the same intervals in the ACE data. In addition to the binned fit for $n_E$, we also run an analysis with a second density, $n_E^{(4.39)}$, in the solar wind density profile model with $p = 4.39$. We chose this index from the fast solar wind model of Equation (7) in order to investigate how a higher order term might effect the binned values of $n_E$ specifically from the few pulsars with very small $\theta_i$. We use a time-constant value of $n_E^{(4.39)}$ in this application for ease of analysis.

This second analysis, with an $n_E^{(4.39)}$ parameter highlights the importance of these higher order terms in the density model when there are observations with lines of sight close to the Sun. The geometric factor for the $\mathrm{DM}_{\odot}$ delay, $(\pi- \theta_i)/\sin\theta_i$, of all TOAs is shown in the bottom panel of Figure 4. There are two occasions in the 11-year dataset where $\theta_i$ is very small for PSR J0030+0451 and PSR J1614−2230. These can be seen as the points where the geometric factor is larger than $\sim100$. The two bins, marked by vertical dots, where these occur correspond to the two bins where the recovered values for $n_E$ differ by many standard deviations from the usual $1/r^2$ model for the density, over estimating $n_E$. While the $p = 4.39$ model seems to ameliorate the effect of these small $\theta_i$ observations, the direct cause of the anomalously high measurements of $n_E$ could also be due to the highly inhomogeneous streaming of the solar wind at these small distances or from individual events, i.e. coronal mass ejections. See Appendix § 9 for more discussion.

Two other bins warrant attention. The first bin has only a few TOAs with narrow frequency coverage, and hence parameters recovered have large error bars. The bin near MJD 54250 has no TOAs due to concurrent down time at both the Arecibo Observatory and the Green Bank Observatory. In this bin, the ACE prior is returned in the posterior. Hence no information beyond the prior is gained.

Continuous Time-Dependent Solar-Wind Perturbations

In addition to the piecewise constant models for $n_E$ discussed above, we also implemented a continuous Fourier-basis model for $n_E(t)$ as a perturbation to a mean $n_E$. The model is reminiscent of the free spectral models used by PTAs to describe achromatic red noise ([10]). It is constructed using the TOAs for each pulsar, parametrized with a separate value for each component of the basis,

Rather than trying to reproduce the fine structure recoverable by the You et al. [76] model, here we are concerned with long to mid scale variations of $n_E$, therefore we use a set of frequencies based on the time span, $T = 11.4$ yr, of NG11 and include 30 linearly spaced frequencies ranging from $[1/T, 30/T]$, using 2 parameters per frequency, $a_j$ and $b_j$.

In various test analyses we recorded the same type of behavior for these continuous models as was seen in Figure 4, i.e., the fit returned large values of $n_E$ during periods of time when J0030+0451 and J1614−2230 observations had small values of $\theta_i$. Here we only report on the continuous model for $n_E$ where we also fitted for the same static, higher order $n_E^{(4.39)}$ term discussed in § 4.2. Figure 5 shows the resulting recovery of $n_E(t) = n_E^{\mathrm{mean}} + n_E^{F}(t)$ from an analysis that uses the model in Equation (13). The use of the Fourier basis gives us immediate access to an analogue of the power spectral density for the perturbations to the solar electron density. Figure 6 shows violin plots for the power in the various frequencies calculated using the posteriors of the Fourier decomposition coefficients, $P = a_i^2 + b_i^2$. Notice that a majority of the frequencies show small amounts of power, while the lowest frequency, $1/(11.4\ \mathrm{yr})$, recovers the most power — very close to an inverse average solar cycle length of 11.75 yr. It will be interesting to see which frequencies have the most power in longer datasets.

A continuous signal like this would be a good candidate for a Gaussian process, however current code bases for PTAs construct different realizations of a process for each pulsar in an array. Since this work reveals a great deal of potential for gleaning information about the solar wind from PTAs we plan on developing per LOS solar wind perturbations for these pulsars in future work by using a simple deterministic inhomogeneous model.

Discussion

This work demonstrates a new set of Bayesian tools that allow for the use of PTA datasets, or any collection of pulsar timing data, as probes of the solar wind. The main innovations include a completely general Bayesian fit that includes timing model marginalization, generic spherically symmetric models for the solar electron density, and time dependent models. These models can be assembled in myriad ways in order to mitigate the solar wind as a noise source in pulsar data, or to use these sensitive datasets to study the Sun’s coronal behavior across the solar system. The NANOGrav data set is not aimed at monitoring the solar wind, nonetheless it is sensitive to interesting solar phenomena. The tools we have developed, if applied to a PTA data set built for the purposes of solar science, have great potential.

One important aspect of this work not yet discussed is how these models can enable more in depth generic ISM studies. These solar wind models allow the removal of the DM variation signal from the local interplanetary medium in a principled way. This would better enable the use of these pulsar timing datasets for studying the ISM unencumbered by local variations in electron density.

The solar wind signal has been known to cause issues [49, 60] when trying to measure some timing effects in pulsar data sets—for instance, timing parallax, critical for constraining pulsar distances. The covariances with measuring parallax imply that improvements in the solar wind model could ultimately improve single source GW searches, as the distance to pulsars is an important aspect of pulsar term searches for GW from supermassive black binaries [1]. The solar wind has also been implicated as one of the systematics in the “triple system” [6], impeding better limits of Einstein’s equivalence principle. Better modeling of the solar wind might allow for more accurate measurements of these limits, in addition to allowing for better noise mitigation when these data are used for GW searches. Improvements to the handling of the solar wind in pulsar timing investigations stand to improve precision pulsar science in a wide variety of ways.

PTA studies of the solar wind have the potential to add a completely independent set of measurements to the myriad space missions studying solar weather. PTA observations survey a much wider swath of the Sun’s environment than any spacecraft would be able to, taking column density information in 70+ lines of sight every month, probing beyond the outermost reaches of the solar environment. The observational campaigns of lower frequency observatories like the Canadian Hydrogen Intensity Mapping Experiment [15], LOFAR [61, 71] and the Long Wavelength Array (LWA, [62]) are continually adding data that is extremely useful for DM variability studies, since these lower frequencies allow for more accurate measurement of the variations. These data can add important tests to the well developed solar weather modeling efforts undertaken by the solar physics community [70]. These large scale measurements would allow for long timescale monitoring of the fast and slow solar wind, giving access to continuous measurements of the solar wind at all solar latitudes.

Acknowledgments

Author contributions. This paper is the result of the work of many people and uses data from over a decade of pulsar timing observations. We list specific contributions below. JSH developed and ran the solar wind analyses and led the paper writing. JS, and DRM contributed substantially to paper writing, discussion and interpretation of results. ZA, KC, MED, PBD, TD, JAE, RDF, ECF, EF, PAG, GJ, MLJ, MTL, LL, DRL, RSL, MAM, CN, DJN, TTP, SMR, PSR, RS, IHS, KS, JKS, and WZ ran observations and developed the 11-year data set.

Acknowledgments. The NANOGrav project receives support from National Science Foundation (NSF) as a Physics Frontier Center award number 1430284. We thank Caterina Tiburzi, Aditya Parsatharathy, Joseph Lazio and Patrick O’Neill for comments on an early draft of this manuscript. TD and MTL are supported by an NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. ECF is supported by NASA under award number 80GSFC17M0002. The Dunlap Institute is funded by an endowment established by the David Dunlap family and the University of Toronto. TTP acknowledges support from the MTA-ELTE Extragalactic Astrophysics Research Group, funded by the Hungarian Academy of Sciences (Magyar Tudományos Akadémia), which was used during the development of this research. SMR is a CIFAR Fellow. Portions of this work performed at NRL were supported by ONR 6.1 basic research funding. Pulsar research at UBC is supported by an NSERC Discovery Grant and by CIFAR. WWZ is supported by the CAS Pioneer Hundred Talents Program and the Strategic Priority Research Program of the Chinese Academy of Sciences, grant No. XDB23000000.

Facilities: AO, GBO

Software: Numpy ([28]), Scipy ([75]), Mat-PlotLib ([35]), Pint ([50]), Enterprise ([23]), enterprise_extensions ([63]), Tempo2 ([33]), PTMCMCSampler ([22])

Appendix

Generic Spherically Symmetric DM for Integer Powers

The expression in Equation (11) simplifies considerably for integer values of $p > 1$. Here we show the first few expressions for powers greater than $p = 2$ for reference.

Coronagraph Images

As discussed in § 4.2, three observation epochs from NG11 have fairly small solar impact angles. In Figure 7 we show images from the Large Angle Spectroscopic Coronagraph (LASCO, Brueckner et al. [14]) on board the Solar and Heliospheric Observatory (SOHO) from the same time as the pulsar observations, with the approximate pulsar positions marked. These white light images from the C2 camera can be related to the solar electron density (Quémérais & Lamy [58]), however such calculations are out of the scope of the current work. These images are provided

to highlight how close to the Sun these observations were taken ($< 6R_{\odot}$) and how inhomogeneous the inner solar system can be. Testing various estimates of 2 and 3-dimensional electron density from these white light images would be an interesting use of dedicated pulsar observations close to the Sun.