© ESO 2022

1 Introduction

Radio pulsars, rapidly rotating and highly magnetized neutron stars, have been a fascinating subject of research since their discovery in 1967 (Hewish et al. 1968). As very compact sources of radio emission, they are also good probes of the interstellar medium (ISM), in particular its ionized component. There are four main propagation effects that occur when the signals from pulsars pass through the ionized interstellar medium (IISM): frequency dispersion, Faraday rotation, interstellar scintillation (ISS), and pulse broadening (see Lorimer & Kramer 2012, for a review). Frequency dispersion and Faraday rotation effects can both be understood by propagation through a homogeneous medium. However, the IISM in general is inhomogeneous and highly turbulent (Rickett 1990), resulting in ISS or pulse broadening (Scheuer 1968). Phase differences in the deflected pulsar signals result in interference observed as intensity fluctuations in frequency and time, similar to the familiar optical “twinkling” of stars caused by the atmosphere of the Earth.

The density inhomogeneities of the IISM at various length scales leads to two different regimes of scintillation, referred to as “weak” and “strong” scintillation. Weak scintillation is weakly modulated, implying that most of the light has traveled almost a single path and that phase perturbations at the observer’s plane are small. Strong scintillation is fully modulated, implying that most of the light arrives via multiple paths and that phase perturbations at the observer’s plane are large. Rickett (1990) provided a scintillation modulation index m that is the standard deviation of the observed flux densities divided by their mean. In weak scintillation, m is much less than unity; in strong scintillation, m is approximately unity or larger than unity. The variations of strong scintillation emerge on two different timescales: the short-term variability (on the order of minutes) caused by changing interference between parts of the scattering disk and long-term variability (on the order of months) with spatial scales on the order of the scattering disk so that observers receive light from a fully different set of scattering points. In this paper, we investigate the short-term variability of the strong scintillation exploiting the long-term observations of pulsars from the European Pulsar Timing Array (EPTA).

Pulsar scintillation properties are best studied using their dynamic spectrum, which is the two-dimensional image of the pulsed intensity as a function of observing time and observing frequency. The interference maxima in the dynamic spectrum are called scintles, where broad and narrow scintles indicate small and large separations of the ray paths, respectively. There are two common approaches that are described in the following to further study scintillation properties. Firstly, employing two parameters, the scintillation bandwidth v_d and the scintillation timescale τ_d, the average characteristics of scintles can be quantified for each observation. τ_d is the half width at 1/e along the time axis and v_d is the half width at half maximum along the frequency axis in the two-dimensional autocorrelation function (ACF) of the dynamic spectrum (Cordes 1986). Secondly, a “criss-cross” pattern can be seen in the dynamic spectrum of some observations (Hewish 1980). In order to investigate such “slopes” or “criss-cross” scintles, Cordes & Wasserman (1984) presented a secondary spectrum by taking the two-dimensional Fourier transform of the dynamic spectrum.

Stinebring et al. (2001) found a faint but clear arc in the secondary spectra of several pulsars at 430 MHz, using high-resolution, high-sensitivity dynamic spectra. Walker et al. (2004); Cordes et al. (2006); Pen & Levin (2014a); Gwinn & Sosenko (2019) elaborated the theory of parabolic arcs in the secondary spectrum. Reardon et al. (2020), Main et al. (2020) and Mall et al. (2022) derived the location and the nature of the scattering screen from fitting arc curvatures for PSRs J0437-4715, J0613-0200 and J1643-1224, respectively. Yao et al. (2021) reported the first evidence for three-dimensional alignment between the spin and velocity vectors based on scintillation arc curvature studies. However, the majority of observations in this work do not show meaningful scintillation arcs. We restrict the analysis to the v_d and the τ_d, which are easier to measure.

In general, v_d is highly frequency dependent. Thus, gathering measurements at multiple frequencies allows one to estimate the scaling index a of the scintillation bandwidth with frequency. By establishing the relation between the scaling index and the spectral index β of the electron density fluctuation spectrum versus wavenumber, the turbulence characteristics of the IISM can be inferred (Goodman & Narayan 1985). Rickett (1977) provided that the scaling index of the scintillation bandwidth with frequency is α = 4.0 in the thin screen model. Another more commonly used scaling index is α = 4.4 for Kolmogorov turbulence with a spectral index β = 11/3 in a uniform medium (Romani et al. 1986). Gupta (2000) summarized which scintillation studies favoured, or disfavored, a pure Kolmogorov spectrum. To explain events that disfavored a pure Kolmogorov spectrum, Lambert & Rickett (1999, 2000) presented a steeper spectrum with a spectral index of β = 4, this type of spectrum could be caused by a medium with abrupt density changes.

However, τ_d is weakly frequency dependent but more modulated by the transverse velocity of the pulsar, Earth and IISM. Lyne & Smith (1982) presented some observations which show that the scintillation characteristics of pulsars are closely related to their velocities, as measured by their proper motions. In turn, measurements of scintillation also can be used to estimate the pulsar’s transverse velocity. Subsequently, Gupta (1995) presented a correlation between pulsar proper motion velocities and scintillation velocities of radio pulsars, and Cordes & Rickett (1998) extended it for different scattering geometry models. In some cases, the time series of τ_d show an annual variation, and for binaries, it is possible to show orbital-annual variations. By employing these variations, the small-scale distribution and inhomogeneities of the IISM, or some orbital parameters can be determined. For examples, Linsky et al. (2008) investigated the transverse velocity and the distribution along the line of sight of plasma clouds using annual variations from three quasar radio sources, Rickett et al. (2014) exploited annual and orbital variations of scintillation timescales in the double pulsar system J0737–3039 to estimate the location, spatial structure and transverse phase gradient of a scattering screen as well as two orbital parameters: the inclination angle i and the position angle of the ascending node Ω. Reardon et al. (2019) also used the annual and the orbit variations to resolve the previous ambiguity in the sense of the i, the Ω and measurements of the spatial structure of the interstellar electron-density fluctuations, including: the spatial scale and anisotropy of the diffraction pattern, the distribution of scattering material along the line of sight, and spatial variation in the strength of turbulence from epoch to epoch.

We describe our observations and analysis methods in Sects. Sect. 2 and 3, respectively. The measurements of scintillation parameters are presented in Sect. 4 and the comparison between our results and previous studies and NE2001 are given in Sect. 5. In Sect. 6, we discuss the significance of our results in the context of PTAs. We present the conclusion and further research in the last section.

2 Data

The EPTA (Desvignes et al. 2016) is a collaboration of European research institutes and radio observatories that was established in 2006 and makes use of the five largest telescopes at decimetric wavelengths in Europe: the Effelsberg 100-m Radio Telescope (EFF), the Lovell Radio Telescope at the Jodrell Bank Observatory (JBO), Nançay radio telescope (NRT), the Westerbork Synthesis Radio Telescope (WSRT) and the Sardinia Radio Telescope (SRT). The EPTA aims to study the astrophysics of millisecond pulsars (MSPs) and to detect cosmological gravitational waves in the nanohertz regime.

In this work, we aim to analyze long-term scintillation observations for EPTA pulsars. To define the sample of pulsars, we considered all EPTA sources included in Desvignes et al. (2016) that were expected to have an easily resolvable scintillation bandwidth, either based on earlier scintillation studies, or based on model predictions given by the NE2001 model (Cordes & Lazio 2002). This resulted in the 13 MSPs listed in Table 1. They span spin periods from 1.56 to 16.45 ms with a DM range of 9–71 pc cm⁻³. The observations used for those 13 pulsars are mainly from continued monitoring of the NRT telescope that has the largest observing bandwidth, and partly from the EFF, JBO and WSRT telescopes. As the SRT has only recently commenced

MSP monitoring, no data from this telescope are included. We describe the observations used in this paper from each telescope in the following subsections, and the relevant details of the observations after removal of the radio-frequency interference (RFI) are included in Table 1.

2.1 Nançay Radio Telescope

The 94-m equivalent NRT is located in the small commune of Nançay, two hours’ drive south of Paris, France. Two receivers: the low-frequency receiver (1.1–1.8 GHz) and the high-frequency receiver (1.7–3.5 GHz) with a total observable bandwidth of 512 MHz centered on a selectable center frequency are used, most of the observations with the low-frequency receiver were carried out at a center frequency of 1484 MHz. All NRT data in this work was recorded with the new wideband NUPPI dedispersion backend (Liu et al. 2014) that has been using since August 2011. Data are recorded with a subintegration length of 15, 31 or 61 s and a channel bandwidth of 4 MHz. The observations used in this work were made between 2011 August and 2020 January. Each observation’s duration is around one hour.

2.2 Effelsberg Radio Telescope

The EFF 100-m radio telescope of the Max-Planck Institute for Radioastronomy (MPIfR) is located about 1.3 km northeast of Effelsberg, a southeastern part of the town of Bad Münstereifel in Germany. Observations used in this paper come from two receivers: a single-pixel receiver with a center frequency of 1347.5 MHz and 200 MHz of bandwidth (observations prior to 2017), and a multibeam receiver with a center frequency of 1397.5 MHz and a bandwidth of 400 MHz (2017 onward). Data are recorded with a channel bandwidth of 1.5625 MHz and a subintegration length of 10 s. The data used in this work have been regularly collected by the PSRIX (200-MHz data) and Automatix (400-MHz data) systems (Lazarus et al. 2016) from 2011 March to 2020 June with typical observation lengths of 20–50 min. As scintillation strength is highly frequency-dependent and the data were not averaged in frequency, in order to produce as homogeneous a data set as possible, all observations observed with the multibeam receiver were downsized to 200 MHz of bandwidth with a center frequency of 1347.5 MHz by cutting out the edges of the band. Additionally, some pulsars are observed at the 11-cm band with a center frequency of 2627 MHz and a bandwidth of 200 MHz from 2011 March to 2015 July, data for those pulsars are recorded with a channel bandwidth of 1.5625 MHz and a subintegration length of 10 s.

2.3 Lovell Radio Telescope

The Lovell telescope is a 76.2-m radio telescope at the JBO, near Goostrey, Cheshire in North West England. Pulsars in this work are regularly observed and we used the data from 2011 April to 2020 January. All observations are recorded with a bandwidth of 400 MHz, a channel bandwidth of 1.5625 MHz and a subintegration time ranging from 10 to 120 s. Observation durations are approximately 30 minutes.

2.4 Westerbork Synthesis Radio Telescope

The WSRT is an aperture synthesis interferometer located north of the village of Westerbork in the northeastern Netherlands. It consists of a 2.7-km east-west array of fourteen 25-m dishes, adding up to a collecting area equivalent to that of a 94-m dish when combined as a tied array. We only use data on three pulsars from this telescope, one (PSR J0636+5128) was observed from 2014 January until 2015 June with an observation length of one hour, a subintegration time of 60 s, a bandwidth of 80 MHz centered on a frequency of 350 MHz and a channel bandwidth of 0.156MHz. The other two pulsars (PSRs J1939+2134 and J1713+0747) were observed from 2007 February until 2015 June for around 30 min per observation, with the same subintegration time of 60s, and a bandwidth of 160MHz centered on a frequency of 1380 MHz and with a channel bandwidth of 0.3125 MHz.

3 Methods

3.1 Scintillation Parameters

If radio-frequency interference (RFI) is impulsive and broadband or narrow-band and persistent, it will show up as an extended feature in the dynamic spectrum – in which case it is likely to have a strong influence on the derived scintillation parameters. Therefore, before getting the dynamic spectrum, RFI rejection is necessary. To distinguish RFI from pulsar signals in each time and frequency bin, we calculated the pulse profile’s Pearson autocorrelation coefficient in the off-pulse area with a phase lag of one bin. The pixels in the dynamic spectrum with an autocorrelation coefficient at least Lσ above the mean (where σ is the standard deviation of the correlation coefficients obtained), were identified as RFI and replaced with zeroes. L is a tunable parameter in the range ~ 1.1–2, which was manually adjusted for optimal performance as it is strongly dependent on the amount of RFI present. In addition the highest and lowest frequency channels were removed since a combination of decreasing telescope sensitivity and increased potential of systematic effects made these channels unreliable. In consequence, the effective center frequency and bandwidths given in Table 1 differ from those of the raw data described in Sect. 2.

After dedispersing, we formed total intensity (Stokes I) profiles for each channel and subintegration by summing polarizations, producing a data cube I (f, t, phase) from each folded archive. A template was obtained by averaging all intensity profiles over the entire observing duration and over the usable range of the observing bandwidth for every observation. The on- and off-pulse regions were determined by the signal-to-noise ratios (S/Ns) of the template bins, with the on-pulse region having S/N > 7, and vice-versa for the off-pulse region with S/N ≤ 7.

To acquire only scintillation patches in the dynamic spectrum, we removed the frequency dependence of our observational sensitivity, which is caused by the combination of instrumental sensitivity and pulsar spectral index. To remove the instrumental sensitivity effects in different frequency channels, we divided the intensity in each channel by the time average of the off-pulse region across the full observation. However, some bright pulsars have emission in the off-pulse region, which contributes to off-pulse emission getting confused with scintillation (Ravi & Deshpande 2018). Thus, we averaged many observations to remove the instrumental bandpass for the bright pulsars. After removing the instrumental bandpass, we averaged the power of the on-pulse region based on as many observations as possible to calculate and then remove the effect of the pulsar’s spectral index. Furthermore, in each time and frequency bin, we subtracted the mean of the off-pulse region to remove variable background flux. Finally, the remaining power in the on-pulse region in each time and frequency bin was used as a point in the dynamic spectrum. Four examples of normalized dynamic spectra are presented in the left panels of Fig. 1, showing the wide range of scintle sizes between our sources.

Following Cordes (1986), in order to quantify the average characteristics of scintles at each observation, τ_d and v_d have been defined, where τd is the half width at 1/e along the time axis and v_d is the half-width at half maximum along the frequency axis in the two-dimensional ACF of the dynamic spectrum. Gupta et al. (1994) defined the ACF as $$\text{ACF}\left(\text{Δ}f,\text{Δ}t\right)=\frac{R\left(\text{Δ}f,\text{Δ}t\right)}{R\left(0,0\right)},$$(1)

where R(Δf, Δt) = <ΔI(f, t) · ΔI(f + Δf, t + Δt) >, ΔI(f, t) is the deviation of the intensity from the mean intensity <I> of the

dynamic spectrum, and <> denotes an average over the range in frequency and time from one observation. In the calculation of the global mean intensity of the dynamic spectrum <I>, we masked all zeros that resulted from RFI removal. Then, we used a convolution algorithm that utilises the fast Fourier transform to do the calculation of the ACF. Due to the finiteness of the data, the mentioned convolution algorithm automatically padded the arrays ΔI(f, t) and ΔI(f + Δf, t + Δt) with zeros, which results in a different number of samples contributing to various pixels of the ACF. To correct this, we normalized the ACF to take into account the different numbers of contributions to each pixel. Some pulsars, especially weaker sources, show a sharp spike at zero frequency and time lag of the ACF, resulting from the noise in the dynamic spectrum. To remove that spike, we replaced it with the mean value of the nearest two points along the time lag axis.

As in previous works, cuts along two axes of the ACF (called 1D ACF) and a Gaussian function (Gupta et al. 1994) were used to obtain the v_d and τ_d, $$\begin{array}{c}\text{ACF}\left(v=0,\tau \right)=\exp \left(-a\cdot {\tau }^2\right),\\ \text{ACF}\left(v,\tau =0\right)=\exp \left(-b\cdot v^2\right).\end{array}$$(2)

We determined the fitting range from the center of the ACF to the point where the uncertainty of the fitting parameters from the least-squares fits is the smallest. For the determination of scintillation timescales, we tried the theoretical form for a Kolmogorov scattering medium ACF(v = 0, τ) = exp(−a · τ^5/3) (Coles et al. 2005, 2010; Reardon et al. 2019). However, we found that the commonly used Gaussian function gives a better fit to the shape along the time axis of the ACF of our pulsars.

Then, the scintillation parameters are derived with $$\begin{array}{l}{\tau }_{\text{d}}=\sqrt{\frac{1}{a},}\hfill \\ v_{\text{d}}=\sqrt{\frac{\ln 2}{b}.}\hfill \end{array}$$(3)

Usually, the uncertainty of the scintillation parameters consists of the formal uncertainty from the fitting procedure and the statistical uncertainty due to the finite number of scintles in the dynamic spectrum. Following Cordes et al. (1985), the statistical fractional uncertainty, σ_est, is given by $${\sigma }_{\text{est}}={\left[f_{\text{d}}\cdot \left(\frac{\text{BW}\cdot {\text{t}}_{\text{obs}}}{v_{\text{d}}{\tau }_{\text{d}}}\right)\right]}^{-0.5},$$(4)

where BW is the bandwidth of the observation, t_obs is the length of the observation and f_d is the filling factor. The filling factor is an estimate of the fraction of the dynamic spectrum that is filled with scintles. For our 13 pulsars, depending to the size of their scintles, we used filling factors from 0.2 to 0.4 for different pulsars, please see the details in Table 1. We plotted the ACF of four representative dynamic spectra in the middle panel of Fig. 1. In the side panels, cuts along the two axes of the ACF are plotted, the black lines are the best Gaussian fits.

However, as can be seen in Fig. 1, since the samples in the ACF are highly correlated, the uncertainties derived from the ACF are heavily underestimated. In order to reflect the true Gaussian variations (as a proxy for the measurement uncertainty) of measurements in Fig. 2, we:

1. Took the Fourier Transform (FT) and then the power spectrum (PS) of the time series of the v_d (or τ_d);

2. Used the PS to identify frequencies with excess power;

3. Removed these frequency channels from the FT, in practice, the remaining FT now describes the white noise (WN) in our data without containing any of the non-white variations;

4. Determined the inverse Fourier Transform of the WN back to the time domain;

5. Calculated the Root Mean Square (RMS) of the whitened data;

6. Added that RMS in quadrature to the error bars.

In all cases, this additional measurement uncertainty dominated the resulting error bar of the measurements.

Fig. 1 Dynamic spectra I(f, t) from four observations with different scintle sizes in the left panels and their normalized autocorrelation functions in the right panels. The horizontal and vertical axes correspond to the time (or time lag) and frequency (or frequency lag), respectively. In the two smaller side plots, the sky blue points are the 1D autocorrelation functions at zero time and frequency lag, the black lines are the best Gaussian fits, the red lines indicate the values for vd and τd. The pulsar name, telescope, epoch of observation, vd and τd are given above the dynamic spectra.

Fig. 2 Long-term variations of scintillation bandwidth vd, scintillation timescale τd and DM for 13 EPTA pulsars. To ensure that the variations shown by the most precise DM values remain visible, for most pulsars, we only show the most precise DM measurements, which are usually derived from 21-cm Nançay data. The different colors indicate measurements at different telescopes, following the legend given in the top-right corner of each sub plot. Following the telescope name, we present the center frequency in units of MHz. Downward-pointing and upward-pointing triangles indicate unreliable estimates (i.e., upper and lower limits) caused by limited frequency resolution and limited observation length, respectively. The vertical gray dotted lines indicate the start of a calendar year.

3.2 The Time Series of DM

Fluctuations in electron density that give rise to scattering should also contribute to variations in DM. If the scattering angle is small, as is typical in the 21-cm band, changes in DM and scintillation are likely to happen at roughly the same time. We therefore expected to find some relations between scintillation parameters and DM. For example, McKee et al. (2018) demonstrated a correlation between scattering and dispersion measure variations in observations of the Crab pulsar (with a correlation coefficient of 0.56 ± 0.01), and Coles et al. (2015) reported two interesting events showing increases in both the intensity of scintillation and dispersion toward PSRs J1603–7202 and J1017–7156, which they referred to as extreme scattering events (ESEs).

Thus, the observations with RFI removed are used to get the time series of DM by using the PSRCHIVE (Hotan et al. 2004) and TEMPO2 (Hobbs et al. 2006) software packages. We first used a fully frequency-averaged high-S/N observation to produce an analytic template by employing the PSRCHIVE program PAAS. Dependent on the mean S/N of each pulsar, then, all observations are scrunched to 10–20 channels. The analytic template and the PSRCHIVE program PAT are used to yield ToAs, whereby each channel yields one ToA. With the TEMPO2 software package and a good-quality timing model (taken from Desvignes et al. 2016), for each pulsar, we removed all ToA outliers and all ToAs with large error bar (several times larger than the median). Finally, the remaining ToAs of each observation are used to fit for the DM at that observing epoch.

The differential time delay across the effective bandwidth, which is naturally related to the DM measurement, $$\text{Δ}t\propto \text{DM}\frac{2\text{BW}}{f_2^3},$$(5)

where BW and f_c are the bandwidth and center frequency in MHz, Δt is expressed in seconds and DM is in units of pc cm⁻³ (Verbiest & Shaifullah 2018). That means the wider BW and the lower center frequency can yield the higher precision DM measurement. For our EPTA data, therefore, a wide range in DM precision can be obtained by fitting DM to individual observations. Consequently, we only use the most precise DM measurements obtained, which are typically from 21-cm Nançay data.

4 Results

According to the method described in Sect. 3, we present the median of v_d and τ_d for 13 pulsars in Cols. (11) and (12) of Table 1, the sub and super scripts are the 5/95 percentiles. To analyze the long-term variations of pulsar scintillation, we plotted the time series of v_d and τ_d of all pulsars in the first and second panels of Fig. 2, respectively. The time series of DM for each pulsar was plotted in lowest panels of Fig. 2. Some pulsars have been observed at both 21-cm and 11-cm bands, in which case the significant difference in scintillation parameters causes the 21-cm values to be compressed at the bottom of the scale, with any variations invisible. To remedy this, for pulsars with measurements at different bands, the scintillation parameters have been plotted on a logarithmic scale. If the number of observations in a given observing band is too small for an analysis for the variation of scintillation parameters, we did not plot the measurements from this observing band in Fig. 2, for example, the measurements of PSR J1857+0943 at 11-cm band. The data for all our measurements are publicly available on zenodo¹.

In addition to the measured scintillation parameters, the relevant observation information is also listed in Table 1, including the effective center frequency f_c, the number of observations N_obs, the channel bandwidth (CHBW), the effective bandwidth (BW), the subintegration length t_sub and the mean observation length ${\overline{t}}_{\text{obs}}$. The effective center frequencies and the effective bandwidths in Table 1 vary slightly as described in Sect. 2 because edge frequency channels were sometimes removed due to RFI.

Observations for some pulsars at some telescopes could not be used to derive scintillation parameters and were therefore not included in the results presented in Table 1 and Fig. 2. We describe the specific details of which observations were not suitable below: (i) due to limited frequency resolution, some data sets were not able to resolve any scintles in the dynamic spectrum, for examples, WSRT observations of PSR J0636+5128 at a center frequency of 346.25 MHz and NRT observations of PSR J1600–3053 at a center frequency of 1484MHz; (ii) data sets with narrow observing bandwidth and weaker scintillation, where one full scintle is larger than the available frequency range, for examples, EFF observations of PSR J2145–0750 at a center frequency of 2627 MHz and NRT observavtions of PSR 1713+0747 at a center frequency of 2539 MHz; and (iii) some observations showed a lot of RFI and the dynamic spectrum was badly damaged after RFI removal, which seriously affects the derived scintillation parameters. If more than half of a set of observations are accompanied by severe RFI, we abandoned all observations from that set, for example, JBO observations of PSR J2317+1439 at a center frequency of 1532MHz. If only a few observations in a set showed extensive RFI, we just removed those observations.

As shown in the first panel of Fig. 1, in some observations, the scintles are not fully resolved by the frequency resolution in the dynamic spectrum. For those observations, usually, the 1D ACF has less than three points within the half-power width, which leads to an unreliable (usually overestimated) measurement of v_d. Besides, we found that almost all v_d of those observations are smaller than 1.5 times the channel bandwidth. Thus, if the value we derived for v_d is smaller than 1.5 times the channel bandwidth, we considered this value to be unreliable, use it as an upper limit and use a downward-pointing triangle to plot it in Fig. 2. If, in a given data set, 5% or more of the v_d measurements are upper limits, then the 5th percentile is unreliable, hence the subscript is shown in red in Table 1; similarly, the median value for v_d is shown in red if 50% or more of measurements are upper limits; and the 95th percentile (superscript) is shown in red if more than 95% of the measurements are upper limits.

Some observations cannot show one full scintle in the time axis of the dynamic spectrum because of the limited observation duration, as shown in the fourth plot of Fig. 1: in this case the minimum value of the ACF along the time axis is larger than 1/e. For those observations, we fit all points of the time axis and extrapolated the result to the delay where the time lag is smaller than or equal to 1/e, and then we used the half-width at 1/e along this virtually extended time axis of the ACF as the value for τ_d. However, the τd estimated in this way is larger than the observing duration t_obs and is unreliable, so we used an upward-pointing triangle to plot them in Fig. 2. Similarly to the treatment of v_d, values for τ_d are shown in red in Table 4, depending on what fraction of measurements were considered to be limited by the observing length of the data. Contrary to the treatment of v_d, we sorted τ_d in descending order and we marked the 5th percentile as the superscript (the 95th percentile as the subscript).

5 Comparison

In this section, we compare the results presented in Table 1 first against previously published observational results and subsequently against the theoretical predictions of the NE2001 Galactic electron density model.

5.1 Comparison with Previous Studies

Many previous studies on the scintillation properties of pulsars have been published. The sources in our analysis overlap most with those of Levin et al. (2016) and Turner et al. (2021). Also, our data span almost matches their data span. It is, therefore, worthwhile to compare our scintillation parameters with theirs. However, their scintillation parameters were determined at different radio frequencies. To compare them with our values, their vd values were rescaled to our observing frequency by using a scaling index α = 4.4. For the scintillation timescale, we assumed that τ_d is proportional to f^1.2 (Cordes & Lazio 2002). The comparison values are listed in Table 2. The sub and super scripts are calculated by first adding and subtracting the uncertainties that were given in literature and then rescaling.

For a visual comparison, we plotted the v_d from previous works and our measurements from Table 2 in the left panel of Fig. 3. Both vertical and horizontal error bars of each point are asymmetric, the smallest and biggest values of error bars are the subscript and superscript in Table 2. Since we used the 5 and 95 percentiles of all measurements as sub and super scripts for each data set in this work in Table 2, a whole horizontal error bar indicates the range of 90% of the probability density. The gray line is the equivalence diagonal, which means that those points where our measurements are in keeping with the previous works lie on or near this line. Comparing data at both of our center frequencies with both literature sources, we created four groups and plotted them with a different color.

To examine the probability density of the difference between our measurements and those from literature for the four groups we mentioned in previous paragraph, we plotted four violin plots of the differences, Delta, in the right panel of Fig. 3, where $$\text{Delta\hspace{0.17em}}=\frac{M_{\text{this\hspace{0.17em}work}}-M_{\text{literature}}}{\sqrt{{\sigma }_{\text{this\hspace{0.17em}work}}^2+{\sigma }_{\text{literature}}^2}}$$(6)

with M the median of the v_d for a given pulsar and at a given observing frequency, and σ is its measurement uncertainty. We used the standard deviation of our measurements as the uncertainty for our own measurements. For the values from literature we determined the upper (lower) uncertainty as the difference between the superscripts (respectively the subscripts) and the most likely value quoted in Table 2. These upper and

lower uncertainties are then combined by taking the root-mean-square (rms) of them; and this rms value is used as σ_literature in Eq. (6).

Comparison of the scintillation bandwidths in both panels of Fig. 3 shows that, at the 21-cm band, our values are typically consistent with those from literature, especially for Turner et al. (2021). In the 11-cm band, the previously published values show a notable bias towards small values compared to the values from our work, but they are still consistent within their error bars. There are some possible reasons to explain those discrepancies:

1. The strength of scintillation changes over time, as the time series of scintillation bandwidth of PSR J0636+5128 in Fig. 2 illustrates, where the values after MJD 57600 increased significantly;

2. Scaling indices can be expected to be different for different pulsars, and may also vary with time (Bansal et al. 2019). A scaling index of α = 4.4 based on the classical Kolmogorov spectrum represents a mean value for all pulsars, but most pulsar lines of sight do not follow a pure Kolmogorov spectrum and show a smaller or larger value of scaling index. So errors can easily be introduced when we rescale the scintillation bandwidth to another frequency using a scaling index of α = 4.4; the size of these errors is proportional to the frequency range scaled. That is why there is a larger bias at the 11-cm band;

3. Unreliable measurements affected by limited frequency resolution and limited bandwidth also contribute to discrepancies. Usually, a higher frequency resolution and a wider bandwidth will produce more accurate measurements;

4. The different data analysis methods can also introduce slight differences.

Usually, the fourth reason only leads to a tiny discrepancy that is within the error bar. For the different data span, the first reason may lead to a significant discrepancy. For example, PSR J0636+5128 show a visible increase round MJD 57600 in the time series of v_d and τ_d, and the previously published works did not cover data after 2017. In this case, the increased monitoring length introduces a discrepancy. Luckily, our data span almost fully covers the data span of Levin et al. (2016) and Turner et al. (2021), and does not extend much beyond theirs, so that in the majority of cases this first reason does not cause discrepancies.

However, there are some significant discrepancies between our work and previous publications, which are highlighted by three ellipses in the left panel of Fig. 3. The points in the red ellipse show that the previous measurements are larger than our own, which is caused by the third reason, the limited frequency resolution. Since we have a higher resolution for the JBO and WSRT, we can get a more accurate measurement when scintles are very narrow in frequency. The points in the blue ellipse are more likely caused by the second reason that is an incorrect scaling index. For the points in the yellow ellipse, the second and third reasons both are responsible for the significant discrepancies.

The finite observation length is the main cause restricting the accuracy of τ_d in both our own and previously published analyzes. All observations of Levin et al. (2016) and Turner et al. (2021) were from NANOGrav and were ~30min in duration, which makes the calculation of τ_d problematic. Consequently, Levin et al. (2016) did not present values for τ_d, while Turner et al. (2021) only gave τ_d for a few pulsars. Most of our observations have a duration of longer than 30 min, which allows us to get reliable τ_d values for more pulsars.

There are four pulsars that have τ_d values derived in both Turner et al. (2021) and our work. Comparison of these values displays slight discrepancies, which are caused by a number of reasons (identical or similar to those mentioned for v_d): (a) the strength of scintillation changes over time, (b) an incorrect scaling index, (c) the limited observation duration, (d) the different data analysis methods used. In addition, τ_d is dependent on the Earth’s velocity and pulsar orbital velocity, both of which vary with time. Therefore, for the limited number of observations, the different observing times also likely cause some discrepancy.

We also tried using another widely used scaling, 4.0 for the v_d and 1.0 for the τ_d, to rescale the previous publish values. Although there are some tiny offsets in the rescale values, the conclusion about the comparison has not been changed.

Fig. 3 The left panel compares the vd from previous works to our measurements. Here, “L” indicates Levin et al. (2016), “T” refers to Turner et al. (2021), “11 cm” compares against measurements we made in the 2-GHz range; “21 cm” refers to measurements we made in the 1-GHz range – in both cases the previously published measurements were rescaled to our observing frequency, as given in Col. (3) of Table 2. Data from Levin et al. (2016) were at 1500 MHz while Turner et al. (2021) used data at 1500 MHz as well as 820 MHz, but we only compare to the former. The three ellipses indicate groups of points that show significant offsets, which are discussed in more detail in the text. The right-hand panel shows violin plot for the distribution of the differences (Delta) of the four groups between the vd of our measurements and the vd from previous studies.

5.2 Comparison with NE2001

We also compared our measurements with the predicted scintillation parameters from NE2001 (Cordes & Lazio 2002) which is a model for the Galactic distribution of free electrons and can be used to estimate the scintillation properties for a given sky position and DM value. Since NE2001 only gives the predictions at 1 GHz, we used a scaling index of α = 4. 4 to rescale vd to our observing frequencies. Note that the value of τ_d depends not only on the frequency but also on the pulsar transverse velocity, as deduced in Eq. (46) of Cordes & Lazio (1991). In the prediction

of τ_d, the NE2001 Electron Density Model assumes that all pulsars have a pulsar transverse velocity of 100 km s⁻¹. In practice the pulsar transverse velocity ranges from a few km s⁻¹ to several thousand km s⁻¹. In addition to using a frequency scaling index of 1.2, the estimated pulsar transverse velocity also has been used to rescale τ_d, as in Eq. (46) of Cordes & Lazio (1991). Here, the pulsar transverse velocity was calculated by $$V_{\mu }=4.74\mu D_{\text{kpc}},$$(7)

where V_µ is in units of km s⁻¹, µ is the proper motion with units of mas yr⁻¹, and D_kpc is the pulsar distance in kpc. All the reference values and reference papers are listed in Table 3.

Then the rescaled values of v_d and τ_d from the predictions of NE2001 are listed in Table 2, the subscripts and superscripts are the predicted values at the bottom and top of our bandwidth, respectively. Similar to the treatment of the comparison with previous works, we plotted the predictions from NE2001 vs. our measurements in Fig. 4 to compare. According to the effective center frequency of observations, we created two groups: 11-cm band and 21-cm band.

As shown in Table 2 and Fig. 4, our measurements for v_d are systematically smaller than the predictions for most pulsars. An incorrect scaling index or limits to our observing bandwidth and resolution could affect these results, although they cannot fully explain this apparent bias. An alternative explanation is a potential sample bias in our work, since the sample was defined by high NE2001 predictions for v_d. Finally, it is possible that NE2001 overestimates scintillation bandwidth for nearby pulsars, although this cannot be rigorously determined based on our limited sample.

For τ_d, some pulsars show that our measurements are smaller than predicted, some pulsars show our measurements are consistent with predictions, some pulsars show our measurements are larger than predicted. Except for an incorrect scaling index and the limited duration of our observations, an incorrect estimate of the pulsar transverse velocity can also contribute to significant discrepancies. In addition, almost all pulsars show that the range of variation from our measurements is larger than those from predictions, which is because, in the prediction of τ_d of the NE2001 model, the Earth’s velocity and pulsar orbital velocity have not been taken into consideration.

Fig. 4 The NE2001 predictions of scintillation parameters vs. our measurements. The left panel shows the comparison of vd, the right panel shows the comparison of τd.

6 Discussion

Based on the measurements presented in Sect. 4, we found some interesting features in the time series of scintillation parameters. In this section, we describe future lines of research based on these features.

6.1 Annual variations of the τ_d

The time series of τ_d shows a strong annual variation for PSRs J0613−0200 and J0636+5128. The scintillation timescale τ_d is given by the ratio of the spatial scale of the diffractive scintillation pattern l_d and the scintillation velocity V_ISS (Cordes & Rickett 1998): $${\tau }_{\text{d}}=\frac{l_{\text{d}}}{V_{\text{ISS}}}.$$(8)

In a thin screen scattering geometry model, the scintillation velocity V_ISS can be determined by the effective transverse velocity V_eff and the fractional distance of the scattering screen s (the pulsar at s = 0 and the Earth at s = 1), V_ISS = V_eff/s (Cordes & Rickett 1998). The V_eff is a combination of the pulsar, Earth, and IISM velocities at scintillation screen position s, $$V_{\text{eff}}\left(s\right)=s{V}_{\text{E}}+\left(1-s\right)\left(V_{\text{p}}+V_{\mu }\right)-V_{\text{IISM}}\left(s\right),$$(9)

where V_E is the Earth velocity, V_p is pulsar orbital transverse velocity and V_IISM is the IISM velocity. Eqs. (8) and (9) clarify that variations in τ_d are determined by six parameters: l_d, s, V_E,V_p, V_µ and V_IISM. Annual variations in τ_d arise from the changes in the transverse velocity of the Earth.

Although the Earth’s motion has contributions to τ_d for all pulsars, not all pulsars showed annual variations in τ_d. Table 4 summarizes which pulsars show signs of annual variations in τ_d. To distinguish under which conditions annual variations can be observed, in Table 4 we also list V_µ (which is the same as in Table 3) and the mean orbital velocity V_o (which was calculated based on the orbital period P_b and the projected semi-major axis A₁, which in turn are available from pulsar timing). The orientation Ω and the inclination angle i of the pulsar orbit are difficult to estimate and are known only for a few pulsars, therefore, we calculated the V_o, not V_p in Table 4. From Table 4, one sees that all pulsars that showed annual variations had relatively small V_µ and V_o. But the reverse is not true: some pulsars with small V_µ

and V_o did not show annual variations, likely because of the location of the scattering screen s reduces the amplitude of annual variations. Indeed, from Eqs. (8) and (9), we can determine that pulsars with a small V_µ, a small V_p and a large s are more likely to show annual variations in τ_d. In case of a large s, the contribution of pulsar velocities on V_eff has been reduced. On the other hand, a large s indicates that the scintillation screen is close to the Earth. Plasma clouds close to the Earth are likely related to two features: the edge of the Local Bubble (LB)³ and the local interstellar medium (LISM)⁴, and their velocities are supposed to be of the order of 10 km s⁻¹.

PSRs J0613−0200 and J0636+5128 both show a very strong annual variation in τ_d, which can be used to estimate the small-scale distribution along the line of sight and inhomogeneities of the IISM, like has been done by Rickett et al. (2014) and Reardon et al. (2019). This analysis will be published in a forthcoming paper.

6.2 Correlations between v_d and τ_d

The time series of scintillation parameters v_d and τ_d show correlated variations in a number of cases:

1. PSR J0636+5128 with a strong annual variation in τ_d shows also a weak annual variation in v_d;

2. The v_d and the τ_d of PSR J0636+5128 both show a visible increase along with periods of persistent highs around MJD 57600;

3. The v_d and the τ_d of PSR J1939+2134 both show dips in the two regions between the vertical red lines of Fig. 5. These two events are most clearly detected in the WSRT data due to their high frequency resolution, so the WRST data are plotted separately in Fig. 5.

The occurrence of correlated variations is not unexpected, since effectively both v_d and τ_d depend on the scale of the diffractive scintillation pattern, l_d, as described by Cordes & Rickett (1998). The consequence is, that if anisotropic scattering causes spatial variations in l_d; and if either the screen is sufficiently close to Earth or the variability of l_d exists on a sufficiently small scale that the line of sight probes parts of the IISM with significantly different l_d over the course of a year, then annual variations would be expected in both τ_d and v_d.

The dependence of τ_d on 1_d was already given in Eq. (8), the dependence of v_d on l_d is given by (Cordes & Rickett 1998): $$v_{\text{d}}\propto l_{\text{d}}^2{f}^2,$$(10)

where f is the observing frequency. Thus, if the annual variations enter l_d, we can also expect an annual variation in the time series of v_d. But the annual variation of v_d is far weaker than that of τ_d because annual variations in l_d are likely small, while those in V_iss can be large.

Equations (8) and (10) also can be used to explain the second and third synchronous variations between v_d and τ_d. When some unpredictable changes happen for l_d, which could happen when particular plasma clouds dominate the ISS at different times, v_d and τ_d would both be affected and show similar changes in a similar time.

6.3 Frequency dependence

After ruling out all limits of v_d and τ_d due to the limited frequency resolution and the limited observing duration, as we expected, Table 1 and Fig. 2 both show that measurements of v_d and τ_d increase with observing frequency in general, and the frequency dependence of v_d is far stronger than that of τ_d. Since the frequency dependence of τ_d is very weak and τ_d is modulated by the Earth’s and the pulsar’s velocities which vary strongly with time, we only investigate the scaling with frequency of v_d. The strength of scintillation fluctuates with time and we do not have synchronous observations at the different observing frequencies. To remedy the impact of that, only measurements of v_d at the different bands and with a time interval to the closest observation of less than 7 days have been used to estimate the scaling index. Then, the frequency scaling index α of v_d is given by: $$\alpha =\frac{\log \left(\raisebox{1ex}{$v_{d1}$}\!\left/ \!\raisebox{-1ex}{$v_{d2}$}\right.\right)}{\log \left(\raisebox{1ex}{$f_1$}\!\left/ \!\raisebox{-1ex}{$f_2$}\right.\right)}$$(11)

where the v_d1, v_d2 and the f₁, f₂ are scintillation bandwidths and observing frequencies for two bands, respectively. Note that if

a set at a given center frequency and a given telescope have 50% or more limits (as opposed to measurements), we disregard the whole data set here. In Fig. 6, we plot the time series of the scaling index for PSRs J1857+0943 and J1939+2134. One sees significant correlated variations, suggesting these variations are not random on these timescales. The mean and standard deviation of the a for PSRs J1857+0943 and J1939+2134 are 2.82 ± 1.98 and 3.18 ± 0.60, respectively.

Fig. 5 Time series of scintillation bandwidth vd, scintillation timescale τd and DM for PSR J1939+2134. There are two events, that show decreases in vd and τd, and increases in DM, which are indicated by the red lines. The first one is between MJDs 55807 and 55845, the second one is between MJDs 56000 and 56090.

Fig. 6 The time series of frequency scaling indexes of the vd for two pulsars. The mean and standard deviation of scaling indexes are printed in the right-bottom corner.

6.4 The correlation between scintillation and DM

In Fig. 2, we see two interesting instances of correlated variations between scintillation parameters and DM. Firstly, during the period of the two dips in v_d that are mentioned above in Fig. 5, the DM shows two increases. Coles et al. (2015) detected similar events in PSRs J1603−7202 and J1017−7156, they exploited those events to make more complete models of the extreme scattering event (ESE), including an estimate of the “outer-scale” of the turbulence in the plasma lens. Secondly, similar to the two events mentioned above, in Fig. 5, PSR J1939+2134 show an inverse relation between τ_d and the DM through much of our data set, when τd increases, the DM decreases, and vice versa.

In the lens model of scintillations (Pen & King 2012; Pen & Levin 2014b), changes in DM and the refracting angle δ of scintillation can both be used to infer the required electron column densities. Lin et al. (2021) have already used this to calculate transient pulse profile changes expected to occur in concert with DM changes, and compared these with what is observed during two “dip” events of DM for PSR J1713+0747. More details will be analyze in a forthcoming paper.

6.5 Impacts on timing

Pulsar Timing Arrays (PTAs) aim to detect gravitational waves through spatially correlated timing residuals between nearby millisecond pulsars (see Verbiest et al. 2021, and references therein). To successfully detect gravitational waves, corrupting effects have to be mitigated down to unprecedented levels. Interstellar scintillation is one effect that has been shown to have the potential to corrupt pulsar-timing experiments at levels significant to PTA science (Lentati et al. 2017). Here we briefly interpret how scintillation impacts pulsar timing.

When pulsar signals travel through the inhomogeneous and highly turbulent IISM, another scintillation result, the scattering delay τ_st, can be calculated from $$2\pi v_{\text{d}}{\tau }_{\text{st}}=C_1$$(12)

where C₁ is a constant and close to unity (its value changes slightly for different geometries and models of the turbulence wavenumber spectrum, for a Kolmogorov spectrum, C₁ = 1.16) (Lambert & Rickett 1999). The τ_st provided two potential ways to affect the timing precision (Coles et al. 2010). Firstly, the pulse shape can be affected by pulse broadening: angular scattering broadens sharp single pulses into quasi-exponential pulses with a scattering delay τ_st, but how this affects the ToAs depends on the overall shape and complexity of the pulse profile. In this case, since fundamentally no distinction can be made between a time-constant amount of interstellar scattering and the intrinsic profile shape, the basic value for τ_st is of no direct relevance to timing experiments. Instead, variations of τ_st determine the impact on timing. Secondly, Hemberger & Stinebring (2008); Main et al. (2020) supposed that the shape of the pulse is effectively unchanged because the τ_st is much less than the width of the pulse profile and pointed out that the leading order effect from scintillation would be a time-shift of the pulse profile.

The scattering delays τ_st of 13 pulsars are calculated exploiting the measurements of v_d and the Eq. (12), we presented them in the last column of Table 1 by the median value and their 5/95 percentiles. All values range from a few nanoseconds to hundreds of nanoseconds, implying that for most pulsars in our sample the median amount of scattering is hard or impossible to determine at these frequencies. From just the lower and upper values of inferred scattering time, one infers that the largest differences are of order 100 ns. This implies that for PTA efforts, scattering variations may be irrelevant in most pulsars, but could be relevant in some.

In order to study potential improvements to the timing precision, however, it would be more meaningful to determine the time variable impact on timing based on the observed variations in scintillation; and while taking into account the actual shape of the pulse profile. This can readily be quantified though simulations, as will be presented in a follow-up paper.

7 Conclusions and further research

In this paper, we presented measurements of v_d and τ_d for 13 EPTA pulsars at multiple epochs by creating dynamic spectra and their ACFs. We plotted the time series of scintillation parameters and DM. We compared our scintillation parameters with previously published values and predictions from the NE2001 electron density model. The comparison showed that our results are typically consistent with those from literature. More specifically, for some pulsars we obtained more precise measurements in v_d because we have higher frequency resolution at WSRT and JBO; since most of our observations have a longer observing length than Levin et al. (2016) and Turner et al. (2021), we were able to determine τ_d for more pulsars. Furthermore, we have a longer time span of data and a larger number of observations for our pulsars, thus we were able to identify more interesting features in the time series of our scintillation parameters, as pointed out in Sect. 6. This will be discussed more fully in subsequent papers. In the comparison with NE2001, we found that the v_d we achieved in this work are mostly smaller than those predicted by NE2001. For PSRs J1857+0943 and J1939+2134, the multiband measurements of v_d allowed us to analyze the frequency scaling index.

Our results could help to improve the timing precision of PTA millisecond pulsars (Wang et al. 2022), by converting the scintillation bandwidth into scattering timescales τ_st and then subtracting these from TOAs. Furthermore, our work helps to determine the required observing length and frequency resolution for PTA observations to allow such schemes to work.

We will further investigate several useful and important applications of these data. Specifically, deeper analysis of this data set is being undertaken along the following lines:

1. Annual variations in the time series of scintillation timescales. It is possible to determine the distance and orientation of IISM clumps, but it is dependent on the physics of scattering and the lens models applied and hence these measurements can help to test the various models;

2. The impact of variable scattering on pulsar timing. In order to go beyond simplified theoretical models and fold in the complexity of the pulse profiles, a simulation of the impact on timing from the variable scattering delays (as derived from the scintillation bandwidth) is being undertaken;

3. Correlations between scintillation parameters and DM, mentioned in Sect. 6.4. Since the DM precision is not high enough to see trends in many pulsars of this work, we did not find many convincing correlations. But, a couple of these pulsars also have LOFAR DM time series (Donner et al. 2020), so a combined analysis could be meaningful.

Our analysis of these aspects with appear in subsequent papers in this series.

Acknowledgments

We thank Bill Coles for useful discussions and suggestions, Julian Donner for allowing us access to the code that is used to obtain the time series of DM, Eleni Graikou for her contribution to the pulsar observing with EFF. Part of this work is based on observations with the 100-m radio telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a consolidated grant from the STFC in the UK. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from “Programme National de Cosmologie et Galaxies” (PNCG) of CNRS/INSU, France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). JPWV acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the Heisenberg programme (Project No. 433075039). YL acknowledges support from the China Scholarship Council (No. 201808510133). J. W. McKee is a CITA Postdoctoral Fellow: This work was supported by Ontario Research Fund-research Excellence Program (ORF-RE) and the Natural Sciences and Engineering Research Council of Canada (NSERC) [funding reference CRD 523638-18].

References

All Tables

All Figures

Fig. 1 Dynamic spectra I(f, t) from four observations with different scintle sizes in the left panels and their normalized autocorrelation functions in the right panels. The horizontal and vertical axes correspond to the time (or time lag) and frequency (or frequency lag), respectively. In the two smaller side plots, the sky blue points are the 1D autocorrelation functions at zero time and frequency lag, the black lines are the best Gaussian fits, the red lines indicate the values for vd and τd. The pulsar name, telescope, epoch of observation, vd and τd are given above the dynamic spectra.

In the text

Fig. 2 Long-term variations of scintillation bandwidth vd, scintillation timescale τd and DM for 13 EPTA pulsars. To ensure that the variations shown by the most precise DM values remain visible, for most pulsars, we only show the most precise DM measurements, which are usually derived from 21-cm Nançay data. The different colors indicate measurements at different telescopes, following the legend given in the top-right corner of each sub plot. Following the telescope name, we present the center frequency in units of MHz. Downward-pointing and upward-pointing triangles indicate unreliable estimates (i.e., upper and lower limits) caused by limited frequency resolution and limited observation length, respectively. The vertical gray dotted lines indicate the start of a calendar year.

In the text

Fig. 3 The left panel compares the vd from previous works to our measurements. Here, “L” indicates Levin et al. (2016), “T” refers to Turner et al. (2021), “11 cm” compares against measurements we made in the 2-GHz range; “21 cm” refers to measurements we made in the 1-GHz range – in both cases the previously published measurements were rescaled to our observing frequency, as given in Col. (3) of Table 2. Data from Levin et al. (2016) were at 1500 MHz while Turner et al. (2021) used data at 1500 MHz as well as 820 MHz, but we only compare to the former. The three ellipses indicate groups of points that show significant offsets, which are discussed in more detail in the text. The right-hand panel shows violin plot for the distribution of the differences (Delta) of the four groups between the vd of our measurements and the vd from previous studies.

In the text

	Fig. 4 The NE2001 predictions of scintillation parameters vs. our measurements. The left panel shows the comparison of vd, the right panel shows the comparison of τd.
In the text

	Fig. 5 Time series of scintillation bandwidth vd, scintillation timescale τd and DM for PSR J1939+2134. There are two events, that show decreases in vd and τd, and increases in DM, which are indicated by the red lines. The first one is between MJDs 55807 and 55845, the second one is between MJDs 56000 and 56090.
In the text

	Fig. 6 The time series of frequency scaling indexes of the vd for two pulsars. The mean and standard deviation of scaling indexes are printed in the right-bottom corner.
In the text