The Thousand-Pulsar-Array programme on MeerKAT -- X. Scintillation arcs of 107 pulsars

We present the detection of 107 pulsars with interstellar scintillation arcs at 856--1712\,MHz, observed with the MeerKAT Thousand Pulsar Array Programme. Scintillation arcs appear to be ubiquitous in clean, high S/N observations, their detection mainly limited by short observing durations and coarse frequency channel resolution. This led the survey to be sensitive to nearby, lightly scattered pulsars with high effective velocity -- from a large proper motion, a screen nearby the pulsar, or a screen near the Earth. We measure the arc curvatures in all of our sources, which can be used to give an estimate of screen distances in pulsars with known proper motion, or an estimate of the proper motion. The short scintillation timescale in J1731$-$4744 implies a scattering screen within 12\,pc of the source, strongly suggesting the association between this pulsar and the supernova remnant RCW 114. We measure multiple parabolic arcs of 5 pulsars, all of which are weakly scintillating with high proper motion. Additionally, several sources show hints of inverted arclets suggesting scattering from anisotropic screens. Building on this work, further targeted MeerKAT observations of many of these pulsars will improve understanding of our local scattering environment and the origins of scintillation; annual scintillation curves would lead to robust screen distance measurements, and the evolution of arclets in time and frequency can constrain models of scintillation.


Introduction
Observations of the radio emission of pulsars exhibits scintillation, an interference pattern owing to electron density fluctuations in the ionized interstellar medium (IISM). Recent studies of several pulsars have revealed "scintillation arcs" -a parabolic distribution of power in the 2D power spectrum of scintillation commonly referred to as the secondary spectrum. Discovered by Stinebring et al. (2001), scintillation arcs have proven a powerful probe of the distribution and geometry of scattering structures in the IISM, as well as a probe of pulsar emission regions (e.g. , and binary orbits (e.g. . In addition to a main parabola, many sources show 'arclets', inverted parabolae of the same curvature stemming from apices on the main parabola, which persist over timescales of months. This is difficult to explain with a picture of turbulent, isotropic, volume filling IISM, instead suggesting compact scattering structures along highly anisotropic thin 'screens' (Walker et al. 2004;Cordes et al. 2006). ★ Email: ramain@mpifr-bonn.mpg.de Some scattering screens have been seen to persist over many years McKee et al. 2022;Mall et al. 2022), suggesting the extent of these screens to be 1000 AU. Moreover, arclets are seen to track between observations and across frequency (Hill et al. 2005;Sprenger et al. 2022), consistent with compact scattering structures at fixed angular positions.
It is currently unknown precisely how and where these screens originate, although screen associations have placed screens in supernova remnants (Cordes et al. 2004;Yao et al. 2021), HII regions , potentially near the edge of the local bubble Sprenger et al. 2022), or in ionized gas surrounding hot stars (Walker et al. 2017). The electron densities required to scatter radio light at such large angles are difficult to reconcile with pressure balance in the IISM, leading to models of scattering in thin over-or under-dense plasma sheets seen edge on (Pen & Levin 2014), or in magnetic filaments (Gwinn 2019).
While most studies of scintillation arcs to date have focused on specific sources of interest, wider surveys are useful in determining the behaviour and occurrence of scintillation arcs, as well as the distribution of screens in the Milky Way. Recent surveys suggest scintillation arcs are common, if not prevalent. Stinebring et al. (2022) found scintillation arcs in 19 out of 22 detected pulsars in a survey of bright, low Dispersion Measure (DM) pulsars using the Green Bank Telescope and the Arecibo Observatory. Wu et al. (2022) detect scintillation arcs in 9 out of 31 sources with LOFAR at 110 − 190 MHz; scintillation becomes finer in both time and frequency, requiring very fine channelization and time sampling, resulting in harsh signal-to-noise limitations in detecting arcs.
The Thousand Pulsar Array (TPA), is part of the Large Survey Project 'MeerTime' on the MeerKAT telescope (Bailes et al. 2020), which has observed 1270 pulsars in the southern sky in a wide band from 856−1712 MHz (details of the TPA in Johnston et al. 2020). With the wide band, and the massive upgrade in sensitivity of MeerKAT compared to Parkes, the TPA can be used as a representative survey of scintillation properties of pulsars in the southern sky.
In this paper, we present the first results of scintillation using the TPA, focusing on detection and measurements of scintillation arcs of 107 pulsars. Section 2 revisits the necessary theory of scintillation arcs, Section 3 describes the data, the reduction, the sample and curvature measurements, Section 4 describes our results, and Section 5 includes ramifications of our results and prospects for future work.

Background
In this section, we briefly review the relevant theory of scintillation arcs, and refer the reader to Walker et al. (2004); Cordes et al. (2006) for more detailed treatment. Additionally, in section 2.2, we outline limiting cases, and the approximations necessary to estimate the screen distance or pulsar proper motion from a single measure of a scintillation arc curvature.

Scintillation Arcs
The dynamic spectrum ( , ) = | ( , )| 2 expresses the pulse intensity as a function of time and frequency. The 2D power spectrum, or secondary spectrum | ( , )| 2 expresses the intensity in terms of its conjugate variables and which are related to the Doppler shift and phase delay. In the limit of strong scattering, the secondary resulting interference pattern can be described through the stationary phase approximation as a discrete set of deflected images. The secondary spectrum then contains the interference of all pairs of images where is the observing wavelength, and the effective distance eff and effective velocity eff are defined as where p , l and p , l are the distances and velocities of the pulsar and the scattering screen, respectively, and ⊕ is the Earth's velocity. In weak scattering, the main pulsar image is much brighter than subimages, resulting in a sharp forward parabola. With the main pulsar image at = 0, the common dependence of and on results in a relation of In strong scattering from a highly anisotropic screen, subimages are of comparable brightness, and the cross-terms result in inverted parabolae ('arclets') of the same curvature with apices located on the main arc. Following Mall et al. (2022) and Sprenger et al. (2022), we choose not to work directly with , but rather with separating the measurements and constants from physical quantities of interest, where has the added benefit of approaching 0 rather than diverging when eff, | | → 0.

Approximations
For an isolated pulsar, contains the degenerate combination of 6 unknown parameters; the pulsar's distance and 2D velocity , the screen's distance , and either the parallel velocity , | | and axis of anisotropy for a 1D screen, or the 2D velocity for an isotropic 2D screen. Even in pulsars where the proper motion and distance are independently well measured, a single measure of is insufficient to measure the screen properties, where this degeneracy is typically either broken by using annual variations of or Mall et al. 2022;McKee et al. 2022;Sprenger et al. 2022), or through VLBI or multi-telescope dynamic spectra (Brisken et al. 2010;Simard et al. 2019). However, one does not typically have annual scintillation curves, or scintillation cross-spectra; under certain limits reviewed here, one can obtain meaningful information from .
Screen is near to the pulsar or the Earth -A value much greater than is likely to imply local scattering, either near to the pulsar or the Earth. When scattering is in the pulsar's local environment (e.g. supernova remnants Cordes et al. 2004;Yao et al. 2021, or eclipsing binaries, Main et al. (2018; Lin et al. (2021) ), then → , and → 0. The Earth's velocity is negligible, and the only relevant terms are the distance and velocity between the source and the screen. Assuming the screen velocity is negligible, and the screen is isotropic or aligned with eff , where ≡ − . Interstellar screens 5 pc of the Earth have been observed from scintillation arcs of B1133+16 (McKee et al. 2022, and scintillation of active galactic nuclei (Wang et al. 2021). Under the same assumptions as above, for a local screen, We caution however that screen velocities or misalignment can greatly influence the estimates of or , as their contributions to either are squared.
Screen is halfway to the pulsarwhen the pulsar proper motion is unknown, it can be estimated given an assumption of the screen distance. If e.g. the screen is halfway to the pulsar, then eff = , = 1/2, and While highly uncertain, this can result in a rough estimate of pulsar proper motions. This argument is equivalent to determining proper motions using scintillation velocities, e.g. Cordes & Rickett (1998).

Data
The data and processing for this paper are identical to those used in the census of Posselt et al. (2022); all of the data were taken under the Thousand Pulsar Array project (Johnston et al. 2020) as part of MeerTime (Bailes et al. 2020). In brief the observations were taken with the MeerKAT L-band receiver (centered at ≈ 1283.6 MHz), and are folded with 8 s integrations and ≈ 0.836 MHz channels, corresponding to a total range of 62.5 mHz in and ≈ 0.6 s in . For one pulsar (see Section 4.3) we used the data available at a high time resolution of 38.28 s to extend the range of . As observation duration is one of the limiting factors to detecting scintillation arcs, we limited this study to the single longest observation of every pulsar, typically ∼ 5 − 20 min, although in a few cases much longer.

Creating Dynamic Spectra
To create the raw dynamic spectra, radio frequency interference (RFI) is flagged and masked by the standard deviation in the offpulse region, using methods identical to Main et al. (2020); Mall et al. (2022). After RFI masking, the off-pulse is subtracted in every time and frequency subintegration. The background-subtracted folded spectrum is then multiplied by the frequency averaged pulse profile, and summed over phase and polarization to form ( , ).
Due to the high sensitivity and wide band of MeerKAT, the dominant source of noise in the dynamic spectrum of many pulsars is not receiver noise, but rather the pulse-to-pulse variations and the window function caused by the RFI mask. The measured dynamic spectrum is meas ( , ) = s ( , ) RFI ( , ) int ( , ). The intrinsic variations are often removed by dividing the frequency-averaged intensity from every subintegration, but this has the effect of upweighting noisy integrations (and vice-versa), and diverges when the pulsar nulls. We try to solve this problem using a Wiener Filter, following the same steps as Lin et al. (2021). The one addition is that we multiply our window function by the frequency-averaged flux in each time bin, int ( , ), in an attempt to filter out the intrinsic pulse variations. For computational reasons, the full band is split into 8 subbands, and the time axis is split in longer observations (> 24 minutes, 180 time bins), and the Wiener Filter computed separately in each. This leads to slight artefacts at the boundaries, resulting in low noise in the secondary spectrum. Characteristic examples of filtered dynamic spectra are shown in Figure 1.

Creating Secondary Spectra
Due to the 2 scaling of the arc curvature , a regular 2D FFT over a wide band will result in a smearing of scintillation arcs. We account for this by using the 'NuT' transform applied in Sprenger et al. (2021), where the Fourier Transform in time is performed on a scaled axis = ref , which fixes arcs to constant curvature. In addition, compact features of fixed angular positions on the sky will be seen at constant positions across frequency in the secondary spectrum. We choose ref = 1400 MHz for all plots in this paper. Before taking the Fourier transform, the mean of the dynamic spectrum was subtracted, and a 10% border in both frequency and time were tapered with a Hanning window.

Arc Detections
For every pulsar, diagnostic plots were created which included the dynamic and secondary spectra, and 2D autocorrelation functions (ACFs) across the full band, as well as zoom ins to the two clean bands of 920-1020 MHz, and 1300-1500 MHz. Any sources with observable modulation were deemed sources with either fully or partially resolved scintillation -for each of these pulsars, secondary spectra were produced with a range of limits in both clean bands. These plots were inspected manually; scintillation arcs were identified as sources with significant cross-power in the secondary spectrum; ie. power away from the = 0, = 0 axes. Characteristic examples showing arc detections, as well as sources with partially resolved scintillation (in either time or frequency) but without an arc, are shown in Figure 1.

Arc Curvature Fitting
To fit for , we used the 'Normalized Secondary Spectra', which remaps the axis to ,norm = √︁ ref / ). This transformation maps parabolae to vertical lines of ,norm ∝ , and in this paper we take ref = 0.6 s, the maximum value of allowed by the channel widths. Summing over the axis results in the power as a function of arc curvature ( ), where arcs can be identified as peaks. We fit peaks in ( ) with an inverted parabola, where the starting values and fitting range of are chosen to vary by source, using an interactive fitter. Arc asymmetries may originate from local phase gradients across the scattering screen (Cordes et al. 2006;Rickett et al. 2014), so we fit separately for positive and negative values of ; final values of are the mean from both sides, while we also record the difference Δ . The error of is computed from the half width at half maximum (HWHM) divided by the S/N of the peak, where the N is the standard deviation of (| | > 2 fit ). For sources with screen distances discussed in the text (all of which are cases of large suggesting local scattering), the two screen solutions are estimated from equations 8 and 9. The screen solution near the pulsar uses previously measured proper motions, while the solution near Earth uses ) to predict Earth's velocity vector in the plane of the sky in the direction of the pulsar.

Results
Of the 1270 pulsars observed, we found 107 to have observable scintillation arcs. The measured values of , Δ and derived estimates of ISS are in Tables 1 and 2, and secondary spectra plots are shown in Figure A1.
The detection of arcs is greatly limited by the short durations of observations and the coarse channels; many sources showed a clear scintillation pattern in frequency, but with too short a duration to resolve the pattern in time, and many more sources have scintles far smaller than the channel bandwidth, and are averaged out in each channel. The channel resolution means we are biased towards detecting arcs in mildly scattered, primarily nearby / low DM sources. The short observing durations mean that only high sources will show an arc; ie. sources with large velocities, or screen near to the pulsar or Earth. This is seen in our results − the inferred proper motion velocities are 100 km/s in most sources, in many cases much larger.
In addition to the constraints from channel bandwidth and duration, the detection of arcs can be limited by S/N, or by unmasked RFI, pulsar moding or nulling, which lead to convolved noise in the secondary spectrum. While we tried to account for the pulsar's variable emission in the dynamic spectra using a Wiener Filter, this is highly suboptimal when pulsars are nulling for 50% of the time, and the RFI flagging and excision may not be perfect in all observations. We plot the DM and S/N of our arc detections compared to the full sample in Fig. 2. The sources with arcs are primarily with DM 100 pc cm −3 , and with S/N 100, where higher DM sources . Dynamic spectra are normalized such that their average intensity is 1, and secondary spectra are plotted in log-scaling, spanning the full range in and . J1136+1551 and J0922+0638 both show arcs; J1136+1551 is weakly scintillating with the intensity slightly modulated around 1, while J0922+0638 strongly scintillates with intensity extending to many times the mean. J2215+1538 shows strong scintillation resolved in frequency, but unresolved in time due to low eff / √ eff and short duration. J0835−4510 shows strong scintillation resolved in time but far unresolved in frequency; due to high S/N, scintillation is still apparent but with greatly reduced modulation. J1818−1422 is scattered to ≈ 12 ms at L-band (Oswald et al. 2021), corresponding to ≈ 10 Hz, entirely unresolved by our channels.
are typically more scattered, with scintillation unresolved due to our coarse channels.

Multiple Screens, and Evolution with Frequency
We detect 5 sources with clear evidence of multiple arcs, shown in Figure 3. These sources all tend to be very low DM, nearby sources, and exhibit weak scintillation (ie. with / 1, see J1136+1551 in Fig 1). The fact that the nearest pulsars show multiple screens suggests that there are many screens along any given line of sight. However, strongly scintillating sources tend to show a single, dominant arc (e.g. J0922+0638 in Fig. 1). This can arise either from physical reasons or from selection -it is possible that there tends to be one dominant screen along the line-of-sight, perhaps due to having the highest ΔDM, and being near to the halfway point which maximises the path length difference (and thus ) for a given deflection angle. In addition, 2-screen effects could quench the scintillation from additional screens once the source appears extended, while in weakly scattered sources the pulsar line-of-sight image dominates and 2-screen effects are subdominant. A possible selection effect is that sources of strong scattering result in inverted arclets, which may drown out additional screens in the secondary spectra. One additional selection effect to consider is that all of the multi-screen sources have large proper motions of 140 − 660 km/s, where the slowest of these sources (J1057−5226) had a comparatively long 87.8 min observation. This suggests that the detection of multiple arcs is greatly influenced by the resolution in , in addition to the effects described above. We discuss future avenues to break these ambiguities in Section 5. The outermost arcs in three of our multi-arc sources imply scattering screens very near the source, with ≈ 35 pc, ≈ 27 pc, ≈ 65 pc in J0536−7543, J1239+2453, and J2048−1616 respectively, or within a parsec of the Earth. If the scattering is near Earth, there would be a strong annual modulation of the arcs, similar to the outermost arc in B1133+16 (McKee et al. 2022, and the corresponding screens would have a large angular extent. These multiscreen sources imply a screen density along the line of sight of 200 − 300 pc/screen, while the 6 screens in B1133+16 shown by McKee et al. (2022) imply 60 pc/screen.

Secondary Spectrum Evolution with Frequency
In Figure 3, in addition to multi-screen sources, we included sources which show interesting behaviour across frequency. This includes sources with discrete features at the same location across frequency, consistent with scattering from compact regions of fixed angular position. PSR J1614+0737 shows an evolution from weak to strong scattering across the band, showing evidence of inverted arclets at 1 GHz indicative of a highly anisotropic screen. PSR J1740−1311 looks like a diffuse arc at 1 GHz, but appears to show structured power at 1.4 GHz. A quantitative comparison of the feature positions across frequency (e.g. Brisken et al. 2010;Rickett et al. 2021), and with time (e.g. Hill et al. 2005;Sprenger et al. 2022) will concretely determine if scattering occurs from compact regions, and may distinguish between scattering from over-or under-dense regions.

PSR J1731-4744
PSR J1731−4744 (B1727−47) was recently determined to be a high proper-motion pulsar with = 151 ± 19 mas/year, and tracing back the pulsar motion suggests an association with supernova remnant RCW 114 (Shternin et al. 2019). In our sample, J1731−4744 stood out as having the highest value of , or equivalently showing an extremely fast ∼several second scintillation timescale with resolvable scintles in frequency, shown in Fig. 4. For this source, we re-reduced the filterbank data with single pulse time integrations (≈ 0.83 s), to fully resolve the arc in . The measured value of = 3830 ± 440 km s −1 kpc −0.5 indicates either a screen near to the source, with estimates of ≈ 0.02, and ≈ 12 pc, or a screen at ≈ 0.03 pc from the Earth. The SNR shell of RCW 114 has a radius of ≈ 2 • , ≈ 20 pc at p ≈ 600 pc, comparable to pl . Our result strongly suggests the association of J1731−4744 with RCW 114, in agreement with Shternin et al. (2019).

Comparison with Literature
Several of the sources in our sample have had previously observed scintillation arcs. We highlight and compare a few of these results.
J0837+0610/B0834+06 is one of the best studied pulsars for scintillation, showing dramatic inverted arclets, being used for phase retrieval (Walker et al. 2008), and imaged through VLBI to be a highly anisotropic screen (Brisken et al. 2010) at a distance of ≈ 420 pc. In our data, we detect only a forward arc without inverted arclets. Our estimate of ≈ 412 pc is quite close to the more robust value from Brisken et al. (2010), suggesting we are still seeing the same scattering screen, where this agreement is likely aided by the fact that the screen is (seemingly coincidentally) closely aligned with the pulsar proper motion. The transition from strong to weak scattering is consistent with the findings of Smirnova et al. (2020), who observed the same change in the source's scintillation between 2012 and 2015.
J1057-5226/B1055-52 was studied with Parkes data at 1400 MHz by Kerr et al. (2018), finding a prolonged period of increased scattering (ie. decreased scintillation bandwidth) from MJD 54725-55725, interpreted as an extreme scattering event (ESE). During the ESE, the secondary spectra were dominated by a single arc with inverted arclets at curvature of outer ≈ 0.02s 3 . In observations at MJD 57506-57509, long past the ESE, a second, interior arc is visible at inner ≈ 0.06s 3 . These arcs correspond to outer ≈ 340 km s −1 kpc −0.5 , inner ≈ 200 km s −1 kpc −0.5 , in close agreement with the two arcs seen in our data at outer = 372 ± 79 km s −1 kpc −0.5 , inner = 166 ± 14 km s −1 kpc −0.5 on MJD 58655. This suggests that the nature of scattering is similar to the state of 'quiescent' scattering observed in the later observations of Kerr et al. (2018); the same arcs being visible over many years, combined with the pulsar's transverse velocity of ≈ 140 km/s implies the transverse extent of the two screens is 100 AU, 300 AU respectively. J1136+1551/B1133+16 was studied by McKee et al. (2022), who measure the variations of 6 arcs across 40 years. In our single observation, we clearly detect three of these arcs, which correspond to their arcs B, C, and E. In longer observations with finer frequency resolution, we could likely further disentangle the multiple screens.
J1239+2453/B1237+25 shows clear evidence of 3 screens in our data, across the full band. In previous LOFAR observations at  MHz, the secondary spectra show only a diffuse blend of power . Additionally, only the inner arc was faintly seen by Fadeev et al. (2018), who observed the source with GBT at 316 − 332 MHz. Wu et al. (2022). Their values of ≈ 4.4 3 at observing frequency of 150 MHz correspond to ≈ 216 km s −1 kpc −0.5 , in agreement with our measured value of = 243 ± 25 km s −1 kpc −0.5 , suggesting we are observing the same screen. Rickett et al. (2021) at 340 and 800 MHz with the GBT, showing clear inverted arclets and being fit well with a 1D screen brightness distribution. In our data, we miss this structure, seeing only the outline of an arc. This shows the limitations of our short observing time (9 min, compared to 60 min observations), and that rich structures can be revealed using longer integrations. It is also possible that our observation is at a low point of eff, | | from annual motion.

Ramifications
In summary, we have presented the discovery of 107 sources with detectable scintillation arcs with the Meerkat TPA. In this section we discuss ramifications and future avenues.
Our findings are consistent with recent suggestions that scintillation arcs are prevalent , their detection primarily limited by frequency resolution and observation duration, provided the dynamic spectrum is high S/N and largely devoid of RFI. The nature of this survey manifests in strong biases, where arcs are easiest to detect in weakly scattered sources with high effective velocity -either through a high proper motion, or local screens, as in the case of J1731−4744. The detection of multi-arc sources include even more selection biases, with multi-arc sources primarily being nearby, weakly scintillating sources with very high transverse velocity. There are likely many screens on any given sightline which cannot be seen in secondary spectra due to insufficient time duration, 2-screen interactions, or inverted arcs.
This survey leaves many avenues for further research. In the TPA  alone, there are hundreds more sources with resolvable scintillation, but without the detection of arcs. Analysis of the full sample can provide measures of scattering time on many more sightlines, which will help improve galactic electron models (e.g. Cordes & Lazio 2002;Yao et al. 2017). The more distant, scattered sources will be best studied at higher frequencies where arcs with be sharper, more easily resolvable. Through wide-band/multi-frequency observations, or through phase-retrieval techniques such as cyclic spectroscopy (Demorest 2011;Walker et al. 2013), one could measure the correspondence between time-domain scattering and scintillation arcs. This may help elucidate the morphology of scattering tails, and investigate if seemingly isotropic scattering could be reproduced through the ensemble of many anisotropic screens (e.g. Oswald et al. 2021).
With only a single arc curvature measurement, only crude esti-  mates of screen distances can be made. Additional observations to obtain annual curves Main et al. 2020;Mc-Kee et al. 2022;Mall et al. 2022), or VLBI (Brisken et al. 2010 can fully solve for screen distances and geometry. In multi-screen sources, VLBI combined with dynamic cross-spectra can be used to simultaneously solve for multiple screens (Simard et al. 2019), and characteristic 2-screen effects may also provide additional constraints on screen geometry, as seen in B1508+55 (Sprenger et al. 2022). With a growing number of scattering screen measurements, it may be possible to map the distribution of screens in the Milky Way, cross-matching with foreground stars, HII regions, supernova remnants, and with the local bubble and other local plasma. With high S/N, one could probe magnetism in scattering screens, and compare scintillation across pulse phase to resolve pulsar emission. The results of this survey will then serve as a valuable stepping stone to studying the distribution of free electrons in the Milky Way, probing the dynamics and emission of pulsars, and elucidating the astrophysical origin of interstellar scintillation. This paper has been typeset from a T E X/L A T E X file prepared by the author. Table 1. Arc measurements in sources with known proper motions. and Δ are computed as described in Section 3.4. Values of d psr with errorbars come from parallax measurements cited below, otherwise they are DM distances from NE2001 (Cordes & Lazio 2002). In Δ , − or + denote cases where the arc could only be measured on the negative or positive side of D,norm . Distance references: (1) Deller et al. (2009)