Testing the accuracy of the ionospheric Faraday rotation corrections through LOFAR observations of bright northern pulsars

TL;DR

This study evaluates the accuracy of ionospheric Faraday rotation corrections using LOFAR pulsar observations, demonstrating that with optimal models, corrections can reach an accuracy of approximately 0.06–0.07 rad m$^{-2}$ over a year.

Contribution

It assesses the performance of different ionospheric models in correcting Faraday rotation in pulsar data, highlighting the potential for high-precision measurements.

Findings

Optimal ionospheric maps achieve correction accuracy of ~0.06–0.07 rad m$^{-2}$ over one year.

Residual unmodelled effects introduce systematics and correlated noise.

Comparison of free-electron density maps informs best practices for ionospheric correction.

Abstract

Faraday rotation of polarized emission from pulsars measured at radio frequencies provides a powerful tool to investigate the interstellar and interplanetary magnetic fields. However, besides being sensitive to the astrophysical media, pulsar observations in radio are affected by the highly time-variable ionosphere. In this article, the amount of ionospheric Faraday rotation has been computed by assuming a thin layer model. For this aim, ionospheric maps of the free electron density (based on Global Positioning System data) and semi-empirical geomagnetic models are needed. Through the data of five highly polarized pulsars observed with the individual German LOw-Frequency ARray stations, we investigate the performances of the ionospheric modelling. In addition, we estimate the parameters of the systematics and the correlated noise generated by the residual unmodelled ionospheric effects, and show the comparison of the different free-electron density maps. For the best ionospheric maps, we have found that the rotation measure corrections on one-year timescales after subtraction of diurnal periodicity are accurate to $\sim$ 0.06--0.07 rad m$^{-2}$.

Figures

Tables

Table 1. Table 2: Estimation of the noise parameters based on the Bayesian analysis of RM residuals using POMME10 geomagnetic model and different ionospheric maps. The results for other two considered in this paper geomagnetic models (EMM and IGRF) are indistinguishable within the uncertainties and are demonstrated in Supplementary material online. The analysis was based on datasets of PSR J0332+5434, PSR J1136+1551 and PSR J0814+7429, respectively (see Table 1 ). The used noise model is the one described in Eq. ( 9 ). We report the median and 1-sigma uncertainty values of the amplitude of the Lorentzian spectrum A L subscript 𝐴 L A_{\text{L}} [rad m -2 ], which effectively quantifies the measurement uncertainties of the ionospheric RM corrections; the turnover frequency of the Lorentzian spectrum f 0 subscript 𝑓 0 f_{0} [day -1 ], the amplitude of the 1-day harmonic signal A d subscript 𝐴 d A_{\text{d}} [rad m -2 ] and maximum likelihood estimation of RM IISM subscript RM IISM \text{RM}_{\text{IISM}} [rad m -2 ]. The latter is assumed to be constant across a timescale of several months. The factor f B = 1.11 subscript 𝑓 𝐵 1.11 f_{B}=1.11 was applied.

		PSR J0332+5434
Model	$A_{\text{L}}^{\text{med}}$	$f_0^{\text{med}}$	$A_{\text{d}}^{\text{med}}$	${\text{RM}}_{\text{IISM}}^{\text{ML}}$
UQRG	${0.045}_{-0.002}^{+0.003}$	${1.5}_{-0.2}^{+0.2}$	${0.012}_{-0.006}^{+0.007}$	$-$64.16
JPLG	${0.050}_{-0.002}^{+0.002}$	${1.7}_{-0.2}^{+0.2}$	${0.025}_{-0.005}^{+0.007}$	$-$64.21
EHRG	${0.054}_{-0.003}^{+0.003}$	${1.2}_{-0.1}^{+0.1}$	${0.012}_{-0.008}^{+0.007}$	$-$64.05
IGSG	${0.060}_{-0.003}^{+0.003}$	${1.2}_{-0.1}^{+0.1}$	${0.02}_{-0.007}^{+0.005}$	$-$64.08
ESAG	${0.068}_{-0.004}^{+0.005}$	${1.1}_{-0.1}^{+0.1}$	${0.025}_{-0.009}^{+0.007}$	$-$64.05
UPCG	${0.073}_{-0.004}^{+0.003}$	${0.9}_{-0.1}^{+0.1}$	${0.025}_{-0.014}^{+0.008}$	$-$64.17
CODG	${0.12}_{-0.01}^{+0.01}$	${0.29}_{-0.06}^{+0.06}$	${0.063}_{-0.009}^{+0.009}$	$-$63.95

Table 2. Table 3: Estimation of the noise parameters based on the Bayesian analysis of RM residuals using EMM geomagnetic model and different ionospheric maps. The results are presented in the same fashion as in Table 2 of the paper. The analysis was based on datasets of PSR J0332+5434, PSR J1136+1551 and of PSR J0814+7429. We report the median and 1-sigma uncertainty values of the amplitude of the Lorentzian spectrum A L subscript 𝐴 L A_{\text{L}} [rad m -2 ], which effectively quantifies the measurement uncertainties of the ionospheric RM corrections; the turnover frequency of the Lorentzian spectrum f 0 subscript 𝑓 0 f_{0} [day -1 ], the amplitude of the 1-day harmonic signal A d subscript 𝐴 d A_{\text{d}} [rad m -2 ] and maximum likelihood estimation of RM IISM subscript RM IISM \text{RM}_{\text{IISM}} [rad m -2 ]. The latter is assumed to be constant across a timescale of several months. The factor f B = 1.11 subscript 𝑓 𝐵 1.11 f_{B}=1.11 was applied.

		PSR J0332+5434
Model	$A_{\text{L}}^{\text{med}}$	$f_0^{\text{med}}$	$A_{\text{d}}^{\text{med}}$	${\text{RM}}_{\text{IISM}}^{\text{ML}}$
UQRG	${0.045}_{-0.003}^{+0.003}$	${1.5}_{-0.1}^{+0.2}$	${0.012}_{-0.006}^{+0.007}$	-64.16
JPLG	${0.051}_{-0.002}^{+0.003}$	${1.8}_{-0.2}^{+0.2}$	${0.025}_{-0.005}^{+0.007}$	-64.21
EHRG	${0.054}_{-0.003}^{+0.003}$	${1.2}_{-0.1}^{+0.1}$	${0.012}_{-0.008}^{+0.007}$	-64.06
IGSG	${0.060}_{-0.003}^{+0.003}$	${1.2}_{-0.1}^{+0.1}$	${0.02}_{-0.006}^{+0.005}$	-64.08
ESAG	${0.068}_{-0.005}^{+0.005}$	${1.2}_{-0.1}^{+0.1}$	${0.025}_{-0.009}^{+0.007}$	-64.07
UPCG	${0.073}_{-0.004}^{+0.003}$	${0.9}_{-0.1}^{+0.1}$	${0.025}_{-0.014}^{+0.008}$	-64.17
CODG	${0.11}_{-0.01}^{+0.01}$	${0.29}_{-0.06}^{+0.06}$	${0.063}_{-0.009}^{+0.009}$	-63.95

Table 3. Table 4: Analogous to Table 2 and 3 , but with IGRF geomagnetic model used for ionospheric RM modeling.

		PSR J0332+5434
Model	$A_{\text{L}}^{\text{med}}$	$f_0^{\text{med}}$	$A_{\text{d}}^{\text{med}}$	${\text{RM}}_{\text{IISM}}^{\text{ML}}$
UQRG	${0.045}_{-0.002}^{+0.002}$	${1.5}_{-0.1}^{+0.2}$	${0.012}_{-0.006}^{+0.007}$	-64.16
JPLG	${0.049}_{-0.002}^{+0.002}$	${1.7}_{-0.3}^{+0.3}$	${0.025}_{-0.005}^{+0.007}$	-64.21
EHRG	${0.054}_{-0.003}^{+0.003}$	${1.2}_{-0.1}^{+0.1}$	${0.012}_{-0.008}^{+0.007}$	-64.05
IGSG	${0.060}_{-0.003}^{+0.003}$	${1.2}_{-0.1}^{+0.1}$	${0.02}_{-0.007}^{+0.005}$	-64.09
ESAG	${0.068}_{-0.004}^{+0.005}$	${1.1}_{-0.1}^{+0.1}$	${0.025}_{-0.009}^{+0.007}$	-64.05
UPCG	${0.073}_{-0.004}^{+0.003}$	${0.9}_{-0.1}^{+0.1}$	${0.025}_{-0.014}^{+0.008}$	-64.17
CODG	${0.11}_{-0.01}^{+0.01}$	${0.29}_{-0.06}^{+0.06}$	${0.063}_{-0.009}^{+0.009}$	-63.95

Equations

$$\psi ={\psi }_0+\text{RM}{\lambda }^2,$$

$$\text{RM}=0.81{\int }_{\text{source}}^{\text{observer}}{n}_{\text{e}}\mathbf{B}\cdot d\mathbf{r}[\text{rad}{\text{m}}^{-2}],$$

$${\text{RM}}_{\text{iono}}=2.6\times 10^{-17}\text{STEC}\times {\text{B}}_{\text{iono}}{\text{ rad m}}^{-2},$$

$${\text{RM}}_{\text{obs}}=\underset{\text{deterministic}}{\underbrace{{\text{RM}}_{\text{iono}}+{\text{RM}}_{\text{1day}}+{\text{RM}}_{\text{IISM}}}}+\underset{\text{stochastic}}{\underbrace{{\text{RM}}_{\text{noise}}+\mathbf{n}}}.$$

$$S(f)=\frac{A_{\text{L}}^2}{f_0\left[1+{\left(\frac{f}{f_0}\right)}^2\right]},$$

$$C_{\text{L}}=A_{\text{L}}^2\exp (-f_0\tau ),$$

$$C_{\text{WN}}={\sigma }_i^2{\delta }_{ij},$$

$$\begin{array}{r}\log P_{\text{pst}}(\mathbf{\Theta })\sim \log P_{\text{pr}}(\mathbf{\Theta })\\ -\sum _{i=1}^N\frac{1}{2}\left({\text{RM}}_{\text{obs}}-{\text{RM}}_{\text{iono}}-{\text{RM}}_{\text{1day}}-{\text{RM}}_{\text{IISM}}\right)\times \\ C^{-1}\times \left({\text{RM}}_{\text{obs}}-{\text{RM}}_{\text{iono}}-{\text{RM}}_{\text{1day}}-{\text{RM}}_{\text{IISM}}\right)-\frac{1}{2}\ln (2\pi \text{det}C),\end{array}$$

$${\text{RM}}_{\text{obs}}=\frac{{\text{RM}}_{\text{iono}}}{f_B}+{\text{RM}}_{\text{1day}}+{\text{RM}}_{\text{IISM}}+{\text{RM}}_{\text{noise}}+\mathbf{n}.$$

$$\begin{array}{r}\Delta \text{RM}\simeq -0.81n_{\text{e}}B\cdot v_{\perp }T\sin \theta +0.81n_{\text{e}}B\cdot v_{\parallel }T\cos \theta \\ \sim 3\times 10^{-6}{\text{rad m}}^{-2}{\left(\frac{L}{1\text{ kpc}}\right)}^{-1}\left(\frac{\text{RM}}{30{\text{ rad m}}^{-2}}\right)\left(\frac{|{\mathbf{v}}_{\perp }|}{100{\text{ km s}}^{-1}}\right)\left(\frac{T}{\text{ yr}}\right)\end{array}$$

$$\begin{array}{r}{\text{rms}}_{\text{RM}}\sim 0.81\sqrt{n_{\text{e}}^2{\sigma }_B^2+B_{\parallel }^2{\sigma }_n^2}L\sim \sqrt{{\int }_{1/T}^{\infty }{\text{PSD}}_{\text{KL}}(f)\text{d}f}\\ =6\times 10^{-5}{\text{rad m}}^{-2}{\left(\frac{L}{1\text{kpc}}\right)}^{\frac{1}{2}}{\left(\frac{|v|}{100{\text{km s}}^{-1}}\right)}^{\frac{5}{6}}{\left(\frac{T}{\text{yr}}\right)}^{\frac{5}{6}}\end{array}$$

$$\begin{array}{r}{Q}_{\text{mod}}=I_{\text{mod}}p\cos 2\chi \cos 2\psi =I_{\text{mod}}p\cos 2\chi \cos (2\text{RM}{\lambda }^2+2{\psi }_0),\\ U_{\text{mod}}=I_{\text{mod}}p\cos 2\chi \sin 2\psi =I_{\text{mod}}p\cos 2\chi \sin (2\text{RM}{\lambda }^2+2{\psi }_0).\end{array}$$

$$\begin{array}{r}{\sigma }_q=\frac{Q}{I}\sqrt{{\left(\frac{{\sigma }_Q}{Q}\right)}^2+{\left(\frac{{\sigma }_I}{I}\right)}^2},\\ {\sigma }_u=\frac{U}{I}\sqrt{{\left(\frac{{\sigma }_U}{U}\right)}^2+{\left(\frac{{\sigma }_I}{I}\right)}^2},\end{array}$$

$$\begin{array}{rl}& q_{\text{obs},i}=A{c}_i+B{s}_i+{\gamma }_q+{ϵ}_{q,i},\text{and}\\ & u_{\text{obs},i}=A{s}_i-B{c}_i+{\gamma }_u+{ϵ}_{u,i},\end{array}$$

$$P_{\text{pst}}(\text{parameters}|\text{data})=\frac{P_{\text{pr}}(\text{parameters})P(\text{data}|\text{parameters})}{P(\text{data})},$$

$$\begin{array}{r}{P}_{\text{posterior}}(A,B,{\gamma }_q,{\gamma }_u,\text{RM},\eta |q_{\text{obs},i},u_{\text{obs},i})\propto \\ {\Pi }_{i=1}^{N_{\text{ch}}}\frac{1}{\sqrt{2\pi }{\sigma }_{q,i}\eta }\exp \left(-\frac{(q_{\text{obs},i}-A{c}_i-B{s}_i-{\gamma }_q{)}^2}{2({\sigma }_{q,i}\eta {)}^2}\right)\times \\ \times {\Pi }_{i=1}^{N_{\text{ch}}}\frac{1}{\sqrt{2\pi }{\sigma }_{u,i}\eta }\exp \left(-\frac{(u_{\text{obs},i}+A{s}_i-B{c}_i-{\gamma }_u{)}^2}{2({\sigma }_{u,i}\eta {)}^2}\right).\end{array}$$

$$\begin{array}{r}\ln P_{\text{posterior}}\left(A,B,{\gamma }_q,{\gamma }_u,\text{RM}|q_{\text{obs},i},u_{\text{obs},i}\right)\propto \\ -\frac{1}{2}\sum _{i=1}^{N_{\text{ch}}}\left[\frac{{\left(q_{\text{obs},i}-A{c}_i-B{s}_i-{\gamma }_q\right)}^2}{{\sigma }_{q,i}^2}+\frac{{\left(u_{\text{obs},i}+A{s}_i-B{c}_i-{\gamma }_u\right)}^2}{{\sigma }_{u,i}^2}\right]=\\ \frac{1}{2}(-\widehat{YY}+2A\widehat{YC}+2B\widehat{YS}+2{\gamma }_q{Y}_q+2{\gamma }_u{Y}_u-A^2\widehat{CC}-\\ -B^2\widehat{SS}-{\gamma }_q^2{W}_q-{\gamma }_u^2{W}_u-2A{\gamma }_q{C}_q-2A{\gamma }_u{C}_u-2B{\gamma }_q{S}_q\\ -2B{\gamma }_u{S}_u).\end{array}$$

$$\begin{array}{rl}{W}_q& =\sum _{i=1}^N{\omega }_{q,i},\text{and}{W}_u=\sum _{i=1}^N{\omega }_{u,i},\\ Y_q& =\sum _{i=1}^N{\omega }_{q,i}{q}_{\text{obs},i},\text{and}{Y}_u=\sum _{i=1}^N{\omega }_{u,i}{u}_{\text{obs},i},\\ \widehat{YY}& =\sum _{i=1}^N{\omega }_{q,i}{q}_{\text{obs},i}^2+{\omega }_{u,i}{u}_{obs,i}^2,\\ \widehat{YC}& =\sum _{i=1}^N{\omega }_{q,i}{q}_{\text{obs},i}{c}_i+{\omega }_{u,i}{u}_{obs,i}{s}_i,\\ \widehat{YS}& =\sum _{i=1}^N{\omega }_{q,i}{q}_{\text{obs},i}{s}_i-{\omega }_{u,i}{u}_{obs,i}{c}_i,\\ \widehat{CC}& =\sum _{i=1}^N{\omega }_{q,i}{c}_i^2+{\omega }_{u,i}{s}_i^2,\widehat{SS}=\sum _{i=1}^N{\omega }_{q,i}{s}_i^2+{\omega }_{u,i}{c}_i^2,\\ C_q& =\sum _{i=1}^N{\omega }_{q,i}{c}_i,\text{and}{C}_u=\sum _{i=1}^N{\omega }_{u,i}{c}_i,\\ S_q& =\sum _{i=1}^N{\omega }_{q,i}{s}_i\text{ and}{S}_u=\sum _{i=1}^N{\omega }_{q,i}{s}_i,\end{array}$$

$$\begin{array}{r}{\omega }_{q,i}=\frac{1}{{\sigma }_{q,i}^2}\text{ and }{\omega }_{u,i}=\frac{1}{{\sigma }_{u,i}^2}.\end{array}$$

$$\begin{array}{r}{P}_{\text{posterior}}\left(\text{RM}|q_{\text{obs}},u_{\text{obs}}\right)\propto \\ \frac{1}{\sqrt{|4DF-E^2|\widehat{CC}\widehat{SS}}}\exp \left(M-\frac{\widehat{YY}}{2}+\frac{D{G}^2-EGJ+F{J}^2}{E^2-4DF}\right),\end{array}$$

$$\begin{array}{rl}D& =\frac{C_q^2\widehat{SS}+S_q^2\widehat{CC}-W_q\widehat{CC}\widehat{SS}}{2\widehat{CC}\widehat{SS}},\\ F& =\frac{C_u^2\widehat{SS}+S_u^2\widehat{CC}-W_u\widehat{CC}\widehat{SS}}{2\widehat{CC}\widehat{SS}},\\ E& =\frac{C_u{C}_q\widehat{SS}+S_u{S}_q\widehat{CC}}{\widehat{CC}\widehat{SS}},\\ J& =\frac{C_q\widehat{YC}\widehat{SS}+S_q\widehat{YS}\widehat{CC}-Y_q\widehat{CC}\widehat{SS}}{\widehat{CC}\widehat{SS}},\\ G& =\frac{C_u\widehat{YC}\widehat{SS}+S_u\widehat{YS}\widehat{CC}-Y_u\widehat{CC}\widehat{SS}}{\widehat{CC}\widehat{SS}}\text{ and }\\ M& =\frac{{\widehat{YC}}^2\widehat{SS}+{\widehat{YS}}^2\widehat{CC}}{2\widehat{CC}\widehat{SS}}.\end{array}$$

$${\hat{\eta }}^2=-\frac{2}{2N_{\text{ch}}-4}\left[M-\frac{\widehat{YY}}{2}+\frac{D{G}^2-EGJ+F{J}^2}{E^2-4DF}\right].$$

Full text

††pubyear: 2018††pagerange: Testing the accuracy of the ionospheric Faraday rotation corrections through LOFAR observations of bright northern pulsars–4

1 Introduction

Since their discovery (Hewish et al., 1968), pulsars have been a powerful tool to probe the magnetoionic plasma. Due to frequency-dependent dispersion delay and scattering of their signals, pulsars can be used to study, e.g., turbulence in the ionised interstellar medium (ISM) on many orders of magnitude (e.g. Rickett, 1977; Armstrong et al., 1995; You et al., 2007a), the distribution of free electrons in the Milky Way and the Local Bubble (e.g. Cordes & Lazio, 2002; Bhat et al., 1998), and the electron content of the Solar wind (e.g. You et al., 2007b; Howard et al., 2016). Magnetised plasma also induces Faraday rotation in linearly polarised radiation, that is, a rotation of the polarization angle $\psi$ depending on the radiation wavelength $\lambda$,

[TABLE]

where RM is the rotation measure defined as:

[TABLE]

with $n_{\text{e}}$ being the electron density in the ionised ISM $[\text{cm}^{-3}]$, $\bm{B}$ the magnetic field $[\mu\text{G}]$ and $d\bm{r}$ [pc] the infinitesimal interval of the distance along the line of sight (LoS). From the above expressions one can see that more accurate RM estimations can be achieved with broad-band instruments operating at longer wavelengths.

Due to the high percentage of linear polarisation, and low levels of magnetospheric Faraday rotation (e.g. Wang et al., 2011), pulsars are useful objects to measure RM induced by the ionised ISM, and hence the Galactic magnetic fields (e.g. Han et al., 2018).

Because the propagation effects are strongly depend on $\lambda$, low frequencies are favoured for studies of these effects in pulsars111For strongly Faraday rotated sources, such as pulsars in the dense regions(e.g. magnetar in the Galactic center) and distant active galactic nuclei, RMs can be as well effectively probed with instruments, operating at 1-2 GHz and higher frequencies.. Moreover, the steep spectra of pulsars (e.g. Bates et al., 2013) and the reduction of the linear polarization fractions at high frequencies in pulsar emission (e.g. Johnston et al., 2008), make the low-frequency band even more preferable for Faraday rotation studies.

Recently, low-frequency pulsar astronomy was revived thanks to a number of cutting-edge facilities such as the Long Wavelength Array (LWA) (Ellingson et al., 2009), the Murchison Widefieled Array (MWA) (Tingay et al., 2013), the Giant Ukrainian Radio Telescope (GURT) (Zakharenko et al., 2016) and the LOw Frequency ARray (LOFAR) (van Haarlem et al., 2013; Stappers et al., 2011).

Nevertheless, polarisation studies at low frequencies are challenging. Besides the effects of the magnetised ionised ISM, linearly polarised radiation can be noticeably rotated by the highly variable terrestial ionosphere. Moreover, it can significantly depolarise observations when averaging over several hours. For a review on the propagation of radio waves through the ionosphere, see e.g. Wilson et al. (2013) or Thompson et al. (2001).

In order to mitigate the ionospheric contribution to FR, numerous techniques have been developed. One very promising approach is based on providing quasi-simultaneous observations of a known background source (e.g. the diffuse polarised background), located within the ionospheric correlation spatial scale with respect to the source of interest, to recover the ionospheric Faraday rotation (Lenc et al., 2016).

Alternatively, the ionospheric Faraday rotation can be estimated by combining models of the ionospheric electron density and of the geomagnetic field. In the majority of the studies that aimed to measure the interstellar Faraday rotation in astronomical sources (e.g. Weisberg et al., 2004; Han et al., 2006; Yan et al., 2011b), the ionospheric electron density was computed through the semi-empirical International Reference Ionosphere (IRI) model (Bilitza et al., 2014), which provides monthly-averaged ionospheric electron density profiles up to 2000 km, as a function of time and location. However, due to the sparsely distributed ground and space observatories that contribute to the IRI model, and the large averaging time, the modelled values of electron densities can significantly deviate from the real ones (Mosert et al., 2007). Higher accuracies can be reached by a technique described in Erickson et al. (2001), where the ionospheric electron densities are obtained through raw dual-frequency GPS data, recorded with a set of local GPS receivers. When applied to PSR J1932+1059, the variance of the differences between the observed RM as obtained at the VLA, and the predicted ionospheric RM as computed with the AIPS APGPS routine222A similar approach is implemented in the ALBUS software https://github.com/twillis449/ALBUS_ionosphere, was found to be 0.2 rad m*-2*.

A handier and less computationally expensive alternative to this approach consists in using global ionospheric maps of electron column densities in the ionosphere, which are based on the available data from all the GPS stations spread around Earth. This technique was implemented and tested on a set of pulsars by Sotomayor-Beltran et al. (2013), showing a qualitatively good agreement between the expected and the observed values of FR. However, Sotomayor-Beltran et al. (2013) have restricted their analysis to probing only two global ionospheric maps (ROBR and CODG), and the research was carried out on a set of observations with timespans of only several hours. The standard deviations of the residuals between the RMs, observed and modelled with CODG and ROBR, varied for different datasets in the ranges 0.12–0.20 rad m*-2* and 0.07–0.20 rad m*-2*, respectively.

In this article, we aim to compare the performance, and estimate the accuracy of different publicly available global ionospheric maps, when applied to correct for ionospheric Faraday rotation in several months of pulsar data. For these goals, we used pulsar observations obtained with the international LOFAR stations in Germany. In Section 2 we describe the instrumental and observational setup and our data reduction, including a first, application of a simple ionospheric modelling. In Section 3 we attempt to model the ionospheric Faraday rotation in our dataset and we analyse the systematics left in the RM residuals after ionospheric mitigation. In Sections 3.1, 3.2 we focus on how to correct for the systematics, and show the results obtained after the implementation of our additional corrections and the comparison of different global ionospheric maps. In Section 4 we then summarise our findings.

2 Observations and Data Reduction

LOFAR, the LOw-Frequency ARray, is an international interferometric telescope operating at very low frequencies, from 10 up to 240 MHz (van Haarlem et al., 2013). LOFAR stations are distributed throughout Europe, with a dense core, the Superterp, located in the Netherlands. Six of the stations are located in Germany: DE601 in Effelsberg, DE602 in Unterweilenbach, DE603 in Tautenburg, DE604 in Bornim, DE605 in Jülich and DE609 in Norderstedt.

Each LOFAR station outside the Netherlands consists of 96 pairs of dipoles of low band antennas (LBAs) and 96 tiles (i.e., groups formed of 16 pairs of dipoles) of high-band antennas (HBAs). Typically three days a week the German stations are used as stand-alone telescopes by the GLOW (German Long Wavelength) consortium333GLOW is an association of German universities and research institutes, which promotes the use of the meter wavelength spectral window for astrophysical purposes. GLOW members operate the work of the German LOFAR stations and GLOW is further involved at the planned Square Kilometer Array (SKA) project. More details can be found in https://www.glowconsortium.de/index.php/en/, to perform an observing campaign of pulsars at low frequencies ($\sim$100–200 MHz) using the HBAs. Commonly, each German station in stand-alone mode observes a unique set of pulsars with $\sim$2-hours integration time per pulsar. The specifics for the dataset used in the presented analysis are summarised in Table 1. As it will be shown in the next sections, we focus our analysis on the characterization of potential short- and long-term trends in the residuals between the ionospheric models and the data. All the selected pulsars have high signal-to-noise ratio (S/N), which varies from $\sim$800 up to $\sim$2000, and a significant fraction of linear polarization (at least 10%), which allows us to measure the RM with high accuracy and precision. Besides this, for the purposes of the short-term analysis, we chose pulsars with a significant fraction of long, continuous observations (from a few hours to entire days). This allows us to properly identify also high-frequency systematics. For the purposes of the long-term analysis, this last requirement is not strictly necessary, and we thus selected pulsars with a long observing baselines.

After digitizing and beamforming, the data have 5.12$\mu$s time resolution and are split into frequency channels of 195 kHz bandwidth. Due to data rate limitations, only 366 (488 for DE601) channels are recorded on machines at the Max-Planck-Institut für Radioastronomie in Bonn and at the Jülich Supercomputing Centre using the LuMP (LOFAR und MPIfR Pulsare) Software444https://github.com/AHorneffer/lump-lofar-und-mpifr-pulsare. The datasets are then, coherently de-dispersed, folded modulo the pulse period and reduced to more manageable 10-second sub-integrations with the DSPSR software555http://dspsr.sourceforge.net/ (van Straten & Bailes, 2011), and stored as PSRFITS archives (Hotan et al., 2004). We then excise the radio-frequency interference with the COASTGUARD’s clean.py surgical algorithm (Lazarus et al., 2016).

In contrast to steerable radio telescopes, the LOFAR antennas are fixed on the ground, which causes a distortion of the polarisation signal, as well as decrease of the intrinsic signal intensity, towards low elevations, due to the projection effects (Noutsos et al., 2015). For instance, such an instrumental response is responsible for the so-called instrumental peak at [math] rad m*-2* in the RM spectrum while performing the RM synthesis analysis (Burn, 1966; Brentjens & de Bruyn, 2005). We mitigate these instrumental effects by applying a Jones calibration matrix based on the Hamaker measurement equations (see Hamaker et al., 1996; Smirnov, 2011). However, Noutsos et al. (2015) showed that across several hours of observations taken with the Superterp, the intrinsic signal intensity of the LOFAR antennas significantly degrades at low elevations ($\lesssim$ $30\degr$) even after the calibration procedure has been applied.

Due to the fact that radio observations in the LOFAR frequency band are quite sensitive to the highly variable ionospheric layer (van Weeren et al., 2016), we split pulsar archives into 15-minute subintegrations with the PSRCHIVE software package666http://psrchive.sourceforge.net/ (van Straten et al., 2012), which corresponds to the minimum time-sampling of ionospheric maps that we have tested (see Section 2.1). This reduces the unresolved contribution of ionospheric RM, while still providing a reasonable S/N.

After this, we estimate the RM for each of the 15-minute subintegrations, building an RM time series for each of the analysed datasets. For this, we use an optimised version of the classical RM synthesis technique, described in Appendix A.

In Fig. 1 we demonstrate the example of the Bayesian Generalised Lomb-Scargle Periodogram (BGLSP) application to one of the 15-min observation of PSR J1136+1551. We clearly see systematic deviation from the modelled $Q$ and $U$, which is reflected in the spectrum as a low-frequency excess of power around 0 rad/m2. The origin of these systematics is not known for certain, but it is highly likely that it is associated with instrumental properties, e.g. non-linearity in the instrumental setup. Because the spurious peak affects a small range of values around 0 rad/m2, we expect sources with significant larger RMs to be uneffected, suggesting little or no influence on our results. However, we point out that the results can be biased when dealing with astronomical sources with low RM values. In order to prove these considerations, we have performed two tests. Firstly, we tested the basic assumption that any discrepancies between the models and the data are induced by an effect that is strongly frequency dependent. Therefore, we have split data into two sub-bands and measured RM values separately for the bottom and upper half of the bandwidth. The results show that both RM values are in excellent agreement within the uncertainties. This suggests that the effect is not strongly depending on frequency. Still, we also tested whether a systematic effect could conspire to mimic a wrong RM value. As a worst case scenario, we have investigated the impact of systematics, in case they had a quadratic dependency on frequency, which would mimic the $\lambda^{2}$ dependency introduced by the physical effect of Faraday rotation. The simulated Stokes Q and U were evenly sampled in frequency with a realistic 20% of data loss due to radio-frequency interference. We run a Monte Carlo simulation with $10^{3}$ realizations of this set-up, for increasing values of RMs from 0 to 20 rad/m2. A range of the systematic amplitudes were tested with reduced $\chi^{2}$ of up to 10, as the reduced $\chi^{2}$ detected in the data did not exceed this value. We found that, starting from an RM value of $\sim\!6$ rad/m2, the mean and variance of the distribution of the recovered RMs are in a good agreement with the results from BGLSP (see Fig. 2). This behaviour is expected, since as soon as the source RM is larger than the width of the systematic feature, the two signals can be separated reliably.

With the reliability of our RM measurements established by these tests, we proceed to do a first attempt to mitigate the Faraday rotation ionospheric contribution.

2.1 On modelling the ionospheric RM variations: thin layer ionospheric model

If not taken into account, the ionosphere introduces noise in the measured RM values. This makes it impossible, for instance, to investigate RM variations caused by the turbulent ionised ISM, which are expected to be $\sim$3–4 orders of magnitude lower than the root-mean-square (rms) of the ionospheric RM fluctuations (see Eq. (11)). We now briefly recap the ionospheric RM behavior and the ways to model it.

The ionospheric layer, partially consisting of free electrons and positively charged ionised molecules and atoms, extends from 50 km to beyond 2000 km above the Earth’s surface (Rishbeth & Garriott, 1969). The ionospheric contribution to RM can be estimated to be of the order of 1–4 rad m*-2*, however, the essential complexity in treating the ionospheric RM comes from its strong variability, which typically changes during the day up to 80$\%$. The ionization fraction of the ionospheric shell, mostly caused by photoionization processes involving the Sun’s extreme ultra-violet and X-ray emission, varies significantly over timescales of minutes (due to Solar flares) up to years (11-year Solar cycle). Besides this, the ionosphere shows diurnal (caused by the relative motion of the Sun on the celestial sphere) and 27-day periodicities (due to the Solar rotation). As the Earth’s atmosphere is not homogeneous and different molecules are dominating at different heights, the ionospheric shell, does not have a homogeneous electron density distribution, and achieves its maximum during the day time in the so-called F sublayer, which implies $\sim$ 50-60 % of all the electrons in the ionosphere (Bilitza et al., 2017).

Because of this, the ionosphere can be reasonably well modelled by a thin shell located at the effective ionospheric height, which is usually estimated to be between 300 and 600 km above the Earth’s surface.

As the projected thickness of the non-uniform ionospheric layer increases out from the zenith, it is common practice to discard data at low elevations. For this work we have used a $30\degr$ elevation cut-off777This number is partially motivated by the limitations of the polarization calibration method used in this work. In the case of the ionosphere, and with the mentioned assumptions, Eq. (2) is reduced to (Sotomayor-Beltran et al., 2013):

[TABLE]

where STEC (Slant TEC, where TEC stands for ‘Total Electron Content’) is equal to the column density of electrons [m*-2*] at the cross-section between the LoS and the ionospheric shell and $B_{\text{iono}}$ is the projection of the magnetic field [G] in the F-layer on the LoS. The thin layer approximation has already been implemented in several codes aimed at the estimation of the ionospheric RM (e.g. in Sotomayor-Beltran et al., 2013). In particular, for the work presented here we use the publicly available RMextract software888https://github.com/lofar-astron/RMextract, that estimates the ionospheric RM along a certain LoS and at a certain point in time making use of a geomagnetic field model and a global ionospheric map. An example of ionospheric RM calibration with RMextract, applied to the RM sequence of PSR J0332+5434, is demonstrated on Fig. 3 (upper panel). From here on in this paper for demonstration purposes we have used JPLG maps, which have showed the second best result in our analysis and are commonly available for the majority of our observing epochs.

The geomagnetic field models are conventionally represented as spherical harmonical expansions of a scalar magnetic potential. Several geomagnetic models are publicly available, among which are the Enhanced Magnetic Model (EMM)999https://www.ngdc.noaa.gov/geomag/EMM/, the International Geomagnetic Reference Field (IGRF) (Thébault et al., 2015), the World Magnetic Model101010https://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml, and POMME10111111http://geomag.org/models/pomme10.html. The lower panel of Fig. 4 shows a comparison of the ionospheric magnetic field given by EMM, POMME10 and IGRF12, for the years 2013 through 2018 for lines of sight from Germany in the direction of 30°in elevation (minimum elevation used in our work). The plot demonstrates clear systematic behaviour, although, on average between 2013 and 2018, there is less than 0.1% difference between different geomagnetic models. The discrepancy seems to be increasing with time. Thus, for the future datasets taken around 2020 geomagnetic models with non-evolving with time geomagnetic parameters will reach few per cent level difference between them and should be used with care. Fig. 5 demonstrates that for low elevation observations this difference can hit 1% from the absolute value.

We have conducted a full analysis by making use of all three geomagnetic models. In order to be concise, we present only the results of POMME10 (Maus et al., 2006) here (see Table 2). In the case one of the other two geomagnetic models the results on parameter estimation and the presence of various systematics in the data remain unchanged, and are presented in Supplementary material online (according to Table 2).

The global ionospheric maps (publicly available121212ftp://cddis.gsfc.nasa.gov/gnss/products/,131313ftp://igs.ensg.ign.fr/) in IONEX141414https://igscb.jpl.nasa.gov/igscb/data/format/ionex1.pdf format, provide estimates of the vertical TEC. Several global ionospheric maps are available: CODG (from the University of Bern), ESAG and EHRG (European Space Agency), JPLG (Jet Propulsion Laboratory), UPCG and UQRG (Technical University of Catalonia, see Orús et al. (2005)), IGSG (International GNSS Service). Although the maps can be based on the same GPS data, the published TEC values can vary from group to group because of different interpolation schemes and different spatial and temporal resolution. In practice, the maps we have used, all have a spatial resolution of $2.5\degr\times 5\degr$(latitude $\times$ longitude). CODG and EHRG have a time resolution of 1 hour, UQRG of 0.25 hours and the remaining maps (ESAG, IGSG, JPLG and UPCG) have a time resolution of 2 hours.

3 Systematics in the RM residuals

The residual RM series, after subtraction of the ionospheric model from the observed RM values, show the presence of correlated structures and strong colored noise. For instance, Fig. 6 shows the residual RMs for three different pulsars across about a 2-month long timespan. As mentioned in the previous section, the ionospheric RM was corrected by using the RMextract software package, the POMME10 geomagnetic field model, and the CODG/JPLG maps. The CODG time series show similar trend (e.g., the jump of the RM around the 15th of May 2016) despite the pulsars being significantly separated on the sky. The magnitude of these RM variations significantly exceeds those expected from astrophysical sources (see Section 4 for details). This implies that the origin of the correlations is not interstellar, but an insufficient modelling of the ionosphere. Moreover, if these maps provide unreliable information about the uncertainties of the TEC values, this will affect the uncertainties of the modelled ionospheric RM and, in turn, our ability to determine the significance of astrophysical RM variations.

In order to solve for these issues, we conduct an independent search of the systematics in the modelled ionospheric RMs, with the aim of obtaining good estimates of the white-noise level and the uncertainties for the ionospheric RM time series.

Some of the observed structures in the residuals can be well explained by a diurnal sinusoid with amplitude $A_{\text{d}}$, the effect of which is demonstrated in Fig. 7, 8 and Fig. 3 B,C. Besides being responsible for a 1-day peak in the power spectrum of the RM residuals, it also creates a 1-year pseudo-periodicity in the data, as the transit time of the source shifts gradually during the day across the year.

After subtracting the 1-day sinusoid, the spectrum shows obvious evidence of red noise at high frequencies, and evolves into a white noise plateau at low frequencies (see Fig. 9). Such a spectrum is described in our model by a Lorentzian function, also known as an Ornstein-Uhlenbeck process (Uhlenbeck & Ornstein, 1930). In the next two paragraphs we provide the reader with the mathematical description of the found systematics and introduce the criteria for the comparison of different ionospheric maps.

3.1 Analysis of RM residuals on timescales up to one year

The observational evidences, discussed in Section 3 allow us to define a mathematical model to describe the contributions to the observed RM time series with timescale shorter than a year. As will be pointed out in Section 4, the interstellar contribution $\text{RM}_{\text{IISM}}$ is very small, and typically only visible on timescales of order of several years (Yan et al., 2011a). We have restricted the dataset considered in this section to only several months (Table 1), so we can assume this parameter to be constant and, thus, the contribution $\text{RM}_{\text{IISM}}$ will not bias the estimates of the parameters of the systematics.

The vector $\textbf{RM}_{\text{obs}}=[\text{RM}_{\text{t}_{1}},\text{RM}_{\text{t}_{2}},...,\text{RM}_{\text{t}_{N}}]$ that contains the RM time series of a certain pulsar observed at $N$ epochs $t_{i}$ can be seen as a combination of deterministic and stochastic contributions:

[TABLE]

$\textbf{RM}_{\text{iono}}$ stands for the semi-empirical ionospheric thin layer model of RM variations described in Section 2.1. $\textbf{RM}_{\text{1day}}$ is the harmonic signal with 1-day period, that can be parametrised as $\textbf{RM}_{\text{1day}}=A_{\text{d}}\sin(2\pi\bm{t}/1\text{day}+\phi)$. $\textbf{RM}_{\text{noise}}$ and $\bm{n}$ are stochastic noise contributions. $\textbf{RM}_{\text{noise}}$ is given by the plateau of the Lorentzian spectrum, whose one-sided spectral density is described as:

[TABLE]

with $A_{\text{L}}$ [rad m*-2*] being the already mentioned amplitude of the stochastic signal and $f_{0}$ [day*-1*] the turnover frequency. From this expression it can be easily shown that, while behaving like red noise on a short time scales, $\textbf{RM}_{\text{noise}}$ reduces to white noise for $f\ll f_{0}$ with a constant variance $A_{\text{L}}^{2}$. By applying the Wiener-Khinchin theorem, the variance-covariance matrix of this process is then given by:

[TABLE]

where $\tau=2\pi|t_{i}-t_{j}|$ with $t_{i}$ and $t_{j}$ are two different epochs.

The uncorrelated white noise component $\bm{n}$ in Eq. (4), coming from the measurement noise of Stokes parameters (see Fig. 13), has a flat power spectral density with variance-covariance of the form:

[TABLE]

with $\delta_{ij}$ being a Kronecker delta and $\bm{\sigma}$ the vector of the formal uncertainties of the observed RMs151515The uncertainties are modified by a factor $\eta$ ,see Eq. (22), determined via the Bayesian Lomb-Scargle Periodogram described at the end of Appendix A.

In order to investigate the properties of the stochastic and deterministic signals that emerge in the RM residuals after the ionospheric correction, we use Bayesian inference in the time domain. Given the model in Eq. (4), and assuming that the stochastic parts are drawn from random Gaussian processes, the posterior probability for the unknown parameters $\mathbf{\Theta}=[A_{\text{L}},f_{0},A_{\text{d}},\phi,\text{RM}_{\text{IISM}}]$ is written as:

[TABLE]

where $C=C_{\text{L}}+C_{\text{WN}}$.

We have applied the model discussed in this paragraph to pulsar datasets that span less than a year (outlined in Table 1). The high S/Ns of our pulsars make us more sensitive to the signals generated by the imperfections of ionospheric RM modelling, described by Eq. (9). The factor $f_{B}$ was fixed to 1.11, as found in Section 3.2, and, thus, was excluded from the set of free parameters of the model.

In order to explore the 5-dimensional parameter space $\mathbf{\Theta}$ of Eq. (4), we ran a Markov Chain Monte Carlo simulation, using the Bayesian inference tool MultiNest (Feroz et al., 2009). The priors of the parameters $P_{\text{pr}}(\Theta)$ were chosen to be uninformative (see Caballero et al., 2016): uniform for $A_{\text{L}}$, $f_{0}$, $\phi$, and $\text{RM}_{\text{IISM}}$ and log-uniform for $A_{\text{d}}$. A representative example of the obtained 2-dimensional posterior probability plot is shown in Fig. 12.

Because we only used data from a limited group of pulsars, and probed a statistically not significant sample of LoS, the values given in Table 2 should not be treated as definitive solutions. However, our results give a qualitative estimation of the accuracy of different ionospheric maps.

3.2 Analyis of RM residuals on timescales beyond one year

After subtracting $\textbf{RM}_{\text{iono}}$ and $\textbf{RM}_{\text{1day}}$ from the observed RM time series, the long term datasets still show a deterministic linear trends. The trend is not visible on a timescale of months, but it becomes obvious across several years (see Fig. 10). We find that such linear trend can be suppressed by scaling the $\textbf{RM}_{\text{iono}}$ time series (resulting from a thin layer model) by a constant factor $f_{B}$. In other words, Eq. (4) is modified as

[TABLE]

A positive trend of the order of $1\text{--}2\times 10^{-4}$ rad m*-2* day*-1* was noticed in four pulsars observed with three GLOW stations. Removal of the linear trend, by applying the factor $f_{B}$, is demonstrated in Fig. 10.

Making use of pulsar datasets that span more than one year (outlined in Table 1) and all six global ionospheric maps, considered in this paper, the least square fit estimate was found to be $f_{B}=1.11^{+0.04}_{-0.04}$. The decrease in root-mean-square (rms), mostly due to the elimination of the linear trend, is illustrated in Fig. 11. The results for all three geomagnetic models are identical.

There are several physical interpretations possible for the factor $f_{B}$. Among them is the possible overestimation of the geomagnetic strength, $B_{\text{LoS}}$, and/or the underestimation of the ionospheric effective height (Birch et al., 2002) due to the poor knowledge about the electron density in Earth’s plasmasphere. For instance, $f_{B}=1.11$ is equivalent to an increase of the effective height from 450 km up to $\sim$ 700 km. An explanation to the trend might be searched in the complex dynamical behavior of the ionospheric effective height. As a matter of fact, it has been shown in multiple investigations (e.g. Hernández-Pajares et al., 2011; Arora et al., 2016), that the ionospheric effective height can vary from 300 up to 800 km, depending on the time of day, season, and level of Solar activity. For instance, the 11-year sunspot cycle, the last maximum of which was in 2014, can cause significant ionization in the ionospheric layer, thus both increasing the ionospheric thickness and ionospheric effective height (Liu et al., 2007). One promising way to improve the model is by using the effective heights determined via the IRI-Plas software (Gulyaeva et al., 2013), which takes into account plasmasphere contribution (Arora et al., 2016).

The nature of the factor $f_{B}$ is still under investigation and is planned to be tested on a larger sample of pulsars in the future.

4 Discussions and Conclusions

In this paper, we have characterised and investigated the deterministic and stochastic RM variations generated by the ionospheric layer through pulsar observations taken with the German LOFAR stations. The main day-to-day variability was modelled by assuming a thin-layer ionosphere, located at 450 km above the Earth’s surface. For this model, the magnetic field was taken from the publicly available geomagnetic maps (POMME10, EMM, IGRF), while the information about electron densities was extracted from the selection of different global ionospheric maps. Besides that, an additional signal peaked at a frequency of 1 day*-1* in the power spectrum, which was significant in almost all processed datasets, and was removed by including in the model a 1-day period sinusoid. The residual noise could be described by a Lorentzian spectrum, which behaves like white noise on long timescales and defines our sensitivity to long-term RM variations. The parameters of the model were estimated by applying a Bayesian framework to the RM time series of three pulsars. The observed RM for each epoch was determined by using an improved RM synthesis technique, based on BGLSP, which accounts for non-regularly sampled data and constant offsets in Stokes $Q$ and $U$ due to instrumental effects. By running a Markov Chain Monte Carlo, we have estimated the amplitude of the Lorentzian spectrum (or variance of white noise) for all the ionospheric maps. An additional linear trend becomes visible on a timescale of several years. To account for this, we have applied a factor $f_{B}=1.11$ to the ionospheric RM contribution that was modelled by RMextract, $\textbf{RM}_{\text{mod}}$, as determined in Section 3.2. This slightly reduces the level of the Lorentzian spectrum plateau for some of the pulsars, determined in Section 3.1.

Our results for the three pulsars are slightly different. Nevertheless, two of them (PSRs J0332+5434, J0814+7429) are consistent within 2 sigma, while J1136+1551 shows slightly higher values. We show that geomagnetic models mostly agree and that consequently the accuracy of ionospheric RM corrections is dominated by the uncertainties and inaccuracies in ionospheric TEC maps, which we have investigated in the paper. On average UQRG and JPLG, combined with one of the geomagnetic models, show better results than the other ionospheric maps. If one is going to use one of these two maps to correct for RM variations, the variance of the white noise can be conservatively set to 0.06–0.07 rad m*-2* for observations taken in Europe after 1-day sin waves and linear trend have been taken into account. This is approximately an order of the magnitude higher than the uncertainties on the observed RM, obtained from BGLSP, for the pulsars considered. As we have used the data of only three pulsars and our observational sites are located only in Germany, this value can vary, e.g., increasing significantly in places with sparse GNSS station coverage. Thus, in order to get reliable estimates of the sensitivity to long-term RM variations for a specific instrument, we recommend to undertake a similar kind of analysis for their sites independently.

Essentially, the determined values along with BGLSP uncertainties define the sensitivity of RM measurements to astrophysical signals. One of the promising signal of interest, when dealing with Faraday rotation studies, is the time-variable interstellar contribution to the RM. Let us assume that the ionised ISM is homogeneous. Then, the relative motion between a pulsar moving with velocity $\bm{v}$ and an observer can cause temporal RM variations induced by the change both in the projection of the magnetic field on the LoS and in the pulsar distance $L$. By differentiating Eq. (2) under the assumption of a small change between the initial and the final position of a pulsar, we have:

[TABLE]

where $\theta$ is the angle between the magnetic field vector and the LoS, and $T$ is the whole timespan.

Besides this deterministic signal, we expect a time-variable stochastic part of the interstellar contribution, as predicted by the Kolmogorov turbulence (Kolmogorov, 1941). As it was shown in Keith et al. (2013) the power spectral density of the stochastic contribution is $\text{PSD}_{\text{KL}}=0.0112\times D(\tau)\tau^{-\frac{5}{3}}f^{-\frac{8}{3}}$, where $D(\tau)$ is the structure function. The estimated rms of RM will increase with time $T$ (Minter & Spangler, 1996; Xu & Zhang, 2016) as:

[TABLE]

where $\sigma_{B}^{2}$ and $\sigma_{n}^{2}$ are the variances of magnetic field and electron density fluctuation, respectively.

These calculations show that the signals of interest are characterised by a very small amplitudes, of the order of $10^{-5}$–$10^{-4}$ rad m*-2* , which is several orders of magnitude lower than the observed RM variations in this work. From the comparison of the power spectral densities161616$\text{PSD}_{\text{KL}}\simeq\text{PSD}_{\text{WN}}=\sigma^{2}/f_{\text{Ny}}$, where $\text{PSD}_{\text{WN}}$ is the power spectral density of the white noise and $f_{\text{Ny}}$ is the Nyquist frequency of our dataset. we can conclude that we need $\sim$40 years of observations with the current sensitivity (mostly limited by the imperfections of the ionospheric modelling) for this kind of signals to become significant.

More promising signals of astrophysical nature could be registered thanks to extreme scattering events (Coles et al., 2015), associated with the passage of a blob of high density plasma through the LoS, extreme magneto-ionic environment of the source (Desvignes et al., 2018), and coronal mass ejections Howard et al. (2016), which may cause more prominent RM perturbations.

A deeper understanding of the physics of ionospheric behaviour and instrumental GPS biases, along with the development of more regular GPS station arrays in the direct vicinity to the radio telescopes will improve the quality of the estimates of TEC in the ionospheric layer, which will, in turn, increase our sensitivity to the astrophysical RM variations.

Acknowledgements

The authors would like to thank William Coles and Nicolas Caballero for fruitful discussions. NP acknowledges the support from IMPRS Bonn/Cologne and the Bonn-Cologne Graduate School (BCGS). S.O. acknowledges the Australian Research Council grant Laureate Fellowship FL150100148.

This paper is based on data obtained with the German stations constructed by ASTRON International LOFAR Telescope (ILT) (van Haarlem et al., 2013) under project codes LC0_014, LC1_048, LC2_011, LC3_029, LC4_025, LT5_001 (run until cycle 8), LC9_039 and during station-owners time. LOFAR (van Haarlem et al., 2013) is the Low Frequency Array designed and constructed by ASTRON. In this work we made use of data from the Effelsberg (DE601) LOFAR station funded by the Max-Planck–Gesellschaft; 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 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 (German LOng-Wavelength array), 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 by the states of Nordrhein-Westfalia and Hamburg.

This research has made extensive use of NASA’s Astrophysics Data System and the ATNF Pulsar Catalog.

Appendix A Classical RM synthesis technique and the ways to improve it

Here we inspect the Faraday rotation induced by a magnetoionic medium on the radiation of a point source, in the case we are not affected by the effects of multibeam propagation such as differential Faraday rotation and wavelength-dependent polarization (Sokoloff et al., 1998). As mentioned in Section 1, the induced variation in the polarization angle is proportional to the square of the observational wavelength $\lambda$. It can be shown (see Eqs. below) that a nonzero RM gives rise to harmonic signals in both of the Stokes parameters $Q$ and $U$ across the bandwidth, which are shifted by $\pi/2$ with respect to each other. This allows us to estimate not only the absolute value of RM, but also its sign through the RM synthesis and related techniques (Burn, 1966; Brentjens & de Bruyn, 2005). These algorithms can recover the harmonic signal in RM ranges from $\pi/[\lambda_{\text{max}}^{2}-\lambda_{\text{min}}^{2}]$ up to approximately average Nyquist boundary $~{}\pi/2\delta(\lambda^{2})$, where $\delta(\lambda^{2})$ is average size of the sample, determined by the size of the frequency channel, and $\lambda_{\text{min}}$ and $\lambda_{\text{max}}$ are the lowest and highest observational wavelengths, respectively. In our case, the range of available RMs is $~{}0.5$–$120\text{ rad m}^{-2}$. As we are working in a relatively low RM regime, the effects of depolarization in the frequency channels are neglected in this work (Schnitzeler & Lee, 2015).

Mathematically, the problem of RM search can be described in the following way. The expected Stokes $Q_{\text{mod}}$ and $U_{\text{mod}}$ can be expressed as functions of the intrinsic intensity $I_{\text{mod}}$, polarization fraction $p$, and angles $\psi$ and $\chi$, which characterize the polarization ellipse,

[TABLE]

In practice, we can only access measured Stokes parameters $I_{\text{obs}}$, $Q_{\text{obs}}$, $U_{\text{obs}}$, corrupted by noise. For sake of simplicity, as we do not know the intrinsic intensity of the source because we do not perform flux calibration, we use $Q_{\text{obs}}/I_{\text{obs}}$ and $U_{\text{obs}}/I_{\text{obs}}$, denoted as $q$ and $u$, respectively. In this work, we assume that $q$ and $u$ are distributed normally around their mean values with variances:

[TABLE]

where $\sigma_{I}$, $\sigma_{Q}$ and $\sigma_{U}$ are the standard deviations of the observed Stokes parameters in the off-pulse region (see Appendix C).

To recover the RM, we apply to the more classical RM synthesis technique the method of the Bayesian Generalised Lomb-Scargle Periodogram (BGLSP), described in (Mortier et al., 2015; Zechmeister & Kürster, 2009).

By writing $c_{i}=\cos(2\text{RM}\lambda_{i}^{2}-\theta)$ and $s_{i}=\sin(2\text{RM}\lambda_{i}^{2}-\theta)$, the normalised Stokes $q_{\text{obs}}$ and $u_{\text{obs}}$ can be expressed as

[TABLE]

where $A$ and $B$ are the amplitudes of oscillation, $\gamma_{q}$ and $\gamma_{u}$ are the constant offsets associated with the instrumental peak, and $\theta$ is an arbitrary phase reference point, which does not affect the RM and is defined in Appendix B. The noise contributions $\epsilon_{q,i}$ and $\epsilon_{u,i}$ are assumed to be normally distributed, with standard deviations $\sigma_{q,i}$ and $\sigma_{u,i}$ and to vary independently across frequency channels. According to Bayes’ theorem, the posterior probability can be written as

[TABLE]

where $P_{pr}(\text{parameters})$ is the prior distribution of the unknown parameters, $P(\text{data}|\text{parameters})$ is the likelihood function, $P(\text{data})$ is the so-called Bayesian evidence, which is a normalization factor in our case and plays an important role in the problem of model selection. Assuming uniform prior distributions of the parameters, the posterior probability is proportional to a likelihood :

[TABLE]

The resultant form for the posterior probability is analytically marginalised over the nuisance parameters [$A$, $B$, $\gamma_{q}$, $\gamma_{u}$] (see Eq. (20)) and is provided in Appendix B.

As the integration time of pulsar observations is not infinitely small, the rate of change of ionospheric RMs during the integration time will introduce an additional ambiguity to measured RMs, which we have taken into account here by introducing the parameter $\eta$ in the denominator of Eq. (16) (see also Schnitzeler & Lee, 2017). It acts effectively as a multiplier for all the $Q$ and $U$ error bars (see Section B) and is correlated with the reduced $\chi^{2}$ value. Particularly, for our case of a 15-min integration time, $\eta$ typically varies between 1.5 and 3. The parameter $\eta$ was estimated separately, using Eq. (22), and fixed to its maximum likelihood value.

By performing a 1D grid search in the RM parameter space, we can successfully recover the RM posterior probability (see Fig. 1). The uncertainty in the RM value is determined as the variance of the normal Gaussian distribution, fit to the resultant shape of the posterior probability.

Appendix B Bayesian Lomb-Scargle Periodogram

Here we provide the derivation of the marginalised posterior probability $P_{\text{posterior}}$ from Eq. (16). Using similar notations to Mortier et al. (2015), we can determine the part of the expression for the posterior probability which depends on unknown parameters (called Sufficient Statistics) as:

[TABLE]

The cross term $AB\sum_{i=1}^{N_{\text{ch}}}\left(\omega_{q,i}-\omega_{u,i}\right)c_{i}s_{i}$ can be suppressed by assuming that $\tan(2\theta)=\sum_{i=1}^{N_{\text{ch}}}\left(\omega_{q,i}-\omega_{u,i}\right)\sin\left(4\text{RM}\lambda^{2}\right)/\sum_{i=1}^{N_{\text{ch}}}\left(\omega_{q,i}-\omega_{u,i}\right)\cos\left(4\text{RM}\lambda^{2}\right)$. In the above expression the following denominations were used:

[TABLE]

and the weights are defined in a traditional way as:

[TABLE]

The resultant expression for the Sufficient Statistics after marginalization over nuisance parameters $A,B,\gamma_{q},\gamma_{u}$ is

[TABLE]

where

[TABLE]

The resultant expression can be easily generalised to the case of underestimated uncertainties in Stokes $Q$ and $U$ by including an extra free parameter $\eta$, such that $\omega_{q}\rightarrow\eta^{-2}\omega_{q}$ and $\omega_{u}\rightarrow\eta^{-2}\omega_{u}$. In this case the resultant marginalised posterior probability will be the function of two parameters, $P_{\text{posterior}}\left(\text{RM},\eta|q_{\text{obs}},u_{\text{obs}}\right)$. In order to account properly for the fuzzy structure around 0 rad m*-2*, a more complex model of the noise should be used, e.g., a power-law red noise, which is not under consideration in this paper.

In order to determine the unknown parameters within the Bayesian framework, one needs to numerically reconstruct the 2D posterior probability, which can be effectively managed by the Markov Chain Monte Carlo. In the frequentist approach, which is less computationally expensive and was used in this paper, we are interested in the maximum likelihood estimation of the unknown parameters, which for $\eta$ can be found analytically:

[TABLE]

Appendix C Notes on the distribution of Q/I and U/I

A full investigation of the Gaussianity of the $q$ and $u$ distributions is beyond the scope of this paper. However, here we will have some general comment on this.

By postulating that the observed Stokes $I_{\text{obs}}$, $Q_{\text{obs}}$ and $U_{\text{obs}}$ are normally distributed, one can derive that $q$ and $u$ will actually follow a Cauchy-like distributions. In Schnitzeler & Lee (2017) it was showed that the non-Gaussianity of $q$ and $u$ can potentially bias the uncertainty of measured RM in the low S/N regime171717 Schnitzeler & Lee (2017) demonstrated that this problem can be avoided by introducing $N_{\text{ch}}$ (number of frequency channels) nuisance parameters $I_{\text{mod},i}$, and found the analytical expression for the likelihood, marginalised over these parameters in the case of weakly polarised sources ($L_{\text{mod}}<<I_{\text{mod}}$).. However, in the case of high S/N (i.e., $\sigma_{I}/I<0.1$, see Hayya et al. (1975); Kuethe et al. (2000)), the resultant Cauchy distribution can be reasonably well approximated by the normal distribution. By selecting frequency channels above the threshold, and simulating the normally distributed Stokes parameters in each of them, we have reconstructed the RM distribution determined with the BGLSP method. For the two pulsars that we have included in the test (PSRs J0332+5434 and J1136+1551), we have found that the resultant RM distribution can be well approximated by a normal one and its parameters (variance and mean) are in good agreement with those determined with the BGLSP. In Fig. 13 we display the RM distribution for PSR J1136+1551. The properties if the reconstructed distribution are given in caption.

The further analysis of non-Gaussianity of $q$ and $u$ and its influence on the distribution of the resultant RMs will be addressed in future work.

Appendix D Results of the analysis of the systematics for the other geomagnetic models