11institutetext:  ASTRON $-$ the Netherlands Institute for Radio Astronomy,
Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands 22institutetext:
Dipartimento di Fisica “G. Occhialini”, Università di Milano-Bicocca, Piazza
della Scienza 3, I-20126, Milano, Italy 33institutetext: Fakultät für Physik,
Universität Bielefeld, Postfach 100131, 33501 Bielefeld, Germany
44institutetext: Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69,
53121 Bonn, Germany 55institutetext: Department of Astrophysics/IMAPP, Radboud
University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands 66institutetext:
Technische Universität Berlin, Institut für Geodäsie und
Geoinformationstechnik, Fakultät VI, Sekr. H 12, Straße des 17. Juni 135,
10623 Berlin, Germany 77institutetext: GFZ German Research Centre for
Geosciences, Telegrafenberg, 14473 Potsdam, Germany 88institutetext: SKA
Organisation, Jodrell Bank, Macclesfield SK11 9FT, UK 99institutetext:
Gravitational Wave Data Centre, Swinburne University of Technology, P.O. Box
218, Hawthorn, VIC 3122, Australia 1010institutetext: Centre for Astrophysics
and Supercomputing, Swinburne University of Technology, P.O. Box 218,
Hawthorn, VIC 3122, Australia 1111institutetext: LPC2E - Université d’Orléans
/ CNRS, 45071 Orléans cedex 2, France 1212institutetext: Station de
Radioastronomie de Nançay, Observatoire de Paris, PSL Research University,
CNRS, Univ. Orléans, OSUC,18330 Nançay, France 1313institutetext: South
African Radio Astronomy Observatory, 2 Fir Street, Black River Park,
Observatory 7925, South Africa 1414institutetext: Department of Physics and
Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South
Africa 1515institutetext: Hamburger Sternwarte, University of Hamburg,
Gojenbergsweg 112, 21029 Hamburg, Germany 1616institutetext: Max-Planck-
Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85748 Garching b.
München, Germany 1717institutetext: Ruhr University Bochum, Faculty of Physics
and Astronomy, Astronomical Institute, 44780 Bochum, Germany
1818institutetext: Thüringer Landessternwarte, Sternwarte 5, 07778 Tautenburg,
Germany 1919institutetext: Jodrell Bank Centre for Astrophysics, University of
Manchester, M13 9PL, UK 2020institutetext: Leibniz-Institut für Astrophysik
Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany

# The impact of Solar wind variability on pulsar timing

C. Tiburzi<EMAIL_ADDRESS>G. M. Shaifullah 1122 C. G. Bassa 11 P. Zucca
11 J. P. W. Verbiest 3344 N. K. Porayko 44 E. van der Wateren 1155 R. A.
Fallows 11 R. A. Main 44 G. H. Janssen 1155 J. M. Anderson 6677 A-.S. Bak
Nielsen 4433 J. Y. Donner 4433 E. F. Keane J. Künsemöller 8833 S. Osłowski
991010 J-.M. Grießmeier 11111212 M. Serylak 13131414 M. Brüggen 1515 B. Ciardi
1616 R.-J. Dettmar 1717 M. Hoeft 1818 M. Kramer 441919 G. Mann 2020 C. Vocks
2020

(Received MM DD, YYYY; accepted MM DD, YYYY)

###### Abstract

Context. High-precision pulsar timing requires accurate corrections for
dispersive delays of radio waves, parametrized by the dispersion measure (DM),
particularly if these delays are variable in time. In a previous paper we
studied the Solar-wind (SW) models used in pulsar timing to mitigate the
excess of DM annually induced by the SW, and found these to be insufficient
for high-precision pulsar timing. Here we analyze additional pulsar datasets
to further investigate which aspects of the SW models currently used in pulsar
timing can be readily improved, and at what levels of timing precision SW
mitigation is possible.

Aims. Our goals are to verify: a) whether the data are better described by a
spherical model of the SW with a time-variable amplitude rather than a time-
invariant one as suggested in literature, b) whether a temporal trend of such
a model’s amplitudes can be detected.

Methods. We use the pulsar-timing technique on low-frequency pulsar
observations to estimate the DM and quantify how this value changes as the
Earth moves around the Sun. Specifically, we monitor the DM in weekly to
monthly observations of 14 pulsars taken with parts of the LOw-Frequency ARray
(LOFAR) across time spans of up to 6 years. We develop an informed algorithm
to separate the interstellar variations in DM from those caused by the SW and
demonstrate the functionality of this algorithm with extensive simulations.
Assuming a spherically symmetric model for the SW density, we derive the
amplitude of this model for each year of observations.

Results. We show that a spherical model with time-variable amplitude models
the observations better than a spherical model with constant amplitude, but
that both approaches leave significant SW induced delays uncorrected in a
number of pulsars in the sample. The amplitude of the spherical model is found
to be variable in time, as opposed to what has been previously suggested.

###### Key Words.:

pulsars:general, solar wind, ISM: general, gravitational waves

##  1 Introduction

High-precision pulsar timing (Lorimer & Kramer, 2004) is a technique used to,
for example, investigate irregularities in the Solar system planetary
ephemerides (e.g. Caballero et al. 2018; Vallisneri et al. 2020), generate
alternative time-scale references (e.g. Hobbs et al. 2020), test general
relativity (e.g. Archibald et al. 2018; Voisin et al. 2020) and alternative
theories of gravity (e.g. Shao et al. 2013) or search for low-frequency
gravitational waves with Pulsar Timing Arrays (PTAs, e.g. Tiburzi 2018; Burke-
Spolaor et al. 2019). In particular, a level of timing residuals below 100 ns
is usually indicated as the white-noise threshold to achieve in PTA
experiments (e.g., Janssen et al. 2015; Siemens et al. 2013).

The sensitivity of these experiments can be significantly degraded by various
noise processes such as those caused by errors in clock standards or
inaccuracies in the planetary ephemerides.

One of the most common sources of noise in pulsar-timing data (Lentati et al.,
2016) is the variable amount of free electrons along the line of sight (LoS).
The radio waves coming from pulsars are dispersed due to the ionized medium,
leading to time delays depending on the observing frequency, following the
relation:

$\Delta t=\frac{e^{2}}{2\pi
m_{e}c}\>\frac{DM}{f^{2}}=\mathcal{D}\>\frac{DM}{f^{2}}$ (1)

where $\Delta t$ is the time-delay (in s) induced at an observing frequency
$f$ (in MHz) with respect to infinite frequency, $e$ is the electron charge,
$m_{e}$ is the electron mass and $c$ is the speed of light ($e^{2}/2\pi
m_{e}c$ is the dispersion constant $\mathcal{D}=1/(2.4\times
10^{-4})\leavevmode\nobreak\ \mathrm{MHz^{2}}\leavevmode\nobreak\
\mathrm{pc}^{-1}\leavevmode\nobreak\ \mathrm{cm}^{3}\leavevmode\nobreak\
\mathrm{s}$, see Manchester & Taylor 1972) and DM is the dispersion measure
(in pc cm-3):

$DM=\int_{LoS}n_{e}(l)\>\mathrm{d}l$ (2)

with $n_{e}$ being the electron density along the LoS.

Variations in $n_{e}$ along the LoS induce DM fluctuations, causing changing
$\Delta t$ contributions to the arrival times of pulsar radiation that need to
be taken into account in pulsar timing experiments. The two main contributions
to the DM along a certain LoS (for a comprehensive review, see Lam et al.
2016) are the ionized interstellar medium (IISM) and the Solar wind (SW). In
particular, the SW contribution to DM depends on the Solar elongation of the
pulsar (i.e., the projected angular separation between the pulsar and the
Sun), whose temporal variations therefore induce DM time-fluctuations.

The SW has been recognized as a noise source that could induce false
detections of gravitational waves in PTA experiments (Tiburzi et al., 2016).

The standard pulsar timing approach to mitigate the SW contribution (e.g., the
International PTA data releases by Verbiest et al. 2016 and Perera et al.
2019a) typically consists in approximating the SW as a spherically symmetric
distribution of electrons $n_{e,sw}$ (Edwards et al., 2006):

$n_{e,sw}=\mathrm{A_{AU}}\left[\frac{1\mathrm{AU}}{r}\right]^{2}$ (3)

where $r$ is the distance between the pulsar and the Sun and $\mathrm{A_{AU}}$
is the free electron density of the Solar wind at 1 AU, reported to be $7.9$
and time-constant (Madison et al., 2019), and which we will henceforth refer
to as the amplitude of the SW density model. The DM contribution of this model
is obtained by integrating Equation 3 along the LoS, and can be expressed as
(Edwards et al., 2006; You et al., 2007b):

$DM_{sw}=4.85\times
10^{-6}\mathrm{A_{AU}}\frac{\rho}{\sin\rho}\leavevmode\nobreak\
\mathrm{pc\leavevmode\nobreak\ cm^{-3}}$ (4)

where $\rho$ is the pulsar-Sun-observer angle. This model implicitly assumes
that the amplitude $\mathrm{A_{AU}}$ is constant with time and independent of
the ecliptic latitude of the pulsar. To account for the fact that this model
may not be an optimal SW approximation at small Solar elongations, it is
common to eliminate data points taken at close ($<5$ degrees) angular
distances from the Sun in pulsar timing experiments (Verbiest et al., 2016).

However, the SW is more complex than implied from this simple model. Under
Solar minimum conditions it is mostly bi-modal, with a fast stream seen above
polar coronal holes, and a slow stream seen above a mostly-equatorial streamer
belt (e.g. Coles, 1996). A polar coronal hole can sometimes extend towards
equatorial latitudes, allowing the fast and slow streams to interact, leading
to denser regions of compression at the leading edge of the fast stream and
rarefied regions following behind (e.g Schwenn, 1990). Coronal Mass Ejections
(CMEs) will further complicate this picture. The picture becomes even more
complex as Solar activity increases towards a maximum, when the bi-modal
structure effectively breaks down allowing coronal streamers to extend to high
latitudes. The reader is referred to Schwenn (2006) for a more detailed
illustration of the SW system.

The shortcomings of the spherical model were neatly illustrated in a recent
publication by Tokumaru et al. (2020) who observed the Crab pulsar in 2018
using the Toyokawa Observatory. This paper also compared their results to SW
conditions assessed from observations of interplanetary scintillation and
coronal white light, and provided a useful discussion on how observations of
pulsars could also be used to assist research into the SW.

The aforementioned shortcomings led You et al. (2007a) to propose a revised SW
model for the pulsar DM which considered the SW as bimodal. The authors used
different free-electron radial distributions for each of the two SW streams,
and used Solar magnetograms to decompose the LoS into parts affected by one or
the other component, the total contribution from the Solar-wind being the sum
of these individual contributions. This was demonstrated in the paper to
better correct the DM for the SW contribution than the basic spherical model,
but it should be noted that the pulsar observations were taken during the
approach to solar minimum, when a bimodal solar wind structure is more
evident.

Tiburzi et al. (2019) compared the performance of these models on highly-
sensitive, low-frequency observations of PSR J0034$-$0534, while also allowing
a time-variable amplitude (following the approach of You et al. 2012) in both
of the models. The authors demonstrated that neither model provided an
adequate description of the SW impact on the dataset, but also that the
spherical one performed better than the other. The observations used in that
paper were, however, taken at solar maximum, which may explain why the bimodal
model did not perform better in that instance. More explanation of the
possible reasons is given in Tiburzi et al. (2019).

In this article we expand on the analysis of Tiburzi et al. (2019), using a
larger sample of pulsars to verify a) whether the SW DM contributions to these
pulsar LoSs are better described by a spherical model with a time-variable
amplitude or a time-constant one as suggested in the literature (Edwards et
al., 2006; Madison et al., 2019), and b) whether a consistent temporal trend
can be detected in the amplitudes of this model which might suggest that the
bimodal approach be revisited in future work.

The article is structured as follows: in Section 2 we present the dataset and
the selected pulsars, while in Section 3 we describe the analysis. The results
are presented in Section 4. In Section 5 we discuss future prospects, and in
Section 6 we draw our conclusions.

##  2 Dataset

The utilized dataset comes from a number of pulsar monitoring campaigns
carried out with the high-band antennas of different subsets of the
International LOFAR (LOw Frequency ARray) telescope (van Haarlem et al., 2013;
Stappers et al., 2011): the six German International LOFAR stations, the
Swedish International LOFAR station, and the LOFAR Core. The observing
bandwidth covers a frequency interval from ${\sim}100$ to ${\sim}190$ MHz,
with a central frequency of about $150$ MHz (variations of a few MHz occur
among the different observing sites). The recorded data were coherently
dedispersed, folded into 10-second-long subintegrations modulo the pulse
period, and divided into frequency channels of 195 kHz with the DSPSR software
suite (van Straten & Bailes, 2011). The integration length ranges from $1$ to
$3$ hours with the international stations, and from $7$ to $20$ minutes with
the LOFAR core (for more details regarding the observational setup, see
Porayko et al. 2019; Donner et al. 2019 and Tiburzi et al. 2019).

Together, the aforementioned observing campaigns monitor more than $100$
pulsars. However, for the scope of this article, we selected pulsars with the
following characteristics: a) ecliptic latitude between $-20$ and $+20$
degrees, b) observing cadence higher than once per month, c) more than one
year of observing time-span, d) without gaps between successive observations
exceeding $100$ consecutive days. This results in a dataset of $43$ pulsars,
whose sky locations are shown in Figure 1 and characteristics are reported in
tables 1 and 2111Part of the reported values come from the ATNF pulsar
catalog, https://www.atnf.csiro.au/research/pulsar/psrcat/ (Manchester et al.,
2005)., and for which we have used all the data available until May $2019$ (up
to August $2019$ in some cases). This source list was further refined to 14
sources that prove useful probes of the Solar wind, as discussed in Section
3.3.

Figure 1: Sky distribution in Galactic coordinates of selected (filled green
circles) and rejected (unfilled red circles) pulsars from the lists in Table 1
and Table 2. Larger markers denote higher (median) precision of the measured
DM. The dashed blue line marks the ecliptic and the dashed gray lines show a
region of $\pm$20 degrees from the ecliptic. The background shows the merged
all-sky H$\alpha$ map from Finkbeiner (2003).

##  3 Data analysis

In the following, we detail the analysis methods applied to the data to obtain
the DM values, disentangle the IISM-induced effects, and identify the pulsars
showing persistent SW signatures. We make use of the pulsar timing technique,
which keeps track of every pulsar rotation to model the evolution of the pulse
period and phase over the timespan of our dataset. This allows us to average
each observation in time and improve the signal-to-noise ratio (S/N).

###  3.1 Calculation of the DM values

For all observations, independent of the observing site, we removed radio-
frequency interference, corrected for azimuth and elevation-dependent gain
using the LOFAR beam model within the dreambeam
package222https://github.com/2baOrNot2ba/dreamBeam, and applied band-
limitations to retain a common frequency range between ${\sim}118$ and
${\sim}188$ MHz333Band-limiting is necessary to avoid biases in our results
because the observing sites record slightly different original bandwidths. The
indicated frequency range is common to all of the observing sites (for more
details, see Donner et al. 2019; Tiburzi et al. 2019).. These operations were
carried out using the psrchive software suite (Hotan et al., 2004; van Straten
et al., 2012) and a modified
version444https://github.com/larskuenkel/iterative_cleaner, see also Kuenkel
(2017). of the Coastguard software suite (Lazarus et al., 2016).

For each pulsar we then computed a DM value per observation through the pulsar
timing technique by proceeding as follows. We first selected the dataset that
covered the longest time-baseline. These observations were weighted by the
square of their S/N, added together and then fully averaged in time and
partially in frequency (usually down to $10$ frequency channels to increase
the S/N555For particularly faint pulsars we applied a larger frequency-
averaging factor.). Finally, the template was smoothed using a wavelet
smoothing scheme first introduced by Demorest et al. (2013). We then collected
the observations obtained by all the observing sites for that pulsar, we fully
averaged them in time and partially in frequency to the same resolution of the
template. For each observation we then generated a set of ToAs associated to
the frequency channels with the psrchive software suite by cross-correlating
the observation with the reference template.

A few exceptions to the aforementioned general template-generation scheme were
adopted in case the pulsar was too faint, or bright but strongly affected by
red noise (due to irregularities in the pulsar rotation or extreme IISM-linked
DM variations, that caused the template to appear broadened). In the first
case, we averaged the longest dataset described earlier over frequency and
time to then produce an analytic template, obtained by approximating the data-
derived pulse profile with a sum of von Mises functions. In the second case,
we used a small subset of phase-aligned observations which were subsequently
averaged and smoothed.

After generating a set of frequency-resolved ToAs per observation, we used the
tempo2 software suite for pulsar timing (Edwards et al., 2006) and the
procedure outlined in Tiburzi et al. (2019) to calculate one DM value per
observation666Note that, as for Tiburzi et al. (2019), the DM time derivatives
included in the original timing model for that pulsar were only applied to
properly dedisperse the template and frequency-average the observations, but
they were not used in the subsequent determination of the DM variations.. We
then combined the DM time series from all the available observing sites, after
subtracting the reference DM value.

###  3.2 Disentangling the IISM contribution

The obtained DM time series show variations due to both the SW and the IISM.
To remove the influence of the IISM, for each pulsar we proceeded as follows.

The DM time series was divided into 460-day long segments centered on the
Solar conjunctions, i.e., adding an additional 1.5 months of baseline to a
6-month time window on either side of the Solar conjunction. Hence, the
segments overlap for about 100 days. Segments that either contained gaps of
more than 55 days between successive observations, or those where the
effective time span is less than 368 days (80% of 460 days) were also
discarded, as they do not provide a long enough baseline to properly define
the IISM effects (this only affects initial and final segments).

The DM time series in each segment was modeled in a Bayesian framework as the
sum of a spherically-symmetric SW model (Equation 4) with $\rm{A_{AU}}$ being
a free parameter, and a polynomial to account for the IISM contribution. By
simulating DM time-series affected by Kolmogorov turbulence (Armstrong et al.,
1995), we found that a cubic polynomial was sufficient to model the DM
variations due to the IISM on 460-day long segments (see the method’s
validation in Appendix A). To account for any common systematic error in the
estimation of the DM uncertainties, we inserted an additional parameter in the
model, summed in quadrature with the DM uncertainties in the likelihood
function. We used a Markov Chain Monte Carlo method (implemented using the
emcee package; Foreman-Mackey et al., 2013) to obtain the parameters of the
cubic polynomial, the uncertainties correction parameter and the SW amplitude
$A_{\mathrm{AU}}$, and account for the covariance of the SW amplitude with the
parameters of the IISM model. We assigned flat priors to the polynomial
coefficients, and a flat and positive prior to the amplitude of the SW model
and the uncertainties correction parameter.

As a final step, the modeled IISM contribution was subtracted from the DM time
series. In the overlapping region between two successive years, data points
that lay within 6 months of the preceding Solar approach were approximated by
the plasma model computed for the first year, while data points within 6
months of the following Solar approach were approximated by the one computed
for the second year777The overlap between adjacent years guarantees a
continuous and smooth IISM model across multiple years.. After estimating the
model, data points collected earlier than 6 months before the first Solar
conjunction, and later than 6 months after the last Solar conjunction were
discarded.

As an example of the final result, Figure 2 shows the IISM disentanglement for
PSR J0030+0451.

Figure 2: IISM-SW disentanglement in PSR J0030+0451. The upper panel shows the
time series of the DM variations and the IISM models (dot-dashed green line)
and of the SW and IISM combined (solid blue line). The used data are
represented as black dots and red diamonds (where the red diamonds indicate
the observations within 50 degrees in Solar elongation and the black dots the
observations beyond 50 degrees), while the discarded data points are in gray.
The middle panel shows the DM residuals after subtraction of the IISM model
only, and the lowest panel shows the DM residuals after subtraction of both
the IISM and the SW models.

We note here that in Tiburzi et al. (2019) the authors demonstrated that the
spherical model is a poor description of the SW, when tested against
sufficiently sensitive data. Nevertheless, we have adopted it in the procedure
described above for two reasons. First, it was proven to be better among the
two available models compared in that work. Secondly, in Tiburzi et al. (2019)
we found that, while the spherical approximation could not model the short-
term SW-induced DM variations, it provided a reasonable description of the
long-term ones.

###  3.3 Final selection

The procedure outlined above allowed us to refine our pulsar selection by
rejecting those sources where the SW signature is not reliably detected.
Specifically, we retained a pulsar if and only if in more than half of the
dataset: a) the model described in the previous section was preferred over an
IISM-only model as evaluated by the Bayesian information criterion (BIC), b)
the posterior distribution of the SW model’s amplitude was significantly
different from zero888In PSR J1300$+$1240 only the first two years (out of a
total five) meet the requirements. However, we included it in the final
selection after visually inspecting the DM time series and manually examining
the results..

Figure 3: DM time series of the selected pulsars, after subtraction of the
IISM component. Observations highlighted in red were taken when the Solar
elongation of the source was less than $50$ degrees, while the gray vertical
lines mark the MJDs of the Solar conjunctions. For the sake of visual clarity,
these plots show only 95% of the most precise measurements although all data
points are used in the analysis. Figure 4: Pulsars included in the final
selection (red circles) versus the excluded pulsars (black stars) as a
function of ecliptic latitude and the logarithm of the median DM error.

A total of $14$ pulsars satisfied the mentioned requirements, as reported in
Table 1 (for the discarded sources, see Table 2 in the Appendix). Among these
there are the PTA-class millisecond pulsars J0030+0451, J0034$-$0534,
J1022+1001, J2145$-$0750 and J2317+1439 (Perera et al., 2019b), which, as
expected, display the best DM precision of the sample. We consider the 14
selected sources as well-suited for studying the electron density in the SW at
low frequencies in the Northern hemisphere. Figure 4 shows that, for equal
ecliptic latitude, the pulsars included in the final selection always present
the best DM precision among the sources at the same Ecliptic latitudes. The
few cases of sources with high DM precision where the SW is not detected, can
be explained typically by a combination of low S/N and a poor sampling cadence
(e.g., PSR J1024$-$0719). For the few pulsars in which these causes are not
applicable (e.g., PSR J0837+0610), we speculate that the reason lies in an
asymmetry of the Solar wind contribution with respect to the heliographic
latitude. A more rigorous test of this hypothesis will be presented in a
future work through comparisons with data-derived magneto-hydrodynamic
simulations of the SW electron density fluctuations, such as EUHFORIA (Poedts
& Pomoell, 2017).

Figure 3 shows the SW component of the DM variations of the final selection.

Table 1: Final source selection, encompassing $14$ pulsars. The table reports the source name, the covered time-span, the number of observing sites that have been monitoring that specific source, the Galactic coordinates, the rotational spin period, ecliptic latitude, dispersion measure (DM) of the pulsar, as measured during the general pulsar timing analysis described in Section 3, the decimal logarithm of the median DM uncertainty and the number of observations used to generate the data-derived template. Name | Time-span | Observing | Galactic | Period | Ecliptic | DM | Log | Observations
---|---|---|---|---|---|---|---|---
|  | Sites | Coo. [deg] | [ms] | Latitude [deg] |  [pc/cm3] | M(eDM) | per template
J0030+0451 | 2013-01 2019-05 | 7 | $113.1$ | $-57.6$ | $4.9$ | $1.45$ | $4.3$ | $-4.02$ | $215$
J0034$-$0534 | 2012-12 2019-05 | 6 | $111.5$ | $-68.1$ | $1.9$ | $-8.53$ | $13.8$ | $-4.66$ | $227$
J0051+0423 | 2013-08 2019-05 | 5 | $123.0$ | $-58.5$ | $354.7$ | $-1.05$ | $13.9$ | $-3.17$ | $208$
J0304+1932 | 2013-08 2019-05 | 6 | $161.1$ | $-33.3$ | $1387.6$ | $2.10$ | $15.7$ | $-3.0$ | $208$
J0407+1607 | 2013-08 2019-05 | 7 | $176.6$ | $-25.7$ | $25.7$ | $-4.74$ | $35.6$ | $-3.37$ | $216$
J0826+2637 | 2013-08 2019-05 | 7 | $197.0$ | $31.7$ | $530.7$ | $7.24$ | $19.5$ | $-3.79$ | $150$
J1022+1001 | 2012-12 2019-08 | 8 | $231.8$ | $51.1$ | $16.5$ | $-0.06$ | $10.3$ | $-3.95$ | $290$
J1136+1551 | 2013-08 2019-05 | 6 | $241.9$ | $69.2$ | $1187.9$ | $12.16$ | $4.8$ | $-3.97$ | $299$
J1300+1240 | 2012-12 2019-08 | 6 | $311.3$ | $75.4$ | $6.2$ | $17.58$ | $10.2$ | $-4.29$ | $239$
J1400$-$1431 | 2015-10 2019-05 | 7 | $327.0$ | $45.1$ | $3.1$ | $-2.11$ | $4.9$ | $-4.33$ | $167$
J1607$-$0032 | 2013-09 2019-05 | 6 | $10.7$ | $35.5$ | $421.8$ | $19.99$ | $10.7$ | $-3.61$ | $26$
J1645$-$0317 | 2013-08 2019-08 | 7 | $14.1$ | $26.1$ | $387.7$ | $18.86$ | $35.8$ | $-4.18$ | $4$
J2145$-$0750 | 2013-01 2019-08 | 7 | $47.8$ | $-42.1$ | $16.1$ | $5.31$ | $9.0$ | $-4.11$ | $232$
J2317+1439 | 2012-12 2019-08 | 8 | $91.4$ | $-42.4$ | $3.4$ | $17.68$ | $21.9$ | $-4.28$ | $7$

##  4 Results

###  4.1 Performance with respect to a SW model with constant amplitude

From the point of view of pulsar timing experiments, it is important to
understand whether a spherical model of the SW performs better (i.e., yields
smaller residuals when subtracted from the observations) when a time-variable
or a static amplitude is assumed. For this aim, we repeated the analysis
described in Section 3.2 on the final pulsar selection by fixing the amplitude
for the spherical SW model to a value of $7.9\leavevmode\nobreak\
\mathrm{cm}^{-3}$ (Madison et al., 2019) and we compared their results. We
stress that, for this analysis, we exclude those segments in the pulsar’s
datasets where the posterior distribution of the SW amplitude was found to be
consistent with zero in the previous Section. Assuming an absence of
frequency-dependence of the DM (cf. Cordes et al. 2016), in Figure 5 we show
the comparison between the two spherical models, with a time-dependent and a
time-invariable amplitude, reported through Equation 1 in terms of residual
time delays at 1400 MHz (the main reference frequency for high-precision
pulsar-timing studies). In particular, we display the rms of the time delays
induced by the residual DM fluctuations in the two different analyses, binned
in Solar elongation. In Figure 5, the black dots and red stars refer
respectively to a constant- and a variable-amplitude SW model.

By assuming the rms of the residual time delays as criterion, a model with a
variable amplitude performs better than one with a constant amplitude in more
than 60% of the cases for Solar elongations up to $20$ degrees. The PTA-class
pulsars of the sample show the most significant improvements depending on the
Solar elongation. For example, the rms of the residual time delays decreases
of a few tens of sigma999Measured as
$(\mathrm{rms}_{\mathrm{C},i}-\mathrm{rms}_{\mathrm{V},i})/\mathrm{erms}_{\mathrm{V},i}$,
with $\mathrm{rms}_{\mathrm{V},i}$ and $\mathrm{rms}_{\mathrm{C},i}$ being the
rms of the residuals left by, respectively, the spherical model with time-
dependent amplitude and the spherical model with time-invariable amplitude at
Solar elongation $i$, and $\mathrm{erms}_{\mathrm{V},i}$ being the uncertainty
to $\mathrm{rms}_{\mathrm{V},i}$. in PSR J0034$-$0534 below 15 degrees in
Solar elongation, and from ${\sim}15$ to ${\sim}10$ sigma in PSRs
J2145$-$0750, J1022+1001, and J0030+0451 at the closest Solar elongations.

Nevertheless, Figure 5 also clearly shows that a simple spherical model, even
when a variable amplitude is applied, does not provide a sufficiently precise
description of the SW contribution to the DM. This is evident from the level
of the residual time-delay’s rms at the smallest elongations (lower than 10 to
20 degrees), which, in the case of the PTA-class pulsars (yielding the best DM
precision), never reaches the noise floor set by the points at the largest
Solar elongations. At the smallest elongations, SW acceleration leads to a
steeper decrease in density than implied by a simple inverse-square. This is
illustrated by Bird et al. (1994) who found that a density model of
$r^{-2.54}$ better fitted data inside of 10∘ elongation. The bimodal model
proposed by You et al. (2007a) used separate density models to account the
fast and slow streams, based on results published by Guhathakurta & Fisher
(1995), Guhathakurta & Fisher (1998) and Muhleman & Anderson (1981), Allen
(1947) respectively. The results presented here further demonstrate the
necessity to account for SW acceleration in measurements taken close to the
Sun.

This confirms the findings of Tiburzi et al. (2019), and extends them to a
larger number of pulsars at different ecliptic latitudes.

Figure 5: Root-mean-square of the residual timing delays in our data, rescaled
to an observing frequency of 1400 MHz. The timing residuals were binned in
solar elongation with a resolution of 5 degrees and a Solar-wind model with
time-constant (black) or variable (red) amplitude was subtracted. The pulsars
are sorted by ecliptic latitude.

###  4.2 Temporal and latitudinal evolution of the SW density

Our results show that the amplitude of the spherical SW model is not constant
with time and ecliptic latitude. This is displayed in the left panel of Figure
6, which reports the temporal evolution of the aforementioned amplitudes, each
obtained from the data of a specific Solar approach for each pulsar. By
averaging the computed SW amplitudes (excluding the upper limits) obtained
from different pulsars in a given year and latitude range (right panel of
Figure 6), a decreasing trend in time becomes evident at ecliptic latitudes
between 5 and 20 degrees, and from $-$5 to $-$10 degrees. The identified
decreasing trend may be related to the cyclical Solar activity, which peaked
in 2014 and is expected to have reached a minimum in 2020. Conversely, the
amplitude of the model does not show an evident decrease in time for pulsars
with an ecliptic latitude included between $-5$ and $+5$ degrees. This may be
due to the smaller distance to the Ecliptic with respect to other pulsars, but
additional studies are needed to identify the cause of these temporal
trends101010Because the spherical model is known to be an imperfect SW
approximation (Tiburzi et al., 2019), we report in Appendix C an analogous
overview of the DM variations..

|
---|---
Figure 6: Temporal trend of the amplitude for the SW spherical model, divided
into bins of ecliptic latitude, for the individual pulsars (left panel) and
averaged year by year (right panel). The gray line marks an amplitude of
$7.9\leavevmode\nobreak\ \mathrm{cm}^{-3}$ (the best-fit value from Madison et
al. 2019). All data points are shown with vertical error bars, while the upper
limits are displayed at the levels of the 95th percentile of the posterior
distribution in the left panel, and not considered in the right panel.

##  5 Future prospects for study of SW-induced DM signatures

Here we comment on possible future studies of this kind, specifically we
identify two possible sources for future studies with LOFAR and comment on the
potential for future radio observatories.

###  5.1 Prospects for LOFAR Studies

Some of the pulsar campaigns whose data we have used in this article are
carried out with the LOFAR core at a monthly cadence. As stated in Section 2,
if a pulsar is observed at a monthly cadence only (i.e., without any coverage
from the international stations, that usually observe at a weekly cadence) it
has not been taken into consideration for this project. However, it is worth
mentioning two particularly promising pulsars for possible SW studies, because
of their low ecliptic latitude and good DM precision: PSRs J1730$-$2304 (with
a median DM precision of $3\times 10^{-4}\leavevmode\nobreak\
\mathrm{pc/cm^{3}}$ and an ecliptic latitude of 0.19 degrees) and J2256$-$1024
(with a median DM precision of $6\times 10^{-5}\leavevmode\nobreak\
\mathrm{pc/cm^{3}}$ and an ecliptic latitude of $-3.41$ degrees). New
observing campaigns with the LOFAR core have been recently carried out with
the aims of increasing the observing coverage of these and other pulsars
during their Solar approaches, and of collecting simultaneous data for
interplanetary scintillation studies to be compared with the pulsar-based
results and white-light observations, as detailed below.

Moreover, as the decreasing trend shown in many pulsars might be related to
the Solar cycle, future LOFAR observations may be able to detect an increase
in the excess DM induced by the SW as the Solar activity is expected to
increase in the upcoming years.

Further research is also planned to investigate more robust methods to
estimate the SW contribution to pulsar DM. This includes more detailed
investigation of the bimodal approach of You et al. (2007a), incorporating a
direct comparison of the line-of-sight with Carrington maps of coronal white
light to ascertain regions of fast and slow SW, as in the approach taken
recently by Tokumaru et al. (2020). Other approaches under investigation
involve the use of 3-D tomographic reconstructions of velocity and density in
the inner heliosphere, obtained from observations of interplanetary
scintillation and coronal white light (e.g Jackson et al., 2020, and
references therein), and the use of space weather models such as EUHFORIA
(Poedts & Pomoell, 2017).

###  5.2 SW studies with next-generation telescopes

While we have demonstrated the capabilities of LOFAR in monitoring the SW with
pulsars, its activity could be complemented by the utilization of other
telescopes.

Observing facilities in the Southern hemisphere with the capability of
covering low-frequency ranges, such as Murchinson Widefield Array (Tingay et
al., 2013) and the upcoming SKA1-Low (Braun et al., 2015) will have access to
a different and more extended set of pulsars in comparison to LOFAR, because
of a more prolonged visibility of the inner Galaxy. Furthermore, with
telescopes such as CHIME (Bandura et al., 2014), MeerKAT (Jonas, 2009),
SKA1-Low, LOFAR2.0 (the upcoming upgrade of LOFAR) and
NenuFAR111111https://nenufar.obs-nancay.fr/ (Bondonneau et al., 2020) there
will be access to a wider frequency bandwidth, with a consequential enabling
of an increased range of scientific studies. As some of these facilities have
the possibility of multi-beaming and sub-arraying, they will also be able to
track more pulsars simultaneously, hence reducing the times of the observing
campaigns.

##  6 Conclusions

We presented a study of the impact of the SW on a large sample of pulsars
observed for up to 6 years with the LOFAR telescope. This study demonstrated
that the spherically symmetric and static SW model, that is commonly used in
pulsar-timing experiments, can be improved by allowing a time-variable
amplitude, but that these improvements do not suffice for correcting the SW
impact at the levels required for high-precision pulsar-timing experiments
(consistent with our earlier findings in Tiburzi et al., 2019). For PTA-class
pulsars, in particular, the residuals between observations taken at small
Solar elongations (lower than 10 to 20 degrees) and the corresponding model
are not able to reach the noise-floor set by residuals corresponding to larger
Solar elongations. As Tiburzi et al. (2016) demonstrated that the SW is a
potential source of false GW detections, it is then advisable to adopt time-
variable SW amplitude and treat carefully observations taken during Solar
conjunction within high-precision pulsar timing experiments.

Moreover, the data indicated that the amplitude of the spherical model has a
dependency on time for pulsars whose ecliptic latitude lies either between
$-10$ and $-5$ degrees or $+5$ and $+20$ degrees. In particular, the amplitude
decreases from the first half of the observed time-span to the second half of
that time-span. On the other hand, the amplitude tends to remain constant for
pulsars with an ecliptic latitude included between $-5$ and $+5$ degrees.
Additional studies will be necessary to recognize the cause of these trends,
and whether they are in connection to the Solar cycle activity.

While the investigation of new SW models for high-precision pulsar timing was
not in the scope of the current work, a series of measures can be adopted to
mitigate the SW-induced noise (assuming the absence of frequency-dependent DM
or scattering):

  * •

Carrying out observations with wideband observing receivers may allow, in
combination with the usage of frequency-resolved templates, a precise
determination of time-dependent DM (due to the IISM, the SW or both) that can
be then used to correct the noise induced by variable dispersion;

  * •

Carrying out simultaneous observations at low and high frequency to allow a
precise DM determination by using the low-frequency data. This can be then
used to correct the high-frequency ToAs;

  * •

If obtaining simultaneous observations is not possible, it may be of aid
carrying out high-cadence, low-frequency observations (ideally once every two
or three days) for two weeks around the Solar conjunction. The DM values
obtained from these observations can be interpolated and used as correction
scheme for high-frequency observations. The high-cadence of the observations
is important in this case because, as shown by Niu et al. (2017), the fast SW
variability may invalidate the DM corrections if the high- and low-frequency
observations are separated by more than one day;

  * •

For pulsars with a flat spectral index, it may be meaningful to carry out
observations at very high frequencies (at S band or more), where the DM-
induced noise is marginal, and apply the spherical SW model by choosing an
amplitude accordingly to, e.g., the results of this article depending on the
pulsar’s ecliptic latitude. However, this approach implies that a certain
amount of timing noise will be left in the data, especially the ones taken
close to the Solar conjunction. Moreover, Lam et al. (2018) demonstrated that
observing frequencies lower than 1 GHz are more optimal to achieve high
precision in pulsar timing experiments.

## Data Availability

The data underlying this article that were collected with the International
LOFAR Stations and with the LOFAR core under still private observing programs,
will be shared on reasonable request to the corresponding author. The data
underlying this article collected with the LOFAR core and under public
observing programs are available at: https://lta.lofar.eu/.
The initial timing models, the ToA files and the DM time series for the
pulsars used in this article are publicly available on Zenodo from the 1st of
January 2021 (DOI: 10.5281/zenodo.4247554).

###### Acknowledgements.

This work is part of the research program Soltrack with project number
016.Veni.192.086, which is partly financed by the Dutch Research Council
(NWO). The authors thank the anonymous referee for their support and useful
comments. CT This paper is partially based on data obtained with: i) the
German stations of the International LOFAR Telescope (ILT), constructed by
ASTRON (van Haarlem et al., 2013) and operated by the German LOng Wavelength
(GLOW) consortium (https://www.glowconsortium.de/) during station-owners time
and proposals LC0_014, LC1_048, LC2_011, LC3_029, LC4_025, LT5_001, LC9_039,
LT10_014; ii) the LOFAR core, during proposals LC0_011, DDT0003, LC1_027,
LC1_042, LC2_010, LT3_001, LC4_004, LT5_003, LC9_041, LT10_004, LPR12_010;
iii) the Swedish station of the ILT during observing proposals carried out
from May 2015 to January 2018. We made use of data from the Effelsberg (DE601)
LOFAR station funded by the Max-Planck-Gesellschaft; the Unterweilenbach
(DE602) LOFAR station funded by the Max-Planck-Institut für Astrophysik,
Garching; the Tautenburg (DE603) LOFAR station funded by the State of
Thuringia, supported by the European Union (EFRE) and the Federal Ministry of
Education and Research (BMBF) Verbundforschung project D-LOFAR I (grant
05A08ST1); the Potsdam (DE604) LOFAR station funded by the Leibniz-Institut
für Astrophysik, Potsdam; the Jülich (DE605) LOFAR station supported by the
BMBF Verbundforschung project D-LOFAR I (grant 05A08LJ1); and the Norderstedt
(DE609) LOFAR station funded by the BMBF Verbundforschung project D-LOFAR II
(grant 05A11LJ1). The observations of the German LOFAR stations were carried
out in the stand-alone GLOW mode, which is technically operated and supported
by the Max-Planck-Institut für Radioastronomie, the Forschungszentrum Jülich
and Bielefeld University. We acknowledge support and operation of the GLOW
network, computing and storage facilities by the FZ-Jülich, the MPIfR and
Bielefeld University and financial support from BMBF D-LOFAR III (grant
05A14PBA) and D-LOFAR IV (grants 05A17PBA and 05A17PC1), and by the states of
Nordrhein-Westfalia and Hamburg. MB acknowledges support from the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy - EXC 2121 ”Quantum
Universe” - 390833306. CT acknowledges support from Onsala Space Observatory
for the provisioning of its facilities/observational support. The Onsala Space
Observatory national research infrastructure is funded through Swedish
Research Council grant No 2017-00648. GS was supported by the Netherlands
Organization for Scientific Research NWO(TOP2.614.001.602). JPWV acknowledges
support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg
program (Project No. 433075039).

##  Appendix A Effectiveness of the IISM-SW disentangling scheme

To confirm that our choice of a third order polynomial is a sufficient
description of the IISM contributions, we have simulated 500 DM time series
affected by Kolmogorov turbulence and white noise as drawn from a zero mean
Gaussian population with a variance of $1\times 10^{-4}\leavevmode\nobreak\ $
pc/cm3, and we have attributed irregular error bars to match the ones of PSR
J1022+1001. To this time-series, we then added the SW contribution from a
spherically symmetric model, with a variable amplitude per year. This final
dataset was then analyzed following the steps detailed in Section 3.2. For
each year and injected value of the spherical SW amplitude, as shown in Figure
7, we were able to recover a statistically comparable value of the amplitude
of spherical SW model.

We further tested whether the third order polynomial model is able to account
for the IISM contribution only. To this end, we repeated the above simulations
without adding the SW contribution, and we analyzed the DM time series using
the MCMC algorithm described in Section 3.2, modified to apply a segment-by-
segment cubic-only model. To determine the success rate of modeling procedure
we used a two-tailed Kolmogorov-Smirnov (KS) test. The results report that in
the 100% of the cases, the distribution of the residuals after applying the
cubic model are identical to a zero-mean Gaussian distribution. This confirms
the ability of our algorithm in removing the IISM contribution from the DM
time series.

To test whether our algorithm is as performing on real, irregular DM time
series as it is on simulated data, we repeated the just-described analysis on
the DM time series of 4 millisecond pulsars presented in Donner et al. (2020),
namely PSRs J0218+4232, J0740+6620, J1125+7819, J1640+2224. These sources have
a high ecliptic latitude ($>30$ degrees); therefore, we can assume that they
are only affected by IISM-induced variations. We used the IISM-only MCMC
algorithm that we described in the previous paragraph (i.e., modified to apply
a segment-by-segment cubic-only model) to model the IISM variations presented
by these pulsars. The two-tailed KS test previously described showed that,
after modeling the DM time series with the results of our algorithm, the
residuals resulted as well compatible with a zero-mean Gaussian distribution
(see Figure 8).

Thus we are able to demonstrate that our modeling method leads to the IISM
contribution being successfully disentangled from the SW component.

Figure 7: Injected (dashed, orange lines) and recovered (filled, black histograms) values of the amplitude of the spherically symmetric SW model. |
---|---
|
Figure 8: DM time series of the four tested millisecond pulsars (upper panel
of each quadrant) and their post-modeling residuals (bottom panels).

##  Appendix B Pulsars without significant SW signatures

Table 2 reports the investigated pulsars that do not show significant SW
signatures.

Table 2: Investigated pulsars that do not show significant SW signatures. The table reports the source name, the covered time-span, the number of observing sites that have been monitoring that specific source, the Galactic coordinates, the rotational spin period, ecliptic latitude, dispersion measure (DM) of the pulsar as measured during the general pulsar timing analysis described in Section 3, the decimal logarithm of the median DM uncertainty and the number of observations used to generate the data-derived template. Name | Time-span | Observing | Galactic | Period | Ecliptic | DM | Log | Observations
---|---|---|---|---|---|---|---|---
|  | Sites | Coo. [deg] | [ms] | Latitude [deg] |  [pc/cm3] | M(eDM) | per template
J0034$-$0721 | 2013-08 2019-08 | 7 | $110.4$ | $-69.8$ | $943.0$ | $-10.15$ | $10.9$ | $-2.57$ | $215$
J0137+1654 | 2013-09 2019-05 | 6 | $138.4$ | $-44.6$ | $414.8$ | $6.31$ | $26.1$ | $-2.64$ | $199$
J0151$-$0635 | 2013-09 2019-05 | 5 | $160.4$ | $-65.0$ | $1464.7$ | $-16.84$ | $25.8$ | $-2.02$ | $234$
J0525+1115 | 2013-09 2019-05 | 5 | $192.7$ | $-13.2$ | $354.4$ | $-11.93$ | $79.4$ | $-2.52$ | $207$
J0528+2200 | 2013-08 2019-08 | 6 | $183.9$ | $-6.9$ | $3745.5$ | $-1.24$ | $50.9$ | $-2.84$ | $181$
J0538+2817 | 2014-02 2019-04 | 6 | $179.7$ | $-1.7$ | $143.2$ | $4.94$ | $40.0$ | $-2.15$ | $190$
J0540+3207 | 2016-03 2019-05 | 5 | $176.7$ | $0.8$ | $524.3$ | $8.76$ | $62.1$ | $-3.13$ | $122$
J0543+2329 | 2013-08 2019-04 | 5 | $184.4$ | $-3.3$ | $246.0$ | $0.10$ | $77.7$ | $-2.39$ | $237$
J0609+2130 | 2013-10 2019-05 | 7 | $189.2$ | $1.0$ | $55.7$ | $-1.92$ | $38.7$ | $-2.91$ | $197$
J0612+3721 | 2013-09 2019-05 | 6 | $175.4$ | $9.1$ | $298.0$ | $13.95$ | $27.2$ | $-2.99$ | $216$
J0614+2229 | 2013-08 2019-08 | 5 | $188.8$ | $2.4$ | $335.0$ | $-0.90$ | $96.9$ | $-2.98$ | $196$
J0629+2415 | 2013-08 2019-05 | 4 | $188.8$ | $6.2$ | $476.6$ | $0.99$ | $84.2$ | $-2.98$ | $189$
J0659+1414 | 2013-09 2019-08 | 5 | $201.1$ | $8.3$ | $384.9$ | $-8.44$ | $14.1$ | $-1.96$ | $251$
J0823+0159 | 2013-08 2019-05 | 4 | $222.0$ | $21.2$ | $864.9$ | $-16.90$ | $23.8$ | $-2.78$ | $220$
J0837+0610 | 2013-08 2019-05 | 6 | $219.7$ | $26.3$ | $1273.8$ | $-11.98$ | $12.9$ | $-3.88$ | $202$
J0943+1631 | 2013-08 2019-05 | 5 | $216.6$ | $45.4$ | $1087.4$ | $2.70$ | $20.3$ | $-1.77$ | $255$
J0953+0755 | 2013-08 2019-05 | 7 | $228.9$ | $43.7$ | $253.1$ | $-4.62$ | $3.0$ | $-3.27$ | $276$
J1024$-$0719 | 2012-12 2019-08 | 1 | $251.7$ | $40.5$ | $5.2$ | $-16.04$ | $6.5$ | $-3.32$ | $68$
J1543$-$0620 | 2013-12 2019-05 | 6 | $0.6$ | $36.6$ | $709.1$ | $13.06$ | $18.4$ | $-3.49$ | $216$
J1705$-$1906 | 2016-03 2019-04 | 4 | $3.2$ | $13.0$ | $299.0$ | $3.72$ | $22.9$ | $-2.3$ | $153$
J1744$-$1134 | 2012-12 2019-08 | 7 | $14.8$ | $9.2$ | $4.1$ | $11.80$ | $3.1$ | $-3.77$ | $394$
J1820$-$0427 | 2013-08 2019-05 | 4 | $25.5$ | $4.7$ | $598.1$ | $18.88$ | $84.3$ | $-2.54$ | $211$
J1825$-$0935 | 2013-12 2019-08 | 5 | $21.4$ | $1.3$ | $769.0$ | $13.71$ | $19.4$ | $-3.29$ | $184$
J1834$-$0426 | 2013-08 2019-05 | 5 | $27.0$ | $1.7$ | $290.1$ | $18.73$ | $79.4$ | $-2.63$ | $156$
J1848$-$1952 | 2016-03 2019-05 | 5 | $14.8$ | $-8.3$ | $4308.2$ | $3.09$ | $18.3$ | $-2.25$ | $103$
J1900$-$2600 | 2013-08 2019-05 | 5 | $10.3$ | $-13.5$ | $612.2$ | $-3.29$ | $37.9$ | $-2.58$ | $218$
J1913$-$0440 | 2013-09 2019-08 | 6 | $31.3$ | $-7.1$ | $825.9$ | $17.53$ | $89.4$ | $-2.83$ | $1$
J2051$-$0827 | 2013-06 2019-03 | 2 | $39.2$ | $-30.4$ | $4.5$ | $8.85$ | $20.7$ | $-3.58$ | $67$
J2222$-$0137 | 2016-03 2019-05 | 6 | $62.0$ | $-46.1$ | $32.8$ | $7.98$ | $3.3$ | $-3.56$ | $37$

##  Appendix C Variations in the SW impact on the DM time series

Section 4.2 reports the amplitude of the spherical SW model as variable across
the years for a number of pulsars. However, Tiburzi et al. (2019) (and this
article) showed the spherical model to be an imperfect SW approximation. As an
additional confirmation of the detected variability of the SW contribution, we
also examined the DM variations themselves, after subtracting the IISM
component. Figure 9 shows the averages of DM variations as a function of the
Solar elongation (in 5 degree bins, up to 50 degrees) from all the 14 pulsars,
organized by ecliptic latitude and time of observation. As for the amplitudes
of the spherical SW model, it is possible to verify that also the magnitude of
the DM variations decreases in time within the ecliptic latitude ranges
spanning from 20 to 5 degrees and $-$5 to $-$10 degrees. For example, for
pulsars with ecliptic latitude ranging from 15 and 20 degrees, the average DM
variations at Solar elongations smaller than 20 degrees decrease from
$4.2\times 10^{-4}$ pc/cm-3 to $5\times 10^{-5}$ pc/cm-3. Similarly, pulsars
whose ecliptic latitudes range between 5 and 10 degrees have DM variations at
Solar elongations lower than 10 degrees that drop from $1.2\times 10^{-3}$
pc/cm-3 to $5.2\times 10^{-4}$ pc/cm-3 across the timespan. On the other hand,
the DM variations for pulsars with ecliptic latitudes included between $-$5
and 5 degrees tend to fluctuate without showing a clear decreasing trend. This
confirms the trend seen in Figure 6.

Figure 9: Average DM variations (after the subtraction of the IISM) as a
function of Solar elongation, divided per time and ecliptic latitude. All data
points are shown with vertical error bars.

## References

  * Allen (1947) Allen, C. W. 1947, 107, 426
  * Archibald et al. (2018) Archibald, A. M., Gusinskaia, N. V., Hessels, J. W. T., et al. 2018, Nature, 559, 73
  * Armstrong et al. (1995) Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995, ApJ, 443, 209
  * Bandura et al. (2014) Bandura, K., Addison, G. E., Amiri, M., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9145, Ground-based and Airborne Telescopes V, ed. L. M. Stepp, R. Gilmozzi, & H. J. Hall, 914522
  * Bird et al. (1994) Bird, M. K., Volland, H., Paetzold, M., et al. 1994, ApJ, 426, 373
  * Bondonneau et al. (2020) Bondonneau, L., Grießmeier, J. M., Theureau, G., et al. 2020, arXiv e-prints, arXiv:2009.02076
  * Braun et al. (2015) Braun, R., Bourke, T., Green, J. A., Keane, E., & Wagg, J. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 174
  * Burke-Spolaor et al. (2019) Burke-Spolaor, S., Taylor, S. R., Charisi, M., et al. 2019, A&A Rev., 27, 5
  * Caballero et al. (2018) Caballero, R. N., Guo, Y. J., Lee, K. J., et al. 2018, MNRAS, 481, 5501
  * Coles (1996) Coles, W. A. 1996, Ap&SS, 243, 87
  * Cordes et al. (2016) Cordes, J. M., Shannon, R. M., & Stinebring, D. R. 2016, ApJ, 817, 16
  * Demorest et al. (2013) Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., et al. 2013, ApJ, 762, 94
  * Donner et al. (2020) Donner, J. Y., Verbiest, J. P. W., Tiburzi, C., et al. 2020, arXiv e-prints, arXiv:2011.13742
  * Donner et al. (2019) Donner, J. Y., Verbiest, J. P. W., Tiburzi, C., et al. 2019, A&A, 624, A22
  * Edwards et al. (2006) Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, 372, 1549
  * Finkbeiner (2003) Finkbeiner, D. P. 2003, 146, 407
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Guhathakurta & Fisher (1998) Guhathakurta, M. & Fisher, R. 1998, 499, L215
  * Guhathakurta & Fisher (1995) Guhathakurta, M. & Fisher, R. R. 1995, 22, 1841
  * Hobbs et al. (2020) Hobbs, G., Guo, L., Caballero, R. N., et al. 2020, MNRAS, 491, 5951
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302
  * Jackson et al. (2020) Jackson, B. V., Buffington, A., Cota, L., et al. 2020, Frontiers in Astronomy and Space Sciences, 7, 76
  * Janssen et al. (2015) Janssen, G., Hobbs, G., McLaughlin, M., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 37
  * Jonas (2009) Jonas, J. L. 2009, IEEE Proceedings, 97, 1522
  * Kuenkel (2017) Kuenkel, L. 2017, LOFAR Studies of Interstellar Scintillation
  * Lam et al. (2016) Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2016, ApJ, 819, 155
  * Lam et al. (2018) Lam, M. T., McLaughlin, M. A., Cordes, J. M., Chatterjee, S., & Lazio, T. J. W. 2018, ApJ, 861, 12
  * Lazarus et al. (2016) Lazarus, P., Karuppusamy, R., Graikou, E., et al. 2016, MNRAS, 458, 868
  * Lentati et al. (2016) Lentati, L., Shannon, R. M., Coles, W. A., et al. 2016, MNRAS, 458, 2161
  * Lorimer & Kramer (2004) Lorimer, D. R. & Kramer, M. 2004, Handbook of Pulsar Astronomy
  * Madison et al. (2019) Madison, D. R., Cordes, J. M., Arzoumanian, Z., et al. 2019, ApJ, 872, 150
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993
  * Manchester & Taylor (1972) Manchester, R. N. & Taylor, J. H. 1972, 10, 67
  * Muhleman & Anderson (1981) Muhleman, D. O. & Anderson, J. D. 1981, 247, 1093
  * Niu et al. (2017) Niu, Z.-X., Hobbs, G., Wang, J.-B., & Dai, S. 2017, Research in Astronomy and Astrophysics, 17, 103
  * Perera et al. (2019a) Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019a, MNRAS, 490, 4666
  * Perera et al. (2019b) Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019b, MNRAS, 490, 4666
  * Poedts & Pomoell (2017) Poedts, S. & Pomoell, J. 2017, in EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, 7396
  * Porayko et al. (2019) Porayko, N. K., Noutsos, A., Tiburzi, C., et al. 2019, MNRAS, 483, 4100
  * Schwenn (1990) Schwenn, R. 1990, Large-Scale Structure of the Interplanetary Medium, ed. R. Schwenn & E. Marsch, 99
  * Schwenn (2006) Schwenn, R. 2006, Living Reviews in Solar Physics, 3, 2
  * Shao et al. (2013) Shao, L., Caballero, R. N., Kramer, M., et al. 2013, Classical and Quantum Gravity, 30, 165019
  * Siemens et al. (2013) Siemens, X., Ellis, J., Jenet, F., & Romano, J. D. 2013, Classical and Quantum Gravity, 30, 224015
  * Stappers et al. (2011) Stappers, B. W., Hessels, J. W. T., Alexov, A., et al. 2011, A&A, 530, A80
  * Tiburzi (2018) Tiburzi, C. 2018, PASA, 35, e013
  * Tiburzi et al. (2016) Tiburzi, C., Hobbs, G., Kerr, M., et al. 2016, MNRAS, 455, 4339
  * Tiburzi et al. (2019) Tiburzi, C., Verbiest, J. P. W., Shaifullah, G. M., et al. 2019, MNRAS, 487, 394
  * Tingay et al. (2013) Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007
  * Tokumaru et al. (2020) Tokumaru, M., Tawara, K., Takefuji, K., Sekido, M., & Terasawa, T. 2020, Sol. Phys., 295, 80
  * Vallisneri et al. (2020) Vallisneri, M., Taylor, S. R., Simon, J., et al. 2020, ApJ, 893, 112
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
  * van Straten & Bailes (2011) van Straten, W. & Bailes, M. 2011, PASA, 28, 1
  * van Straten et al. (2012) van Straten, W., Demorest, P., & Osłowski, S. 2012, Astronomical Research and Technology, 9, 237
  * Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267
  * Voisin et al. (2020) Voisin, G., Cognard, I., Freire, P. C. C., et al. 2020, A&A, 638, A24
  * You et al. (2012) You, X. P., Coles, W. A., Hobbs, G. B., & Manchester, R. N. 2012, MNRAS, 422, 1160
  * You et al. (2007a) You, X. P., Hobbs, G., Coles, W. A., et al. 2007a, 378, 493
  * You et al. (2007b) You, X. P., Hobbs, G. B., Coles, W. A., Manchester, R. N., & Han, J. L. 2007b, 671, 907


Prospects for diffractive dijet photoproduction at the EIC

V. Guzey1 and M. Klasen2$\star$

1 National Research Center “Kurchatov Institute”, Petersburg Nuclear Physics
Institute (PNPI), Gatchina, 188300, Russia

2 Institut für Theoretische Physik, Westfälische Wilhelms-Universität Münster,
Wilhelm-Klemm-Straße 9, 48149 Münster, Germany

*<EMAIL_ADDRESS>
Proceedings for the XXVIII International Workshop
on Deep-Inelastic Scattering and Related Subjects,
Stony Brook University, New York, USA, 12-16 April 2021
10.21468/SciPostPhysProc.?

## Abstract

We discuss the prospects of diffractive dijet photoproduction at the EIC to
distinguish different fits of diffractive proton PDFs, different schemes of
factorization breaking, to determine diffractive nuclear PDFs and pion PDFs
from leading neutron production.

##  1 Introduction

The HERA electron-proton collider, which operated at DESY Hamburg from 1992 to
2007, has provided us with a wealth of experimental information on the strong
interaction and the quark and gluon structure of the proton [1]. In
particular, parton density functions (PDFs) and the strong coupling constant
$\alpha_{s}$ were determined with unprecedented precision, not only in
inclusive deep-inelastic scattering (DIS), but also from jet production in DIS
and photoproduction. It came as a surprise that a large fraction of events
(e.g. 10-15% in DIS) were diffractive, leaving the scattered proton intact or
in an excited state. Using QCD [2] and assuming Regge factorization [3], these
processes could be described in terms of diffractive PDFs.

In photoproduction [4], factorization for diffractive processes turned out to
be broken [5], similarly to observations at hadron colliders [6]. The question
whether this occurred only for photons resolved into their hadronic structure
or also for direct photons (and thus globally) could not be answered, also
because the two processes become related at next-to-leading order (NLO) of QCD
and beyond [7].

The advent of a new electron-ion collider (EIC) at BNL raises hopes that a
solution for this long-standing problem might be found, although we will see
that this will not be easy due to the lower beam energies of the EIC compared
to HERA and the ensuing dominance of direct photon processes [8]. However, the
EIC will also provide access to so far completely unknown nuclear diffractive
PDFs, which would then allow one to test theoretical models of nuclear
shadowing and its relation to diffraction [9]. Finally, a leading neutron
spectrometer would permit to resume studies of the virtual pion cloud in the
nucleon and thus of the pion PDFs at low momentum fractions $x$ [10].

##  2 Diffractive dijet photoproduction and HERA results

At HERA, electrons or positrons of 27.5 GeV energy collided with protons of
first 820, then 920 GeV. Instrumentation with forward taggers allowed the
study of both photoproduction and diffraction. Dijet events accounted for a
large fraction of the cross section (e.g. 10-20% in DIS) [11]. The ratio of
diffractive over inclusive dijets in photoproduction was 1-4% [12].

Figure 1: Diffractive production of dijets with invariant mass $M_{12}$ in
direct (left) and resolved (right) photon-pomeron collisions, leading to the
production of one or two additional remnant jets. The hadronic systems $X$ and
$Y$ are separated by the largest rapidity gap in the final state.

The two mechanisms of direct and resolved diffractive dijet photoproduction
are illustrated in Fig. 1. In inclusive DIS (i.e. with virtual photons
$\gamma^{*}$), diffractive PDFs were extracted by both the H1 [13] and ZEUS
[14] collaborations and tested successfully in DIS dijet and heavy quark
production. When only a rapidity gap between the hadronic systems $X$ and $Y$
was imposed and the intact proton was not measured as in the H1 2006 Fits A
and B, the normalization was larger by a factor of $1.23$ than when the
diffractive PDFs were determined with a leading proton as in ZEUS 2009 Fit SJ.
The H1 2006 Fits A and B differed mostly in their singular vs. regular gluon
behavior at large observed parton momentum fraction $z_{I\\!\\!P}^{\rm obs.}$
in the exchanged pomeron $I\\!\\!P$.

##  3 Electron-proton collisions at the EIC

At the EIC, electrons of energy 21…18 GeV will collide with protons of energy
100…275 GeV, leading to center-of-mass energies of $\sqrt{s}=92...141$ GeV.
Since only longitudinal momentum fractions $y$ and $x_{I\\!\\!P}$,
$x_{\gamma}^{\rm obs.}$ and $z_{I\\!\\!P}^{\rm obs.}$ are transmitted from the
electron and proton to the photon and pomeron and their partonic constituents,
the invariant mass $M_{12}$ of the hard dijet system is typically small and
requires low kinematic cuts. From experience at HERA, the (asymmetric) cuts on
the transverse momenta $p_{T1,(2)}$ can be as low as 5 (4.5) GeV [15]. The
$K$-factor of NLO/LO cross sections is about two for an anti-$k_{T}$ jet
radius of $R=1$ except at the kinematic edges, where it can be larger [16].

As one can see in Fig. 2, with 21 GeV $\times$ 100 GeV beams and the expected
luminosity the two H1 2006 Fits A and B (which is similar to ZEUS 2009 Fit SJ)
could be distinguished and the evolution of the diffractive PDFs tested,
although only over relatively small ranges in $z_{I\\!\\!P}^{\rm obs.}$ and
$Q^{2}\sim\bar{p}_{T}^{2}$.

With 18 GeV $\times$ 275 GeV beams, the range in jet $p_{T}$ then extends to
12 rather than only 8 GeV and in $x_{\gamma}^{\rm obs.}$ and
$z_{I\\!\\!P}^{\rm obs.}$ down to 0.2 rather than only 0.5, giving at least
some access to the resolved photon contribution as shown in Fig. 3. As also
shown there, a similar effect can be obtained by increasing the cut on the
pomeron momentum fraction $x_{I\\!\\!P}$ from 0.03 to 0.10, although this also
increases the otherwise subleading reggeon contribution. These conditions give
at least some hope to distinguish global from resolved-photon-only
factorization breaking at the EIC [16].

Figure 2: Left: Distributions in the observed pomeron momentum fraction
$z_{I\\!\\!P}^{\rm obs.}$ for three different sets of diffractive PDFs at NLO
QCD. Right: Evolution in $Q^{2}\sim\bar{p}_{T}^{2}$ of the contribution from
gluons in the pomeron. Here, $x_{I\\!\\!P}<0.03$ is imposed.

Figure 3: Distributions in the observed photon momentum fraction
$x_{\gamma}^{\rm obs.}$ in the 18 GeV $\times$ 275 GeV and 21 GeV $\times$ 100
GeV beam configurations (left) and with a higher cut on $x_{I\\!\\!P}<0.10$
(right).

##  4 Electron-ion collisions at the EIC

Coherent diffraction on nuclei can be selected with a large rapidity gap and
by requiring that no neutrons are produced in the zero-degree calorimeter. At
small $x$ nuclear diffractive PDFs are expected to be suppressed compared to
the coherent sum of free nucleon diffractive PDFs due to nuclear shadowing. In
the model of leading twist nuclear shadowing, it can be shown that [9]

$f_{i/A}^{D}(z_{I\\!\\!P},Q^{2},x_{I\\!\\!P})\approx A\cdot 0.65\cdot
f_{i/p}^{D}(z_{I\\!\\!P},Q^{2},x_{I\\!\\!P}).$ (1)

In addition, since the nuclei stay intact, the diffractive free nucleon PDFs
obtained from processes including diffractive dissociation like H1 2006 Fits A
and B must be divided by a factor of 1.23.

Figure 4: NLO QCD cross sections for coherent diffractive dijet
photoproduction $eA\to e^{\prime}+2\ {\rm jets}+X+A$ with various nuclear
beams and a center-of-mass energy per nucleon of $\sqrt{s}=92$ GeV at the EIC.
The cross sections are shown as functions of the jet average transverse
momentum (top left) and rapidity difference (bottom left) as well as the
longitudinal momentum fractions in the photon (top right) and pomeron (bottom
right).

Our NLO predictions for diffractive dijet photoproduction on different nuclei
are shown in Fig. 4. Scaling these predictions with a global suppression
factor of 0.5 or with a factor 0.04 for resolved-only factorization breaking
as expected from the commonly used two-state eikonal model[17] leads to very
similar results, making a distinction of these two schemes again difficult
[16].

##  5 Conclusion

Identifying the mechanism of factorization breaking in diffractive
photoproduction is a desideratum of the HERA program which should be addressed
at the EIC. However, we have seen that this will not be easy due to the lower
EIC center-of-mass energy and the ensuing dominance of direct photon
processes. Access to nuclear diffractive PDFs is, however, a novel unique
opportunity at the EIC, where in particular the connection of diffraction and
nuclear shadowing can be addressed experimentally.

## Acknowledgements

Financial support by the DFG through the grant KL 1266/9-1 within the
framework of the joint German-Russian project “New constraints on nuclear
parton distribution functions at small $x$ from dijet production in $\gamma A$
collisions at the LHC" is gratefully acknowledged. VG’s research is also
supported in part by RFBR, research project 17-52-12070.

## References

  * [1] P. Newman and M. Wing,  _The Hadronic Final State at HERA_ ,  Rev. Mod. Phys. 86(3), 1037 (2014),  10.1103/RevModPhys.86.1037,  1308.3368.
  * [2] J. C. Collins,  _Proof of factorization for diffractive hard scattering_ ,  Phys. Rev. D 57, 3051 (1998),  10.1103/PhysRevD.61.019902,  [Erratum: Phys.Rev.D 61, 019902 (2000)],  hep-ph/9709499.
  * [3] G. Ingelman and P. E. Schlein,  _Jet Structure in High Mass Diffractive Scattering_ ,  Phys. Lett. B 152, 256 (1985),  10.1016/0370-2693(85)91181-5.
  * [4] M. Klasen,  _Theory of hard photoproduction_ ,  Rev. Mod. Phys. 74, 1221 (2002),  10.1103/RevModPhys.74.1221,  hep-ph/0206169.
  * [5] M. Klasen and G. Kramer,  _Review of factorization breaking in diffractive photoproduction of dijets_ ,  Mod. Phys. Lett. A 23, 1885 (2008),  10.1142/S0217732308027461,  0806.2269.
  * [6] M. Klasen and G. Kramer,  _Survival probability for diffractive dijet production in p anti-p collisions from next-to-leading order calculations_ ,  Phys. Rev. D 80, 074006 (2009),  10.1103/PhysRevD.80.074006,  0908.2531.
  * [7] V. Guzey and M. Klasen,  _A fresh look at factorization breaking in diffractive photoproduction of dijets at HERA at next-to-leading order QCD_ ,  Eur. Phys. J. C 76(8), 467 (2016),  10.1140/epjc/s10052-016-4304-5,  1606.01350.
  * [8] R. Abdul Khalek _et al._ ,  _Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report_ (2021),  2103.05419.
  * [9] L. Frankfurt, V. Guzey and M. Strikman,  _Leading Twist Nuclear Shadowing Phenomena in Hard Processes with Nuclei_ ,  Phys. Rept. 512, 255 (2012),  10.1016/j.physrep.2011.12.002,  1106.2091.
  * [10] M. Klasen and G. Kramer,  _Photoproduction of jets on a virtual pion target in next-to-leading order QCD_ ,  Phys. Lett. B 508, 259 (2001),  10.1016/S0370-2693(01)00513-5,  hep-ph/0103056.
  * [11] F. D. Aaron _et al._ ,  _Jet Production in ep Collisions at High $Q^{2}$ and Determination of $\alpha_{s}$_,  Eur. Phys. J. C 65, 363 (2010),  10.1140/epjc/s10052-009-1208-7,  0904.3870.
  * [12] F. D. Aaron _et al._ ,  _Diffractive Dijet Photoproduction in ep Collisions at HERA_ ,  Eur. Phys. J. C 70, 15 (2010),  10.1140/epjc/s10052-010-1448-6,  1006.0946.
  * [13] A. Aktas _et al._ ,  _Measurement and QCD analysis of the diffractive deep-inelastic scattering cross-section at HERA_ ,  Eur. Phys. J. C 48, 715 (2006),  10.1140/epjc/s10052-006-0035-3,  hep-ex/0606004.
  * [14] S. Chekanov _et al._ ,  _A QCD analysis of ZEUS diffractive data_ ,  Nucl. Phys. B 831, 1 (2010),  10.1016/j.nuclphysb.2010.01.014,  0911.4119.
  * [15] M. Klasen and G. Kramer,  _Dijet cross-sections at O( $\alpha\alpha_{s}^{2}$) in photon-proton collisions_,  Phys. Lett. B 366, 385 (1996),  10.1016/0370-2693(95)01352-0,  hep-ph/9508337.
  * [16] V. Guzey and M. Klasen,  _Diffractive dijet photoproduction at the EIC_ ,  JHEP 05, 074 (2020),  10.1007/JHEP05(2020)074,  2004.06972.
  * [17] V. A. Khoze, A. D. Martin and M. G. Ryskin,  _Soft diffraction and the elastic slope at Tevatron and LHC energies: A MultiPomeron approach_ ,  Eur. Phys. J. C 18, 167 (2000),  10.1007/s100520000494,  hep-ph/0007359.