<figure><img src="image_1.png" /></figure>

HAL
open science

# Variable Scintillation Arcs of Millisecond Pulsars observed with the Large European Array for Pulsars

R.A Main, J Antoniadis, S Chen, I Cognard, H Hu, J Jang, R Karuppusamy,

M Kramer, K Liu, Y Liu, et al.

► To cite this version:

R.A Main, J Antoniadis, S Chen, I Cognard, H Hu, et al. Variable Scintillation Arcs of Millisecond Pulsars observed with the Large European Array for Pulsars. Monthly Notices of the Royal Astronomical Society, 2023, 525 (1), pp.1079-1096. 10.1093/mnras/stad1980. hal-04160212

HAL Id: hal-04160212

https://hal.science/hal-04160212

Submitted on 21 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Variable scintillation arcs of millisecond pulsars observed with the Large European Array for Pulsars

R. A. Main<sup>1,★</sup>, J. Antoniadis<sup>1,2,3</sup>, S. Chen<sup>4</sup>, I. Cognard,<sup>5,6</sup> H. Hu<sup>1</sup>, J. Jang,<sup>1</sup> R. Karuppusamy,<sup>1</sup>

M. Kramer<sup>1,7</sup>, K. Liu,<sup>1</sup> Y. Liu,<sup>8,9</sup> G. Mall,<sup>1,10,11</sup> J. W. McKee,<sup>12,13</sup> M. B. Mickaliger,<sup>7</sup> D. Perrodin,<sup>14</sup>

S. A. Sanidas,<sup>7</sup> B. W. Stappers,<sup>7</sup> T. Sprenger,<sup>1</sup> O. Wucknitz,<sup>1</sup> C. G. Bassa,<sup>15</sup> M. Burgay,<sup>14</sup> R. Concu,<sup>14</sup>

M. Gaikwad,<sup>1</sup> G. H. Janssen,<sup>15,16</sup> K. J. Lee,<sup>4</sup> A. Melis,<sup>14</sup> M. Pilia,<sup>14</sup> A. Possenti,<sup>14,17</sup> L. Wang<sup>4</sup>

and W. W. Zhu<sup>18</sup>

*Affiliations are listed at the end of the paper*

Accepted 2023 June 26. Received 2023 June 18; in original form 2023 February 16

## ABSTRACT

We present the first large sample of scintillation arcs in millisecond pulsars (MSPs), analysing 12 sources observed with the Large European Array for Pulsars (LEAP), and the Effelsberg 100-m telescope. We estimate the delays from multipath propagation, measuring significant correlated changes in scattering time-scales over a 10 yr time span. Many sources show compact concentrations of power in the secondary spectrum, which in PSRs J0613−0200 and J1600−3053 can be tracked between observations, and are consistent with compact scattering at fixed angular positions. Other sources such as PSRs J1643−1224 and J0621+1002 show diffuse, asymmetric arcs which are likely related to phase-gradients across the scattering screen. PSR B1937+21 shows at least three distinct screens which dominate at different times and evidence of varying screen axes or multiscreen interactions. We model annual and orbital arc curvature variations in PSR J0613−0200, providing a measurement of the longitude of ascending node, resolving the sense of the orbital inclination, where our best-fit model is of a screen with variable axis of anisotropy over time, corresponding to changes in the scattering of the source. Unmodelled variations of the screen’s axis of anisotropy are likely to be a limiting factor in determining orbital parameters with scintillation, requiring careful consideration of variable screen properties, or independent very long baseline interferometry (VLBI) measurements. Long-term scintillation studies such as this serve as a complementary tool to pulsar timing, to measure a source of correlated noise for pulsar timing arrays, solve pulsar orbits, and to understand the astrophysical origin of scattering screens.

**Key words:** pulsars: general – ISM: general.

## 1 INTRODUCTION

Pulsar timing arrays (PTAs) involve timing an ensemble of millisecond pulsars (MSPs) at different sky positions to detect nHz gravitational waves (GWs) from coalescing supermassive black holes. Recently, PTAs have detected a common red noise signal, which is a time-correlated signal of similar amplitude and spectrum shared among pulsars in the array (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov et al. 2021b; Antoniadis et al. 2022). While it is possible that a gravitational wave signature is responsible for this effect, there is yet no detection of a spatial correlation that would be a smoking-gun of the gravitational wave background (Hellings & Downs 1983).

The ionized interstellar medium (IISM) is one of the largest contributors of correlated noise to PTAs (see Verbiest & Shaifullah 2018 for a review), and understanding all of its effects is crucial, especially as a GW detection may be imminent. The total column density of electrons induces a $\lambda^2$ dispersive delay (where $\lambda$ is the observing wavelength), where variations are seen prominently in

low-frequency observations (Donner et al. 2020; Tarafdar et al. 2022). Spatial variations in the electron column density results in multipath propagation, resulting in delays scaling as $\lambda^\alpha$, with $\alpha \sim 4.0 \pm 0.6$ (Oswald et al. 2021). In the time domain, this effect can be seen through the broadening of pulses by a characteristic scattering tail [e.g. Bhat et al. 2004, although with the additional complication that different scattered paths can encounter a different electron column (Cordes, Shannon & Stinebring 2016; Donner et al. 2019)]. In the Fourier domain, this is observed as scintillation, where temporal and spectral variations of flux density arise from interference between deflected, coherent images of the pulsar.

Pulsar scintillation is now commonly studied through the secondary spectrum, which is the 2D power spectrum of the scintillation pattern. In this space, many sources have been seen to have ‘scintillation arcs’ (Stinebring et al. 2001), parabolic distributions of power which indicate scattering being dominated by highly localized regions, or ‘thin screens’ (Walker et al. 2004; Cordes et al. 2006). Furthermore, the presence of sharp inverted parabolic ‘arclets’ stemming from the main parabola are seen in some cases, implying strong anisotropy (seen in $\sim 20$ per cent of sources in Stinebring et al. 2022).

<footnote>* E-mail: ramain@mpifr-bonn.mpg.de</footnote>
# Variable scintillation arcs of millisecond pulsars observed with the Large European Array for Pulsars

R. A. Main<sup>1,★</sup>, J. Antoniadis<sup>1,2,3</sup>, S. Chen<sup>4</sup>, I. Cognard,<sup>5,6</sup> H. Hu<sup>1</sup>, J. Jang,<sup>1</sup> R. Karuppusamy,<sup>1</sup>

M. Kramer<sup>1,7</sup>, K. Liu,<sup>1</sup> Y. Liu,<sup>8,9</sup> G. Mall,<sup>1,10,11</sup> J. W. McKee,<sup>12,13</sup> M. B. Mickaliger,<sup>7</sup> D. Perrodin,<sup>14</sup>

S. A. Sanidas,<sup>7</sup> B. W. Stappers,<sup>7</sup> T. Sprenger,<sup>1</sup> O. Wucknitz,<sup>1</sup> C. G. Bassa,<sup>15</sup> M. Burgay,<sup>14</sup> R. Concu,<sup>14</sup>

M. Gaikwad,<sup>1</sup> G. H. Janssen,<sup>15,16</sup> K. J. Lee,<sup>4</sup> A. Melis,<sup>14</sup> M. Pilia,<sup>14</sup> A. Possenti,<sup>14,17</sup> L. Wang<sup>4</sup>

and W. W. Zhu<sup>18</sup>

*Affiliations are listed at the end of the paper*

Accepted 2023 June 26. Received 2023 June 18; in original form 2023 February 16

## ABSTRACT

We present the first large sample of scintillation arcs in millisecond pulsars (MSPs), analysing 12 sources observed with the Large European Array for Pulsars (LEAP), and the Effelsberg 100-m telescope. We estimate the delays from multipath propagation, measuring significant correlated changes in scattering time-scales over a 10 yr time span. Many sources show compact concentrations of power in the secondary spectrum, which in PSRs J0613−0200 and J1600−3053 can be tracked between observations, and are consistent with compact scattering at fixed angular positions. Other sources such as PSRs J1643−1224 and J0621+1002 show diffuse, asymmetric arcs which are likely related to phase-gradients across the scattering screen. PSR B1937+21 shows at least three distinct screens which dominate at different times and evidence of varying screen axes or multiscreen interactions. We model annual and orbital arc curvature variations in PSR J0613−0200, providing a measurement of the longitude of ascending node, resolving the sense of the orbital inclination, where our best-fit model is of a screen with variable axis of anisotropy over time, corresponding to changes in the scattering of the source. Unmodelled variations of the screen’s axis of anisotropy are likely to be a limiting factor in determining orbital parameters with scintillation, requiring careful consideration of variable screen properties, or independent very long baseline interferometry (VLBI) measurements. Long-term scintillation studies such as this serve as a complementary tool to pulsar timing, to measure a source of correlated noise for pulsar timing arrays, solve pulsar orbits, and to understand the astrophysical origin of scattering screens.

**Key words:** pulsars: general – ISM: general.

## 1 INTRODUCTION

Pulsar timing arrays (PTAs) involve timing an ensemble of millisecond pulsars (MSPs) at different sky positions to detect nHz gravitational waves (GWs) from coalescing supermassive black holes. Recently, PTAs have detected a common red noise signal, which is a time-correlated signal of similar amplitude and spectrum shared among pulsars in the array (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov et al. 2021b; Antoniadis et al. 2022). While it is possible that a gravitational wave signature is responsible for this effect, there is yet no detection of a spatial correlation that would be a smoking-gun of the gravitational wave background (Hellings & Downs 1983).

The ionized interstellar medium (IISM) is one of the largest contributors of correlated noise to PTAs (see Verbiest & Shaifullah 2018 for a review), and understanding all of its effects is crucial, especially as a GW detection may be imminent. The total column density of electrons induces a $\lambda^2$ dispersive delay (where $\lambda$ is the observing wavelength), where variations are seen prominently in

low-frequency observations (Donner et al. 2020; Tarafdar et al. 2022). Spatial variations in the electron column density results in multipath propagation, resulting in delays scaling as $\lambda^\alpha$, with $\alpha \sim 4.0 \pm 0.6$ (Oswald et al. 2021). In the time domain, this effect can be seen through the broadening of pulses by a characteristic scattering tail [e.g. Bhat et al. 2004, although with the additional complication that different scattered paths can encounter a different electron column (Cordes, Shannon & Stinebring 2016; Donner et al. 2019)]. In the Fourier domain, this is observed as scintillation, where temporal and spectral variations of flux density arise from interference between deflected, coherent images of the pulsar.

Pulsar scintillation is now commonly studied through the secondary spectrum, which is the 2D power spectrum of the scintillation pattern. In this space, many sources have been seen to have ‘scintillation arcs’ (Stinebring et al. 2001), parabolic distributions of power which indicate scattering being dominated by highly localized regions, or ‘thin screens’ (Walker et al. 2004; Cordes et al. 2006). Furthermore, the presence of sharp inverted parabolic ‘arclets’ stemming from the main parabola are seen in some cases, implying strong anisotropy (seen in $\sim 20$ per cent of sources in Stinebring et al. 2022).

<footnote>* E-mail: ramain@mpifr-bonn.mpg.de</footnote>
The majority of the brightest known radio pulsars are isolated and slowly rotating (Manchester et al. 2005), which due to high signal-to-noise (S/N) requirements have been the focus of the widest-ranging studies of scintillation arcs to date (Stinebring et al. 2022; Wu et al. 2022; Main et al. 2023). However, recent studies have begun to show the power of studying scintillation arcs in MSPs. For the precision timing of PTA pulsars, scattering variations may be a source of uncorrected correlated noise (Goncharov et al. 2021a; Chalumeau et al. 2022), which can be estimated through scintillation arcs (Hemberger & Stinebring 2008), or through the frequency scale of scintillation (for applications to PTA data, see e.g. Levin et al. 2016; Liu et al. 2022). Additionally, scintillation arcs encode the relative velocity and distance of the pulsar, scattering screen, and the Earth, so modelling of their annual and orbital variations can be used to precisely measure pulsar orbital parameters as well as screen distances (Reardon et al. 2020; McKee et al. 2022; Walker et al. 2022).

The Large European Pulsar Array (LEAP) is a 195-m tied-array beam telescope comprised of many of the largest telescopes in Europe, and has been observing >20 MSPs at monthly cadence since 2012 as part of the European Pulsar Timing Array (EPTA; Stappers & Kramer 2011; Kramer & Champion 2013; Bassa et al. 2016). Owing to its sensitivity and data products which can be reduced to any time and frequency resolution, LEAP is well suited to study MSP scintillation. In studies to date, secondary spectra have been used to measure the total time delays from scattering in PSR J0613−0200 (Main et al. 2020), and to associate the scattering screen of PSR J1643−1224 with a known H II region (Mall et al. 2022).

In this paper, we present the first large sample of scintillation arcs of MSPs, observed over the last 10 yr with LEAP, and the Effelsberg 100-m telescope. The paper is organized as follows: in Section 2, we revisit the necessary theory of scintillation arcs for this work, in Section 3, we describe our observations and data, and in Section 4, we describe our analysis. We discuss the results of particular pulsars in Section 5, and the whole sample in Section 6. Section 7 contains the ramifications of our findings and prospects for the future.

## 2 BACKGROUND OF SCINTILLATION

### 2.1 Arc curvature and scintillation velocity

Here, we briefly review the relevant theory of scintillation arcs, which we detailed in Main et al. (2020) (originally developed, and explained in more detail in Walker et al. 2004; Cordes et al. 2006).

The dynamic spectrum $I(t, v)$ is the measured flux density as a function of time and frequency (typically averaged over many pulses), showing variations owing to interstellar scintillation. In the 2D power spectrum of the dynamic spectrum $S(f_t, f_v) = |\tilde{I}(f_t, f_v)|^2$, referred to as the 'secondary spectrum', the conjugate variable of time $f_t = -f_D$ is related to the Doppler shift between deflected paths, and depends on the angles $\theta$ of two deflected paths (ij) as

$$f_{D,ij} = - \frac{(\theta_i - \theta_j) \cdot \mathbf{v}_{\text{eff}}}{\lambda}, \qquad (1)$$

and the conjugate variable of frequency $f_v \equiv \tau$ is related to the delay between image pairs, described as

$$\tau_{ij} = \frac{d_{\text{eff}}(\theta_i^2 - \theta_j^2)}{2c}. \qquad (2)$$

The effective distance $d_{\text{eff}}$ and effective velocity $\mathbf{v}_{\text{eff}}$ depend on the relative distances and velocities of the pulsar ($d_{\text{psr}}$, $\mathbf{v}_{\text{psr}}$), screen ($d_{\text{scr}}$, $\mathbf{v}_{\text{scr}}$), and Earth ($\mathbf{v}_{\oplus}$), as

$$d_{\text{eff}} = (1/s - 1)d_{\text{psr}}, \qquad (3)$$

$$\mathbf{v}_{\text{eff}} = (1/s - 1)\mathbf{v}_{\text{psr}} + \mathbf{v}_{\oplus} - \mathbf{v}_{\text{scr}}/s, \qquad (4)$$

where the fractional screen distance from the pulsar is defined as $s \equiv 1 - d_{\text{scr}}/d_{\text{psr}}$.

When one of the two scattered angles is 0 (i.e. the theoretical undeflected line-of-sight), the common dependence of $f_D$ and $\tau$ on the observed angle $\theta$ results in a parabolic distribution of power

$$\tau = \eta f_D^2. \qquad (5)$$

The proportionality constant, or 'arc-curvature' $\eta$ depends on the relative distances and velocities of the pulsar, screen, and Earth, as

$$\eta = d_{\text{eff}} \lambda^2 / (2c \mathbf{v}_{\text{eff}}^2 \cos(\alpha)^2), \qquad (6)$$

where $\alpha$ is the angle between $\mathbf{v}_{\text{eff}}$ and $\psi$.

Throughout this paper, we work with the distance weighted effective velocity $W$, which rearranges $\eta$ to separate the unknown values, and is proportional to $|v_{\text{eff}}, \parallel|$,

$$W \equiv \frac{|v_{\text{eff},\parallel}|}{\sqrt{d_{\text{eff}}}} = \frac{\lambda}{\sqrt{2c\eta}}. \qquad (7)$$

This is the same approach as in Main et al. (2023), and proportional to the quantity used in Mall et al. (2022); the main benefits of $W$ is that it does not diverge as $|v_{\text{eff}}, \parallel| \to 0$ and is independent of the observing frequency. In the absence of arcs, a characteristic curvature $\eta$ can be estimated from the time and frequency scale from the 2D autocorrelation function (ACF) of scintillation. The scintillation bandwidth $v_s$ is defined as the half-width at half-maximum of the ACF in frequency, and the scintillation time-scale $t_s$ is defined as the $1/e$ point of the ACF in time (Cordes 1986).

Using thin screen relations (with a phase structure function with index of 2, details in Cordes & Rickett 1998), the corresponding time delay and Doppler shifts can be inferred from the scintillation bandwidth and time-scale as

$$\tau_s \approx 1/2\pi v_s, \quad f_{D,s} \approx 1/\sqrt{2\pi t_s}. \qquad (8)$$

These relations are approximate, as the prefactors are model dependent. Then we can estimate $W$ as

$$W \approx \frac{\lambda}{t_s} \sqrt{\frac{v_s}{2c\pi}}. \qquad (9)$$

A measure of either the arc curvature $\eta$, or the scintillation time-scale and bandwidth, are then a direct measure of $v_{\text{eff},\parallel}/\sqrt{d_{\text{eff}}}$. In this paper, we measure $W$ from the secondary spectrum wherever possible; scintillation velocities derived from $t_s$ and $v_s$ are dependent on the distribution of power across the screen, and will systematically vary as different regions of the same scattering screen are seen (Cordes & Rickett 1998; Rickett et al. 2014; Reardon et al. 2019). Measurements of scintillation arc curvatures are more robust, demonstrated to be stable to changes in screen's anisotropy (Reardon et al. 2020), or in the presence of significant substructure (Sprenger et al. 2022).

### 2.2 Time-scale for feature movement

Secondary spectra often show compact features at fixed angular positions, which over time are seen to travel through the secondary spectrum due to their relative velocity (Hill et al. 2005; Sprenger et al. 2022). An important quantity is the time-scale for features to pass through the secondary spectrum, which sets a time-scale for
The majority of the brightest known radio pulsars are isolated and slowly rotating (Manchester et al. 2005), which due to high signal-to-noise (S/N) requirements have been the focus of the widest-ranging studies of scintillation arcs to date (Stinebring et al. 2022; Wu et al. 2022; Main et al. 2023). However, recent studies have begun to show the power of studying scintillation arcs in MSPs. For the precision timing of PTA pulsars, scattering variations may be a source of uncorrected correlated noise (Goncharov et al. 2021a; Chalumeau et al. 2022), which can be estimated through scintillation arcs (Hemberger & Stinebring 2008), or through the frequency scale of scintillation (for applications to PTA data, see e.g. Levin et al. 2016; Liu et al. 2022). Additionally, scintillation arcs encode the relative velocity and distance of the pulsar, scattering screen, and the Earth, so modelling of their annual and orbital variations can be used to precisely measure pulsar orbital parameters as well as screen distances (Reardon et al. 2020; McKee et al. 2022; Walker et al. 2022).

The Large European Pulsar Array (LEAP) is a 195-m tied-array beam telescope comprised of many of the largest telescopes in Europe, and has been observing >20 MSPs at monthly cadence since 2012 as part of the European Pulsar Timing Array (EPTA; Stappers & Kramer 2011; Kramer & Champion 2013; Bassa et al. 2016). Owing to its sensitivity and data products which can be reduced to any time and frequency resolution, LEAP is well suited to study MSP scintillation. In studies to date, secondary spectra have been used to measure the total time delays from scattering in PSR J0613−0200 (Main et al. 2020), and to associate the scattering screen of PSR J1643−1224 with a known H II region (Mall et al. 2022).

In this paper, we present the first large sample of scintillation arcs of MSPs, observed over the last 10 yr with LEAP, and the Effelsberg 100-m telescope. The paper is organized as follows: in Section 2, we revisit the necessary theory of scintillation arcs for this work, in Section 3, we describe our observations and data, and in Section 4, we describe our analysis. We discuss the results of particular pulsars in Section 5, and the whole sample in Section 6. Section 7 contains the ramifications of our findings and prospects for the future.

## 2 BACKGROUND OF SCINTILLATION

### 2.1 Arc curvature and scintillation velocity

Here, we briefly review the relevant theory of scintillation arcs, which we detailed in Main et al. (2020) (originally developed, and explained in more detail in Walker et al. 2004; Cordes et al. 2006).

The dynamic spectrum $I(t, v)$ is the measured flux density as a function of time and frequency (typically averaged over many pulses), showing variations owing to interstellar scintillation. In the 2D power spectrum of the dynamic spectrum $S(f_t, f_v) = |\tilde{I}(f_t, f_v)|^2$, referred to as the 'secondary spectrum', the conjugate variable of time $f_t = -f_D$ is related to the Doppler shift between deflected paths, and depends on the angles $\theta$ of two deflected paths (ij) as

$$f_{D,ij} = - \frac{(\theta_i - \theta_j) \cdot \mathbf{v}_{\text{eff}}}{\lambda}, \qquad (1)$$

and the conjugate variable of frequency $f_v \equiv \tau$ is related to the delay between image pairs, described as

$$\tau_{ij} = \frac{d_{\text{eff}}(\theta_i^2 - \theta_j^2)}{2c}. \qquad (2)$$

The effective distance $d_{\text{eff}}$ and effective velocity $\mathbf{v}_{\text{eff}}$ depend on the relative distances and velocities of the pulsar ($d_{\text{psr}}$, $\mathbf{v}_{\text{psr}}$), screen ($d_{\text{scr}}$, $\mathbf{v}_{\text{scr}}$), and Earth ($\mathbf{v}_{\oplus}$), as

$$d_{\text{eff}} = (1/s - 1)d_{\text{psr}}, \qquad (3)$$

$$\mathbf{v}_{\text{eff}} = (1/s - 1)\mathbf{v}_{\text{psr}} + \mathbf{v}_{\oplus} - \mathbf{v}_{\text{scr}}/s, \qquad (4)$$

where the fractional screen distance from the pulsar is defined as $s \equiv 1 - d_{\text{scr}}/d_{\text{psr}}$.

When one of the two scattered angles is 0 (i.e. the theoretical undeflected line-of-sight), the common dependence of $f_D$ and $\tau$ on the observed angle $\theta$ results in a parabolic distribution of power

$$\tau = \eta f_D^2. \qquad (5)$$

The proportionality constant, or 'arc-curvature' $\eta$ depends on the relative distances and velocities of the pulsar, screen, and Earth, as

$$\eta = d_{\text{eff}} \lambda^2 / (2c \mathbf{v}_{\text{eff}}^2 \cos(\alpha)^2), \qquad (6)$$

where $\alpha$ is the angle between $\mathbf{v}_{\text{eff}}$ and $\psi$.

Throughout this paper, we work with the distance weighted effective velocity $W$, which rearranges $\eta$ to separate the unknown values, and is proportional to $|v_{\text{eff}}, \parallel|$,

$$W \equiv \frac{|v_{\text{eff},\parallel}|}{\sqrt{d_{\text{eff}}}} = \frac{\lambda}{\sqrt{2c\eta}}. \qquad (7)$$

This is the same approach as in Main et al. (2023), and proportional to the quantity used in Mall et al. (2022); the main benefits of $W$ is that it does not diverge as $|v_{\text{eff}}, \parallel| \to 0$ and is independent of the observing frequency. In the absence of arcs, a characteristic curvature $\eta$ can be estimated from the time and frequency scale from the 2D autocorrelation function (ACF) of scintillation. The scintillation bandwidth $v_s$ is defined as the half-width at half-maximum of the ACF in frequency, and the scintillation time-scale $t_s$ is defined as the $1/e$ point of the ACF in time (Cordes 1986).

Using thin screen relations (with a phase structure function with index of 2, details in Cordes & Rickett 1998), the corresponding time delay and Doppler shifts can be inferred from the scintillation bandwidth and time-scale as

$$\tau_s \approx 1/2\pi v_s, \qquad f_{D,s} \approx 1/\sqrt{2\pi t_s}. \qquad (8)$$

These relations are approximate, as the prefactors are model dependent. Then we can estimate $W$ as

$$W \approx \frac{\lambda}{t_s} \sqrt{\frac{v_s}{2c\pi}}. \qquad (9)$$

A measure of either the arc curvature $\eta$, or the scintillation time-scale and bandwidth, are then a direct measure of $v_{\text{eff},\parallel}/\sqrt{d_{\text{eff}}}$. In this paper, we measure $W$ from the secondary spectrum wherever possible; scintillation velocities derived from $t_s$ and $v_s$ are dependent on the distribution of power across the screen, and will systematically vary as different regions of the same scattering screen are seen (Cordes & Rickett 1998; Rickett et al. 2014; Reardon et al. 2019). Measurements of scintillation arc curvatures are more robust, demonstrated to be stable to changes in screen's anisotropy (Reardon et al. 2020), or in the presence of significant substructure (Sprenger et al. 2022).

### 2.2 Time-scale for feature movement

Secondary spectra often show compact features at fixed angular positions, which over time are seen to travel through the secondary spectrum due to their relative velocity (Hill et al. 2005; Sprenger et al. 2022). An important quantity is the time-scale for features to pass through the secondary spectrum, which sets a time-scale for
scattering delays to correlate. This time-scale (effectively the same as the traditional 'refractive time-scale' of scintillation) of the screen is the time it takes to traverse to a new section of the screen. For an observable portion of the screen with angle $\theta_0$ with corresponding size length $l_r$,

$$l_r = 2\theta_0 d_{scr} = \sqrt{\frac{8c\tau_0}{d_{eff}}} d_{scr}, \qquad (10)$$

then the time-scale is related to the arc curvature and the maximum delay $\tau_0$ as

$$t_r \approx \frac{l_r}{s v_{eff}} = \frac{\sqrt{8c\tau_0}}{W}. \qquad (11)$$

In Sprenger et al. (2022), all of the observables were expressed as a function of the feature movement in $\sqrt{\tau}$, related to the distance weighted effective velocity as

$$\partial_t \sqrt{\tau} = \frac{W}{\sqrt{2c}} = \frac{1}{2\nu\sqrt{\eta}}, \qquad (12)$$

which can also be used to obtain equation (10). The relevant $\tau_0$ value can either be the highest $\tau$ that features are visible, which gives the time-scale for features to completely pass through the secondary spectrum, or the $1/e$ value of $\tau$.

# 3 OBSERVATIONS AND SAMPLE

## 3.1 LEAP

LEAP is a tied-array telescope comprised of the Effelsberg 100-m telescope, the Lovell telescope at Jodrell Bank Observatory, the Westerbork Synthesis Radio Telescope, the Nançay Radio Telescope, and the Sardinia Radio Telescope, simultaneously observing MSPs with monthly cadence. Each observing run has some subset of these telescopes; the voltage data are recorded at each site, then shipped to Jodrell Bank to be correlated and coherently added using the pipeline developed by the LEAP team (details of LEAP in Bassa et al. 2016, and details of the correlator in Smits et al. 2017). When all telescopes are available, this results in an effective 195-m diameter dish. The scans are typically 30–60 min. The data are recorded in contiguous 16-MHz sub-bands covering 1332–1460 MHz, where the total bandwidth per observation varies between 80–128 MHz depending on the telescopes used; Jodrell Bank and Sardinia never use the full 128-MHz bandwidth. While the standard folding pipeline produces archives with 10 s integrations, 1 MHz channels, the coherently added voltages are stored on tape, allowing us to later reduce the data with much higher spectral resolution to study fine scintillation features.

## 3.2 Effelsberg

### 3.2.1 Long targeted observations

Several of the sources have scintillation which is barely resolved in a ~30–60 min observation, and with scintiles comparable to LEAP’s 128 MHz band. This motivated us to take tailored observations with a wider bandwidth and a longer duration. In addition to the LEAP observations, we obtained 2–3 hr observations of several LEAP sources where the scintillation was not quite resolved in frequency or time. The data were taken using the PSRIX backend (Lazarus et al. 2016), using the central feed of the seven beam receiver, recorded in 25 MHz sub-bands with a

usable bandwidth of 1250–1450 MHz and saved in PSRDADA<fnref n="1" /> format.

### 3.2.2 EPTA observations

Since 2021 March, Effelsberg observations for the EPTA record a separate parallel data stream suitable for scintillation studies, using the Effelsberg Direct Digitization system, folding with 64 000 channels across 1200–1600 MHz (sensitive to delays of $\tau < 80 \mu s$). These are now a regular data product in EPTA observations, and along with LEAP will double the cadence of observations suitable for studying scintillation arcs.

## 3.3 Sample

In our sample, there are six sources with resolvable, variable arcs within the 30–60 min LEAP observations: PSRs J0613–0200, J0621+1002, J1600–3053, J1643–1224, J1918–0642, J1939+2134 (B1937+21). In these sources, we study the time variability of the arc curvature, and time delays. A scintillation arc is also faintly seen in the highest S/N observation of PSR J1824–2452A (B1821–24A).

In addition, we also investigate sources where scintillation is marginally resolved by LEAP, supplementing this study with longer 2–3 hr Effelsberg observations. These include single scans on PSRs J0751+1807, J1713+0747, J1832–0836, J1857+0943 (B1855+09), J2010–1323, as well as a high-cadence approximately bi-weekly campaign on PSR J0613–0200 from 2020 March–June, totalling 19 observations, which was previously included in Main et al. (2020). The first year of EPTA fine-channel scintillation products are included for PSRs J0613–0200 and J1643–1224, with 12 and 9 observations respectively, to demonstrate the value of these data products to increase the cadence and time span of scintillation data products. A summary of the samples is given in Table 1. Our sample partially overlaps with the scintillation study of EPTA pulsars in Liu et al. (2022), where the larger bandwidth, but lower frequency resolution of the Nançay data set allowed for studies of slightly less scattered sources.

# 4 SCINTILLATION DATA PRODUCTS AND MEASUREMENTS

## 4.1 Dynamic and secondary spectra

The creation of dynamic and secondary spectra are almost identical to the methods described in Main et al. (2020); Mall et al. (2022), however we briefly review and describe differences here.

Data were folded using the `dspsr` software package (van Straten & Bailes 2011), beginning from the base band data. The time, phase, and frequency bins were different for each specific pulsar, chosen to fully resolve the scintillation in time and frequency, or equivalently, to extend sufficiently far in $f_D$ and $\tau$ to capture the full extent of arcs in the secondary spectrum. The sub-bands were combined in frequency using `psradd` in the `psrarchive` software package (Hotan, van Straten & Manchester 2004), and polarizations were summed to form total intensity. Radio frequency interference (RFI) flagging and treatment of masked pixels is identical to Main et al. (2020). The outer 10 per cent of the dynamic spectra are tapered with a Hanning window, before forming the secondary spectrum.

<footnote n="1">¹http://psrdada.sourceforge.net/</footnote>
scattering delays to correlate. This time-scale (effectively the same as the traditional 'refractive time-scale' of scintillation) of the screen is the time it takes to traverse to a new section of the screen. For an observable portion of the screen with angle $\theta_0$ with corresponding size length $l_r$,

$$l_r = 2\theta_0 d_{scr} = \sqrt{\frac{8c\tau_0}{d_{eff}}} d_{scr}, \qquad (10)$$

then the time-scale is related to the arc curvature and the maximum delay $\tau_0$ as

$$t_r \approx \frac{l_r}{s v_{eff}} = \frac{\sqrt{8c\tau_0}}{W}. \qquad (11)$$

In Sprenger et al. (2022), all of the observables were expressed as a function of the feature movement in $\sqrt{\tau}$, related to the distance weighted effective velocity as

$$\partial_t \sqrt{\tau} = \frac{W}{\sqrt{2c}} = \frac{1}{2\nu\sqrt{\eta}}, \qquad (12)$$

which can also be used to obtain equation (10). The relevant $\tau_0$ value can either be the highest $\tau$ that features are visible, which gives the time-scale for features to completely pass through the secondary spectrum, or the $1/e$ value of $\tau$.

# 3 OBSERVATIONS AND SAMPLE

## 3.1 LEAP

LEAP is a tied-array telescope comprised of the Effelsberg 100-m telescope, the Lovell telescope at Jodrell Bank Observatory, the Westerbork Synthesis Radio Telescope, the Nançay Radio Telescope, and the Sardinia Radio Telescope, simultaneously observing MSPs with monthly cadence. Each observing run has some subset of these telescopes; the voltage data are recorded at each site, then shipped to Jodrell Bank to be correlated and coherently added using the pipeline developed by the LEAP team (details of LEAP in Bassa et al. 2016, and details of the correlator in Smits et al. 2017). When all telescopes are available, this results in an effective 195-m diameter dish. The scans are typically 30–60 min. The data are recorded in contiguous 16-MHz sub-bands covering 1332–1460 MHz, where the total bandwidth per observation varies between 80–128 MHz depending on the telescopes used; Jodrell Bank and Sardinia never use the full 128-MHz bandwidth. While the standard folding pipeline produces archives with 10 s integrations, 1 MHz channels, the coherently added voltages are stored on tape, allowing us to later reduce the data with much higher spectral resolution to study fine scintillation features.

## 3.2 Effelsberg

### 3.2.1 Long targeted observations

Several of the sources have scintillation which is barely resolved in a ~30–60 min observation, and with scintiles comparable to LEAP’s 128 MHz band. This motivated us to take tailored observations with a wider bandwidth and a longer duration. In addition to the LEAP observations, we obtained 2–3 hr observations of several LEAP sources where the scintillation was not quite resolved in frequency or time. The data were taken using the PSRIX backend (Lazarus et al. 2016), using the central feed of the seven beam receiver, recorded in 25 MHz sub-bands with a

usable bandwidth of 1250–1450 MHz and saved in PSRDADA<fnref n="1" /> format.

### 3.2.2 EPTA observations

Since 2021 March, Effelsberg observations for the EPTA record a separate parallel data stream suitable for scintillation studies, using the Effelsberg Direct Digitization system, folding with 64 000 channels across 1200–1600 MHz (sensitive to delays of $\tau < 80 \mu s$). These are now a regular data product in EPTA observations, and along with LEAP will double the cadence of observations suitable for studying scintillation arcs.

## 3.3 Sample

In our sample, there are six sources with resolvable, variable arcs within the 30–60 min LEAP observations: PSRs J0613–0200, J0621+1002, J1600–3053, J1643–1224, J1918–0642, J1939+2134 (B1937+21). In these sources, we study the time variability of the arc curvature, and time delays. A scintillation arc is also faintly seen in the highest S/N observation of PSR J1824–2452A (B1821–24A).

In addition, we also investigate sources where scintillation is marginally resolved by LEAP, supplementing this study with longer 2–3 hr Effelsberg observations. These include single scans on PSRs J0751+1807, J1713+0747, J1832–0836, J1857+0943 (B1855+09), J2010–1323, as well as a high-cadence approximately bi-weekly campaign on PSR J0613–0200 from 2020 March–June, totalling 19 observations, which was previously included in Main et al. (2020). The first year of EPTA fine-channel scintillation products are included for PSRs J0613–0200 and J1643–1224, with 12 and 9 observations respectively, to demonstrate the value of these data products to increase the cadence and time span of scintillation data products. A summary of the samples is given in Table 1. Our sample partially overlaps with the scintillation study of EPTA pulsars in Liu et al. (2022), where the larger bandwidth, but lower frequency resolution of the Nançay data set allowed for studies of slightly less scattered sources.

# 4 SCINTILLATION DATA PRODUCTS AND MEASUREMENTS

## 4.1 Dynamic and secondary spectra

The creation of dynamic and secondary spectra are almost identical to the methods described in Main et al. (2020); Mall et al. (2022), however we briefly review and describe differences here.

Data were folded using the `dspsr` software package (van Straten & Bailes 2011), beginning from the base band data. The time, phase, and frequency bins were different for each specific pulsar, chosen to fully resolve the scintillation in time and frequency, or equivalently, to extend sufficiently far in $f_D$ and $\tau$ to capture the full extent of arcs in the secondary spectrum. The sub-bands were combined in frequency using `psradd` in the `psrarchive` software package (Hotan, van Straten & Manchester 2004), and polarizations were summed to form total intensity. Radio frequency interference (RFI) flagging and treatment of masked pixels is identical to Main et al. (2020). The outer 10 per cent of the dynamic spectra are tapered with a Hanning window, before forming the secondary spectrum.

<footnote n="1">¹http://psrdada.sourceforge.net/</footnote>
<table><caption>**Table 1.** Summary of the pulsars in our sample. The pulsar distances shown with 1σ errorbars are VLBI measurements using timing constraints as priors from table 6 of Ding et al. (2023) (marked with †), or, where not available, timing parallax measurements compiled by the IPTA in Perera et al. (2019). Distances shown without errors are DM distance estimates from YMW16 model (Yao, Manchester & Wang 2017). The DM values are taken from PSRCAT (Manchester et al. 2005). While LEAP has heightened sensitivity, the targeted Effelsberg observations have a larger bandwidth, and longer durations as described in Section 3.</caption><thead><tr><th>Pulsar name</th><th>DM<br>(pc cm<sup>-3</sup>)</th><th>P<sub>b</sub><br>(d)</th><th>d<sub>psr</sub><br>(kpc)</th><th>μ<sub>α</sub><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub><br>(mas yr<sup>-1</sup>)</th><th>⟨S/N⟩</th><th>⟨t<sub>obs</sub>⟩<br>(min)</th><th>ν<sub>chan</sub><br>(kHz)</th><th>N<sub>phase</sub></th><th>Res.<br>Scint</th></tr></thead><tbody><tr><td>J0613-0200</td><td>38.78</td><td>1.20</td><td>1.11 ± 0.05</td><td>1.828(5)</td><td>-10.35(1)</td><td>134</td><td>51</td><td>62.5</td><td>128</td><td>LEAP</td></tr><tr><td>J0621+1002</td><td>36.47</td><td>8.32</td><td>1.6<sup>+0.5</sup><sub>-0.3</sub>†</td><td>3.2(1)</td><td>0.6(5)</td><td>71</td><td>39</td><td>62.5</td><td>128</td><td>LEAP</td></tr><tr><td>J0751+1807</td><td>30.25</td><td>0.26</td><td>1.4<sup>+0.4</sup><sub>-0.3</sub></td><td>-2.72(6)</td><td>-13.4(3)</td><td>140</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>J1600-3053</td><td>52.33</td><td>14.35</td><td>2.0<sup>+0.3</sup><sub>-0.2</sub></td><td>-0.97(1)</td><td>-7.04(5)</td><td>139</td><td>54</td><td>31.125</td><td>64</td><td>LEAP</td></tr><tr><td>J1643-1224</td><td>62.41</td><td>147.02</td><td>0.95<sup>+0.15</sup><sub>-0.11</sub>†</td><td>6.03(3)</td><td>4.1(1)</td><td>231</td><td>33</td><td>7.8125</td><td>16</td><td>LEAP</td></tr><tr><td>J1713+0747</td><td>15.92</td><td>67.83</td><td>1.20 ± 0.03</td><td>4.924(1)</td><td>-3.913(2)</td><td>408</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>B1821-24A</td><td>119.89</td><td>–</td><td>3.7</td><td>-0.2(2)</td><td>-6(4)</td><td>132</td><td>44</td><td>31.125</td><td>64</td><td>LEAP</td></tr><tr><td>J1832-0836</td><td>28.19</td><td>–</td><td>0.8</td><td>-7.97(5)</td><td>-21.2(2)</td><td>55</td><td>120</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>B1855+09</td><td>13.31</td><td>12.33</td><td>1.1 ± 0.1</td><td>-2.652(4)</td><td>-5.423(6)</td><td>138</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>J1918-0642</td><td>26.46</td><td>10.91</td><td>1.48<sup>+0.19</sup><sub>-0.14</sub>†</td><td>-7.15(2)</td><td>-5.94(5)</td><td>53</td><td>30</td><td>125</td><td>128</td><td>LEAP</td></tr><tr><td>B1937+21</td><td>71.02</td><td>–</td><td>2.9<sup>+0.3</sup><sub>-0.2</sub>†</td><td>0.074(2)</td><td>-0.410(3)</td><td>485</td><td>45</td><td>31.125</td><td>32</td><td>LEAP</td></tr><tr><td>J2010-1323</td><td>22.18</td><td>–</td><td>1.9<sup>+0.8</sup><sub>-0.5</sub></td><td>2.56(6)</td><td>-5.9(2)</td><td>82</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr></tbody></table>

As the Effelsberg observations cover a larger bandwidth, we use the 'NuT transform' instead of a direct Fast Fourier Transform in time to form secondary spectra, as described in Sprenger et al. (2021). This transformation is a direct Fourier Transform over a scaled time axis of $t' = t\nu/\nu_{\text{ref}}$ in every channel with frequency $\nu$. This prevents smearing of scintillation arcs owing to the $\eta \propto \lambda^2$ dependence, and ensures the contribution to scintillation from sources of fixed angular position are at fixed position in the secondary spectrum of a single observation. We apply this in all of our observations, referencing to $\nu_{\text{ref}} = 1400$ MHz. We note that this transform can lead to some artefacts in the secondary spectrum, as power on the $\tau = 0$ axis (arising from e.g. RFI, pulse-to-pulse variations) is spread to diagonal lines. We do not see these artefacts prominently in our observations, as they are most prominent in observations with long-durations, large fractional bandwidths, and large pulse-to-pulse flux variations.

Representative dynamic spectra of the sources in our sample are shown in Fig. 1, and their associated secondary spectra are shown in Fig. 2. The secondary spectra of every observation used in this work are included in the Appendix.

## 4.2 Arc curvature measurement

We have measured the arc curvature using the 'normalized secondary spectrum' as done in Reardon et al. (2020); Walker et al. (2022), in which the $f_D$ axis of the secondary spectrum is mapped to $f_{D,\text{norm}} = f_D\sqrt{\tau_{\text{ref}}/\bar{\tau}}$ (where we set the arbitrary reference time delay $\tau_{\text{ref}} = \tau_{\text{max}}$ throughout). This transformation effectively stretches the $f_D$ axis, mapping parabolas to vertical lines of constant $f_{D,\text{norm}}$. Then $W$ can be identified by finding peaks in $S(f_{D,\text{norm}})$ after performing a weighted sum over $\tau$.

Arcs blend together at low values of $\tau$, becoming more clearly demarcated at high values of $\tau$. As the optimal range of $\tau$ to sum over and $f_{D,\text{norm}}$ to fit vary between pulsars and between epochs, we applied the arc curvature fitting algorithm interactively. Alongside the fitting, the dynamic and secondary spectrum of each observation are verified for corruption by RFI or phasing artefacts. The range in $\tau$ to sum to form $S(f_{D,\text{norm}})$ is adjusted per source, and peaks in $f_{D,\text{norm}}$ can be given as initial guesses for a least-squares fit of a parabola to the local region of $S(f_{D,\text{norm}})$, and the arc curvature is given as $\eta = f_{D,\text{norm}}^2/\tau_{\text{ref}}$. An example is shown in Fig. 3.

## 4.3 Time delays from secondary spectra

Secondary spectra express the power of scintillation in terms of the conjugate variables of time and frequency respectively, $f_D$ and $\tau$. In the strong scattering regime, the secondary spectrum contains contributions from all pairs of interfering images. Hemberger & Stinebring (2008) showed how one can estimate the averaged time delay $\langle\tau\rangle$ from the secondary spectrum, and in Main et al. (2020) it was argued that this technique is valid in the limit of a strong, central image arising from a single thin screen, or when the response function is close to an exponential. The total geometric time delay $\langle\tau\rangle$ can then be estimated as

$$ \langle\tau\rangle = \frac{\int_{0}^{T} \tau |I(\tau)|^2 d\tau}{\int_{0}^{T} |I(\tau)|^2 d\tau}. \qquad (13) $$

For each source, a range in $f_D$ was chosen to fully encompass the power of the visible scintillation arc, and the background noise was estimated from a region of the same size offset by $f_D \pm 30$ mHz. The integrated profile of $\langle\tau\rangle$ against maximum $\tau_{\text{max}} \equiv T$ was examined in all cases; the value and the error on $\langle\tau\rangle$ were estimated by the mean and standard deviation of the profile once it plateaus, taken in the range of $3T/4 < \tau < T$. Examples of secondary spectra and their associated profiles of $I(\tau)$ are shown in Fig. 4.

## 4.4 Scintillation parameters from the ACF

In some of the sources with diffuse arcs, it is difficult to determine parameters from the secondary spectra. In these cases, we measure the scintillation time-scale and bandwidth in a more traditional way, through the 2D ACF of the dynamic spectrum $R(\Delta\nu) = (I*I)(\Delta\nu)$. The ACF in frequency is fit with a Lorentzian, and with a Gaussian in time (Cordes 1986). The scintillation bandwidth is inversely proportional to the bulk scattering delay as $\langle\tau\rangle = C/2\pi\nu_s$, where the model-dependent constant $C$, depending on the distribution of scattered power, is assumed to be 1, which is the value for a thin screen with a square-law phase structure function (Cordes & Rickett 1998). The 2D ACFs of sources without clear arcs is shown in Fig. 5.

The precision on the measurements on the ACF is high, and the true error is dominated by the fact that there are finite scintils within
<table><caption>**Table 1.** Summary of the pulsars in our sample. The pulsar distances shown with 1σ errorbars are VLBI measurements using timing constraints as priors from table 6 of Ding et al. (2023) (marked with †), or, where not available, timing parallax measurements compiled by the IPTA in Perera et al. (2019). Distances shown without errors are DM distance estimates from YMW16 model (Yao, Manchester & Wang 2017). The DM values are taken from PSRCAT (Manchester et al. 2005). While LEAP has heightened sensitivity, the targeted Effelsberg observations have a larger bandwidth, and longer durations as described in Section 3.</caption><thead><tr><th>Pulsar name</th><th>DM<br>(pc cm<sup>-3</sup>)</th><th>P<sub>b</sub><br>(d)</th><th>d<sub>psr</sub><br>(kpc)</th><th>μ<sub>α</sub><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub><br>(mas yr<sup>-1</sup>)</th><th>⟨S/N⟩</th><th>⟨t<sub>obs</sub>⟩<br>(min)</th><th>ν<sub>chan</sub><br>(kHz)</th><th>N<sub>phase</sub></th><th>Res.<br>Scint</th></tr></thead><tbody><tr><td>J0613-0200</td><td>38.78</td><td>1.20</td><td>1.11 ± 0.05</td><td>1.828(5)</td><td>-10.35(1)</td><td>134</td><td>51</td><td>62.5</td><td>128</td><td>LEAP</td></tr><tr><td>J0621+1002</td><td>36.47</td><td>8.32</td><td>1.6<sup>+0.5</sup><sub>-0.3</sub>†</td><td>3.2(1)</td><td>0.6(5)</td><td>71</td><td>39</td><td>62.5</td><td>128</td><td>LEAP</td></tr><tr><td>J0751+1807</td><td>30.25</td><td>0.26</td><td>1.4<sup>+0.4</sup><sub>-0.3</sub></td><td>-2.72(6)</td><td>-13.4(3)</td><td>140</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>J1600-3053</td><td>52.33</td><td>14.35</td><td>2.0<sup>+0.3</sup><sub>-0.2</sub></td><td>-0.97(1)</td><td>-7.04(5)</td><td>139</td><td>54</td><td>31.125</td><td>64</td><td>LEAP</td></tr><tr><td>J1643-1224</td><td>62.41</td><td>147.02</td><td>0.95<sup>+0.15</sup><sub>-0.11</sub>†</td><td>6.03(3)</td><td>4.1(1)</td><td>231</td><td>33</td><td>7.8125</td><td>16</td><td>LEAP</td></tr><tr><td>J1713+0747</td><td>15.92</td><td>67.83</td><td>1.20 ± 0.03</td><td>4.924(1)</td><td>-3.913(2)</td><td>408</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>B1821-24A</td><td>119.89</td><td>—</td><td>3.7</td><td>-0.2(2)</td><td>-6(4)</td><td>132</td><td>44</td><td>31.125</td><td>64</td><td>LEAP</td></tr><tr><td>J1832-0836</td><td>28.19</td><td>—</td><td>0.8</td><td>-7.97(5)</td><td>-21.2(2)</td><td>55</td><td>120</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>B1855+09</td><td>13.31</td><td>12.33</td><td>1.1 ± 0.1</td><td>-2.652(4)</td><td>-5.423(6)</td><td>138</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr><tr><td>J1918-0642</td><td>26.46</td><td>10.91</td><td>1.48<sup>+0.19</sup><sub>-0.14</sub>†</td><td>-7.15(2)</td><td>-5.94(5)</td><td>53</td><td>30</td><td>125</td><td>128</td><td>LEAP</td></tr><tr><td>B1937+21</td><td>71.02</td><td>—</td><td>2.9<sup>+0.3</sup><sub>-0.2</sub>†</td><td>0.074(2)</td><td>-0.410(3)</td><td>485</td><td>45</td><td>31.125</td><td>32</td><td>LEAP</td></tr><tr><td>J2010-1323</td><td>22.18</td><td>—</td><td>1.9<sup>+0.8</sup><sub>-0.5</sub></td><td>2.56(6)</td><td>-5.9(2)</td><td>82</td><td>180</td><td>100</td><td>128</td><td>Eff</td></tr></tbody></table>

As the Effelsberg observations cover a larger bandwidth, we use the 'NuT transform' instead of a direct Fast Fourier Transform in time to form secondary spectra, as described in Sprenger et al. (2021). This transformation is a direct Fourier Transform over a scaled time axis of $t' = t\nu/\nu_{\text{ref}}$ in every channel with frequency $\nu$. This prevents smearing of scintillation arcs owing to the $\eta \propto \lambda^2$ dependence, and ensures the contribution to scintillation from sources of fixed angular position are at fixed position in the secondary spectrum of a single observation. We apply this in all of our observations, referencing to $\nu_{\text{ref}} = 1400$ MHz. We note that this transform can lead to some artefacts in the secondary spectrum, as power on the $\tau = 0$ axis (arising from e.g. RFI, pulse-to-pulse variations) is spread to diagonal lines. We do not see these artefacts prominently in our observations, as they are most prominent in observations with long-durations, large fractional bandwidths, and large pulse-to-pulse flux variations.

Representative dynamic spectra of the sources in our sample are shown in Fig. 1, and their associated secondary spectra are shown in Fig. 2. The secondary spectra of every observation used in this work are included in the Appendix.

## 4.2 Arc curvature measurement

We have measured the arc curvature using the 'normalized secondary spectrum' as done in Reardon et al. (2020); Walker et al. (2022), in which the $f_D$ axis of the secondary spectrum is mapped to $f_{D,\text{norm}} = f_D\sqrt{\tau_{\text{ref}}/\bar{\tau}}$ (where we set the arbitrary reference time delay $\tau_{\text{ref}} = \tau_{\text{max}}$ throughout). This transformation effectively stretches the $f_D$ axis, mapping parabolas to vertical lines of constant $f_{D,\text{norm}}$. Then $W$ can be identified by finding peaks in $S(f_{D,\text{norm}})$ after performing a weighted sum over $\tau$.

Arcs blend together at low values of $\tau$, becoming more clearly demarcated at high values of $\tau$. As the optimal range of $\tau$ to sum over and $f_{D,\text{norm}}$ to fit vary between pulsars and between epochs, we applied the arc curvature fitting algorithm interactively. Alongside the fitting, the dynamic and secondary spectrum of each observation are verified for corruption by RFI or phasing artefacts. The range in $\tau$ to sum to form $S(f_{D,\text{norm}})$ is adjusted per source, and peaks in $f_{D,\text{norm}}$ can be given as initial guesses for a least-squares fit of a parabola to the local region of $S(f_{D,\text{norm}})$, and the arc curvature is given as $\eta = f_{D,\text{norm}}^2/\tau_{\text{ref}}$. An example is shown in Fig. 3.

## 4.3 Time delays from secondary spectra

Secondary spectra express the power of scintillation in terms of the conjugate variables of time and frequency respectively, $f_D$ and $\tau$. In the strong scattering regime, the secondary spectrum contains contributions from all pairs of interfering images. Hemberger & Stinebring (2008) showed how one can estimate the averaged time delay $\langle\tau\rangle$ from the secondary spectrum, and in Main et al. (2020) it was argued that this technique is valid in the limit of a strong, central image arising from a single thin screen, or when the response function is close to an exponential. The total geometric time delay $\langle\tau\rangle$ can then be estimated as

$$ \langle\tau\rangle = \frac{\int_{0}^{T} \tau |I(\tau)|^2 d\tau}{\int_{0}^{T} |I(\tau)|^2 d\tau}. \quad (13) $$

For each source, a range in $f_D$ was chosen to fully encompass the power of the visible scintillation arc, and the background noise was estimated from a region of the same size offset by $f_D \pm 30$ mHz. The integrated profile of $\langle\tau\rangle$ against maximum $\tau_{\text{max}} \equiv T$ was examined in all cases; the value and the error on $\langle\tau\rangle$ were estimated by the mean and standard deviation of the profile once it plateaus, taken in the range of $3T/4 < \tau < T$. Examples of secondary spectra and their associated profiles of $I(\tau)$ are shown in Fig. 4.

## 4.4 Scintillation parameters from the ACF

In some of the sources with diffuse arcs, it is difficult to determine parameters from the secondary spectra. In these cases, we measure the scintillation time-scale and bandwidth in a more traditional way, through the 2D ACF of the dynamic spectrum $R(\Delta\nu) = (I*I)(\Delta\nu)$. The ACF in frequency is fit with a Lorentzian, and with a Gaussian in time (Cordes 1986). The scintillation bandwidth is inversely proportional to the bulk scattering delay as $\langle\tau\rangle = C/2\pi\nu_s$, where the model-dependent constant $C$, depending on the distribution of scattered power, is assumed to be 1, which is the value for a thin screen with a square-law phase structure function (Cordes & Rickett 1998). The 2D ACFs of sources without clear arcs is shown in Fig. 5.

The precision on the measurements on the ACF is high, and the true error is dominated by the fact that there are finite scintils within
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Dynamic spectra of pulsars in the sample. Several of the LEAP sources, only a subset of the band is shown to display the frequency scale of the scintles. The mean of each dynamic spectrum is normalized to 1, and the intensity map covers the range −1–6.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Dynamic spectra of pulsars in the sample. Several of the LEAP sources, only a subset of the band is shown to display the frequency scale of the scintles. The mean of each dynamic spectrum is normalized to 1, and the intensity map covers the range −1–6.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Panorama of scintillation arcs, from the corresponding dynamic spectra in Fig. 1. Scintillation arcs are seen in many sources, with highly varied extents in τ and f_D, owing to differences in the scattering screens, and relative velocities. The arcs are qualitatively different between sources, showing compact regions of power in some (e.g. PSRs J0613-0200, J1600-3053, J0751+1807), clear parabola in others (e.g. PSRs J0621+1002, J1643-1224), and multiple parabolas in PSR B1937+21. The results are described in more detail in Section 5. The intensity maps are logarithmic, extending three orders of magnitude in most cases, or one in B1821-24A.</figcaption></figure>

the observation. For an observation of duration Tobs with bandwidth BW, the fractional error is given by

$$ \sigma_s \approx \left( f_d \frac{T_{\text{obs}}}{t_s} \frac{\text{BW}}{\nu_s} \right)^{-1/2}, \qquad (14) $$

where the filling factor is assumed to be $f_d \approx 0.2$ (Cordes 1986).

# 5 RESULTS FOR INDIVIDUAL SOURCES

The variable time delays and arc curvatures of all of our sources are shown in Fig. 6, and a compilation of derived results are in Table 2. In this section, we describe the results for specific sources.

## 5.1 Isolated pulsars

There are two isolated pulsars in the sample showing scintillation arcs, PSR B1937+21, and PSR B1821-24A. They are in principle useful control sources, where for a fixed screen, variations in W should arise only from Earth's motion. Despite this, they show a range of interesting behaviour owing to dynamic screens.

### 5.1.1 PSR B1937+21

In low-frequency observations of PSR B1937+21, there is previous evidence of a broad scintillation arc (Walker, Demorest & van Straten 2013). In our observations, we detect between one and three screens in a given observation. The secondary spectra are largely devoid of structure in most observations, possibly due to convolution of multiple interacting screens, although an exception is shown in Fig. 4, with distinct structures along the primary arc at delays > 8μs. During certain ranges of time, there does appear to be a dominant screen showing annual variation. However, the phase of the annual curve is not consistent over time, suggesting that the screen orientation is not fixed, or that different screens are varying at different times, shown in Fig. 6. The curvature of the secondary and tertiary screens appear roughly consistent whenever they do reappear, suggesting stability over the eight years of observations. However, due to the difficulty of unambiguously identifying scattering screens, annual fits were not performed.

PSR B1937+21 often emits intrinsically narrow and bright ‘giant pulses’, used in previous LEAP data to directly measure the time variable scattering time-scale (McKee et al. 2019).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Panorama of scintillation arcs, from the corresponding dynamic spectra in Fig. 1. Scintillation arcs are seen in many sources, with highly varied extents in τ and f_D, owing to differences in the scattering screens, and relative velocities. The arcs are qualitatively different between sources, showing compact regions of power in some (e.g. PSRs J0613-0200, J1600-3053, J0751+1807), clear parabola in others (e.g. PSRs J0621+1002, J1643-1224), and multiple parabolas in PSR B1937+21. The results are described in more detail in Section 5. The intensity maps are logarithmic, extending three orders of magnitude in most cases, or one in B1821-24A.</figcaption></figure>

the observation. For an observation of duration Tobs with bandwidth BW, the fractional error is given by

$$ \sigma_s \approx \left( f_d \frac{T_{\text{obs}}}{t_s} \frac{\text{BW}}{\nu_s} \right)^{-1/2}, \qquad (14) $$

where the filling factor is assumed to be $f_d \approx 0.2$ (Cordes 1986).

# 5 RESULTS FOR INDIVIDUAL SOURCES

The variable time delays and arc curvatures of all of our sources are shown in Fig. 6, and a compilation of derived results are in Table 2. In this section, we describe the results for specific sources.

## 5.1 Isolated pulsars

There are two isolated pulsars in the sample showing scintillation arcs, PSR B1937+21, and PSR B1821-24A. They are in principle useful control sources, where for a fixed screen, variations in W should arise only from Earth's motion. Despite this, they show a range of interesting behaviour owing to dynamic screens.

### 5.1.1 PSR B1937+21

In low-frequency observations of PSR B1937+21, there is previous evidence of a broad scintillation arc (Walker, Demorest & van Straten 2013). In our observations, we detect between one and three screens in a given observation. The secondary spectra are largely devoid of structure in most observations, possibly due to convolution of multiple interacting screens, although an exception is shown in Fig. 4, with distinct structures along the primary arc at delays > 8μs. During certain ranges of time, there does appear to be a dominant screen showing annual variation. However, the phase of the annual curve is not consistent over time, suggesting that the screen orientation is not fixed, or that different screens are varying at different times, shown in Fig. 6. The curvature of the secondary and tertiary screens appear roughly consistent whenever they do reappear, suggesting stability over the eight years of observations. However, due to the difficulty of unambiguously identifying scattering screens, annual fits were not performed.

PSR B1937+21 often emits intrinsically narrow and bright ‘giant pulses’, used in previous LEAP data to directly measure the time variable scattering time-scale (McKee et al. 2019).
<figure><img src="image_3.png" /><figcaption>**Figure 3.** An example of the arc curvature fitting algorithm from a LEAP observation of PSR J0613−0200. Left: secondary spectrum, with best-fitting parabola overlaid as the blue dashed line. Bottom: normalized secondary spectrum. Top: weighted sum of the normalized secondary spectrum over τ. The centre and range of $f_D$, norm being fit is shown by the dashed and dotted grey lines, respectively, and the red parabola shows the fit.</figcaption></figure>

Measurements of giant pulse scattering are direct, and unaffected by multiscreen effects. Reassuringly, the trends of $\langle \tau \rangle$ over time is broadly in agreement between what was observed with giant pulses and the values derived from scintillation arcs in this work, with scattering rising from ~0.2 µs to ~0.5 µs in late 2014. This serves as a useful cross-check of both methods. The scattering seems to significantly vary from observation to observation implying variations on less than a month; variations in the DM of PSR B1937+21 have been observed on this time-scale through high-cadence DM measurements (CHIME/Pulsar Collaboration 2021).

### 5.1.2 PSR B1821−24A

PSR B1821−24A is a MSP in the globular cluster M28. The source emits giant pulses, and has been seen to have variable scattering times at L band, with values as large as 25 ± 8µs (Bilous et al. 2015). While this source has a low average S/N, we were able to detect a faint arc in our highest S/N observation, with time delays extending to ~12 µs. Similarly to PSR B1937+21, this source could be a useful control for measuring time delays through scintillation or directly using giant pulse scattering.

The scattering appears to be consistent with the Milky Way ISM, rather than the intracluster gas. Scattering in the intracluster gas would result in very large values of W; approximating $v_{pl}$ with the core radius of $R_c \approx 0.37$ pc, and $v_{pl}$ with the velocity dispersion of $\approx 11 \text{ km s}^{-1}$ (Baumgardt & Hilker 2018; Oliveira et al. 2022), would result in $W \approx 600 \text{ km s kpc}^{-0.5}$, much greater than our measured $W = 36 \pm 16 \text{ km s kpc}^{-0.5}$. Comparably, our measured value of W is easily compatible with typical Earth and screen velocities, at a screen of ~1 kpc. Additionally, NE2001 predicts of $\tau = 2.1$ µs (Cordes & Lazio 2002), comparable to our measured $\langle \tau \rangle = 1.0 \pm 0.5$ µs.

## 5.2 Binary MSPs

### 5.2.1 PSR J0613−0200

The scintillation of this pulsar was previously studied in Main et al. (2020), using LEAP and Effelsberg data from 2013 to 2020. The time delays in 2013 were much larger than in subsequent years, corresponding also to a different screen orientation. In the last two years of observations from 2020 to 2022, including EPTA scintillation data in the past year, the source has experienced heightened scattering fluctuations, with $\langle \tau \rangle$ decreasing to its lowest state in late 2020, and rising to its highest state in 2021 and beyond. This change in scattering properties likely corresponds a change in the observed scattering screen, where either a different screen becomes dominant, or the screen's properties are changing. The modelling of the variable arc curvature of this source will be covered in Section 6.2. The secondary spectra often show distinct, compact features of $\lesssim 0.2$ µs in extent along the main parabola, which can be tracked between observations during our high-cadence Effelsberg campaign in 2020 March–June, covered in Section 6.3.

### 5.2.2 PSR J0621+1002

This pulsar shows well-resolved, low curvature (i.e. large W) scintillation arcs, with faint indications of inverted arclets suggesting an anisotropic screen (Walker et al. 2004; Cordes et al. 2006). The time delays are of order ~0.5 µs, but with large measurement uncertainties due to low S/N per pixel in the secondary spectra arising from the diffuse arcs. The arcs are often featureless but highly anisotropic, changing on the time-scale of months. This likely reflects large, time-variable DM gradients across the screen, discussed in Section 6.4.

This source is an ideal target for annual and orbital fitting of arc curvature (e.g. Reardon et al. 2020; Mall et al. 2022, see Section 6.2). The advance of periastron $\dot{\omega}$ is significantly detected in timing (Perera et al. 2019), which will allow for the component masses to be disentangled when combined with an inclination measurement from scintillation. Additionally, the well-resolved arcs also may enable novel techniques such as the $\theta-\theta$ transformation (Sprenger et al. 2021), which can be used for precise arc curvature measurements (Baker et al. 2022; Sprenger et al. 2022). This will be left to future work.

### 5.2.3 PSR J1600−3053

Similar to PSR J0613−0200, this source shows compact features, with power extending at times to ~16 µs in the τ-axis of the secondary spectrum. The qualitative behaviour of the arcs changes over the course of our observations. In 2016, arcs appear rather broad and diffuse compared to other years, suggesting the combined contributions of multiple screens, or a larger degree of isotropy of the primary screen. In 2019–2020, the secondary spectra are dominated by a small number of discrete moving features, and in observations from 2020 July 25 onwards, there is only a small, featureless concentration of power at low $\tau \lesssim 1$ µs. The total time delays were variable around a mean value of ~400 ns, decreasing to <200 ns between 2019 and 2021. Variable scattering of this source at L band was also found in analyses of the PTA pulsar noise contributions from timing (Goncharov et al. 2021a; Alam et al. 2021b; Chalumeau et al. 2022) The nature of the secondary spectra meant that precise arc curvature measurements were difficult for most observations;
<figure><img src="image_4.png" /><figcaption>**Figure 3.** An example of the arc curvature fitting algorithm from a LEAP observation of PSR J0613−0200. Left: secondary spectrum, with best-fitting parabola overlaid as the blue dashed line. Bottom: normalized secondary spectrum. Top: weighted sum of the normalized secondary spectrum over $τ$. The centre and range of $f_D$, norm being fit is shown by the dashed and dotted grey lines, respectively, and the red parabola shows the fit.</figcaption></figure>

Measurements of giant pulse scattering are direct, and unaffected by multiscreen effects. Reassuringly, the trends of $\langle \tau \rangle$ over time is broadly in agreement between what was observed with giant pulses and the values derived from scintillation arcs in this work, with scattering rising from ~0.2 µs to ~0.5 µs in late 2014. This serves as a useful cross-check of both methods. The scattering seems to significantly vary from observation to observation implying variations on less than a month; variations in the DM of PSR B1937+21 have been observed on this time-scale through high-cadence DM measurements (CHIME/Pulsar Collaboration 2021).

### 5.1.2 PSR B1821−24A

PSR B1821−24A is a MSP in the globular cluster M28. The source emits giant pulses, and has been seen to have variable scattering times at L band, with values as large as $25 \pm 8$ µs (Bilous et al. 2015). While this source has a low average S/N, we were able to detect a faint arc in our highest S/N observation, with time delays extending to ~12 µs. Similarly to PSR B1937+21, this source could be a useful control for measuring time delays through scintillation or directly using giant pulse scattering.

The scattering appears to be consistent with the Milky Way ISM, rather than the intracluster gas. Scattering in the intracluster gas would result in very large values of W; approximating $v_{pl}$ with the core radius of $R_c \approx 0.37$ pc, and $v_{pl}$ with the velocity dispersion of $\approx 11$ km s$^{-1}$ (Baumgardt & Hilker 2018; Oliveira et al. 2022), would result in $W \approx 600$ km s kpc$^{-0.5}$, much greater than our measured $W = 36 \pm 16$ km s kpc$^{-0.5}$. Comparably, our measured value of W is easily compatible with typical Earth and screen velocities, at a screen of ~1 kpc. Additionally, NE2001 predicts of $\tau = 2.1$ µs (Cordes & Lazio 2002), comparable to our measured $\langle \tau \rangle = 1.0 \pm 0.5$ µs.

## 5.2 Binary MSPs

### 5.2.1 PSR J0613−0200

The scintillation of this pulsar was previously studied in Main et al. (2020), using LEAP and Effelsberg data from 2013 to 2020. The time delays in 2013 were much larger than in subsequent years, corresponding also to a different screen orientation. In the last two years of observations from 2020 to 2022, including EPTA scintillation data in the past year, the source has experienced heightened scattering fluctuations, with $\langle \tau \rangle$ decreasing to its lowest state in late 2020, and rising to its highest state in 2021 and beyond. This change in scattering properties likely corresponds a change in the observed scattering screen, where either a different screen becomes dominant, or the screen's properties are changing. The modelling of the variable arc curvature of this source will be covered in Section 6.2. The secondary spectra often show distinct, compact features of $\lesssim 0.2$ µs in extent along the main parabola, which can be tracked between observations during our high-cadence Effelsberg campaign in 2020 March–June, covered in Section 6.3.

### 5.2.2 PSR J0621+1002

This pulsar shows well-resolved, low curvature (i.e. large W) scintillation arcs, with faint indications of inverted arclets suggesting an anisotropic screen (Walker et al. 2004; Cordes et al. 2006). The time delays are of order ~0.5 µs, but with large measurement uncertainties due to low S/N per pixel in the secondary spectra arising from the diffuse arcs. The arcs are often featureless but highly anisotropic, changing on the time-scale of months. This likely reflects large, time-variable DM gradients across the screen, discussed in Section 6.4.

This source is an ideal target for annual and orbital fitting of arc curvature (e.g. Reardon et al. 2020; Mall et al. 2022, see Section 6.2). The advance of periastron $\dot{\omega}$ is significantly detected in timing (Perera et al. 2019), which will allow for the component masses to be disentangled when combined with an inclination measurement from scintillation. Additionally, the well-resolved arcs also may enable novel techniques such as the $\theta-\theta$ transformation (Sprenger et al. 2021), which can be used for precise arc curvature measurements (Baker et al. 2022; Sprenger et al. 2022). This will be left to future work.

### 5.2.3 PSR J1600−3053

Similar to PSR J0613−0200, this source shows compact features, with power extending at times to ~16 µs in the $\tau$-axis of the secondary spectrum. The qualitative behaviour of the arcs changes over the course of our observations. In 2016, arcs appear rather broad and diffuse compared to other years, suggesting the combined contributions of multiple screens, or a larger degree of isotropy of the primary screen. In 2019–2020, the secondary spectra are dominated by a small number of discrete moving features, and in observations from 2020 July 25 onwards, there is only a small, featureless concentration of power at low $\tau \lesssim 1$ µs. The total time delays were variable around a mean value of ~400 ns, decreasing to <200 ns between 2019 and 2021. Variable scattering of this source at L band was also found in analyses of the PTA pulsar noise contributions from timing (Goncharov et al. 2021a; Alam et al. 2021b; Chalumeau et al. 2022) The nature of the secondary spectra meant that precise arc curvature measurements were difficult for most observations;
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Secondary spectra (images), and the associated estimates of the scattering tail $I(\tau)$ by summing over $f_D$ along the arc (top panels). Several of the sources show features at large time delays, beyond a simple exponential tail.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5.** *Images:* 2D ACFs $R(\Delta v, \Delta t)$ of the four pulsars without clear arcs. *Side panels:* cuts of $R(\Delta v, 0)$, $R(0, \Delta t)$ through the 2D ACF (solid black line), and best-fitting models from which the scintillation bandwidth and time-scale are derived (dotted red line).</figcaption></figure>

along with complications from potential changes in the properties of the screen, we leave modelling of the variable arc curvatures of this source to future work.

5.2.4 PSR J1643−1224

The annual and orbital variations of this pulsar were previously studied in Mall et al. (2022), placing the dominant screen distance coincident with Sh 2–27, a large diameter foreground H II region.

PSR J1643−1224 was regularly observed with LEAP from 2012 to 2018, and is observed with Effelsberg as part of EPTA observations, extending our data set from 2012 to 2022. The properties of scintillation arcs in the last year of EPTA observations (i.e. both $\eta$ and their extent in $\tau$) are still consistent with the previous trends, suggesting long-term stability of the screen(s). While the annual variations of the arcs and thus the screen geometry have been stable for ~10 yr, the time delays $\langle\tau\rangle$ are seen to be variable. The scattering measured by scintillation decreases from ~5–2.5 $\mu$s, and roughly correlates with the decreasing DM of the source (Alam et al. 2021a). The scintillation arcs show a persistent asymmetry which is likely related to the DM gradient in this system, which we discuss in Section 6.4.

We note that this source has the finest scintiles of all of the LEAP sources, with $v_s \sim 100$ kHz. At the highest time delays, the scale of the scintillation pattern on Earth is ~2500 km, nearing the length of LEAP’s longest baselines. Equivalently, the angle corresponding to the largest delays is $\theta \sim 18$ mas, while the resolution of LEAP’s
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Secondary spectra (images), and the associated estimates of the scattering tail $I(\tau)$ by summing over $f_D$ along the arc (top panels). Several of the sources show features at large time delays, beyond a simple exponential tail.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** *Images:* 2D ACFs $R(\Delta v, \Delta t)$ of the four pulsars without clear arcs. *Side panels:* cuts of $R(\Delta v, 0)$, $R(0, \Delta t)$ through the 2D ACF (solid black line), and best-fitting models from which the scintillation bandwidth and time-scale are derived (dotted red line).</figcaption></figure>

along with complications from potential changes in the properties of the screen, we leave modelling of the variable arc curvatures of this source to future work.

5.2.4 PSR J1643−1224

The annual and orbital variations of this pulsar were previously studied in Mall et al. (2022), placing the dominant screen distance coincident with Sh 2–27, a large diameter foreground H II region.

PSR J1643−1224 was regularly observed with LEAP from 2012 to 2018, and is observed with Effelsberg as part of EPTA observations, extending our data set from 2012 to 2022. The properties of scintillation arcs in the last year of EPTA observations (i.e. both $\eta$ and their extent in $\tau$) are still consistent with the previous trends, suggesting long-term stability of the screen(s). While the annual variations of the arcs and thus the screen geometry have been stable for ~10 yr, the time delays $\langle\tau\rangle$ are seen to be variable. The scattering measured by scintillation decreases from ~5–2.5 $\mu$s, and roughly correlates with the decreasing DM of the source (Alam et al. 2021a). The scintillation arcs show a persistent asymmetry which is likely related to the DM gradient in this system, which we discuss in Section 6.4.

We note that this source has the finest scintiles of all of the LEAP sources, with $v_s \sim 100$ kHz. At the highest time delays, the scale of the scintillation pattern on Earth is ~2500 km, nearing the length of LEAP’s longest baselines. Equivalently, the angle corresponding to the largest delays is $\theta \sim 18$ mas, while the resolution of LEAP’s
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Time delays (Top), and measurements of *W* (Bottom) for the six sources with fully resolvable scintillation at LEAP. For PSR J0613−0200, the fit of the annual and orbital variations of *W* is shown in red, while the jumps of $\psi$ are shown in blue (fitting described in Section 6.2). For PSR B1937+21, the different colours denote when 1, 2, or 3 parabolas can be identified.</figcaption></figure>

longest baselines is $\lambda/D \approx 21\text{cm}/1200\text{km} \approx 35$ mas. For a pulsar scattered much beyond this, the angle $\theta$ of the furthest images on the sky will be outside of the LEAP beam and be resolved out during coherent addition, complicating the use of LEAP as a single

effective telescope. By the same effect, pulsars with scattering comparable to PSR J1643−1224 can have their screens imaged through VLBI, providing an independent way to determine scattering screen parameters (Brisken et al. 2010). Indeed, Ding et al. (2023) measure
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Time delays (Top), and measurements of *W* (Bottom) for the six sources with fully resolvable scintillation at LEAP. For PSR J0613−0200, the fit of the annual and orbital variations of *W* is shown in red, while the jumps of ψ are shown in blue (fitting described in Section 6.2). For PSR B1937+21, the different colours denote when 1, 2, or 3 parabolas can be identified.</figcaption></figure>

longest baselines is $\lambda/D \approx 21\text{cm}/1200\text{km} \approx 35$ mas. For a pulsar scattered much beyond this, the angle $\theta$ of the furthest images on the sky will be outside of the LEAP beam and be resolved out during coherent addition, complicating the use of LEAP as a single

effective telescope. By the same effect, pulsars with scattering comparable to PSR J1643−1224 can have their screens imaged through VLBI, providing an independent way to determine scattering screen parameters (Brisken et al. 2010). Indeed, Ding et al. (2023) measure
<table><caption>**Table 2.** Summary of quantities derived from the arc curvatures, and from ACF fits. The values of $W \equiv v_{\text{eff},ii}/\sqrt{d_{\text{eff}}}$ and $t_r$ quoted here are an average, representative value, and both will be time-variable owing to the annual motion, pulsar binary motion, and changing screen properties. The values of $t_r$ are computed as in Section 2, and represent the time for a fixed feature to pass from $\langle\tau\rangle$ to $\langle\tau\rangle$ given the average value of $W$.</caption><thead><tr><th>Pulsar name</th><th>⟨τ⟩<br>(μs)</th><th>⟨W⟩<br>(km s<sup>-1</sup> kpc<sup>-0.5</sup>)</th><th>⟨t<sub>r</sub>⟩<br>(d)</th></tr></thead><tbody><tr><td>J0613-0200</td><td>0.13 ± 0.02</td><td>34 ± 3</td><td>33 ± 3</td></tr><tr><td>J0621+1002</td><td>0.47 ± 0.02</td><td>111 ± 6</td><td>19 ± 1</td></tr><tr><td>J0751+1807</td><td>0.62 ± 0.04</td><td>35 ± 2</td><td>70 ± 5</td></tr><tr><td>J1600-3053</td><td>0.34 ± 0.02</td><td>20 ± 1</td><td>93 ± 7</td></tr><tr><td>J1643-1224</td><td>2.8 ± 0.1</td><td>33 ± 3</td><td>157 ± 15</td></tr><tr><td>J1713+0747</td><td>0.010 ± 0.003</td><td>74 ± 21</td><td>4.2 ± 1.3</td></tr><tr><td>B1821-24A</td><td>1.0 ± 0.5</td><td>36 ± 16</td><td>88 ± 45</td></tr><tr><td>J1832-0836</td><td>0.031 ± 0.002</td><td>240 ± 16</td><td>2.3 ± 0.2</td></tr><tr><td>B1855+09</td><td>0.040 ± 0.004</td><td>45 ± 5</td><td>14 ± 2</td></tr><tr><td>J1918-0642</td><td>0.05 ± 0.01</td><td>16 ± 3</td><td>45 ± 9</td></tr><tr><td>B1937+21</td><td>0.18 ± 0.01</td><td>57 ± 5</td><td>23 ± 2</td></tr><tr><td>J2010-1323</td><td>0.049 ± 0.006</td><td>28 ± 3</td><td>25 ± 3</td></tr></tbody></table>

PSR J1643−1224 to be angularly broadened to $\theta = 3.65 \pm 0.43$ mas using the Very Large Baseline Array, and confirm the association of the dominant scattering screen with Sh 2–27.

### 5.2.5 PSR J1918−0642

The scintillation time-scale is ~10 min at its shortest, sufficiently short to reveal a faint arc in LEAP observations. At its slowest, the scintillation time-scale is greater than the observation length, appearing in the secondary spectrum as power along the $f_D = 0$ axis. Moreover, the source changed significantly around 2016, as scintles beforehand were tens of MHz, and transitioned to ≲1 MHz afterwards, indicating a large rise in scattering time (see Fig. A5). Coincidentally, LEAP observations taken until 2016 January were only 20 min long, insufficient to resolve scintillation in time. As the arcs were difficult to resolve, we have too few curvatures measurements to fit for annual and orbital variations but we do note that the arc curvature is clearly variable with contribution from both, as it is not always the same at a given time of year or orbital phase.

# 6 DISCUSSION

## 6.1 Time delays

The measured values of $\langle\tau\rangle$ are shown in Fig. 6, with a summary of results in Table 2. The precision of pulse times of arrival is ~1μs in the most precisely-timed EPTA sources (Chen et al. 2021). While the scattering time-scales are less than this for most of the sources in our sample, uncorrected scattering variations could be a significant source of red noise, as they are correlated in time (Goncharov et al. 2021b). Moreover, large-scale variations in $\langle\tau\rangle$ are correlated on the time-scale of years in several sources, and variations on similar time-scales could masquerade as a common signal shared between pulsars. The subset of pulsars with resolvable scintillation with LEAP have the longest scattering times by design, but this subset also includes some of the pulsars with the highest timing precision in the EPTA. In particular, PSRs J0613−0200, J1600−3053, and J1918−0642 are among the top five most significant contributors to the common red noise process detected by the IPTA (Antoniadis et al. 2022), and all show variable scattering on ≳ 100 ns level in this work. As we

approach a potential GW detection, and as PTA sensitivity increases, the variable time delays from scattering will be important to consider in GW searches.

## 6.2 Annual and orbital fitting of PSR J0613−0200

In Main et al. (2020), a strong annual trend in the arc curvatures of PSR J0613−0200 was seen. Fitting the annual variations resulted in a fractional screen distance of $s = 0.58 \pm 0.10$ during the period of increased scattering in 2013, and $s = 0.62 \pm 0.06$ afterwards, with screen axis $\psi$ changing by ~50°. It was argued that the scattering could originate in the same screen, with orientation changing over time. The orbital variations were ignored, and were an additional source of scatter in individual measurements.

Here, we revisit the curvature variations of PSR J0613−0200, including the most recent data, improved measurement of arc curvatures, and including fitting orbital variations using the same framework as in Mall et al. (2022). We use Gaussian priors on the proper motion and distance from the most recent IPTA values, $\mu_\alpha = 1.828(5)$ mas yr$^{-1}$, $\mu_\delta = -10.35(1)$ mas yr$^{-1}$, $d_{\text{psr}} = 1.11 \pm 0.05$ kpc (Perera et al. 2019). We fit the distance weighted effective velocities $W$ with a model of an anisotropic scattering screen,

$$W = \frac{1}{\sqrt{d_{\text{eff}}}} \left| \left( \frac{1-s}{s} v_{\text{psr},\alpha} + v_{\oplus,\alpha} - \frac{1}{s} v_{\text{scr},\alpha} \right) \sin(\psi) + \left( \frac{1-s}{s} v_{\text{psr},\delta} + v_{\oplus,\delta} - \frac{1}{s} v_{\text{scr},\delta} \right) \cos(\psi) \right|. \quad (15)$$

We compare several different models, detailed in the following sections.

### 6.2.1 Variations in screen properties

From Main et al. (2020), we know a single screen is a poor fit to the full 2013–2019 dataspan. Here, we try two models to account for the variable screen. For each model, we allow for $N$ jumps of the screen parameters, where the times of the jumps are free parameters, bounded between the time of the first and last observation. In the first model, we assume that the scattering originates from a single screen, which only changes in orientation over time [similar to the 1D scattering screens of PSR B0834+06, and B1508 + 55 (Simard et al. 2019; Sprenger et al. 2022)]. In this case, the fractional screen distance $s$, and 2D velocity $v_{\text{scr},\alpha}$, $v_{\text{scr},\delta}$ are free parameters and constant over the full duration. In the second model, we allow $s$, $\psi$, and $v_\psi$ to vary in each jump. This is highly similar to the approach in Walker et al. (2022), who model the arc curvatures of PSR J1603−7202. They allow all screen properties to change between jumps, and find moderately strong support for two jumps, one corresponding to a region of enhanced DM and scattering.

In addition, we fit for the orbital inclination $i$ and angle of nodes $\Omega$, and we include white noise parameters $F$ and $Q$ as free parameters of the fit, such that the scaled errors are $\delta W_{\text{corr}} = \sqrt{(F \times \delta W)^2 + Q^2}$. These parameters can account for biases and underestimated errors, but could also arise physically from variations in screen axis and velocity, which could vary on the refractive time-scale $t_r$ (e.g. Askew, Reardon & Shannon 2023).

The first model, allowing for all screen parameters to vary in each jump, results in $BIC = 483.4$, while the second model with a screen at a fixed distance, changing only the orientation has $BIC = 442.5$. Both models have almost identical white noise parameters $F$ and $Q$, suggesting that they fit the data comparably well, but the first model
<table><caption>**Table 2.** Summary of quantities derived from the arc curvatures, and from ACF fits. The values of $W \equiv v_{\text{eff},ii}/\sqrt{d_{\text{eff}}}$ and $t_r$ quoted here are an average, representative value, and both will be time-variable owing to the annual motion, pulsar binary motion, and changing screen properties. The values of $t_r$ are computed as in Section 2, and represent the time for a fixed feature to pass from $-\langle\tau\rangle$ to $\langle\tau\rangle$ given the average value of $W$.</caption><thead><tr><th>Pulsar name</th><th>⟨τ⟩<br>(μs)</th><th>⟨W⟩<br>(km s<sup>-1</sup> kpc<sup>-0.5</sup>)</th><th>⟨t<sub>r</sub>⟩<br>(d)</th></tr></thead><tbody><tr><td>J0613-0200</td><td>0.13 ± 0.02</td><td>34 ± 3</td><td>33 ± 3</td></tr><tr><td>J0621+1002</td><td>0.47 ± 0.02</td><td>111 ± 6</td><td>19 ± 1</td></tr><tr><td>J0751+1807</td><td>0.62 ± 0.04</td><td>35 ± 2</td><td>70 ± 5</td></tr><tr><td>J1600-3053</td><td>0.34 ± 0.02</td><td>20 ± 1</td><td>93 ± 7</td></tr><tr><td>J1643-1224</td><td>2.8 ± 0.1</td><td>33 ± 3</td><td>157 ± 15</td></tr><tr><td>J1713+0747</td><td>0.010 ± 0.003</td><td>74 ± 21</td><td>4.2 ± 1.3</td></tr><tr><td>B1821-24A</td><td>1.0 ± 0.5</td><td>36 ± 16</td><td>88 ± 45</td></tr><tr><td>J1832-0836</td><td>0.031 ± 0.002</td><td>240 ± 16</td><td>2.3 ± 0.2</td></tr><tr><td>B1855+09</td><td>0.040 ± 0.004</td><td>45 ± 5</td><td>14 ± 2</td></tr><tr><td>J1918-0642</td><td>0.05 ± 0.01</td><td>16 ± 3</td><td>45 ± 9</td></tr><tr><td>B1937+21</td><td>0.18 ± 0.01</td><td>57 ± 5</td><td>23 ± 2</td></tr><tr><td>J2010-1323</td><td>0.049 ± 0.006</td><td>28 ± 3</td><td>25 ± 3</td></tr></tbody></table>

PSR J1643−1224 to be angularly broadened to $\theta = 3.65 \pm 0.43$ mas using the Very Large Baseline Array, and confirm the association of the dominant scattering screen with Sh 2–27.

### 5.2.5 PSR J1918−0642

The scintillation time-scale is ~10 min at its shortest, sufficiently short to reveal a faint arc in LEAP observations. At its slowest, the scintillation time-scale is greater than the observation length, appearing in the secondary spectrum as power along the $f_D = 0$ axis. Moreover, the source changed significantly around 2016, as scintles beforehand were tens of MHz, and transitioned to ≲1 MHz afterwards, indicating a large rise in scattering time (see Fig. A5). Coincidentally, LEAP observations taken until 2016 January were only 20 min long, insufficient to resolve scintillation in time. As the arcs were difficult to resolve, we have too few curvatures measurements to fit for annual and orbital variations but we do note that the arc curvature is clearly variable with contribution from both, as it is not always the same at a given time of year or orbital phase.

# 6 DISCUSSION

## 6.1 Time delays

The measured values of $\langle\tau\rangle$ are shown in Fig. 6, with a summary of results in Table 2. The precision of pulse times of arrival is ~1μs in the most precisely-timed EPTA sources (Chen et al. 2021). While the scattering time-scales are less than this for most of the sources in our sample, uncorrected scattering variations could be a significant source of red noise, as they are correlated in time (Goncharov et al. 2021b). Moreover, large-scale variations in $\langle\tau\rangle$ are correlated on the time-scale of years in several sources, and variations on similar time-scales could masquerade as a common signal shared between pulsars. The subset of pulsars with resolvable scintillation with LEAP have the longest scattering times by design, but this subset also includes some of the pulsars with the highest timing precision in the EPTA. In particular, PSRs J0613−0200, J1600−3053, and J1918−0642 are among the top five most significant contributors to the common red noise process detected by the IPTA (Antoniadis et al. 2022), and all show variable scattering on ≳ 100 ns level in this work. As we

approach a potential GW detection, and as PTA sensitivity increases, the variable time delays from scattering will be important to consider in GW searches.

## 6.2 Annual and orbital fitting of PSR J0613−0200

In Main et al. (2020), a strong annual trend in the arc curvatures of PSR J0613−0200 was seen. Fitting the annual variations resulted in a fractional screen distance of $s = 0.58 \pm 0.10$ during the period of increased scattering in 2013, and $s = 0.62 \pm 0.06$ afterwards, with screen axis $\psi$ changing by ~50°. It was argued that the scattering could originate in the same screen, with orientation changing over time. The orbital variations were ignored, and were an additional source of scatter in individual measurements.

Here, we revisit the curvature variations of PSR J0613−0200, including the most recent data, improved measurement of arc curvatures, and including fitting orbital variations using the same framework as in Mall et al. (2022). We use Gaussian priors on the proper motion and distance from the most recent IPTA values, $\mu_\alpha = 1.828(5)$ mas yr$^{-1}$, $\mu_\delta = -10.35(1)$ mas yr$^{-1}$, $d_{psr} = 1.11 \pm 0.05$ kpc (Perera et al. 2019). We fit the distance weighted effective velocities $W$ with a model of an anisotropic scattering screen,

$$W = \frac{1}{\sqrt{d_{\text{eff}}}} \left| \left( \frac{1-s}{s} v_{\text{psr},\alpha} + v_{\oplus,\alpha} - \frac{1}{s} v_{\text{scr},\alpha} \right) \sin(\psi) + \left( \frac{1-s}{s} v_{\text{psr},\delta} + v_{\oplus,\delta} - \frac{1}{s} v_{\text{scr},\delta} \right) \cos(\psi) \right|. \quad (15)$$

We compare several different models, detailed in the following sections.

### 6.2.1 Variations in screen properties

From Main et al. (2020), we know a single screen is a poor fit to the full 2013–2019 dataspan. Here, we try two models to account for the variable screen. For each model, we allow for $N$ jumps of the screen parameters, where the times of the jumps are free parameters, bounded between the time of the first and last observation. In the first model, we assume that the scattering originates from a single screen, which only changes in orientation over time [similar to the 1D scattering screens of PSR B0834+06, and B1508 + 55 (Simard et al. 2019; Sprenger et al. 2022)]. In this case, the fractional screen distance $s$, and 2D velocity $v_{\text{scr},\alpha}$, $v_{\text{scr},\delta}$ are free parameters and constant over the full duration. In the second model, we allow $s$, $\psi$, and $v_\psi$ to vary in each jump. This is highly similar to the approach in Walker et al. (2022), who model the arc curvatures of PSR J1603−7202. They allow all screen properties to change between jumps, and find moderately strong support for two jumps, one corresponding to a region of enhanced DM and scattering.

In addition, we fit for the orbital inclination $i$ and angle of nodes $\Omega$, and we include white noise parameters $F$ and $Q$ as free parameters of the fit, such that the scaled errors are $\delta W_{\text{corr}} = \sqrt{(F \times \delta W)^2 + Q^2}$. These parameters can account for biases and underestimated errors, but could also arise physically from variations in screen axis and velocity, which could vary on the refractive time-scale $t_r$ (e.g. Askew, Reardon & Shannon 2023).

The first model, allowing for all screen parameters to vary in each jump, results in $BIC = 483.4$, while the second model with a screen at a fixed distance, changing only the orientation has $BIC = 442.5$. Both models have almost identical white noise parameters $F$ and $Q$, suggesting that they fit the data comparably well, but the first model
<table><caption>**Table 3.** Results of modelling arc curvature variations of PSR J0613−0200, and comparison to Main et al. (2020). The model assumes an anisotropic scattering screen at a fixed distance, but allowing for two jumps in the screen orientation $\psi$, as described in Section 6.2. The top set of parameters are the free parameters of the model, while the bottom set includes the pulsar distance prior (distances from*: Perera et al. 2019, †: Desvignes et al. 2016), and derived quantities.</caption><thead><tr><th>Parameters</th><th>Model 1</th><th>Model 2</th><th>Main et al. 2020</th></tr></thead><tbody><tr><td>s</td><td>0.71* ± 0.02</td><td>—</td><td>—</td></tr><tr><td>s0</td><td>—</td><td>0.40* ± 0.19</td><td>0.58† ± 0.10</td></tr><tr><td>s1</td><td>—</td><td>0.73* ± 0.03</td><td>0.62† ± 0.06</td></tr><tr><td>s2</td><td>—</td><td>0.50* ± 0.12</td><td>—</td></tr><tr><td>v<sub>scr, α</sub> (km s<sup>-1</sup>)</td><td>17.2 ± 0.9</td><td>—</td><td>—</td></tr><tr><td>v<sub>scr, δ</sub> (km s<sup>-1</sup>)</td><td>-4.4 ± 0.8</td><td>—</td><td>—</td></tr><tr><td>ψ<sub>0</sub>(°)</td><td>-30 ± 3</td><td>-31 ± 3</td><td>-36 ± 9</td></tr><tr><td>ψ<sub>1</sub>(°)</td><td>16 ± 2</td><td>16 ± 2</td><td>16 ± 2</td></tr><tr><td>ψ<sub>2</sub>(°)</td><td>40 ± 3</td><td>42 ± 3</td><td>—</td></tr><tr><td>i(°)</td><td>58 ± 4</td><td>55 ± 5</td><td>—</td></tr><tr><td>Ω(°)</td><td>124 ± 4</td><td>126 ± 5</td><td>—</td></tr><tr><td>T<sub>1,jump</sub> (mjd)</td><td>56670 ± 30</td><td>56667 ± 30</td><td>56 658 (fixed)</td></tr><tr><td>T<sub>2,jump</sub> (mjd)</td><td>58820 ± 30</td><td>58820 ± 25</td><td>—</td></tr><tr><td>Q</td><td>4.4 ± 0.9</td><td>4.4 ± 1.0</td><td>—</td></tr><tr><td>F</td><td>3.1 ± 0.9</td><td>3.3 ± 0.9</td><td>—</td></tr><tr><td>d<sub>psr</sub> (kpc)</td><td>1.11* ± 0.05</td><td>1.12* ± 0.05</td><td>0.78† ± 0.08</td></tr><tr><td>d<sub>scr</sub> (kpc)</td><td>0.32 ± 0.04</td><td>0.30 ± 0.03</td><td>0.30 ± 0.07</td></tr><tr><td>v<sub>ψ<sub>0</sub></sub> (km s<sup>-1</sup>)</td><td>0.5 ± 0.8</td><td>2.1 ± 1.9</td><td>-1.2 ± 2.5</td></tr><tr><td>v<sub>ψ<sub>1</sub></sub> (km s<sup>-1</sup>)</td><td>12.4 ± 0.8</td><td>11.3 ± 1.2</td><td>12.8 ± 2.8</td></tr><tr><td>v<sub>ψ<sub>2</sub></sub> (km s<sup>-1</sup>)</td><td>7.7 ± 0.9</td><td>10.7 ± 1.3</td><td>—</td></tr><tr><td>BIC</td><td>455.8</td><td>470.5</td><td>—</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Bayesian information criteria (BIC) in fitting W variations of PSR J0613−0200, for models with N jumps in the screen orientation $\psi$. Models with fewer than two jumps result in a poor fit, while models with more than two jumps result in no additional improvement and result in increased BIC due to having more parameters.</figcaption></figure>

is penalized for having more free parameters. We suggest that a single screen can reproduce our observations, but we cannot rule out the possibility that the variable screen properties correspond to different screens dominating at different times. The best-fitting values of these models are tabulated in Table 3, and the BIC of the first model as a function of the number of jumps is shown in Fig. 7.

### 6.2.2 Orbital constraints

Scintillation arcs provide a way to measure resolve the ambiguity in the sense of the inclination, i.e. $i < 90^{\circ}$, or $i > 90^{\circ}$. We obtain a fit orbital inclination of $i = 58 \pm 4^{\circ}$ is consistent with one value from timing of $i_{\text{timing}} = 68^{\circ}_{-10} (\text{Fonseca et al. 2016})$, but mildly inconsistent with the IPTA value of $i_{\text{timing}} = 70 \pm 3^{\circ}$ by $<2\sigma$. We obtain the first measurement of $\Omega = 124 \pm 4^{\circ}$. Additionally, we compare to the fit restricting $i > 90^{\circ}$. This results in values of $i = 116.4 \pm 8.4$, and $\Omega = 274.0 \pm 12.8$, but is disfavoured, with BIC = 523

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Results of modelling annual and orbital scintillation arc variations of PSR J0613−0200, shown during two periods of heightened scattering where arcs could be measured precisely. The orbital amplitude, and the phase of the annual curves are clearly different between the two, indicating a change of screen geometry; while $\Omega$ and $\psi$ are almost perpendicular in 2021–2022, the changing screen orientation ensures that the orbital modulation is seen.</figcaption></figure>

The best fit model is overlaid on the time series of W measurements in Fig. 6, and the decomposition to annual and orbital velocity is shown in Fig. 8. The variation of the properties of the observed screen can be clearly seen as the phase of the annual maxima and minima change over time, and the amplitude of the orbital curve changes due to the changing alignment of $\psi$ and $\Omega$. The times of the jumps both correspond to regions of changing $\langle\tau\rangle$, suggesting that both of these effects trace physical changes of the screen.

### 6.2.3 Comparison to previous work

Our results are largely consistent with, yet more precise than Main et al. (2020). This is unsurprising, as both analysis contain much of the same data, but the measurements of W are made differently, and our present analysis includes data beyond 2020, orbital variations, and times of screen jumps as free parameters rather than fixed. However, our value of $s = 0.71 \pm 0.02$ is smaller than the value of $s = 0.62 \pm 0.06$. This is a result of the different pulsar distance used, $d_{\text{psr}} = 780 \pm 80$ pc from Desvignes et al. (2016), and both measurements result in a consistent screen distance.

## 6.3 Movement of features

Several of the sources in our sample show discrete, compact regions of power in their secondary spectra, as described in Section 2.2, the movement of compact features through the secondary spectrum can be predicted through the arc curvature. We investigate the feature movement in PSR J0613−0200 and PSR J1600−3053.

Following the techniques of Sprenger et al. (2022), we remap the secondary spectra in terms of $\sqrt{\tau} \propto \theta$. Surrounding the best-fit parabola $\eta$ of each secondary spectrum, we take a slice of $I(f_D, \tau_i)$ for each bin in $\tau$, remapping every $\tau$ bin to the closest value of $\sqrt{\tau}$ to form a profile of $I(f_D, \sqrt{\tau})$. The shift of features between observations is predicted by the annual model of $\eta$ (assuming the orbital variations average out over time). The results of a subset of the data in PSR J0613−0200 and PSR J1600−3053 is shown in Figs 9 and 10, respectively. In both cases, features can be seen to persist at a fixed location over the time they traverse the secondary spectra, indicating scattering from compact regions of fixed $\theta$. This can only be seen if there is significant power near the undeflected image of the pulsar; we see that there is persistently significant power
<table><caption>**Table 3.** Results of modelling arc curvature variations of PSR J0613−0200, and comparison to Main et al. (2020). The model assumes an anisotropic scattering screen at a fixed distance, but allowing for two jumps in the screen orientation $\psi$, as described in Section 6.2. The top set of parameters are the free parameters of the model, while the bottom set includes the pulsar distance prior (distances from*: Perera et al. 2019, †: Desvignes et al. 2016), and derived quantities.</caption><thead><tr><th>Parameters</th><th>Model 1</th><th>Model 2</th><th>Main et al. 2020</th></tr></thead><tbody><tr><td>s</td><td>0.71* ± 0.02</td><td>—</td><td>—</td></tr><tr><td>s0</td><td>—</td><td>0.40* ± 0.19</td><td>0.58† ± 0.10</td></tr><tr><td>s1</td><td>—</td><td>0.73* ± 0.03</td><td>0.62† ± 0.06</td></tr><tr><td>s2</td><td>—</td><td>0.50* ± 0.12</td><td>—</td></tr><tr><td>v<sub>scr, α</sub> (km s<sup>-1</sup>)</td><td>17.2 ± 0.9</td><td>—</td><td>—</td></tr><tr><td>v<sub>scr, δ</sub> (km s<sup>-1</sup>)</td><td>-4.4 ± 0.8</td><td>—</td><td>—</td></tr><tr><td>ψ<sub>0</sub>(°)</td><td>-30 ± 3</td><td>-31 ± 3</td><td>-36 ± 9</td></tr><tr><td>ψ<sub>1</sub>(°)</td><td>16 ± 2</td><td>16 ± 2</td><td>16 ± 2</td></tr><tr><td>ψ<sub>2</sub>(°)</td><td>40 ± 3</td><td>42 ± 3</td><td>—</td></tr><tr><td>i(°)</td><td>58 ± 4</td><td>55 ± 5</td><td>—</td></tr><tr><td>Ω(°)</td><td>124 ± 4</td><td>126 ± 5</td><td>—</td></tr><tr><td>T<sub>1, jump</sub> (mjd)</td><td>56670 ± 30</td><td>56667 ± 30</td><td>56 658 (fixed)</td></tr><tr><td>T<sub>2, jump</sub> (mjd)</td><td>58820 ± 30</td><td>58820 ± 25</td><td>—</td></tr><tr><td>Q</td><td>4.4 ± 0.9</td><td>4.4 ± 1.0</td><td>—</td></tr><tr><td>F</td><td>3.1 ± 0.9</td><td>3.3 ± 0.9</td><td>—</td></tr><tr><td>d<sub>psr</sub> (kpc)</td><td>1.11* ± 0.05</td><td>1.12* ± 0.05</td><td>0.78† ± 0.08</td></tr><tr><td>d<sub>scr</sub> (kpc)</td><td>0.32 ± 0.04</td><td>0.30 ± 0.03</td><td>0.30 ± 0.07</td></tr><tr><td>v<sub>ψ<sub>0</sub></sub> (km s<sup>-1</sup>)</td><td>0.5 ± 0.8</td><td>2.1 ± 1.9</td><td>-1.2 ± 2.5</td></tr><tr><td>v<sub>ψ<sub>1</sub></sub> (km s<sup>-1</sup>)</td><td>12.4 ± 0.8</td><td>11.3 ± 1.2</td><td>12.8 ± 2.8</td></tr><tr><td>v<sub>ψ<sub>2</sub></sub> (km s<sup>-1</sup>)</td><td>7.7 ± 0.9</td><td>10.7 ± 1.3</td><td>—</td></tr><tr><td>BIC</td><td>455.8</td><td>470.5</td><td>—</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Bayesian information criteria (BIC) in fitting W variations of PSR J0613−0200, for models with N jumps in the screen orientation $\psi$. Models with fewer than two jumps result in a poor fit, while models with more than two jumps result in no additional improvement and result in increased BIC due to having more parameters.</figcaption></figure>

is penalized for having more free parameters. We suggest that a single screen can reproduce our observations, but we cannot rule out the possibility that the variable screen properties correspond to different screens dominating at different times. The best-fitting values of these models are tabulated in Table 3, and the BIC of the first model as a function of the number of jumps is shown in Fig. 7.

### 6.2.2 Orbital constraints

Scintillation arcs provide a way to measure resolve the ambiguity in the sense of the inclination, i.e. $i < 90^{\circ}$, or $i > 90^{\circ}$. We obtain a fit orbital inclination of $i = 58 \pm 4^{\circ}$ is consistent with one value from timing of $i_{\text{timing}} = 68^{\circ}_{-10} (\text{Fonseca et al. 2016})$, but mildly inconsistent with the IPTA value of $i_{\text{timing}} = 70 \pm 3^{\circ}$ by $<2\sigma$. We obtain the first measurement of $\Omega = 124 \pm 4^{\circ}$. Additionally, we compare to the fit restricting $i > 90$. This results in values of $i = 116.4 \pm 8.4$, and $\Omega = 274.0 \pm 12.8$, but is disfavoured, with BIC = 523

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Results of modelling annual and orbital scintillation arc variations of PSR J0613−0200, shown during two periods of heightened scattering where arcs could be measured precisely. The orbital amplitude, and the phase of the annual curves are clearly different between the two, indicating a change of screen geometry; while $\Omega$ and $\psi$ are almost perpendicular in 2021–2022, the changing screen orientation ensures that the orbital modulation is seen.</figcaption></figure>

The best fit model is overlaid on the time series of W measurements in Fig. 6, and the decomposition to annual and orbital velocity is shown in Fig. 8. The variation of the properties of the observed screen can be clearly seen as the phase of the annual maxima and minima change over time, and the amplitude of the orbital curve changes due to the changing alignment of $\psi$ and $\Omega$. The times of the jumps both correspond to regions of changing $\langle\tau\rangle$, suggesting that both of these effects trace physical changes of the screen.

### 6.2.3 Comparison to previous work

Our results are largely consistent with, yet more precise than Main et al. (2020). This is unsurprising, as both analysis contain much of the same data, but the measurements of W are made differently, and our present analysis includes data beyond 2020, orbital variations, and times of screen jumps as free parameters rather than fixed. However, our value of $s = 0.71 \pm 0.02$ is smaller than the value of $s = 0.62 \pm 0.06$. This is a result of the different pulsar distance used, $d_{\text{psr}} = 780 \pm 80$ pc from Desvignes et al. (2016), and both measurements result in a consistent screen distance.

## 6.3 Movement of features

Several of the sources in our sample show discrete, compact regions of power in their secondary spectra, as described in Section 2.2, the movement of compact features through the secondary spectrum can be predicted through the arc curvature. We investigate the feature movement in PSR J0613−0200 and PSR J1600−3053.

Following the techniques of Sprenger et al. (2022), we remap the secondary spectra in terms of $\sqrt{\tau} \propto \theta$. Surrounding the best-fit parabola $\eta$ of each secondary spectrum, we take a slice of $I(f_D, \tau_i)$ for each bin in $\tau$, remapping every $\tau$ bin to the closest value of $\sqrt{\tau}$ to form a profile of $I(f_D, \sqrt{\tau})$. The shift of features between observations is predicted by the annual model of $\eta$ (assuming the orbital variations average out over time). The results of a subset of the data in PSR J0613−0200 and PSR J1600−3053 is shown in Figs 9 and 10, respectively. In both cases, features can be seen to persist at a fixed location over the time they traverse the secondary spectra, indicating scattering from compact regions of fixed $\theta$. This can only be seen if there is significant power near the undeflected image of the pulsar; we see that there is persistently significant power
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Feature alignment in PSR J0613−0200, during the dense Effelsberg observing campaign. Top: Secondary spectra, where several discrete features can be seen to move throughout. Bottom: Profiles of $\sqrt{\tau}$ versus $f_D$ from the corresponding secondary spectra, made as described in Section 6.3, and shifted by the predicted movement between observations. The value $\sqrt{\tau} \propto \theta$ is a proxy for the image positions; features connected vertically between observations suggest persistent scattering at regions of fixed $\theta$.</figcaption></figure>

surrounding $\theta = 0$ which is difficult to track between observations. This contains the bulk of the power, dominating the changing values of $\langle \tau \rangle$.

## 6.4 DM gradients and asymmetric arcs

Scintillation arcs often show a clear asymmetry in power, related to the phase structure across the scattering screen. A local linear DM slope along the direction of $v_{eff}$ creates a refractive shift (Cordes et al. 2006; Rickett et al. 2014),

$$ \theta_r = \frac{\lambda^2 r_e}{2\pi v_{eff,11}} \partial_t DM \quad (16) $$

<figure><img src="image_10.png" /><figcaption>**Figure 10.** Same as 9, but for a series of LEAP PSR J1600−3053 observations showing moving features.</figcaption></figure>

leading to a new zero point of a secondary spectrum offset by $\tau_r$ as:

$$ \partial_t DM = \frac{2\pi v_{eff,11}}{\lambda^2 r_e} \sqrt{\frac{2c\tau_r}{d_{eff}}} \quad (17) $$

Under these assumptions, the gradient in DM within the screen can be estimated from the asymmetry in scintillation arcs, and vice versa (shown in practice in Reardon & Coles 2023 using scintillation ACFs). Additionally, the relation between the two depends on the distance weighted effective velocity; connecting all related observables will allow for the maximum amount of information to be extracted about intervening scattering screens.

In our sample, PSRs J0621+1002 and J1643−1224 are the clearest examples showing diffuse, highly asymmetric scintillation arcs, likely reflecting significant variations in DM (example secondary spectra shown in Fig. 11). The DM curve of PSR J1643−1224 from the NANOGrav 12.5 yr data release shows a persistent downwards trend in DM of $\Delta DM \sim 10^{-3}$ pc cm$^{-3}$ yr$^{-1}$ from ~2013 to 2016 (Alam et al. 2021a), during which time the scintillation arcs showed persistent asymmetric power to the right quadrant of the secondary spectrum. The sign of asymmetry in PSR J0621+1002 changes on the time-scale of months, which may suggest rapidly varying DM. This explanation is plausible as observations with the Low Frequency Array (LOFAR) at frequencies of about 140 MHz the DM of J0621+1002 has been seen to vary by ~$10^{-2}$ pc cm$^{-3}$ on several month time-scales (Donner et al. 2020). However, we note
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Feature alignment in PSR J0613−0200, during the dense Effelsberg observing campaign. Top: Secondary spectra, where several discrete features can be seen to move throughout. Bottom: Profiles of $\sqrt{\tau}$ versus $f_D$ from the corresponding secondary spectra, made as described in Section 6.3, and shifted by the predicted movement between observations. The value $\sqrt{\tau} \propto \theta$ is a proxy for the image positions; features connected vertically between observations suggest persistent scattering at regions of fixed $\theta$.</figcaption></figure>

surrounding $\theta = 0$ which is difficult to track between observations. This contains the bulk of the power, dominating the changing values of $\langle \tau \rangle$.

## 6.4 DM gradients and asymmetric arcs

Scintillation arcs often show a clear asymmetry in power, related to the phase structure across the scattering screen. A local linear DM slope along the direction of $v_{eff}$ creates a refractive shift (Cordes et al. 2006; Rickett et al. 2014),

$$ \theta_r = \frac{\lambda^2 r_e}{2\pi v_{eff,11}} \partial_t DM \qquad (16) $$

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Same as 9, but for a series of LEAP PSR J1600−3053 observations showing moving features.</figcaption></figure>

leading to a new zero point of a secondary spectrum offset by $\tau_r$ as:

$$ \partial_t DM = \frac{2\pi v_{eff,11}}{\lambda^2 r_e} \sqrt{\frac{2c\tau_r}{d_{eff}}} \qquad (17) $$

Under these assumptions, the gradient in DM within the screen can be estimated from the asymmetry in scintillation arcs, and vice versa (shown in practice in Reardon & Coles 2023 using scintillation ACFs). Additionally, the relation between the two depends on the distance weighted effective velocity; connecting all related observables will allow for the maximum amount of information to be extracted about intervening scattering screens.

In our sample, PSRs J0621+1002 and J1643−1224 are the clearest examples showing diffuse, highly asymmetric scintillation arcs, likely reflecting significant variations in DM (example secondary spectra shown in Fig. 11). The DM curve of PSR J1643−1224 from the NANOGrav 12.5 yr data release shows a persistent downwards trend in DM of $\Delta DM \sim 10^{-3}$ pc cm$^{-3}$ yr$^{-1}$ from ~2013 to 2016 (Alam et al. 2021a), during which time the scintillation arcs showed persistent asymmetric power to the right quadrant of the secondary spectrum. The sign of asymmetry in PSR J0621+1002 changes on the time-scale of months, which may suggest rapidly varying DM. This explanation is plausible as observations with the Low Frequency Array (LOFAR) at frequencies of about 140 MHz the DM of J0621+1002 has been seen to vary by ~$10^{-2}$ pc cm$^{-3}$ on several month time-scales (Donner et al. 2020). However, we note
<figure><img src="image_11.png" /><figcaption>**Figure 11.** *Top:* PSR J1643−1224 secondary spectra at a similar time of year, showing a very similar distribution of power, with an asymmetric arc with power at positive $f_D$. The persistent asymmetric power distribution suggests decreasing DM along $v_{eff}$, and thus decreasing DM with time. *Bottom:* PSR J0621+1002 secondary spectra, showing diffuse, highly asymmetric power of changing signs. This likely suggests variable DM in PSR J0621+1002, but could also arise from the sign of $v_{eff}$ changing from orbital motion; the orbit of PSR J0621+1002 will be investigated in future work.</figcaption></figure>

that the scintillation arcs are sensitive to DM gradients within the scattering screen, not necessarily the total changing electron column which is measured by timing. A detailed analysis comparing high-cadence DM and scintillation arc asymmetries will be valuable, but is beyond the scope of this paper.

# 7 CONCLUSIONS AND FUTURE PROSPECTS

In this paper, we performed the first large sample study of scintillation arcs in MSPs, where of 22 sources regularly observed at LEAP, we observed scintillation arcs in 12. We are able to measure the time-variable arc curvature and scattering in six of these sources, with ~monthly cadence over 5–10 yr.

The scintillation arcs reveal the structure along the dominant scattering screens in these sources, revealing varying phenomena, including compact sources of scattering in PSRs J0613−0200 and J1600−3053, asymmetric distributions of power likely reflecting DM gradients in PSRs J0621+1002 and PSR J1643−1224, and multiple arcs indicating scattering by multiple thin screens along the line of sight in PSR B1937+21. In fitting of the variable scintillation arc curvatures of PSR J0613−0200, we were able to measure $\Omega$, and resolve the sense of *i*, finding a value of *i* consistent with pulsar timing. The screen axis of PSR J0613−0200 changes by tens of degrees over 10 yr (~100 AU), corresponding to visible changes in the extent of scattering.

The time delays measured through scintillation can be compared and combined with other methods, including scattering measured at lower frequencies, measured through sharp features such as giant pulses (Bilous et al. 2015; Main et al. 2017; McKee et al. 2019) or microstructure (Liu et al. 2022). The effects of correlated, variable scattering, as well as correction methods will be assessed using simulations and applied to PTA data in future work. Orbital studies using scintillation can be improved with better understanding of scattering screens, and with more precise measurements of the arc curvature. Studies to date have all been incoherent, attempting to measure the primary scintillation arc without full information of inverted arclets results from interfering pairs of images which arise in highly anisotropic screens. Phase retrieval

techniques such as holography (Walker et al. 2008; Osłowski & Walker 2023), cyclic spectroscopy (Demorest 2011; Walker et al. 2013), and the $\theta-\theta$ transformation (Sprenger et al. 2021; Baker et al. 2022) can greatly increase the precision. In sources with discrete features, the movement of features between observations gives another constraint on the average arc curvature in the time between observations, and can be used as an additional precise constraint (Sprenger et al. 2022). Even without these advanced techniques, improved cadence of observations offer a great improvement, to better fill the annual and orbital planes, to track features between observations, and to track screen changes. Using all available data, including measurements of scintillation velocities in conjunction with arcs, will result in better constraints on pulsar orbits and screens.

Much larger than the effect of scattering variations are the changes in DM, for which there is significant effort to measure (e.g. Jones et al. 2017; Donner et al. 2020; CHIME/Pulsar Collaboration 2021; Tarafdar et al. 2022). Scintillation, scattering, and refractive flux variations are all physically linked, and related to the changing column density of electrons. Detailed mappings of these quantities, as is now being attempted in eclipsing binaries (Lin et al. 2021, 2023), will be valuable, and lead to a more complete physical understanding of the effects of the IISM on pulsar signals.

# ACKNOWLEDGEMENTS

We thank Daniel Reardon for useful comments, and for advice on modelling arc curvature variations. This work was supported by the ERC Advanced Grant ‘LEAP’, Grant Agreement Number 227947 (PI M. Kramer). The European Pulsar Timing Array (EPTA) is a collaboration between European Institutes, namely ASTRON (NL), INAF/Osservatorio Astronomico di Cagliari (IT), the Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), The University of Manchester (UK), The University of Birmingham (UK), The University of Cambridge (UK), and The University of Bielefeld (GER), with an aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves.

The Effelsberg 100-m telescope is operated by the Max-Planck-Institut für Radioastronomie. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a consolidated grant from the STFC in the UK. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Foundation for Radio Astronomy, ASTRON, with support from NWO. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique. The Sardinia Radio Telescope (SRT) is funded by the Department of Universities and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS), and is operated as a National Facility by the National Institute for Astrophysics (INAF). From 2014 March–2016 January, the SRT data were acquired as part of the Astronomical Validation of the SRT. We thus thank the SRT Astronomical Validation Team, and in particular: S. Casu, E. Egron, N. Iacolina, A. Pellizzoni, and A. Trois.

JA acknowledges support from the European Commission (ARGOS–CDS, Grant Agreement number: 101094354), the Niarchos Foundation (SNF), and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of the ‘Science and Society’-action always strive for excellence—Theodoros Papazoglou’ (Project Number: 01431). SC acknowledges the support by the ANR Programme d’Investissement d’Avenir (PIA) under
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Top: PSR J1643−1224 secondary spectra at a similar time of year, showing a very similar distribution of power, with an asymmetric arc with power at positive $f_D$. The persistent asymmetric power distribution suggests decreasing DM along $v_{eff}$, and thus decreasing DM with time. Bottom: PSR J0621+1002 secondary spectra, showing diffuse, highly asymmetric power of changing signs. This likely suggests variable DM in PSR J0621+1002, but could also arise from the sign of $v_{eff}$ changing from orbital motion; the orbit of PSR J0621+1002 will be investigated in future work.</figcaption></figure>

that the scintillation arcs are sensitive to DM gradients within the scattering screen, not necessarily the total changing electron column which is measured by timing. A detailed analysis comparing high-cadence DM and scintillation arc asymmetries will be valuable, but is beyond the scope of this paper.

## 7 CONCLUSIONS AND FUTURE PROSPECTS

In this paper, we performed the first large sample study of scintillation arcs in MSPs, where of 22 sources regularly observed at LEAP, we observed scintillation arcs in 12. We are able to measure the time-variable arc curvature and scattering in six of these sources, with ~monthly cadence over 5–10 yr.

The scintillation arcs reveal the structure along the dominant scattering screens in these sources, revealing varying phenomena, including compact sources of scattering in PSRs J0613−0200 and J1600−3053, asymmetric distributions of power likely reflecting DM gradients in PSRs J0621+1002 and PSR J1643−1224, and multiple arcs indicating scattering by multiple thin screens along the line of sight in PSR B1937+21. In fitting of the variable scintillation arc curvatures of PSR J0613−0200, we were able to measure $\Omega$, and resolve the sense of *i*, finding a value of *i* consistent with pulsar timing. The screen axis of PSR J0613−0200 changes by tens of degrees over 10 yr (~100 AU), corresponding to visible changes in the extent of scattering.

The time delays measured through scintillation can be compared and combined with other methods, including scattering measured at lower frequencies, measured through sharp features such as giant pulses (Bilous et al. 2015; Main et al. 2017; McKee et al. 2019) or microstructure (Liu et al. 2022). The effects of correlated, variable scattering, as well as correction methods will be assessed using simulations and applied to PTA data in future work. Orbital studies using scintillation can be improved with better understanding of scattering screens, and with more precise measurements of the arc curvature. Studies to date have all been incoherent, attempting to measure the primary scintillation arc without full information of inverted arclets results from interfering pairs of images which arise in highly anisotropic screens. Phase retrieval

techniques such as holography (Walker et al. 2008; Osłowski & Walker 2023), cyclic spectroscopy (Demorest 2011; Walker et al. 2013), and the $\theta-\theta$ transformation (Sprenger et al. 2021; Baker et al. 2022) can greatly increase the precision. In sources with discrete features, the movement of features between observations gives another constraint on the average arc curvature in the time between observations, and can be used as an additional precise constraint (Sprenger et al. 2022). Even without these advanced techniques, improved cadence of observations offer a great improvement, to better fill the annual and orbital planes, to track features between observations, and to track screen changes. Using all available data, including measurements of scintillation velocities in conjunction with arcs, will result in better constraints on pulsar orbits and screens.

Much larger than the effect of scattering variations are the changes in DM, for which there is significant effort to measure (e.g. Jones et al. 2017; Donner et al. 2020; CHIME/Pulsar Collaboration 2021; Tarafdar et al. 2022). Scintillation, scattering, and refractive flux variations are all physically linked, and related to the changing column density of electrons. Detailed mappings of these quantities, as is now being attempted in eclipsing binaries (Lin et al. 2021, 2023), will be valuable, and lead to a more complete physical understanding of the effects of the IISM on pulsar signals.

## ACKNOWLEDGEMENTS

We thank Daniel Reardon for useful comments, and for advice on modelling arc curvature variations. This work was supported by the ERC Advanced Grant ‘LEAP’, Grant Agreement Number 227947 (PI M. Kramer). The European Pulsar Timing Array (EPTA) is a collaboration between European Institutes, namely ASTRON (NL), INAF/Osservatorio Astronomico di Cagliari (IT), the Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), The University of Manchester (UK), The University of Birmingham (UK), The University of Cambridge (UK), and The University of Bielefeld (GER), with an aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves.

The Effelsberg 100-m telescope is operated by the Max-Planck-Institut für Radioastronomie. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a consolidated grant from the STFC in the UK. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Foundation for Radio Astronomy, ASTRON, with support from NWO. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique. The Sardinia Radio Telescope (SRT) is funded by the Department of Universities and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS), and is operated as a National Facility by the National Institute for Astrophysics (INAF). From 2014 March–2016 January, the SRT data were acquired as part of the Astronomical Validation of the SRT. We thus thank the SRT Astronomical Validation Team, and in particular: S. Casu, E. Egron, N. Iacolina, A. Pellizzoni, and A. Trois.

JA acknowledges support from the European Commission (ARGOS–CDS, Grant Agreement number: 101094354), the Niarchos Foundation (SNF), and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of the ‘Science and Society’-action always strive for excellence—Theodoros Papazoglou’ (Project Number: 01431). SC acknowledges the support by the ANR Programme d’Investissement d’Avenir (PIA) under
the FIRST–TF network (ANR-10-LABX-48-01) project and the Oscillator IMP project (ANR-11-EQPX-0033-OSC-IMP), and by grants from the Région Bourgogne Franche Comté intended to support the PIA. SC and IC acknowledge financial support from Programme National de Cosmologie and Galaxies (PNCG) and Programme National Hautes Energies (PNHE) funded by CNRS, CEA, and CNES, France. HH acknowledges the support by the Max-Planck Society as part of the ‘LEGACY’ collaboration with the Chinese Academy of Sciences on low-frequency gravitational wave astronomy. JWM gratefully acknowledges support by the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference #CITA 490888–16]. KL is supported by the European Research Council for the ERC Synergy Grant BlackHoleCam under contract number 610058. KJLee gratefully acknowledges support from National Basic Research Program of China, 973 Program, 2015CB857101 and NSFC 11373011. DP gratefully acknowledges financial support from the research grant ‘iPeska’ (P.I. Andrea Possenti) funded under the INAF national call PRIN-SKA/CTA with Presidential Decree 70/2016. WWZ was supported by the National Natural Science Foundation of China Grant No. 11873067 the CAS-MPG LEGACY project and the Strategic Priority Research Program of the Chinese Academy of Sciences Grant No. XDB23000000.

DATA AVAILABILITY

The dynamic spectra used in this work are published on zenodo at
https://doi.org/10.5281/zenodo.7415215.

REFERENCES

Alam M. F. et al., 2021a, ApJS, 252, 4

Alam M. F. et al., 2021b, ApJS, 252, 5

Antoniadis J. et al., 2022, MNRAS, 510, 4873

Arzoumanian Z. et al., 2020, ApJ, 905, L34

Askew J., Reardon D. J., Shannon R. M., 2023, MNRAS, 519, 5086

Baker D., Brisken W., van Kerkwijk M. H., Main R., Pen U.-L., Sprenger T., Wucknitz O., 2022, MNRAS, 510, 4573

Bassa C. G. et al., 2016, MNRAS, 456, 2196

Baumgardt H., Hilker M., 2018, MNRAS, 478, 1520

Bhat N. D. R., Cordes J. M., Camilo F., Nice D. J., Lorimer D. R., 2004, ApJ, 605, 759

Bilous A. V., Pennucci T. T., Demorest P., Ransom S. M., 2015, ApJ, 803, 83

Brisken W. F., Macquart J. P., Gao J. J., Rickett B. J., Coles W. A., Deller A. T., Tingay S. J., West C. J., 2010, ApJ, 708, 232

Chalumeau A. et al., 2022, MNRAS, 509, 5538

Chen S. et al., 2021, MNRAS, 508, 4970

CHIME/Pulsar Collaboration et al., 2021, ApJS, 255, 5

Cordes J. M., 1986, ApJ, 311, 183

Cordes J. M., Lazio T. J. W., 2002, preprint (arXiv:astro-ph/0207156)

Cordes J. M., Rickett B. J., 1998, ApJ, 507, 846

Cordes J. M., Rickett B. J., Stinebring D. R., Coles W. A., 2006, ApJ, 637, 346

Cordes J. M., Shannon R. M., Stinebring D. R., 2016, ApJ, 817, 16

Demorest P. B., 2011, MNRAS, 416, 2821

Desvignes G. et al., 2016, MNRAS, 458, 3341

Ding H. et al., 2023, MNRAS, 519, 4982

Donner J. Y. et al., 2019, A&A, 624, A22

Donner J. Y. et al., 2020, A&A, 644, A153

Fonseca E. et al., 2016, ApJ, 832, 167

Goncharov B. et al., 2021a, MNRAS, 502, 478

Goncharov B. et al., 2021b, ApJ, 917, L19

Hellings R. W., Downs G. S., 1983, ApJ, 265, L39

Hemberger D. A., Stinebring D. R., 2008, ApJ, 674, L37

Hill A. S., Stinebring D. R., Asplund C. T., Berwick D. E., Everett W. B., Hinkel N. R., 2005, ApJ, 619, L171

Hotan A. W., van Straten W., Manchester R. N., 2004, Publ. Astron. Soc. Aust., 21, 302

Jones M. L. et al., 2017, ApJ, 841, 125

Kramer M., Champion D. J., 2013, Class. Quantum Gravity, 30, 224009

Lazarus P., Karuppusamy R., Graikou E., Caballero R. N., Champion D. J., Lee K. J., Verbiest J. P. W., Kramer M., 2016, MNRAS, 458, 868

Levin L. et al., 2016, ApJ, 818, 166

Lin F. X., Lin H.-H., Luo J., Main R., McKee J., Pen U.-L., Simard D., van Kerkwijk M. H., 2021, MNRAS, 508, 1115

Lin F. X., Main R. A., Jow D., Li D. Z., Pen U. L., van Kerkwijk M. H., 2023, MNRAS, 519, 121

Liu Y. et al., 2022, A&A, 664, A116

Main R., van Kerkwijk M., Pen U.-L., Mahajan N., Vanderlinde K., 2017, ApJ, 840, L15

Main R. A. et al., 2020, MNRAS, 499, 1468

Main R. A. et al., 2023, MNRAS, 518, 1086

Mall G. et al., 2022, MNRAS, 511, 1104

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

McKee J. W. et al., 2019, MNRAS, 483, 4784

McKee J. W., Zhu H., Stinebring D. R., Cordes J. M., 2022, ApJ, 927, 99

Oliveira R. A. P. et al., 2022, A&A, 657, A123

Osłowski S., Walker M. A., 2023, MNRAS, 519, 1261

Oswald L. S. et al., 2021, MNRAS, 504, 1115

Perera B. B. P. et al., 2019, MNRAS, 490, 4666

Reardon D. J., Coles W. A., 2023, MNRAS, 521, 6392

Reardon D. J., Coles W. A., Hobbs G., Ord S., Kerr M., Bailes M., Bhat N. D. R., Venkatraman Krishnan V., 2019, MNRAS, 485, 4389

Reardon D. J. et al., 2020, ApJ, 904, 104

Rickett B. J. et al., 2014, ApJ, 787, 161

Simard D., Pen U. L., Marthi V. R., Brisken W., 2019, MNRAS, 488, 4963

Smits R. et al., 2017, Astron. Comput., 19, 66

Sprenger T., Wucknitz O., Main R., Baker D., Brisken W., 2021, MNRAS, 500, 1114

Sprenger T., Main R., Wucknitz O., Mall G., Wu J., 2022, MNRAS, 515, 6198

Stappers B., Kramer M., 2011, American Astronomical Society, AAS Meeting #217, id.234.05; Bulletin of the American Astronomical Society, Vol. 43, 2011

Stinebring D. R., McLaughlin M. A., Cordes J. M., Becker K. M., Goodman J. E. E., Kramer M. A., Sheckard J. L., Smith C. T., 2001, ApJ, 549, L97

Stinebring D. R. et al., 2022, ApJ, 941, 34

Tarafdar P. et al., 2022, Publ. Astron. Soc. Aust., 39, e053

van Straten W., Bailes M., 2011, Publ. Astron. Soc. Aust., 28, 1

Verbiest J. P. W., Shaifullah G. M., 2018, Class. Quantum Gravity, 35, 133001

Walker M. A., Demorest P. B., van Straten W., 2013, ApJ, 779, 99

Walker M. A., Koopmans L. V. E., Stinebring D. R., van Straten W., 2008, MNRAS, 388, 1214

Walker M. A., Melrose D. B., Stinebring D. R., Zhang C. M., 2004, MNRAS, 354, 43

Walker K., Reardon D. J., Thrane E., Smith R., 2022, ApJ, 933, 16

Wu Z. et al., 2022, A&A, 663, A116

Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29

APPENDIX: SECONDARY SPECTRA

Here, we show the secondary spectra of all of the sources with resolv-
able scintillation with LEAP, as well as the Effelsberg observations
of PSRs J0613–0200 and J1643–1224.
the FIRST-TF network (ANR-10-LABX-48-01) project and the Oscillator IMP project (ANR-11-EQPX-0033-OSC-IMP), and by grants from the Région Bourgogne Franche Comté intended to support the PIA. SC and IC acknowledge financial support from Programme National de Cosmologie and Galaxies (PNCG) and Programme National Hautes Energies (PNHE) funded by CNRS, CEA, and CNES, France. HH acknowledges the support by the Max-Planck Society as part of the 'LEGACY' collaboration with the Chinese Academy of Sciences on low-frequency gravitational wave astronomy. JWM gratefully acknowledges support by the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference #CITA 490888-16]. KL is supported by the European Research Council for the ERC Synergy Grant BlackHoleCam under contract number 610058. KJLee gratefully acknowledges support from National Basic Research Program of China, 973 Program, 2015CB857101 and NSFC 11373011. DP gratefully acknowledges financial support from the research grant 'iPeska' (P.I. Andrea Possenti) funded under the INAF national call PRIN-SKA/CTA with Presidential Decree 70/2016. WWZ was supported by the National Natural Science Foundation of China Grant No. 11873067 the CAS-MPG LEGACY project and the Strategic Priority Research Program of the Chinese Academy of Sciences Grant No. XDB23000000.

## DATA AVAILABILITY

The dynamic spectra used in this work are published on zenodo at https://doi.org/10.5281/zenodo.7415215.

## REFERENCES

* Alam M. F. et al., 2021a, ApJS, 252, 4

* Alam M. F. et al., 2021b, ApJS, 252, 5

* Antoniadis J. et al., 2022, MNRAS, 510, 4873

* Arzoumanian Z. et al., 2020, ApJ, 905, L34

* Askew J., Reardon D. J., Shannon R. M., 2023, MNRAS, 519, 5086

* Baker D., Brisken W., van Kerkwijk M. H., Main R., Pen U.-L., Sprenger T., Wucknitz O., 2022, MNRAS, 510, 4573

* Bassa C. G. et al., 2016, MNRAS, 456, 2196

* Baumgardt H., Hilker M., 2018, MNRAS, 478, 1520

* Bhat N. D. R., Cordes J. M., Camilo F., Nice D. J., Lorimer D. R., 2004, ApJ, 605, 759

* Bilous A. V., Pennucci T. T., Demorest P., Ransom S. M., 2015, ApJ, 803, 83

* Brisken W. F., Macquart J. P., Gao J. J., Rickett B. J., Coles W. A., Deller A. T., Tingay S. J., West C. J., 2010, ApJ, 708, 232

* Chalumeau A. et al., 2022, MNRAS, 509, 5538

* Chen S. et al., 2021, MNRAS, 508, 4970

* CHIME/Pulsar Collaboration et al., 2021, ApJS, 255, 5

* Cordes J. M., 1986, ApJ, 311, 183

* Cordes J. M., Lazio T. J. W., 2002, preprint (arXiv:astro-ph/0207156)

* Cordes J. M., Rickett B. J., 1998, ApJ, 507, 846

* Cordes J. M., Rickett B. J., Stinebring D. R., Coles W. A., 2006, ApJ, 637, 346

* Cordes J. M., Shannon R. M., Stinebring D. R., 2016, ApJ, 817, 16

* Demorest P. B., 2011, MNRAS, 416, 2821

* Desvignes G. et al., 2016, MNRAS, 458, 3341

* Ding H. et al., 2023, MNRAS, 519, 4982

* Donner J. Y. et al., 2019, A&A, 624, A22

* Donner J. Y. et al., 2020, A&A, 644, A153

* Fonseca E. et al., 2016, ApJ, 832, 167

* Goncharov B. et al., 2021a, MNRAS, 502, 478

* Goncharov B. et al., 2021b, ApJ, 917, L19

* Hellings R. W., Downs G. S., 1983, ApJ, 265, L39

* Hemberger D. A., Stinebring D. R., 2008, ApJ, 674, L37

* Hill A. S., Stinebring D. R., Asplund C. T., Berwick D. E., Everett W. B., Hinkel N. R., 2005, ApJ, 619, L171

* Hotan A. W., van Straten W., Manchester R. N., 2004, Publ. Astron. Soc. Aust., 21, 302

* Jones M. L. et al., 2017, ApJ, 841, 125

* Kramer M., Champion D. J., 2013, Class. Quantum Gravity, 30, 224009

* Lazarus P., Karuppusamy R., Graikou E., Caballero R. N., Champion D. J., Lee K. J., Verbiest J. P. W., Kramer M., 2016, MNRAS, 458, 868

* Levin L. et al., 2016, ApJ, 818, 166

* Lin F. X., Lin H.-H., Luo J., Main R., McKee J., Pen U.-L., Simard D., van Kerkwijk M. H., 2021, MNRAS, 508, 1115

* Lin F. X., Main R. A., Jow D., Li D. Z., Pen U. L., van Kerkwijk M. H., 2023, MNRAS, 519, 121

* Liu Y. et al., 2022, A&A, 664, A116

* Main R., van Kerkwijk M., Pen U.-L., Mahajan N., Vanderlinde K., 2017, ApJ, 840, L15

* Main R. A. et al., 2020, MNRAS, 499, 1468

* Main R. A. et al., 2023, MNRAS, 518, 1086

* Mall G. et al., 2022, MNRAS, 511, 1104

* Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

* McKee J. W. et al., 2019, MNRAS, 483, 4784

* McKee J. W., Zhu H., Stinebring D. R., Cordes J. M., 2022, ApJ, 927, 99

* Oliveira R. A. P. et al., 2022, A&A, 657, A123

* Osłowski S., Walker M. A., 2023, MNRAS, 519, 1261

* Oswald L. S. et al., 2021, MNRAS, 504, 1115

* Perera B. B. P. et al., 2019, MNRAS, 490, 4666

* Reardon D. J., Coles W. A., 2023, MNRAS, 521, 6392

* Reardon D. J., Coles W. A., Hobbs G., Ord S., Kerr M., Bailes M., Bhat N. D. R., Venkatraman Krishnan V., 2019, MNRAS, 485, 4389

* Reardon D. J. et al., 2020, ApJ, 904, 104

* Rickett B. J. et al., 2014, ApJ, 787, 161

* Simard D., Pen U. L., Marthi V. R., Brisken W., 2019, MNRAS, 488, 4963

* Smits R. et al., 2017, Astron. Comput., 19, 66

* Sprenger T., Wucknitz O., Main R., Baker D., Brisken W., 2021, MNRAS, 500, 1114

* Sprenger T., Main R., Wucknitz O., Mall G., Wu J., 2022, MNRAS, 515, 6198

* Stappers B., Kramer M., 2011, American Astronomical Society, AAS Meeting #217, id.234.05; Bulletin of the American Astronomical Society, Vol. 43, 2011

* Stinebring D. R., McLaughlin M. A., Cordes J. M., Becker K. M., Goodman J. E. E., Kramer M. A., Scheckard J. L., Smith C. T., 2001, ApJ, 549, L97

* Stinebring D. R. et al., 2022, ApJ, 941, 34

* Tarafdar P. et al., 2022, Publ. Astron. Soc. Aust., 39, e053

* van Straten W., Bailes M., 2011, Publ. Astron. Soc. Aust., 28, 1

* Verbiest J. P. W., Shaifullah G. M., 2018, Class. Quantum Gravity, 35, 133001

* Walker M. A., Demorest P. B., van Straten W., 2013, ApJ, 779, 99

* Walker M. A., Koopmans L. V. E., Stinebring D. R., van Straten W., 2008, MNRAS, 388, 1214

* Walker M. A., Melrose D. B., Stinebring D. R., Zhang C. M., 2004, MNRAS, 354, 43

* Walker K., Reardon D. J., Thrane E., Smith R., 2022, ApJ, 933, 16

* Wu Z. et al., 2022, A&A, 663, A116

* Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29

## APPENDIX: SECONDARY SPECTRA

Here, we show the secondary spectra of all of the sources with resolvable scintillation with LEAP, as well as the Effelsberg observations of PSRs J0613–0200 and J1643–1224.
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Secondary spectra of all observations of PSR J0613−0200, including targeted Effelsberg observations and recent scintillation-stream EPTA data.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure A1.** Secondary spectra of all observations of PSR J0613−0200, including targeted Effelsberg observations and recent scintillation-stream EPTA data.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure A2.** Secondary spectra of all observations of PSR J1643−1224, including recent scintillation-stream EPTA data.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure A3.** Secondary spectra of all observations of PSR J0621+1002.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure A2.** Secondary spectra of all observations of PSR J1643−1224, including recent scintillation-stream EPTA data.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure A3.** Secondary spectra of all observations of PSR J0621+1002.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure A4. Secondary spectra of all observations of PSR J1600−3053.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure A5. Secondary spectra of all observations of PSR J1918−0642.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure A4. Secondary spectra of all observations of PSR J1600−3053.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure A5. Secondary spectra of all observations of PSR J1918−0642.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure A6.** Secondary spectra of all observations of PSR B1937+21.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.

<footnote n="1">¹Max Planck Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="2">²Department of Physics, Institute of Astrophysics, FORTH, University of Crete, Voutes, University Campus, GR-71003 Heraklion, Greece</footnote>

<footnote n="3">³Argelander Institut für Astronomie, Auf dem Hügel 71, D-53121 Bonn, Germany</footnote>

<footnote n="4">⁴Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, P. R. China</footnote>

<footnote n="5">⁵Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E CNRS-Université d'Orléans, F-45071 Orléans, France</footnote>

<footnote n="6">⁶Station de radioastronomie de Nançay, Observatoire de Paris, PSL Research University, CNRS/INSU F-18330 Nançay, France</footnote>

<footnote n="7">⁷Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="8">⁸CAS Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China</footnote>

<footnote n="9">⁹Fakultät für Physik, Universität Bielefeld, Postfach 100131, D-33501 Bielefeld, Germany</footnote>

<footnote n="10">¹⁰Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada</footnote>

<footnote n="11">¹¹Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada</footnote>

<footnote n="12">¹²E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK</footnote>

<footnote n="13">¹³Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK</footnote>

<footnote n="14">¹⁴INAF – Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy</footnote>

<footnote n="15">¹⁵ASTRON, the Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, NL-7991 PD Dwingeloo, the Netherlands</footnote>

<footnote n="16">¹⁶Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, NL-6500 GL Nijmegen, the Netherlands</footnote>

<footnote n="17">¹⁷Dipartimento di Fisica, Università di Cagliari, S.P. Monserrato-Sestu Km 0,700, I-09042 Monserrato (CA), Italy</footnote>

<footnote n="18">¹⁸National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Rd, Chaoyang District, Beijing 100012, P. R. China</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure A6.** Secondary spectra of all observations of PSR B1937+21.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.

<footnote n="1">¹Max Planck Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="2">²Department of Physics, Institute of Astrophysics, FORTH, University of Crete, Voutes, University Campus, GR-71003 Heraklion, Greece</footnote>

<footnote n="3">³Argelander Institut für Astronomie, Auf dem Hügel 71, D-53121 Bonn, Germany</footnote>

<footnote n="4">⁴Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, P. R. China</footnote>

<footnote n="5">⁵Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E CNRS-Université d'Orléans, F-45071 Orléans, France</footnote>

<footnote n="6">⁶Station de radioastronomie de Nançay, Observatoire de Paris, PSL Research University, CNRS/INSU F-18330 Nançay, France</footnote>

<footnote n="7">⁷Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="8">⁸CAS Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China</footnote>

<footnote n="9">⁹Fakultät für Physik, Universität Bielefeld, Postfach 100131, D-33501 Bielefeld, Germany</footnote>

<footnote n="10">¹⁰Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada</footnote>

<footnote n="11">¹¹Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada</footnote>

<footnote n="12">¹²E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK</footnote>

<footnote n="13">¹³Centre of Excellence for Data Science, Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK</footnote>

<footnote n="14">¹⁴INAF – Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy</footnote>

<footnote n="15">¹⁵ASTRON, the Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, NL-7991 PD Dwingeloo, the Netherlands</footnote>

<footnote n="16">¹⁶Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, NL-6500 GL Nijmegen, the Netherlands</footnote>

<footnote n="17">¹⁷Dipartimento di Fisica, Università di Cagliari, S.P. Monserrato-Sestu Km 0,700, I-09042 Monserrato (CA), Italy</footnote>

<footnote n="18">¹⁸National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Rd, Chaoyang District, Beijing 100012, P. R. China</footnote>