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

We present the first large sample of scintillation arcs in millisecond pulsars, 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 timescales over a 10-year timespan. 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 multi-screen 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. Unmodeled 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 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.


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;Goncharov et al. 2021b;Chen et al. 2021;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  2 dispersive delay (where  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   , with  ∼ 4.0 ± 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 et al. 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 ∼ 20% of sources in Stinebring et al. 2022).
The majority of the brightest known radio pulsars are isolated and slowly rotating (Manchester et al. 2005), which due to high 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 millisecond pulsars.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;Walker et al. 2022;McKee 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 re-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 HII region (Mall et al. 2022).
In this paper, we present the first large sample of scintillation arcs of MSPs, observed over the last 10 years 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 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

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  (, ) 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 (   ,   ) = | Ĩ (   ,   )| 2 , referred to as the 'secondary spectrum', the conjugate variable of time   ≡  D is related to the Doppler shift between deflected paths, and depends on the angles  of two deflected paths ( ) as and the conjugate variable of frequency   ≡  is related to the delay between image pairs, described as . (2) The effective distance  eff and effective velocity  eff depend on the relative distances and velocities of the pulsar ( psr ,  psr ), screen ( scr ,  scr ), and Earth ( ⊕ ), as where the fractional screen distance from the pulsar is defined as  ≡ 1 −  scr / psr .When one of the two scattered angles is 0 (i.e. the theoretical undeflected line of sight), the common dependence of  D and  on the observed angle  results in a parabolic distribution of power The proportionality constant, or 'arc-curvature'  depends on the relative distances and velocities of the pulsar, screen, and Earth, as where  is the angle between  eff and .
Throughout this paper, we work with the distance weighted effective velocity , which rearranges  to separate the unknown values, and is proportional to 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  is that it does not diverge as | eff, | → 0 and is independent of the observing frequency.
In the absence of arcs, a characteristic curvature  can be estimated from the time and frequency scale from the 2D autocorrelation function (ACF) of scintillation.The scintillation bandwidth  s is defined as the half-width at half-maximum of the ACF in frequency, and the scintillation timescale  s is defined as the 1/ 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 timescale as These relations are approximate, as the prefactors are model dependent.Then we can estimate  as A measure of either the arc curvature , or the scintillation timescale and bandwidth, are then a direct measure of  eff, / √︁  eff .In this paper, we measure  from the secondary spectrum wherever possible; scintillation velocities derived from  s and  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).

Timescale 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 timescale for features to pass through the secondary spectrum, which sets a timescale for scattering delays to correlate.This timescale (effectively the same as the traditional "refractive timescale" 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  0 with corresponding size length  r , then the timescale is related to the arc curvature and the maximum delay  0 as In Sprenger et al. (2022), all of the observables were expressed as a function of the feature movement in √ , related to the distance weighted effective velocity as which can also be used to obtain equation (10).The relevant  0 value can either be the highest  that features are visible, which gives the timescale for features to completely pass through the secondary spectrum, or the 1/ value of .

LEAP
The Large European Array for Pulsars (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 minutes.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 bandwdith.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.

Long Targeted Observations
Several of the sources have scintillation which is barely resolved in a ∼ 30 − 60 minute observation, and with scintles 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 hour 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 7-beam receiver, recorded in 25 MHz subbands with a usable bandwidth of 1250-1450 MHz and saved in † format.

EPTA observations
Since March 2021, Effelsberg observations for the EPTA record a separate parallel data stream suitable for scintillation studies, using the Effelsberg Direct Digitisation system, folding with 64000 channels across 1200-1600 MHz (sensitive to delays of  < 80 µ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.
In addition, we also investigate sources where scintillation is marginally resolved by LEAP, supplementing this study with longer 2 − 3 hour 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 March−June 2020, 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 timespan 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 dataset allowed for studies of slightly less scattered sources.

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 baseband 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  D and  to capture the full extent of arcs in the secondary spectrum.The sub-bands were combined in frequency using psradd in the psrchive software package (Hotan et al. 2004), and polarisations 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% of the dynamic spectra are tapered with a Hanning window, before forming the secondary spectrum.
As the Effelsberg observations cover a larger bandwidth, we use the 'NuT transform' instead of a direct FFT 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  = / ref in every channel with frequency .This prevents smearing of scintillation arcs owing to the  ∝  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  ref = 1400 MHz.We note that this transform can lead to some artefacts in the secondary spectrum, as power on the  = 0 axis (arising from e.g.RFI, pulseto-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 Figure 1, and their associated secondary spectra are shown in Figure 2. The secondary spectra of every observation used in this work are included in the Appendix.

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  D axis of the secondary spectrum is mapped to  D,norm =  D √︁  ref / (where we set the arbitrary reference time delay  ref =  max throughout).This transformation effectively stretches the  D axis, mapping parabolae to vertical lines of constant  D,norm .Then  can be identified by finding peaks in (  D,norm ) after performing a weighted sum over .
Arcs blend together at low values of , becoming more clearly demarcated at high values of .As the optimal range of  to sum over and  D,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  to sum to form (  D,norm ) is adjusted per source, and peaks in  D,norm can be given as initial guesses for a least-squares fit of a parabola to the local region of (  D,norm ), and the arc curvature is given as An example is shown in Figure 3.

Time delays from secondary spectra
Secondary spectra express the power of scintillation in terms of the conjugate variables of time and frequency respectively,  D and .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  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  can then be estimated as For each source, a range in  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  D ± 30 mHz.The integrated profile of  against maximum  max ≡  was examined in all cases; the value and the error on  were estimated by the mean and standard deviation of the profile once it plateaus, taken in the range of 3/4 <  < .Examples of secondary spectra and their associated profiles of  () are shown in Figure 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 timescale and bandwidth in a more traditional way, through the 2D autocorrelation function (ACF) of the dynamic spectrum (Δ) = ( * )(Δ).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  = /2 s , where the model-dependent constant , 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 Figure 5.
The precision on the measurements on the ACF is high, and the true error is dominated by the fact that there are finite scintles within the observation.For an observation of duration  obs with bandwidth BW, the fractional error is given by where the filling factor is assumed to be   ≈ 0.2 (Cordes 1986).

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

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  should arise only from Earth's motion.Despite this, they show a range of interesting behaviour owing to dynamic screens.

PSR B1937+21
In low frequency observations of PSR B1937+21, there is previous evidence of a broad scintillation arc (Walker et al. 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 Figure 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 Figure 6.The curvature of the secondary and tertiary screens appear roughly consistent whenever they do reappear, suggesting stability over the 8 years of observations.However, due to the difficulty of unambiguously identifying scattering screens, as well as the large distance uncertainty of the source, 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 timevariable scattering timescale (McKee et al. 2019).Measurements of giant pulse scattering are direct, and unaffected by multi-screen effects.Reassuringly, the trends of  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 timescale through high-cadence DM measurements (CHIME/Pulsar Collaboration et al. 2021).

PSR B1821−24A
PSR B1821−24A is a millisecond pulsar 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 signal-to-noise, 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 intra-cluster gas.Scattering in the intracluster gas would result in very large values of ; approximating  pl with the core radius of   ≈ 0.37 pc, and  pl with the velocity dispersion of ≈ 11 km s −1 (Oliveira et al. 2022;Baumgardt & Hilker 2018), would result in  ≈ 600 km s kpc −0.5 , much greater than our measured  = 36 ± 16 km s kpc −0.5 .Comparably, our measured value of  is easily compatible with typical Earth and screen velocities, at a screen of ∼ 1 kpc.Additionally, NE2001 predicts of  = 2.1 µs (Cordes & Lazio 2002), comparable to our measured  = 1.0 ± 0.5 µs.

PSR J0613−0200
The scintillation of this pulsar was previously studied in Main et al. (2020), using LEAP and Effelsberg data from 2013−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−2022, including EPTA scintillation data in the past year, the source has experienced heightened scattering fluctuations, with  decreasing to its lowest state in late 2020, and rising to its highest yet 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 0.2µs in extent along the main parabola, which can be tracked between observations during our high-cadence Effelsberg campaign in March−June 2020, covered in section 6.3.

PSR J0621+1002
This pulsar shows well resolved, low curvature (i.e.large ) 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 timescale of months.This likely reflects large, timevariable 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 Sec. 6.2).The advance of periastron  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 - 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.

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-07-25 onwards, there is only a small, featureless concentration of power at low  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 .Time delays (Top), and measurements of  (Bottom) for the 6 sources with fully resolvable scintillation at LEAP.For PSR J0613−0200, the fit of the annual and orbital variations of  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 parabolae can be identified.(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; 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.

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 HII region.PSR J1643−1224 was regularly observed with LEAP from 2012−2018, and is observed with Effelsberg as part of EPTA observations, extending our dataset from 2012-2022.The properties of scintillation arcs in the last year of EPTA observations (i.e. both  and their extent in ) 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 years, the time delays  are seen to be variable.The scattering measured by scintillation decreases from ∼ 5 − 2.5 µ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 scintles of all of the LEAP sources, with   ∼ 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  ∼ 18 mas, while the resolution of LEAP's longest baselines is / ≈ 21cm/1200 km ≈ 35 mas.For a pulsar scattered much beyond this, the angle  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 PSR J1643−1224 to be angularly broadened to  = 3.65 ± 0.43 mas using the Very Large Baseline Array, and confirm the association of the dominant scattering screen with Sh 2-27.

PSR J1918−0642
The scintillation timescale is ∼10 minutes at its shortest, sufficiently short to reveal a faint arc in LEAP observations.At its slowest, the scintillation timescale is greater than the observation length, appearing in the secondary spectrum as power along the  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 A5). Coincidentally, LEAP observations taken until January 2016 were only 20 minutes 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.

Time Delays
The measured values of  are shown in Figure 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 timescales 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  are correlated on the timescale of years in several sources, and variations on similar timescales 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, pulsars J0613−0200, J1600−3053, and J1918−0642 are among the top 5 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.

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  = 0.58 ± 0.10 during the period of increased scattering in 2013, and  = 0.62 ± 0.06 afterwards, with screen axis  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,   = 1.828(5)mas/yr,   = −10.35(1)mas/yr,  psr = 1.11 ± 0.05 kpc (Perera et al. 2019).We fit the distance weighted effective velocities  with a model of an anisotropic scattering screen, We compare several different models, detailed in the following sections.

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  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 , and 2D velocity  scr,  ,  scr,  are free parameters and constant over the full duration.In the second model, we allow , , and   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 2 jumps, one corresponding to a region of enhanced DM and scattering.
In addition, we fit for the orbital inclination  and angle of nodes Ω, and we include white noise parameters F and Q as free parameters of .Bayesian Information Criteria in fitting  variations of PSR J0613−0200, for models with  jumps in the screen orientation .Models with fewer than 2 jumps result in a poor fit, while models with more than than 2 jumps result in no additional improvement and result in increased BIC due to having more parameters.
the fit, such that the scaled errors are  corr = √︁ (F × 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 timescale   (e.g.Askew et al. 2023).
The first model, allowing for all screen parameters to vary in each jump, results in  = 483.4,while the second model with a screen at a fixed distance, changing only the orientation has  = 442.5.Both models have almost identical white noise parameters  and , suggesting that they fit the data comparably well, but the first model 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 fit values of these models are tabulated in Table 3.

Orbital Constraints
Scintillation arcs provide a way to measure resolve the ambiguity in the sense of the inclination, i.e.  < 90 • , or  > 90 • .We obtain a fit orbital inclination of  = 58 ± 4 • is consistent with one value from timing of  timing = 68 • +7 −10 (Fonseca et al. 2016), but mildly inconsistent with the IPTA value of  timing = 70 ± 3 • by < 2.We obtain the first measurement of Ω = 124 ± 4 • .Additionally, we compare to the fit restricting  > 90.This results in values of  = 116.4± 8.4, and Ω = 274.0± 12.8, but is disfavoured, with  = 523 The best fit model is overlaid on the timeseries of  measurements in Figure 6, and the decomposition to annual and orbital velocity is shown in Figure 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  and Ω.The times of the jumps both correspond to regions of changing  , suggesting that both of these effects trace physical changes of the screen.

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  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  = 0.71 ± 0.02 is smaller than the value of  = 0.62 ± 0.06.This is a result of the different pulsar distance used,  psr = 780 ± 80 pc from Desvignes et al. (2016), and both measurements result in a consistent screen distance.

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 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 Figure 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 .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 surrounding  = 0 which is difficult to track between observations.This contains the bulk of the power, dominating the changing values of  .

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  eff creates a refractive shift (Cordes et al. 2006;Rickett et al. 2014), DM ( 16) leading to a new zero-point of a secondary spectrum offset by   as: 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 Figure 11).The DM curve of PSR J1643−1224 from the NANOGrav 12.5 year data release shows a persistent downwards trend in DM of ΔDM ∼ 10 −3 pc cm −3 year −1 from ∼ 2013 − 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 timescale of months, which may suggest rapidly varying DM.This explanation is plausible as observations with 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 timescales (Donner et al. 2020).However, we note 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.

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 6 of these sources, with ∼monthly cadence over 5 − 10 years.
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 Ω, and resolve the sense of , finding a value of  consistent with pulsar timing.The screen axis of PSR J0613−0200 changes by tens of degrees over 10 years (∼ 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  −  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 et al. 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(Lin et al. , 2022)), will be valuable, and lead to a more complete physical understanding of the effects of the IISM on pulsar signals.

A 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 1 .
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.

Figure 2 .Figure 3 .
Figure 2. Panorama of scintillation arcs, from the corresponding dynamic spectra in Figure 1.Scintillation arcs are seen in many sources, with highly varied extents in  and  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 parabolae in others (e.g.PSRs J0621+1002, J1643−1224), and multiple parabolae in PSR B1937+21.The results are described in more detail in section 5.The intensity maps are logarithmic, extending 3 orders of magnitude in most cases, or 1 in B1821−24A

Figure 4 .Figure 5 .
Figure 4. Secondary Spectra (images), and the associated estimates of the scattering tail  ( ) by summing over  D along the arc (top panels).Several of the sources show features at large time delays, beyond a simple exponential tail.
Figure6.Time delays (Top), and measurements of  (Bottom) for the 6 sources with fully resolvable scintillation at LEAP.For PSR J0613−0200, the fit of the annual and orbital variations of  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 parabolae can be identified.
Figure7.Bayesian Information Criteria in fitting  variations of PSR J0613−0200, for models with  jumps in the screen orientation .Models with fewer than 2 jumps result in a poor fit, while models with more than than 2 jumps result in no additional improvement and result in increased BIC due to having more parameters.

Figure 8 .
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 Ω and  are almost perpendicular in 2021−2022, the changing screen orientation ensures that the orbital modulation is seen.

√Figure 9 .
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 √ vs.  D from the corresponding secondary spectra, made as described in Section 6.3, and shifted by the predicted movement between observations.The value √  ∝  is a proxy for the image positions; features connected vertically between observations suggest persistent scattering at regions of fixed .

Figure 10 .
Figure 10.Same as 9, but for a series of LEAP PSR J1600−3053 observations showing moving features.

Figure 11 .
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  D .The persistent asymmetric power distribution suggests decreasing DM along  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  eff changing from orbital motion; the orbit of PSR J0621+1002 will be investigated in future work.
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 no.610058.K.J.Lee 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.

Table 1 .
(Manchester et al. 2005)n our sample.The pulsar distances shown with 1 errorbars are all timing parallax measurements compiled by the IPTA inPerera et al. (2019), while the distances shown without errors are DM distance estimates from YMW16 modelYao et al. (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.

Table 2 .
Summary will be time-variable owing to the annual motion, pulsar binary motion, and changing screen properties.The values of   are computed as in Section 2, and represent the time for a fixed feature to pass from −  to  given the average value of  .
of quantities derived from the arc curvatures, and from ACF fits.The values of  ≡  eff, / √  eff and  r quoted here are an average, representative value, and both

Table 3 .
Main et al. (2020)ng arc curvature variations of PSR J0613−0200, and comparison toMain et al. (2020).The model assumes an anisotropic scattering screen at a fixed distance, but allowing for 2 jumps in the screen orientation , 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.