The Thousand-Pulsar-Array programme on MeerKAT -- V. Scattering analysis of single-component pulsars

We have measured the scattering timescale, $\tau$, and the scattering spectral index, $\alpha$, for 84 single-component pulsars. Observations were carried out with the MeerKAT telescope as part of the Thousand-Pulsar-Array programme in the MeerTime project at frequencies between 0.895 and 1.670 GHz. Our results give a distribution of values for $\alpha$ (defined in terms of $\tau$ and frequency $\nu$ as $\tau\propto\nu^{-\alpha}$) for which, upon fitting a Gaussian, we obtain a mean and standard deviation of $\langle\alpha\rangle = 4.0 \pm 0.6$. This is due to our identification of possible causes of inaccurate measurement of $\tau$, which, if not filtered out of modelling results, tend to lead to underestimation of $\alpha$. The pulsars in our sample have large dispersion measures and are therefore likely to be distant. We find that a model using an isotropic scatter broadening function is consistent with the data, likely due to the averaging effect of multiple scattering screens along the line of sight. Our sample of scattering parameters provides a strong data set upon which we can build to test more complex and time-dependent scattering phenomena, such as extreme scattering events.


INTRODUCTION
A commonly observed signature of the effect of the interstellar medium (ISM) upon the radio emission from pulsars is the scattering of the radio flux. Regions of cold (< 10, 000K), dense plasma along the line of sight, which can be approximated as one or multiple thin screens, cause some of the radio flux arriving at the observer to be delayed with respect to the direct line of sight. Observationally, this results in a scattered pulse profile with a characteristic exponential scattering tail (Cronyn 1970), which can be explained assuming isotropic scattering screens, for which the scattering angles have no preferred direction. Measurements of the scattering properties of pulsars reveal information about the structure of the ISM and allow ★ E-mail: lucy.oswald@physics.ox.ac.uk (LSO) us to disentangle these effects from the intrinsic properties of the pulse shape.
In the context of isotropic scattering with an exponential transfer function the key parameter is the scattering timescale, . It is observed that evolves with frequency according to the power law relationship where is the scattering spectral index. The simplest thin screen scattering model predicts = 4 (e.g. Cronyn 1970;Lang 1971), whilst a medium exhibiting Kolmogorov turbulence would have a value of 4.4 (e.g. Lee & Jokipii 1976;Rickett 1977). The simplified spectrum of turbulence is described by a power law for wavenumbers , for values of −1 that lie well within inner and outer fluctuation scales ( −1 and −1 respectively). This is written as ( ) = 2 − (Rickett 1977), where 2 is the proportionality constant, dependent on the electron density. The fluctuation spectral index cannot exceed 4 (Romani et al. 1986). It is related to by = 2 /( − 2), which leads to = 4 being a theoretical upper limit. Observationally however, lower values of , i.e. flatter spectra, are commonly seen (see for example Krishnakumar et al. 2019). A power law break at low frequencies with the scattering timescale approaching frequency independence, accompanied by a loss of flux, may be attributed to a truncated scattering screen (Cordes & Lazio 2001). Where −1 drops below the inner scale of the plasma turbulence, −1 , the power law description of the turbulence spectrum is no longer a valid simplification. Lewandowski et al. (2013Lewandowski et al. ( , 2015a write that, for a single pulsar, this may lead to a flattening of at lower frequencies. The effect of this inner scale, which corresponds to the shortest scale length in the scattering material, should also become apparent in the scattered pulse shape, particularly at large delays in the profile tail (Rickett et al. 2009). Furthermore, the exact shape of the scattered profile will depend upon the degree of anisotropy in the scattering material. Evidence of such anisotropy has been seen in the parabolic arcs observed in the secondary spectra of scattered pulsar observations (e.g. Stinebring et al. 2001;Walker et al. 2004) and high resolution mapping of the locations of scintils (Brisken et al. 2010;Pen et al. 2014). Geyer & Karastergiou (2016) showed that fitting an isotropic model to simulated anisotropically scattered data could lead to inferring smaller values of .
In the literature, values of are usually calculated by performing a power law fit to against frequency, where is measured through one of two ways. In the cases where the scatter broadening is evident in the pulse profile, we measure in the time domain. This is done either through forward modelling (e.g. Geyer et al. 2017) or through deconvolution analysis, such as the algorithm (Högbom 1974;Bhat et al. 2003). Scattering surveys that have employed forward modelling, and are thus directly comparable to the work presented here, are as follows: Cordes et al. (1985); Löhmer et al. (2001Löhmer et al. ( , 2004; Kuzmin & Losovsky (2007); Lewandowski et al. (2011Lewandowski et al. ( , 2013Lewandowski et al. ( , 2015a; Krishnakumar et al. (2017); Geyer et al. (2017) and Krishnakumar et al. (2019). The algorithm has been employed by, for example, Bhat et al. (2004) and Kirsten et al. (2019). A complementary technique is to infer from the scintillation bandwidth (Cordes et al. 1985). This is generally done using the equation 2π = 1 , where 1 is a factor with a value that depends on the electron density wavenumber spectrum and distribution of scattering material. When transforming observational values it is usually set to 1 = 1, the solution for a thin screen, and other values are presented in Cordes & Rickett (1998). This technique is relevant for observations where is too small to be measured directly, generally at higher frequencies. The inverse is also true: where is large, will be too small to be measurable.
The scattering strength, which determines the choice of methodology used to measure and , is expected to be correlated with distance for a given pulsar, added to which is a further stochastic element due to there generally being a small number of scattering screens along the line of sight. Modelling the scattering material as a power-law electron density spectrum with a Kolmogorov distribution of irregularities, the relationship between scattering timescale and DM is expected to be ∝ 2 −4.4 DM 2.2 (Romani et al. 1986;Cordes & Rickett 1998;Krishnakumar et al. 2015).
The new MeerKAT telescope (Bailes et al. 2020) provides high sensitivity in the observing band 856-1712 MHz. The high quality pulse profiles and the broad bandwidth are suitable for a survey of scattered pulsars where measurements of and can be made using a single instrumental setup. This largely eliminates systematic errors arising from measurements of made using multiple receivers and backends at different frequencies. Our survey of scattered pulsars is part of the work of the Thousand-Pulsar-Array (TPA) (Johnston et al. 2020). This is an observing project that is being carried out as part of MeerTime, a large scale project to observe known pulsars using the MeerKAT telescope (Bailes et al. 2020).
The TPA programme will obtain an overview of the properties of a large sample of the observed pulsar population. The structure of this paper is as follows: section 2 describes the observations and data processing. Section 3 describes our new time domain methodology to determine , which uses a Markov Chain Monte Carlo (MCMC) fit. The results of the scattering analysis, both a table of measured scattering parameters and a description of the distributions of these parameters across the population of observed pulsars, are presented in section 4. We also describe how our observational constraints affect the sample of pulsars for which we are able to accurately measure scattering parameters, and we investigate how the level of covariance between the parameters affects our modelling outcomes, before discussing our results in the context of previous work in section 5. Conclusions are summarized in section 6.

OBSERVATIONS AND DATA
We selected 205 pulsars from the TPA programme with simple profiles that resemble a single Gaussian component convolved with an exponential function. We chose these by eye out of 1164 pulsars observed up until 1st June 2020. The initial sample was deliberately selected to be as broad as possible within these constraints, so that the sample was not biased by visual selection. This allowed us to investigate the limitations of our modelling choices, but required us to discard at later stages those pulsars whose profiles were not well modelled with single Gaussians modified by scattering tails. The observing and data processing was carried out as described in Johnston et al. (2020) and Serylak et al. (2020). The pulsars were observed in fold mode over multiple epochs using MeerKAT at L-band. These observations were de-dispersed, cleaned of Radio Frequency Interference (RFI) using C G (Lazarus et al. 2016(Lazarus et al. , see ascl.net/2003) and time-integrated to a single Stokes I profile per observation with a resolution of 1024 bins per pulse period . Observations from the early commissioning phase of the MeerKAT telescope had only the band between 895 and 1670 MHz available out of the ultimate full band of 856-1712 MHz. In order to make all our observations directly comparable with each other, we therefore kept the band between 895 and 1670 MHz, divided into 8 sub-bands, in the final data products. As originally shown by Geyer & Karastergiou (2016) 1 , assuming a power-law relationship between scattering timescale and frequency, the appropriate frequency associated with a given measurement of is given in terms of the centre frequency, , and bandwidth, , of the profile as = 10 [log 10 ( + /2)+log 10 ( − /2)]/2 . (2) Where a pulsar has multiple good quality observations taken at different epochs, we aligned and added these observations to produce a single profile with increased signal-to-noise ratio (S/N). We did 1 Equations 5 and 6 in Geyer & Karastergiou (2016) have typos, with and interchanged. As of 4th February 2021, this has been corrected in the arXiv version of the paper, found at https://arxiv.org/pdf/1607.04994.pdf this through the following method. Within the TPA-project, pulse templates were obtained. Ephemerides were updated when the data showed apparent deviations from the known ones. Employing these ephemerides, we used 2 Edwards et al. 2006) and (van Straten et al. 2012) to obtain phase shifts between individual observations of a pulsar, which we used to align them in phase before adding.

Method
We modelled the scattered pulse profile as a single Gaussian component representing the intrinsic emission, convolved with an exponential scattering function. We assumed that the scattering material took the form of a thin screen with isotropic scattering properties. The temporal broadening function was therefore described in terms of the scattering timescale as where the unit step function ( ) constrains the equation to time > 0 (Cordes & Lazio 2001). We fit for 5 components: the amplitude ( ), mean ( ) and standard deviation ( ) of the Gaussian; the scattering timescale ( ); and any DC offset of the profile baseline. Early modelling also investigated the applicability of an extreme anisotropic model, as given by equation 3 in Geyer et al. (2017). However, since this was found to be generally less successful at fitting the profiles than the isotropic model, we did not take this further.
We fit each of the 8 sub-bands per pulsar independently. We used the train + DC method as described in Geyer & Karastergiou (2016), updating the fitting process to sample the log-likelihood function over parameter space using the MCMC algorithm (Foreman-Mackey et al. 2013). We imposed flat priors, with the only constraint that the parameters had to be physical, in order to allow the fit to converge freely. This meant we constrained and to be positive, and 0 < < 1024 and 0 < < 1024, where 1024 is the number of bins across the pulse period. We ran the MCMC for either 20,000 steps per frequency, or until the chain autocorrelation time estimate was changing by less than 1% and the chain was longer than 100 times the autocorrelation time, whichever was shorter. The large number of model fits that must be run necessitates the upper limit on the number of steps, in order to complete the modelling in a reasonable time frame. We implemented a burn-in of 50%, discarding the first half of the total number of chain steps, i.e. 10,000 steps discarded if the MCMC has run for 20,000 steps. This large burn-in ensures that the results are not biased by the initial conditions. We kept only those sub-bands for which the chains converged onto a single set of parameters for the model. Cases for which the chains failed to converge, or converged on two or more sets of parameters simultaneously, we concluded to be unsuccessful and excluded from the subsequent analysis. We calculated the best fit parameters to be the 50th percentiles of the MCMC samples, with the errors taken as the average of the differences between the 16th and 50th, and 50th and 84th percentiles respectively. We made use of GNU Parallel (Tange 2011) to run our modelling on multiple pulsars simultaneously. Given the values of obtained through our MCMC fit, we calculated . Initial calculations of were done with least squares fitting of a power law to vs frequency. However, since it was found that different least squares algorithms tended to produce different answers, we decided instead to use a further MCMC to map out the probability distribution and gain a better understanding of convergence on the solution. We kept only those 142 pulsars which had at least 4 channels for which the fit was successful, and then also converged on a successful fit for . The value of and its uncertainty we took from the same percentiles of the MCMC samples as described above.

Applicability of methodology
When drawing conclusions on the properties of a population of scattered pulsars, we require reliable measurements of the scattering timescale as a function of frequency. The first test is to check the diagnostic output of our methodology showing as a function of frequency and reject all 21 sources which showed no significant evolution of across the band. Even though some pulsars matched the selection criterion of a Gaussian profile convolved with an exponential, they showed no frequency evolution consistent with scattering. There is an expectation that time domain fits will break down when is too small to measure. If the underlying profile is not Gaussian in shape, this too will affect the fitting, particularly when is small. In Appendix A we show that attempting measurements of when / ≥ 1 is likely to lead to inaccuracy, where is the standard deviation of the Gaussian describing the intrinsic pulse profile. The second test therefore is to inspect our results for and and reject 31 sources with unreliable measurements. We describe how we identify this unreliability with an example in section 4.2. A further two pulsars (J1320-3512 and J1738-3211) we discarded due to the model being visibly inaccurate at the leading edge of the profile. The rising edges of the model and profile respectively did not overlap for these profiles. PSR J1738-3211 has an asymmetric profile with the trailing edge steeper than the rising edge, implying that it is not scattered and should be discarded regardless. PSR J1320-3512 has a leading edge that is shallower than the model leading edge. This suggests either that the intrinsic profile shape is not well modelled by a Gaussian, or that the scattering function for this pulsar might be better described by a thick screen model (Williamson 1972;Kirsten et al. 2019), which has a more gradual leading edge shape than the thin screen model applied in this work. We rejected PSR J1655-3844 after closer inspection at high frequencies revealed it not to fit our single-component criterion. Of the pulsars remaining, there are 15 for which we have rejected a section of the band, for reasons detailed in section 4.3 and Appendix B. The selection criterion of requiring 4 successful fits for then removes a further 3 pulsars from the sample. Having completed the processing of the modelling output as described above, from our original sample of 205 pulsars we have a final subset of 84 for which we have scattering fits in which we can be confident, with at least 4 measured values of .

RESULTS
Before presenting the results for our full set of pulsars, we show an example of a successful fit, and cases where the results of the fit are unreliable.

Modelling results: example PSR J1818-1422
Fig. 1 shows the eight observed pulse profiles of PSR J1818-1422 with the best fit model at each frequency overlaid. We show how the parameters for the model at 0.94 GHz were obtained by representing the MCMC chain as a corner plot in Fig. 2. Fig. 3 shows the corner plot for the power law fit of = − , which gives = 3.787 ± 0.008. The figure shows that and are covariant, but that is nevertheless tightly constrained. The top subplot of Fig. 4 shows how both and evolve with frequency for this pulsar, with the power law fit indicated with a straight line. In addition, the modelling indicates intrinsic profile evolution, as increases with decreasing frequency below ∼1 GHz, evolution comparable to that described for many other pulsar observations (see for example Thorsett 1991).

Covariant model parameters: an example
We find that the results for several pulsars, particularly those that are less strongly scattered, exhibit an inter-relation between the measured values of and that is not inherent to the pulsar and must result from the fitting process. An example is PSR J1743-3153, for which we show a plot of both and against frequency in the middle subplot of Fig. 4. The decrease in with frequency seen for the first three measurements is mirrored by an increase in . The converse is then true for the next measurement: a drop in is counteracted by an increase in . Similar mirroring behaviour can also be seen in the measurements for the highest frequencies. This suggests that the modelling process has a covariance in and : multiple different pairings of these parameters all generate similarly shaped profiles and are therefore indistinguishable. Each sub-band is modelled individually, and the inaccuracy of model parameters converged upon for a given sub-band is then only apparent when compared to the results for the other sub-bands.
It seems likely that PSR J1743-3153 is scattered, but the mirroring behaviour of and means we cannot trust the accuracy of the parameter measurements. We therefore exclude this pulsar, and 32 more sources showing similar behaviour, as identified by visual inspection of similar figures, from further analysis.

Power law deviations and non-Gaussian pulse shapes
There are 15 pulsars which show a visible deviation from a single power law at the higher frequencies, with larger than expected. We hypothesize that this may be caused by the intrinsic profile, and not the scatter broadening, being the dominant factor in the overall pulse shape at high frequencies for these pulsars. As shown in Appendix A, time domain fitting is unlikely to be accurate when / ≥ 1 and the intrinsic profile dominates over the scattering timescale. It is logical that there is a subset of pulsars for which scattering is measurable at lower frequencies, but the pulse shape is dominated by the intrinsic profile at higher frequencies. If so, our assumption of a single Gaussian component is insufficient to replicate the profile shape, resulting in over-estimates of . We address as an example the case of PSR J1653-4249, for which and are plotted against frequency in the bottom subplot of Fig. 4.
It is possible to simulate this pulsar using a profile with an additional weak trailing edge component. At high frequencies the absence of scattering is compensated by the presence of the trailing edge component, such that the measured values of do not continue to decrease with increasing frequency (see Appendix B for details). Although the simulation is very specific, it demonstrates that it is possible to obtain such a result under these conditions and therefore attributing the anomalous behaviour seen in the bottom subplot of Fig. 4 to properties of the ISM requires caution.
In general, for this sample of 15 pulsars, we conclude that the intrinsic profile shape is the dominant factor in determining pulse shape at higher frequencies, rather than the scattering timescale, and that as a result any non-Gaussianity in the intrinsic profile shape may be mis-modelled as contributing to the scattering tail. A possible alternative method for modelling pulse profiles with intrinsic shapes more complex than a Gaussian would be to apply the algorithm. However, this would require its own set of assumptions and would result in modelling results that are not directly comparable to the rest of the results presented in this paper. For these 15 pulsars therefore, we kept only that frequency range where follows a power law, selecting the values through visual inspection, and inferred from those.

Scattering parameters for 84 pulsars
We present the results of the modelling in Table 1. For each pulsar we list the scattering timescale at 1 GHz, along with and the corrected dispersion measure (DM) that we measure through our calculation of the position of the intrinsic pulse profile at each frequency. We also indicate which of the 8 frequencies have been used to calculate these values, indicating which channels were excluded due either to failing to reach convergence on a solution, or to deviating from a power law at high frequencies, as explained in section 4.3. We use our MCMC power law fit = − to compute at = 1 GHz. We do this by calculating the value of at 1 GHz associated with every pair of values of and explored by the MCMC chain, and taking the 50th quantile of these values as and the difference between the 16th and 84th quantiles for its error. We discuss the origin of our corrections to the DMs of these scattered pulsars, along with the relevance of using these corrected DMs, in section 5.2.

The scattering spectral index distribution
Fig. 5 shows histograms of the distributions of obtained for our survey. We show the final filtered sample of 84 pulsars (hatched histogram) and the 33 pulsars which were filtered out due to problems in the modelling (histogram with circle pattern). These have primarily been removed due to covariance between the values of and , as described in sections 3.2, 4.2 and appendix A, but this subset also includes the two pulsars rejected due to the leading edge not being correctly modelled (PSRs J1320-3512 and J1738-3211). The removal of pulsars is done without reference to the value of obtained, so that we are not biased in favour of those pulsars with values of close to what might be expected from theory. The result of the filtering is that pulsars with small values of are largely removed from the sample, and those remaining show a nearsymmetric distribution. We have shown, in Section 4.3 and in the appendices, that poor model fits have a tendency to overestimate , and that poor model fits are more likely for smaller vales of , which are seen more often at higher frequencies. It is therefore expected that poorly modelled pulsars will tend to underestimate , and so it is unsurprising that the majority of the poorly modelled pulsars returned small values of . Fig. 5 also shows a histogram of the values of found in the literature for cases where pulsar profiles observed at ≥ 400 MHz have been fitted with an isotropic time domain model (Cordes et al. 1985;Löhmer et al. 2001Löhmer et al. , 2004Lewandowski et al. 2011Lewandowski et al. , 2013Lewandowski et al. , 2015aKrishnakumar et al. 2017Krishnakumar et al. , 2019. The histogram of literature values has been weighted so that it represents the same total count of pulsars as those shown for this work. We apply the ≥ 400 MHz cut because we observe different behaviour in the scattering surveys at lower frequencies, which we discuss in section 5.3. Comparing our distribution of with that of previously published values, we see that, when we include our poorly modelled pulsars, the distributions are similar. The population of pulsars around 0 < < 3 shown in the published values has been filtered out from our own distribution. Performing a Kolmogorov-Smirnov (KS) test to compare our sample to that from the literature, we obtain the following values. When we use just our final filtered sample in the KS test, we obtain a KS statistic of 0.247 with a corresponding p-value of 0.001. Adding in our poorly modelled pulsars as well gives KS statistic 0.109 and p-value 0.335. This indicates that, whereas our final filtered sample is clearly not drawn from the same distribution as that of the combined literature values, when we include our poorly modelled pulsars in our distribution, we cannot reject the null hypothesis. Fitting a Gaussian to our filtered distribution of we find that it has a mean of 4.0 and a standard deviation of 0.6. The theoretical values of both = 4 and = 4.4 fall within one standard deviation. Fig. 6 shows our measurements of plotted against 0 / , where 0 is measured for our lowest observing frequency of 950 MHz. It is interesting to note that if we split our population into sources where 0 is less or greater than 10% of the period, we find a mean value for each group of 3.7 ± 0.6 and 4.1 ± 0.4 respectively. For = 4, a pulsar with 0 / = 0.1 at the lowest frequency will have / ∼ 0.01 at the highest frequency. We have shown in section 4.3 and in the appendices that our methodology overestimates small values of . This suggests that we measure systematically smaller values of for such cases, which is a limitation of the time domain methodology. Fig. 6 shows that the sources for which 0 / < 0.1 are comprised of pulsars that we have filtered out as being either poorly modelled (orange diagonal crosses), or not scattered (blue Table 1. Values for at 1 GHz, and DM for our filtered sample of scattered pulsars. The values and their uncertainties are calculated as described in the text. The DM is the value that best aligns the modelled intrinsic profile at each frequency. We indicate which of the 8 sub-bands (ordered in increasing frequency) were used to generate these values with ones, and those not included with zeros. The minimum allowed number of sub-bands is 4. Reasons for exclusion of sub-bands (failure of model convergence or deliberate exclusion) are explained in the text.    vertical crosses), or they are not filtered out (black points) and yet still show systematically lower values of , as described above.

PSRJ
We scaled the scattering timescale by the period for the analysis above because all of our observations have the same number of bins across the pulse period (1024 bins). The number of bins containing information about the profile shape will affect the accuracy of the modelling, and a profile that takes up a smaller fraction of the pulse period (smaller / ) will have fewer bins spanning the profile. It is therefore the ratio / < 0.1 that is important in terms of likelihood of the model being able to correctly capture the scattering behaviour of a given pulsar. However, the time resolution of the bins is also relevant to consider. The larger a pulsar's period, the lower the time resolution of a single bin, meaning that less information is captured in each bin. We find that the distribution of periods for the / < 0.1 cases is shifted slightly towards longer periods than the distribution for the / > 0.1 cases, which reflects the time resolution limitation.

Correcting for dispersion
As part of our scatter modelling we identify the position of the underlying Gaussian that describes the intrinsic profile in our model. By fitting for a correction to the DM on these values, we obtain the DM that best aligns the intrinsic pulse profile independent of the effects of scattering. Our new values for the DMs of the pulsars in our sample are presented in Table 1. An alternate approach would be to simultaneously model the data in all sub-bands, which would enable fitting for DM and the pulse scattering index. This approach was taken in modelling the pulse profiles of fast radio bursts (FRBs) in Qiu et al. (2020).
We show a histogram of the measured corrections to the DM in Fig. 7. The majority of these are small and negative, as expected (e.g. Geyer et al. 2017). The tail of large negative ΔDM values is partially attributable to the small subset of pulsars that are more strongly scattered and therefore require larger corrections to the DM to align intrinsic profiles. However, it also reflects the distribution of periods in our pulsar sample. This can be understood as follows. We can assume very roughly that pulsars have the same duty cycle. For profiles at two discrete frequencies that are misaligned by an amount on the order of the pulse width, the shift required to align the profiles will be the same fraction of the pulse period. For pulsars with very different periods, this shift will comprise different absolute lengths of time; correspondingly the DM correction to perform this shift will be larger for the pulsar with the larger period.
Since the choice of DM affects the profile shape within each sub-band, a question arises as to whether correcting the DM would change the measurements of . We tested this on PSR J1630-4733, chosen because its ΔDM of 11.3 cm −3 pc results in the largest relative shift of profiles at different frequencies in our sample set. PSR J1850-0006 has a larger absolute ΔDM of −21.0 cm −3 pc, however, since its period is also larger, the corresponding time shift between profiles at different frequencies is a smaller fraction of the period and so the relative misalignment of profiles is smaller. We identified the magnitude of the DM correction required through scatter modelling, then re-processed the data with the new DM and then repeated the scatter modelling. We compare the results for obtained before and after the reprocessing in Fig. 8. The difference is largest at low frequencies, as expected. It is also uniformly negative, meaning that the scattering timescales in the uncorrected DM case are larger. However, even for this pulsar, where the change in DM is most extreme (shifting the top of the observing band with respect to the bottom of the band by 6% of the pulse period), the measured values of are still equivalent to within 1 in all but the third channel (counting from lowest to highest frequency), and equivalent to within 3 for all channels.

Scattering parameters in context
Cross-comparing our sample of pulsars with the literature, we find 9 pulsars in our sample that have previous measurements for . We list these measurements in Table 2. Some of these have been measured multiple times, giving a total of 13 measurements of to compare to our results, found in Löhmer et al. (2001), Lewandowski et al. (2013) and Lewandowski et al. (2015b). We note that 6 of the 7 estimates of given in Lewandowski et al. (2013) are larger than our corresponding estimates, whilst 4 of the 5 estimates given in Löhmer et al. (2001) are smaller than ours. This suggests that differing modelling approaches may tend to give systematic differences in the Figure 7. Histogram of the values of ΔDM calculated for our filtered scattered pulsar sample, where ΔDM is the difference between the original DM used to dedisperse the pulsar, and the best fit DM generated by our scattering model. Figure 8. Plot of Δ vs. frequency for PSR J1630-4733, where Δ is the difference between the scattering timescales calculated for the data processed with the initial, visibly incorrect, DM, and when reprocessed with the best fit DM produced by the scattering model. resultant parameters estimated. Of the 9 pulsars, 5 have at least one measurement consistent with our own, to within the uncertainties. The 4 that are entirely inconsistent are all found in Lewandowski et al. (2013) and all have larger values of quoted there than those presented here. Fig. 9 shows calculated at 1 GHz plotted against DM. The values for the final sample we plot as black circles, and for compar- ison we also plot those values obtained from pulsars subsequently determined to be either not scattered or poorly modelled. The horizontal lines mark the smallest and largest time resolutions for a single bin of a pulse profile. Since all observations have 1024 bins across the profile, these values correspond to 1/1024 th of the smallest and largest periods in the sample. We take the equation relating and DM fit by Krishnakumar et al. (2015) at 327 MHz, scale it to 1 GHz using our best fit of 4.0 and plot it as a black line. We also plot the scaled fits for ± about the mean, = 3.4 and = 4.6, as dashed lines. The functional form of this equation is s = 3.6 × 10 −6 DM 2.2 (1.0 + 0.00194DM 2 )(1000/327) −4.0 , (4) where has units of ms and DM has units cm −3 pc. We see from Fig. 9 that the pulsars in our final sample have high DMs and that they are clustered close to the best fit line of Krishnakumar et al. (2015). There is some spread in the values, meaning that our results are also compatible with previous fits in the literature made by Ramachandran et al. (1997);Löhmer et al. (2004); Bhat et al. (2004) and Lewandowski et al. (2015a). Since our sample is limited to a narrow DM range, it is not possible to perform an equivalent model fit for the -DM relationship using our own results. We note that our results are not symmetrically distributed about the Krishnakumar et al. (2015) model. This may be related to a difference in modelling approach: whereas we model the intrinsic pulse profile as a Gaussian without constraining its parameters, Krishnakumar et al. (2015) use a high frequency unscattered profile as the template for the intrinsic pulse profile.
The choice of how to model the intrinsic pulse profile necessarily affects the results for the measured scattering timescales and spectral indices. For example, the slower rise time of the thick screen model may absorb the intrinsic pulse shape into its measurement, or conversely the choice of a Gaussian intrinsic pulse may obscure evidence for a thick screen scattering function. Another example is the comparison of our scattering timescale result for PSR J1316-6232 with the estimate published by Crawford et al. (2001) of ∼ 150 ms at 1.35 GHz. Scaling our result to this frequency gives a scattering timescale almost twice as large. This difference follows directly from the difference in modelling choices. Whereas Crawford et al. (2001) provide an estimate of the scattering timescale based only on the exponential decay of the profile intensity, our model involves a convolution of that exponential decay with a Gaussian representing the intrinsic profile shape. For this reason, we caution against direct comparison of individual values measured using different methods. A systematic shift of values caused by the choice of method is likely to have less of an impact on , provided all the values of used to calculate it were estimated using the same method. Nevertheless, this example highlights the importance of applying the same modelling approach to a large-scale sample of pulsars, such as this one, in order to be able to compare the scattering results for different pulsars.
Of the values of associated with poor modelling, there are several at smaller DMs for which is larger than would be expected from the Krishnakumar et al. (2015) model. This is easily explained based on the results of our investigations into the causes of poor model fits. It is expected that pulsars with smaller DMs will, in general, be less strongly scattered. Our work shows that small scattering timescales are likely to be over-estimated, particularly when the timescale size is similar to the width of the intrinsic pulse profile. These poorly modelled pulsars also have estimates close to the temporal resolutions of the pulsars. Indeed, some of these pulsars were concluded not to be scattered at the observing band, which takes the concept of small values of being overestimated to the logical extreme. We caution therefore that there may be a wider tendency to overestimate at the lower limits of both DM and attainable at a given frequency.
This may go some way to explaining what we observe in Fig.  10, which is a plot of against DM for our values (black circles) and those given in the literature as described in section 5.1 and shown in Fig. 5. Our best fit of = 4.0, shown as a horizontal black line, is consistent with the results for the literature. On calculating mean values for 5 bins across the DM range we saw no significant evolution of with DM, save for a slight increase in at the highest DM bin, encompassing the three points in Fig. 10 at around 1000 cm −3 pc). More data is required to determine whether this increase is statistically significant. We also note that weighting the averages by the uncertainties in tends to favour higher values in comparison to the unweighted averages. Inspecting the literature values of , we note a greater spread in the literature values at lower DMs. As we have described, the over-estimation of smaller values of tends to lead to under-estimation of . The reduced strength of scattering at lower DMs should cause this to have a greater effect at low DMs, which is what we see in Fig. 10.
There is a further consideration to which attention must be brought: whereas our modelling explains the distributions of and obtained for pulsar observations at frequencies ≥ 400 MHz, other behaviour is seen in the low frequency studies performed by Kuzmin & Losovsky (2007) and Geyer et al. (2017). Kuzmin & Losovsky (2007) identified a -DM relation of = 60(DM/100) 2.2 ms at 100 MHz, and Geyer et al. (2017), whose results for also corresponded well to this equation, measured a distribution of values of that is systematically shifted to lower values in comparison to ours. An explanation for this may be that low frequency pulsar observations are probing a different scattering environment. One aspect of this is that lower frequency pulsar emission probes a wider region of space than that at higher frequencies since it is scattered more strongly, as shown in fig. 1 of Cordes et al. (2016). Applying the same isotropic model may therefore result in different measurements of at different frequencies, something not measurable for the frequency range of the data presented here but meriting further exploration. A further consideration is the effect of distance. Scattering analyses performed at low frequencies focus on nearby  Krishnakumar et al. (2015), scaled to 1 GHz using our best fit of 4.0. Dashed lines: the same equation, now scaled using = 4.0 ± 0.6. Pale grey crosses: pulsars rejected from the sample due either to inaccurate modelling or the pulsars not being scattered. Horizontal dot-dash lines: these mark the minimum and maximum time resolutions of a single bin for the pulsar sample, where all observations had a resolution of 1024 bins across the pulse period.
pulsars, due to the extreme scatter broadening that occurs in high DM pulsars. It is interesting to consider the effect of anisotropy in this context. First, having compared anisotropic model fits with the isotropic fits presented in this paper and found them to be generally unsuccessful, we see no evidence of anistropic scattering in the sample presented in this paper. Secondly, for nearby pulsars, it is likely that the pulsar emission only passes through a single or a very small number of screens, so that the scattering behaviour observed may bear hallmarks of anisotropy. There is evidence of anistropic scattering in low frequency surveys, e.g. Geyer et al. (2017), and using an isotropic model to identify the scattering parameters where such anisotropy is evident then results in under-estimation of , as described by Geyer & Karastergiou (2016). This may be another factor contributing to the smaller values of at lower DMs seen in Fig. 10: pulsars with smaller DMs are likely to be nearby and hence an anisotropic model may be more appropriate. By contrast, the emission from our high DM pulsars is likely to have passed through multiple scattering screens, all of which may have different levels of anisotropy oriented along different axes. The net scattering effect observed will therefore approximate isotropy, as observed in our data.

CONCLUSIONS
As part of the TPA project we have identified 205 single-component profiles that appear to be scattered, out of a total of 1164 pulsars observed with the MeerKAT telescope. Of these we have obtained good scattering model parameters, and , for 84 pulsars, the largest sample of scattering parameters observed with a single telescope over a continuous bandwidth. Through this we have also obtained estimates of these pulsars' dispersion measures. The preferred choice of DM depends on the purpose for which dedispersion is being applied, and these DMs may be used to obtain accurate alignment of intrinsic pulse profiles for these pulsars, independent of scattering.
Our investigations have identified and highlighted a key cause of inaccurate scatter modelling as covariance of model parameters for the intrinsic pulse profile and the scattering timescale. We have investigated the regimes under which such covariance is likely to cause problems and identified the effects of poor modelling as a tendency to over-estimate and under-estimate . Putting our results into the context of previous work, we propose that high-frequency, high-DM scattered pulsar observations are in general well modelled by the simple assumption of isotropy. Our measured average scattering spectral index, = 4.0 ± 0.6, agrees with both of the simplest theoretical models, 4 (thin screen) or 4.4 (Kolmogorov continuum), to within 1-uncertainty. More complex cases, in particular pulsars where anisotropic scattering is measurable, are observable at lower frequencies, whereas the large distances to the pulsars observable at our frequency band will result in an average scattering effect that is net isotropic. Our results show no evidence of flatter power laws at lower frequencies that might indicate a truncated scattering screen or the effect of the inner scale.
This large body of measured scattering parameters provides a strong basis from which small variations can be determined accurately. For example, extreme scattering events (Walker 2007) could be detected and studied by investigating how varies in observations of the same pulsars at different times. The software we developed to perform the scattering analysis is called SCAMP-I (SCAtter Modelling for Pulsars -I). It is suitable for time domain modelling of single-component scattered pulse profiles and is available for use at https://github.com/pulsarise/SCAMP-I Foundation, an agency of the Department of Science and Innovation. MeerTime data is housed and processed on the OzSTAR supercomputer at Swinburne University of Technology with the support of ADACS and the gravitational wave data centre via AAL. LO acknowledges funding from the UK Science and Technology Facilities Council (STFC) Grant Code ST/R505006/1. AK also acknowledges funding from the STFC consolidated grant to Oxford Astrophysics, code ST/000488. RMS acknowledges support through Australian Research Council Future Fellowship FT190100155.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.
fit at all frequencies, we have calculated / at the highest observing frequency, and plotted the values as a shaded grey histogram. For the pulsars showing the spectral flattening behaviour, we have taken / at the highest frequency for which the values are still following power law behaviour: we treat this as the cut-off point. The histogram of these values is plotted in Fig. A2 as a stepped histogram. We see that the majority of pulsars with no spectral flattening have / < 0.7, in accordance with our theoretical estimates that the parameters can be measured accurately for smaller values of / . By contrast, the histogram of cut-off values peaks at / ∼ 0.8, implying that the flattening we see at higher frequencies is indeed due to the loss of ability to separate the intrinsic profile shape from the exponential scattering. Our theoretical estimate of / = 1 being the cut-off for successful modelling is intended only to be indicative, since it takes no account of the variety of signalto-noise ratios or intrinsic profile shapes of our observations. It is therefore unsurprising that we see large spreads in both histograms in Fig. A2. In particular, the left-hand side of the stepped histogram is consistent with explaining those pulsars which, like PSR J1653-4249, may have extra components that are altering the scattering behaviour of the pulsar in comparison to what we might expect.
These results indicate the extent to which individual measurements of of single profiles are vulnerable to many sources of bias and error. Scattering properties of pulsars can only be characterized reliably in cases where scattering measurements can be performed across large frequency bands. Further, the scattering results for a single pulsar, and the science that can be inferred from them, are best understood in the wider context of all of the other pulsars observed and analysed with the same method.

APPENDIX B: SIMULATION: HOW EXTRA PROFILE COMPONENTS AFFECT SCATTERING TIMESCALE MEASUREMENT
We attempted to replicate the output of the scatter modelling process of PSR J1653-4249 with simulated data, to test whether the presence of a hidden second component in the pulse profile could be responsible for the spectral flattening behaviour seen in this pulsar's modelling results. We simulated a pulse profile made up of two Gaussian components: a large main component and a smaller, narrow secondary component that sits to the right of the main peak. On top of this, we introduced profile evolution: the width of the main component decreases with increasing frequency according to a power law plus constant relationship (Thorsett 1991), and the flux spectral index of the secondary component is flatter than that of the main one, reflecting what is commonly seen in observations. We defined the scattering timescale at the lowest frequency to be the same as that measured for the data. We then applied a spectral index of = 4.4 to to obtain its value at other frequencies. We defined the height, width and position the second component such that the overall profile shapes of the simulation appear similar to those of the real pulsar. Fig. B1 shows how the two components combine to make the simulation profile shapes at the lowest and highest frequencies.
Performing the MCMC fit on this simulation, we obtained results for and that we compare to those of the data in Fig.  B2. The effect of the second component is as expected: its presence does little to alter the fit parameters at the low frequency end, where they are recovered well, but at higher frequencies we see a flattening off in the spectral behaviour of both and that mimics what we witnessed in our measurements for PSR J1653-4249. This lends Figure A1. Plots indicating how the Gaussian and exponential components of our scattering model contribute to the overall model shape, and how that changes with increasing / . Dashed line: the exponential decay function described by equation 3. Dotted line: a Gaussian function. Solid line: the convolution of these two, described by equation A1. Grey vertical line: this marks the phase point at which the combined model intensity (solid line) drops to 10% of its maximum. Grey shading: this indicates the region where the Gaussian (dotted line) is contributing strongly to the overall profile shape (solid line). The definition of a strong contribution is given in the text. Each subplot shows the same set of functions: the changed curve shapes results from the changed ratio of / , which is marked in the top right corner of each subplot. credence to our choice, for pulsars like PSR J1653-4249, to keep only those values of where a power law is still being followed, and infer an from those. Figure A2. Histograms of the ratios of measurement results / . Grey: values for the highest frequency profile fits for the pulsars with a standard power-law -frequency relationship. Transparent, black-edge: values for the highest frequency profile fit that is still following the power-law relationship for those pulsars exhibiting tau-flattening.