© The Authors 2026

1 Introduction

Scintillation is an interference pattern manifested over time and frequency. It is a common phenomenon in pulsars whose radiation in the radio regime is often scattered into multiple paths by the free electrons in the ionized interstellar medium (IISM). These paths are spatially coherent if the source is very compact, and hence interfere with each other. This understanding of pulsar scintillation has been accepted since the work by Scheuer (1968). Nevertheless, the physical nature of the structures causing scintillation remains largely unknown, with competing geometries proposed thus far, such as turbulence (Rickett 1977; Cordes et al. 1985), magnetotails of clouds (Walker 2007), reconnection sheets (Pen & Levin 2014), noodle-like filaments (Gwinn 2019; Gwinn & Sosenko 2019), or discrete scatterers (Kramer et al. 2024). Moreover, the connection to larger-scale associates that are observable with other methods than scintillation is unclear. Practically all known objects in the IISM have been hypothesized, whereby evidence has been collected for H II regions (e.g., Mall et al. 2022), supernova remnants (e.g., Rankin & Counselman 1973), and bow shocks in pulsar wind nebulae (e.g., Reardon et al. 2025). In addition, associations with H I filaments have been found by Stock & van Kerkwijk (2025). However, the role of other structures such as the local bubble (e.g., Ocker et al. 2024) is still unclear and many sources of scattering have no known associates.

Scintillometry is the science of using scintillation to infer properties of astrophysical objects. It is very sensitive to tiny structures because scintillation preserves some phase information of the incoming radiation, enabling an approach akin to interferometry. Scintillometry has been transformed with the introduction of the Fourier transform of the dynamic spectrum as a standard tool, which was named the secondary spectrum by Rickett et al. (1997). Secondary spectra have revealed the scintillation arc phenomenon as a new scintillation observable, which was first noted by Stinebring et al. (2001) to be ubiquitous in pulsar scintillation. Their parabolic structure was derived by Walker et al. (2004) and Cordes et al. (2006) as a result of scattering by a thin plasma screen. Hill et al. (2003, 2005) first described the remarkable stability of these arcs and their substructure over frequency and time, which can be directly mapped onto the distribution of propagation paths or images. The distribution of these images often seems to be very anisotropic and can even form a straight line of a large number of separated images. This was observed by Brisken et al. (2010), who first combined scintillometry and very long baseline interferometry (VLBI).

While most scintillometry studies focused on single-dish observations, methods to combine it with VLBI to get superior and instantaneous astrometry of scattering screens have been improved recently: Simard et al. (2019) developed a method to make use of single-dish intensities in addition to visibilities and Baker et al. (2023) combined VLBI with the θ-θ transform. Both methods were demonstrated on the observations of PSR B0834+06 by Brisken et al. (2010). Applications to other observations are at the point of writing limited to PSR B1508+55 by Marthi et al. (2021) and PSR B1133+16 by Stock et al. (2025).

The pulsar B1508+55 is a bright and isolated pulsar discovered by Huguenin et al. (1968) with a very high proper motion and is located at a distance of ${2.10}_{-0.14}^{+0.13}$ kpc, as measured by Chatterjee et al. (2009). Prior to late 2020, its scintillation arcs did not show arclets, but flat stripes (Stinebring 2007; Marthi et al. 2021; Stinebring et al. 2022). Its current scintillation features prominent inverted arclets, which makes it a promising target for another multitelescope study. The transition to the new state happened during observations by Sprenger et al. (2022), who interpreted the stripes as an effect of a second screen. The first screen lies at a distance within the first 10% from the Earth to the pulsar while the second one lies in the last 10%. The presence of the farther screen can still be observed in the new state of scintillation, even though the second arc ended up disappearing. The characteristic stripe-like blurring of secondary spectra as well as the temporal modulation of the amplitude of individual features along the scintillation arc remained visible in 2023 and 2024 at the time of the new observations presented here. However, the transition seems to arise solely from changes in the screen closer to Earth. It is correlated in time with the line-of-sight crossing some bright scattered images of the pulsar as observed at much lower frequencies by Wucknitz (2018). Using the annual variation of scintillation arcs over several years, the movement of tracked features along the arcs, and the consistent modulation by the second screen, a model of the location, velocity, and orientation of anisotropy of the scattering screens was constructed by Sprenger et al. (2022).

The aim of this work is to test whether two telescopes observing scintillation simultaneously can measure scattering screens with comparable or even superior sensitivity compared to long-time scintillation arc monitoring, without having to rely on the measurement of visibilities and additional telescopes. This is not only important as a proof of concept, but it would also deliver instantaneous 2D maps of the screen that cannot be obtained from long-time monitoring with a single dish.

This paper is organized as follows. The observations are described in Section 2 and the theoretical background is given in Section 3. The cross-spectra and images of the scattering screen are reported in Section 4, while the analysis methods used to constrain scattering screen models are presented in Section 5. The implications are discussed in Section 6.

2 Observations and data processing

Observations of PSR B1508+55 were made simultaneously with the Five-hundred-meter Aperture Spherical Telescope (FAST) and the Effelsberg 100-m telescope. To cover a range of different velocity vectors of Earth and to observe a possible screen evolution, three observations were spread throughout the course of one year. One observation failed due to technical problems at Effelsberg. The two successful observations were performed between December 4 and 5 in 2023 and on May 20 in 2024. The pulsar was observed for the full 4.5 hours of its elevation being high enough to be observed at FAST. The modified Julian dates (MJD) at the time of highest elevation at FAST were 60 283 and 60 450. The baselines covered are shown in Fig. 1. The positions of the telescopes and the declination of the pulsar allow for a wide range of baseline orientations during a single observation.

To maintain flexibility in the channelization, baseband data were recorded at Effelsberg using the Effelsberg Direct Digitisation (EDD) system (Barr et al. 2023). To cross-correlate the data, we aimed for both datasets to have the same bin widths in terms of time and frequency. Furthermore, it is beneficial to have the same edge coordinates of these bins. These goals are challenging to achieve because of the different bandwidth of the receivers as well as the properties of DSPSR (van Straten & Bailes 2011) which was used for folding and channelization, introducing small variations in the start time and temporal bin width. While folding at the period of the pulsar as given by the ephemeris obtained from the ATNF pulsar catalogue (Manchester et al. 2005), the data were divided into subintegrations of the same manually fixed duration at both telescopes. For this number, approximately seven rotations of the pulsar were chosen, whose period is ~0.74 s. This number is a balance between the temporal resolution and signal-to-noise ratio (S/N). The observed band at FAST with the 19-beam receiver ranged over 500 MHz from 1000 to 1500 MHz and was divided into 8192 channels. At Effelsberg, the P217mm receiver was used, which ranges over 400 MHz from 1200 to 1600 MHz, although not the whole band is usable. The different bandwidth makes it impossible to match the FAST channels with an integer number of channels of same width. Thus, the data were divided into 32 768 channels, which were five times finer than the FAST channels. The time bins were aligned as well as possible with the Lepoch option of DSPSR that lets the first time bin start at a given time corresponding to the first time bin of the FAST data.

The functions dedisperse(), remove_baseline(), and convert_state(‘Intensity’) of the Python wrapper of PSRCHIVE (van Straten et al. 2012) were used to remove the baseline, convert the data to Stokes I, and to dedisperse them. At this point, the data consisted of data cubes over time, frequency, and rotational phase separated into 64 phase bins.

Radio frequency interference (RFI) is an important limitation for this method because of its demand for high S/N within short subintegrations and small channels at both telescopes. Using a custom implementation in Python, strong known RFI bands at Effelsberg were manually masked while the remaining data were cleaned by investigating statistics in the manually identified off-pulse region of each subintegration and channel. The investigated statistics were median, mean, modulation index (standard deviation divided by mean), kurtosis, skewness, range (difference between maximum and minimum), maximal Fourier coefficient divided by mean, and Pearson cross-correlation coefficient between neighboring phase bins. A function, P, predicting the likelihood of a sample not being RFI based on its vector, x, of off-pulse statistics was calibrated on manually labeled samples. We chose a multivariate Gaussian for this function and fit for the mean and standard deviation of each statistic by minimizing $$⟨P(\mathit{x}){⟩}_{\mathrm{R}\mathrm{F}\mathrm{I}}+⟨1-P(\mathit{x}){⟩}_{\mathrm{n}\mathrm{o}\mathrm{t}\mathrm{R}\mathrm{F}\mathrm{I}}.$$(1)

This fit converged after collecting ~100 of randomly chosen samples of each category, which can be done in short time by using intermediate fits to search for candidates of the category of lower sample number. The statistic with the most prediction power for RFI turned out to be the range followed with some distance by the kurtosis. Channels where more than 75% of subintegrations had a likelihood of less than 0.5 to not be RFI were masked for the whole observation. All completely masked data were treated as if their P was zero.

Next, template pulse profiles were formed by averaging over all subintegrations, subtracting the median of the off-pulse region for each channel, and taking the mean over all channels, in this order. Masked values were ignored. Finally, all off-pulse phase bins of the obtained profile were set to zero, along with any remaining values smaller than zero at the edges of the pulse.

The dynamic spectrum was then formed from data that had not been masked in advance. Also, the median of the off-pulse region was subtracted for each subintegration and channel. Afterward, the data were multiplied by the template profile and then summed over all phase bins. The described method was used to dynamically remove background and noise level, while measuring the flux of the pulsar while taking into account the whole pulse, with phases of high S/N being given more weight.

All dynamic spectra were now converted to the same frequency channels that were chosen to be equal to the channels of the FAST data, but extended to the band covered by Effelsberg. This was done by binning the channels of the Effelsberg data into these common frequency grids, where the average was weighted with P². The remaining RFI (often caused by contamination solely in the on-pulse region) was suppressed by decreasing the weights of data values larger than seven times the weighted standard deviation of all data and masking all data values greater than ten times above it. These thresholds were determined from past observations of this pulsar in Sprenger et al. (2022) that showed that the real signal almost never exceeds them. Even so, such contamination could not be completely removed. In addition, the square root of the mean weight of each bin was used to obtain the P value for the whole bin. In the converted data, all bins with P < 0.5 were masked. In summary, almost 50% of data needed to be masked in both Effelsberg observations and around 10% of the data in both FAST observations. Hence, only data from 1300 MHz to 1425 MHz were used. In this subband, both telescopes were sensitive enough to observe with less than 3% of the data masked as RFI.

Finally, the dynamic spectra were normalized to remove intrinsic pulse-to-pulse variations and remaining frequency-dependent variations caused by the bandpass applied at the telescopes. This was done iteratively. First, each time bin was divided by its mean over frequency, not counting the masked data. Here, time bins whose mean was below 10% of the average of the whole dynamic spectrum were also masked because these would contaminate the subsequent analysis with their low S/N. Then, the same normalization was done for every channel by dividing by the mean over all time bins and it was repeated for every time bin to ensure that the mean of each time bin is exactly unity. The resulting dynamic spectra are shown in Fig. 2.

Using ASTROPY (Astropy Collaboration 2022), the telescope positions corresponding to the observation times are obtained in the Geocentric Celestial Reference System (GCRS), while the Earth’s velocity is obtained in the International Celestial Reference System (ICRS). These were then projected to 2D vectors in right ascension (RA) and declination (Dec) using the position of the pulsar.

Fig. 1 Telescope positions as viewed from the pulsar throughout the observations. Color represents the time since the start of the observation. Position of the land masses at the beginning (blue) and end (red) of the observations and that of the baseline between the telescopes are also shown.

3 Scintillation arcs from multiple telescopes

PSR B1508+55’s scintillation contains dominant 1D scattering and has shown strong signs of two scattering screens as reported in Marthi et al. (2021) and Sprenger et al. (2022). In 2020, it underwent a phenomenological change, which also affected the annual variation of its scintillation arcs that is indicative of a change in orientation of the axis of anisotropy of the scattering medium (Sprenger et al. 2022) within the closer screen. Thus, we did not treat this angle as a constant in our modeling.

We applied both a one-screen and a two-screen model here. Although past studies have strongly implied a second screen, this screen’s properties are much less constrained; whereas a one-screen model strengthens the comparability to other studies that apply similar models. Furthermore, the second screen is likely to be much farther away and hence smaller in angular size. Thus, a one-screen model is expected to be able to describe any signal arising from the difference between two observer locations on Earth to a good approximation.

All symbols used for geometric variables are defined in Fig. 3. Scintillation is insensitive to small variations in positions along the line of sight. Thus, all positions and velocities can be treated as 2D vectors. The coordinates are defined along the positive directions of RA and Dec. The angles α_x and α_y of the orientation of straight lines of scattered images are defined such that they are zero in the direction of positive RA and grow toward the direction of positive Dec. This notation is the same as used in Sprenger et al. (2022) and can be converted to angles measured east of north via α_{East of North} = 90° − α.

3.1 One-screen theory

The intrinsic temporal and spectral variation of the pulsar is removed by the normalization. The remaining variations are dominated by interference of different paths of propagation whose phase differs due to the different path lengths. This geometric phase difference Φ is simply given by the number of wavelengths, λ = c/ν, that fit into the length difference, ΔL, of the scattered path of propagation and the pulsar’s distance along the direct line of sight, such that Φ = 2πΔL/λ. In the small-angle approximation, this length depends quadratically on the scattering angle. Since at least the work by Brisken et al. (2010), it has been customary to express ΔL = D_effθ²/2, using the effective distance $$D_{\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{d_{\mathrm{s}\mathrm{c}\mathrm{r}}{d}_{\mathrm{p}\mathrm{s}\mathrm{r}}}{d_{\mathrm{p}\mathrm{s}\mathrm{r}}-d_{\mathrm{s}\mathrm{c}\mathrm{r}}},$$(2)

where d_scr is the distance from the observer to the screen and d_psr is the distance from the observer to the pulsar. Temporal evolution happens because of changes of the observed scattering angle, θ, due to the relative velocities of telescopes, the ISM, and the pulsar. An effective velocity term can be defined as $${\mathit{V}}_{\mathrm{e}\mathrm{f}\mathrm{f}}={\mathit{V}}_{\oplus }-\frac{d_{\mathrm{p}\mathrm{s}\mathrm{r}}}{d_{\mathrm{p}\mathrm{s}\mathrm{r}}-d_{\mathrm{s}\mathrm{c}\mathrm{r}}}{\mathit{V}}_{\mathrm{s}\mathrm{c}\mathrm{r}}+\frac{d_{\mathrm{s}\mathrm{c}\mathrm{r}}}{d_{\mathrm{p}\mathrm{s}\mathrm{r}}-d_{\mathrm{s}\mathrm{c}\mathrm{r}}}{\mathit{V}}_{\mathrm{p}\mathrm{s}\mathrm{r}},$$(3)

such that the angular velocity of locations fixed on the screen is dθ/dt = −V_eff/D_eff. This definition of the effective velocity is not unique in the current literature. It is a convenient definition in the case of multiple telescopes because it is equal to the velocity of an observer with respect to a constant point in the scintillation pattern. Hence, the offset, p, of a telescope from the center of the Earth can simply be added to the cumulative shift by this velocity to obtain the phase of a scattered ray exactly at the position, p(t), of a telescope in the plane perpendicular to the line of sight at a point, t, in time, $$\mathrm{\Phi }=\frac{\pi \nu }{c}{D}_{\mathrm{e}\mathrm{f}\mathrm{f}}{(\mathit{\theta }-\frac{\mathit{p}(t)+{\mathit{V}}_{\mathrm{e}\mathrm{f}\mathrm{f}}t}{D_{\mathrm{e}\mathrm{f}\mathrm{f}}})}^2,$$(4)

where we have chosen to define V_⊕ as the velocity of the center of the Earth while any movement of the observer on its surface relative to the center is encoded in p(t). All angles are defined with respect to the line of sight to the pulsar at t = 0 which is the beginning of the observation. In the case of a 1D screen, this expression can be shortened because only contributions that depend on the angular position θ along one axis of the screen vary between rays. Removing all constant contributions and keeping only the components of vectors that are parallel to the straight line of images, we obtain $$\mathrm{\Phi }=\frac{2\pi \nu }{c}[\frac{D_{\mathrm{e}\mathrm{f}\mathrm{f}}}{2}{\theta }^2-(p_{\parallel }(t)+V_{\mathrm{e}\mathrm{f}\mathrm{f},\parallel }t)\theta ].$$(5)

The baseline parallel to the screen between two telescopes A and B can be defined as $$\mathit{b}(t_A,t_B)={\mathit{p}}_A\left(t_A\right)-{\mathit{p}}_B\left(t_B\right).$$(6)

It can be easily read off now that telescope A will observe exactly the same situation as telescope B after a time, Δt, that is given by $${\mathrm{\Delta }}_t(t_A,t_B)=-b_{\parallel }(t_A,t_B)/V_{\mathrm{e}\mathrm{f}\mathrm{f},\parallel }.$$(7)

The times, t_A and t_B, at which the data were taken do not need to be equal. The parallel baseline is given by $$b_{\parallel }(t_A,t_B)=\mathit{b}(t_A,t_B)\cdot (\genfrac{}{}{0ex}{}{\cos {\alpha }_x}{\sin {\alpha }_x}).$$(8)

Pairs of individual scattered paths can be made visible by forming secondary spectra, which are the square modulus of the 2D Fourier transform of the intensity dynamic spectrum, where the Doppler rate, f_D, and delay, τ, are the Fourier conjugates of time and frequency, respectively. Walker et al. (2004) and Cordes et al. (2006) first derived that secondary spectra represent a map of the scattering screen in differential coordinates, where $$f_{\mathrm{D}}({\mathit{\theta }}_1,{\mathit{\theta }}_2,t,\nu )=\frac{\mathrm{d}}{\mathrm{d}t}\frac{\mathrm{\Phi }({\mathit{\theta }}_1,t,\nu )-\mathrm{\Phi }({\mathit{\theta }}_2,t,\nu )}{2\pi },$$(9) $$\tau ({\mathit{\theta }}_1,{\mathit{\theta }}_2,t,\nu )=\frac{\mathrm{d}}{\mathrm{d}\nu }\frac{\mathrm{\Phi }({\mathit{\theta }}_1,t,\nu )-\mathrm{\Phi }({\mathit{\theta }}_2,t,\nu )}{2\pi },$$(10)

for all pairs (θ₁, θ₂) of scattering angles. Even for fixed image positions θ₁ and θ₂, the corresponding delay and Doppler rate evolve with time and frequency. Since secondary spectra cannot be formed from a single data point, this results in a smearing effect for which different mitigation techniques have been proposed by Fallows et al. (2014) and Sprenger et al. (2021). The secondary spectrum coordinates of the 1D one-screen case are $$f_{\mathrm{D}}=-\frac{\nu }{c}{V}_{\mathrm{e}\mathrm{f}\mathrm{f},\parallel }({\theta }_1-{\theta }_2)-\frac{\nu }{c}\frac{\mathrm{d}}{\mathrm{d}t}{p}_{\parallel }(t)({\theta }_1-{\theta }_2),$$(11) $$\tau =\frac{D_{\mathrm{e}\mathrm{f}\mathrm{f}}}{2c}({\theta }_1^2-{\theta }_2^2)-\frac{p_{\parallel }(t)+V_{\mathrm{e}\mathrm{f}\mathrm{f},\parallel }t}{c}({\theta }_1-{\theta }_2).$$(12)

Usually, the second terms in Eqs. (11) and (12) are neglected. They cause a periodic variation of the scintillation arc over the course of a day. We confirmed that our scintillation arc measurements are not sensitive enough to identify such a variation and, hence, we chose to neglect these terms in the following. A double-parabolic appearance of scintillation arcs and inverted arclets that both follow $\tau =\eta f_{\mathrm{D}}^2$ can be read off from these relations. Thus, secondary spectra containing such features are strong evidence for 1D scattering, as is the case for PSR B1508+55. Adopting the approach by Sprenger et al. (2022), we reformulated the arc curvature, η, in terms of the effective drift rate, ζ, to avoid frequency dependencies and divergences, $$\zeta =\frac{1}{2\nu \sqrt{\eta }}=\frac{\left|V_{\mathrm{e}\mathrm{f}\mathrm{f},\parallel }\right|}{\sqrt{2c{D}_{\mathrm{e}\mathrm{f}\mathrm{f}}}}.$$(13)

Simard et al. (2019) developed generalizations for the secondary spectrum in the case of multistation observations. The intensity cross-spectrum, S_I, can be formed without having to measure the visibilities and is defined as $$S_I(f_{\mathrm{D}},\tau )={\tilde{I}}_A(f_{\mathrm{D}},\tau )\times {\tilde{I}}_B(-f_{\mathrm{D}},-\tau ),$$(14)

where $\tilde{I}$ is the Fourier transform of the dynamic spectrum, I. The complex phase of S_I is given by $$\mathrm{\Phi }(f_{\mathrm{D}},\tau )=-\frac{2\pi \nu }{c}\mathit{b}\cdot ({\mathit{\theta }}_1-{\mathit{\theta }}_2)\stackrel{1\mathrm{D}}{=}-\frac{2\pi \nu }{c}{b}_{\parallel }({\theta }_1-{\theta }_2).$$(15)

The intensity cross-spectrum is the Fourier transform of the cross-correlation function, just as the secondary spectrum is the Fourier transform of the autocorrelation function. For 1D screens, the cross-correlation function is equal to the autocorrelation function shifted by Δ_t, defined in Eq. (7), $$\mathrm{CCF}(\mathrm{\Delta }t)=\mathrm{ACF}(\mathrm{\Delta }t-{\mathrm{\Delta }}_t)$$(16)

where Δt is the time lag for which the correlation is computed. This is a result of the scintillation being identical at the two telescopes. This situation arises necessarily for a 1D screen but can only be observed in 2D screens for the unlikely case of the effective velocity and the baseline being exactly aligned.

The temporal ACF of scintillation can be modeled as a Gaussian and is defined by the scintillation time, t_s, $$\mathrm{ACF}(\mathrm{\Delta }t)=m^2\exp (-\frac{\mathrm{\Delta }{t}^2}{2t_s^2}).$$(17)

Here, m is the modulation index which is equal to the standard deviation of the data divided by its mean. An exponential including the exponent 5/3 instead of 2 has been proposed as an alternative to a Gaussian to fit the temporal ACF (e.g., Reardon et al. 2019). However, testing a different function is not necessary here because only the peak of the model function is relevant for a measurement of the shift, Δ_t.

Fig. 2 Dynamic spectra of the two simultaneous observations. The data have been normalized to a mean of unity in every time bin and the band has been cut to the optimal subband. The time axis is synchronous and measured from the beginning of simultaneous observations.

Fig. 3 Definitions relevant to this study. The directions of the arrows are examples and do not define the positive direction of signed quantities. Alternative notations in the one-screen case are given in parentheses.

3.2 Two-screen theory

An analytical solution for the wave phase in a model of two 1D screens, where each ray is scattered by both screens has been developed in Sprenger et al. (2022) and written down in Eqs. A8 and A15 therein. Instead of setting p = 0, we are able to obtain effective positions, p_eff, that are different for each screen, $$\begin{array}{rl}\mathrm{\Phi }=& \frac{\pi \nu }{c}(D_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{X}}{\theta }_x^2+D_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{y}}{\theta }_y^2-2D_{\mathrm{m}\mathrm{i}\mathrm{x}}{\theta }_x{\theta }_y\\ & -2(p_{\mathrm{e}\mathrm{f}\mathrm{f},x}+V_{\mathrm{e}\mathrm{f}\mathrm{f},x,\parallel }t){\theta }_x-2(p_{\mathrm{e}\mathrm{f}\mathrm{f},y}+V_{\mathrm{e}\mathrm{f}\mathrm{f},y,\parallel }t){\theta }_y).\end{array}$$(18)

The distance and velocity terms are given in Sprenger et al. (2022). The effective positions result from the different projections of the telescope position on the two 1D screens and can be obtained with the same logic. They are given by $$p_{\mathrm{e}\mathrm{f}\mathrm{f},x}=p_{\parallel }+\frac{\gamma \delta D_x{D}_{y,s}{p}_{\perp }}{D_y{D}_{x,s}-D_x{D}_{y,s}{\delta }^2},$$(19) $$p_{\mathrm{e}\mathrm{f}\mathrm{f},y}=-\frac{\gamma D_y{D}_{x,s}{p}_{\perp }}{D_y{D}_{x,s}-D_x{D}_{y,s}{\delta }^2},$$(20)

where γ = sin(α_x − α_y) and δ = cos (α_x − α_y).

Each screen produces a different time shift that is obtained by dividing the effective baseline by the effective velocity of that screen, analogously to Eq. (7). These time shifts could be observed separately and are not added up, such that the dynamic spectrum at the two telescopes becomes more different than just being shifted. In particular, the temporal shift between two telescopes A and B due to the first screen is given by $${\mathrm{\Delta }}_t(t_A,t_B)=-\frac{p_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{x},\mathrm{A}}\left(t_A\right)-p_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{x},\mathrm{B}}\left(t_B\right)}{V_{\mathrm{e}\mathrm{f}\mathrm{f},x,\parallel }}$$(21)

4 Imaging the scattered pulsar

The secondary spectrum is closely connected to the scattered image of the pulsar. If it is sufficiently sparse, each point in it represents a single pair of subimages. Its power is their multiplied intensity and its location is determined by their relative phase. As such, the secondary spectrum is a holograph (Walker et al. 2008), preserving phase information within the intensity data. The idea behind the work by Brisken et al. (2010) was to combine this phase information with the phase information obtained via VLBI. Brisken et al. (2010) and Simard et al. (2019) used the fact that at the apex of each arclet one of the two interfering images lies at the origin to create images of the scattering screen of PSR B0834+06. In summary, the location along the scintillation arc separates subimages from each other whose VLBI phase can then be used to infer their angular positions.

Without the VLBI information, the pulsar can be imaged from the secondary spectrum information alone. One example for doing so is given in Sprenger et al. (2021). However, very sharp arclets are required. The orientation of the obtained image on the sky is not constrained and the locations are ambiguous with two solutions on each side of the effective velocity vector. Since the arclets of PSR B1508+55 are blurrier than those of PSR B0834+06 and contaminated by the stripe effect, such a single-dish imaging method was not performed.

4.1 Intensity cross-spectra

To form the intensity cross-spectra, both dynamic spectra need to be located on an identical and regular grid in time and frequency. The procedure to get identical frequency channels and approximately identical time bins was described in Section 2. In addition, small shifts between the coordinate grids can be corrected by considering the Fourier transform, ${\tilde{I}}^{\mathrm{\prime }}$, of a dynamic spectrum whose coordinates are shifted by δt in time and by δν in frequency, $$\begin{array}{rl}{\tilde{I}}^{\mathrm{\prime }}(f_{\mathrm{D}},\tau )& =\int I(t+\delta t,\nu +\delta \nu ){\mathrm{e}}^{-2\pi i(f_{\mathrm{D}}t+\tau \nu )}\mathrm{d}t\mathrm{d}\nu \\ & ={\mathrm{e}}^{2\pi i(f_{\mathrm{D}}\delta t+\tau \delta \nu )}\tilde{I}(f_{\mathrm{D}},\tau ).\end{array}$$(22)

This factor was divided out before the computation of any cross-spectrum. Evolving shifts between time bins due to corrections of the time of arrival of the pulses could not be corrected with this method.

In contrast to the case discussed in Brisken et al. (2010), the baseline evolved significantly throughout the observation and the phase of the cross-spectrum as given in Eq. (15) is expected to vary strongly as a result. Thus, meaningful cross-spectra can only be obtained from subsets of the data. In Figs. 4 and 5, cross-spectra are shown that were computed from subsets each consisting of a fourth of the observation time. Clear linear phase evolutions are visible in these spectra in addition to the double parabolic shape of the scintillation arc.

A phase evolution that is linear in Doppler rate and independent of delay corresponds to a simple temporal shift of the dynamic spectra at the two telescopes with respect to each other. Hence, the corresponding cross-correlation functions that are also reported in Figs. 4 and 5 show the correlation to be maximal at an offset from the origin in time lag.

4.2 Producing the image

The apexes of inverted arclets lie in the vicinity of the main arc. These locations result from pairs of angles where one angle is the origin. According to Eq. (15), the phase of the intensity cross-spectrum at these points becomes $${\mathrm{\Phi }}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{x}}(f_{\mathrm{D}}(\mathit{\theta }),\tau (\mathit{\theta }))=-\frac{2\pi \nu }{c}\mathit{b}\cdot \mathit{\theta }.$$(23)

Even though we only have a single pair of telescopes, we can leverage the evolution of this baseline and divide the observation into subsets to obtain multiple baselines. Thus, Eq. (23) can be solved for θ to form an image. We developed an imaging method closely related to the one described by Simard et al. (2019).

First, the dynamic spectra were divided into ten adjacent subsets of equal length in time. To obtain smaller bins in f_D, they were each padded in time with zeros by twice the subset length at each side. The obtained cross-spectra where then smoothed by a top-hat function of 0.1 μs in τ and 0.5 mHz or 0.25 mHz in f_D where the former was used for the observation in 2023 and the latter for the one in 2024. As a result of the smoothing, the quality of the phases is sufficient to be unambiguously corrected for wrapping at ±π.

Next, points are selected which lie close enough to the main arc such that Eq. (23) can be applied. The scintillation arc’s coordinates are known from measurements of the curvature that are described in Section 5.1 below. We estimated by trial and error that the same f_D windows as those used for the smoothing define a region were our approximations are valid. For each point of positive delay in this region, Eq. (23) can be solved by a least-squares optimization for θ, inserting the baselines of the central time of each subset. The solutions were binned into a 300 × 300 grid spanning over two astronomical units (AU) according to the screen distance at the two epochs taken from the second column of Table 3 below. Each point was weighted by the square root of the mean cross-spectrum amplitude multiplied by |f_D|, which corresponds to its amplitude in angular space according to Sprenger et al. (2021). The amplitude at the origin cannot be measured.

The resulting images are shown in Figs. 6 and 7. This method contains several parameters that were optimized manually and affect the achieved resolution. Still, the dominant structures within these images are stable with respect to some variation of these parameters. These contain the expected straight line but also some substructure deviating from it, as well as brighter regions that relate to individual arclets. The method described here assumes that all images lie close enough to a straight line for their corresponding arclets to be close to the main arc. Offset arclets would need to be imaged separately. We also neglected contributions of the second screen to the phase in this method. The result confirms that they are indeed negligible.

After the solutions for each point are found, the amplitudes in the weighting can be replaced. This can be used to obtain the temporal evolution of the image by using secondary spectra created from a moving window of the dynamic spectrum. It confirms the modulation of the second screen moving along the image that was described in Sprenger et al. (2022).

Fig. 4 Cross-secondary spectra (top row) and cross-correlations (bottom row) from December 4–5, 2023. The observation has been split into four parts of equal length in time, shown in chronological order. The colors denote the complex phase of the cross-secondary spectra, while their brightness denotes the power. The cross-correlation functions are obtained as the Fourier transform of the cross-secondary spectra. They illustrate the temporal shift corresponding to the linear phase evolution in the cross-spectra.

Fig. 5 Cross-secondary spectra (top row) and cross-correlations (bottom row) from 20 May 2024. The setup is identical to Fig. 4.

Fig. 6 Image of the scattering screen at December 4–5, 2023. The central point cannot be imaged and is instead shown as fully saturated to highlight the position of the pulsar. The color scale is linear and normalized to the highest value. A resolution of about 0.1 mas is achieved. The evolution of the image over the course of the observation can be found as an online movie.

Fig. 7 Image of the scattering screen at 21 May 2024. See Fig. 6 for more details. An online movie is available.

5 Measurement of screen parameters

The cross-spectrum analysis revealed that the most prominent difference between the scintillation at the two telescopes is an evolving temporal shift. It can also be seen by eye in Fig. 2. This shift has been theoretically predicted in Eq. (7) and Eq. (21) for a one-screen and a two-screen model, respectively. In a one-screen model, measuring this shift between the telescopes while the baseline rotates with the Earth provides a direct measure of the screen orientation and the parallel effective velocity. At the same time, the arc curvature in Eq. (13) is given by the parallel effective velocity and the effective distance. Thus, measuring the arc curvature in combination with the varying interferometric shift unambiguously constrains the parameters α_x, D_eff, and V_eff,∥. We used these measurements to constrain one-screen and two-screen models.

5.1 Arc curvature measurement

Despite being the defining parameter of a scintillation arc, the arc curvature can be difficult to obtain precisely. Even under ideal observational conditions, the signal itself may significantly depart from the ideal double-parabolic structure because of deviations from a 1D alignment of images. The images might even evolve in time and frequency, be further scattered by additional screens, or follow a distribution not centered at the origin due to large-scale refraction.

Multiple methods have been proposed to improve or automate the measurement of arc curvatures. Bhat et al. (2016) and Reardon et al. (2020) introduced methods based on the Hough transform where the power along different trial arcs is summed and the best curvature is obtained by maximizing this sum. This method works better the more concentrated arcs are. For the opposite case of an arc that is rather double-parabolic rather than a single parabola, Baker et al. (2022) developed a method that explicitly uses the 1D nature of the screen corresponding to this case. This is the method adopted here and we refer to it as the θ-θ method. Its idea is to use the θ-θ transform introduced in Sprenger et al. (2021) that transforms arclets into straight lines for different trial curvatures. For the correct curvature, these lines are horizontal and vertical. This condition is equivalent to maximizing the highest eigenvalue of an eigenvector decomposition of the θ-θ diagram.

As discussed in Sprenger et al. (2022), θ-θ method is sometimes contaminated by additional peaks that are likely non-physical. There, a manual method of measurement has been used in addition to cross-check the measured values. As can be seen in Fig. 8, the θ-θ method was also affected here but their effect never shifts the best fit to or over the edges of the manually identified possible range of arc curvatures. Thus, the same methodology as in Sprenger et al. (2022) was used without the combination with manual measurements; namely, the dynamic spectrum was divided into subsets of size Δt in time and Δν in frequency that were set to $$\begin{array}{rl}\mathrm{\Delta }t& =60\mathrm{s}\times \frac{0.006\sqrt{\mu \mathrm{s}}/\mathrm{h}}{{\zeta }_{\text{est.}}},\\ \mathrm{\Delta }\nu & =2.5\mathrm{M}\mathrm{H}\mathrm{z}.\end{array}$$(24)

The value for ζ_est. used here was estimated by eye using the graphical method described in Sprenger et al. (2022). For each of these subsets, θ-θ transformations were performed for 600 trial values of ζ that ranged from 0.00036 to $0.01008\sqrt{\mu \mathrm{s}}/\mathrm{h}$. For each of these, the highest eigenvalue of an eigenvector decomposition was saved. These eigenvalues were then summed over all frequencies to reduce scatter and a parabolic fit was performed on a region of width $0.0006\sqrt{\mu \mathrm{s}}/\mathrm{h}$ centered around the maximum value. The mean and standard deviation over all fit results was then used as the measured value and its uncertainty, respectively. The measured values are shown in Table 1. For subsequent parameter inference, the combined value and standard deviation was computed. The measurements are in agreement, but still greater than the values predicted from the best-fit model by Sprenger et al. (2022).

Fig. 8 Arc curvature measurement. Bottom row: secondary spectra with equal color scale gauged by the 99.95 percentile. Mid row: sum of eigenvalues for all subsets in frequency as a function of subsets in time and trial ζ. To row: corresponding mean over all subsets in time. The best-fit value for ζ or the scintillation arc is shown in red. In the top and bottom plots its 1σ region is shown in red and the range of possible values identified by eye is shown in yellow as a cross-check.

5.2 Correlation analysis

Since the baseline is continuously evolving, correlating the temporal scintillation is problematic. Instead, we cross-correlated the single time bins. This method is sufficient if there is enough structure in the spectral scintillation to differentiate time bins from each other. It has several advantages. First, assumptions on a regular grid in time are no longer required, which removes a potential source of systematic errors. More importantly, it completely removes the correlation peak at zero shift that arises from remaining intrinsic pulse-to-pulse variations because the temporal intensity information is dismissed. This effect would have tended to bias the measured shift toward smaller absolute values. In addition, baselines from telescope positions evaluated at different times can now be correctly taken into account. In summary, the cross-correlation function (CCF) is computed as $$\mathrm{CCF}(t_{\mathrm{E}\mathrm{f}\mathrm{f}},t_{\mathrm{F}\mathrm{A}\mathrm{S}\mathrm{T}})=\frac{1}{N_{\nu }}\sum _n{I}_{\mathrm{E}\mathrm{f}\mathrm{f}}(t_{\mathrm{E}\mathrm{f}\mathrm{f}},{\nu }_n)I_{\mathrm{F}\mathrm{A}\mathrm{S}\mathrm{T}}(t_{\mathrm{F}\mathrm{A}\mathrm{S}\mathrm{T}},{\nu }_n),$$(25)

where N_ν is the number of channels with center frequencies, ν_n, and the intensities, I, were normalized as described in Section 2 and their mean was subtracted.

In Fig. 9, we show the cross-correlation of each time bin with each other translated to a temporal shift as a function of observation time. This relation can be well described by a linear function with two free parameters. To investigate more complicated fit functions would require a solid uncertainty estimation of the CCF, which has not yet been accomplished because of strong correlations that depend on the amplitude of the signal. Since nearby data points are not independent, the cross-correlation of shifted (but otherwise identical) data are nonzero, even at lags that deviate from the true temporal shift between the two telescopes. Fortunately, the amount of data is large enough to allow for the creation of subsets. Thus, the uncertainty of the two fitted parameters can be estimated as follows. A fit to the full dataset produces values for the two parameters describing the temporal shift as well as values for the modulation index, m, and the scintillation time, t_s, as defined in Eq. (17). Then, the CCF was separated into subsets, where the only time bins at Effelsberg that were considered were those separated by 4t_s. The fit was repeated on these subsets and the standard deviation of the results for each fit parameter was taken as an uncertainty estimate.

For a one-screen model, a natural choice for the two fit parameters are the effective velocity, V_eff,∥, along the screen and the orientation, α_x, of the screen. Then, the temporal shift is given by Eq. (7). Results for this parametrization are shown in Table 2. The two parameters can be independently constrained because the orientation sets the time of zero shift and the effective velocity sets the slope of the time shift. For a two-screen model, the parameters obtained by the fit need to be translated according to Eq. (21) as $$\alpha \mapsto \arctan \left(\frac{\sin \left({\alpha }_x\right)+\beta \cos \left({\alpha }_x\right)}{\cos \left({\alpha }_x\right)-\beta \sin \left({\alpha }_x\right)}\right),$$(26) $$V_{\mathrm{e}\mathrm{f}\mathrm{f},\parallel }\mapsto \frac{V_{\mathrm{e}\mathrm{f}\mathrm{f},x,\parallel }}{\sqrt{1+{\beta }^2}},$$(27)

where $$\beta =\frac{\gamma \delta D_x{D}_{y,s}}{D_y{D}_{x,s}-D_x{D}_{y,s}{\delta }^2}.$$(28)

For comparison, the same analysis can be repeated by correlating the dataset from each telescope to itself instead, which yields the autocor-relation function (ACF). Then, the modulation index and scintillation time scale of the expected Gaussian ACF defined in Eq. (17) can directly be fitted without any shift. However, this means that the evolution of the telescope’s position has been neglected. Different positions on Earth view the scintillation pattern at different times, which affects the observation time at which a certain value of autocorrelation is to be expected. This effect can be accounted for by computing the expected shift along the baseline between the positions of the telescope at the two different times correlated, just as the shift due to the baseline given by Eq. (6) between two different telescopes was computed. Since the ACF drops to zero quickly, these single-station baselines are too small for α_x and V_eff,∥ to be independently measured. Instead, we fix them to the values obtained fitting to the CCF.

Table 2 shows the results for the four fit parameters for the CCF and for the ACF; for the latter, the results are shown both in the case where we are accounting for the additional effect of single-telescope baselines as well as neglecting it. The difference between both is statistically significant, even though the absolute difference is tiny. The absolute values of m are quite similar reflecting the visible similarity of the observations. The more sensitive FAST data are in statistical tension with both the CCF and the ACF results from Effelsberg. This can also be explained by contributions of data gaps, RFI, and noise because small differences can also be seen in the inferred scintillation time. Thus, the result allows only for small differences of the signal arriving at the two telescopes, confirming the dominant one-dimensionality of the screen. Marthi et al. (2021) reported the opposite result of a low cross-correlation between simultaneous observations with a very long baseline. However, the corresponding observations took place before the mentioned transition in 2020, when the farther screen strongly affected the dynamic spectra, which means that the results are not directly comparable. For each simultaneous observation, the sensitivity on the screen orientations is at least four times higher than that obtained from long-time monitoring in Sprenger et al. (2022).

Fig. 9 Evolution of the cross-correlation at December 4–5, 2023 (top) and from May 20, 2024 (bottom). Assuming a regular grid, the corresponding time shift between the telescopes is shown. Positive values mean that the scintles arrive at FAST first. The horizontal white line marks vanishing shift between the scintillation.

5.3 Model fitting

To infer the screen’s orientation, distance, and velocity, we performed a least-squares fit and a Markov chain Monte Carlo (MCMC) analysis using EMCEE (Foreman-Mackey et al. 2013). This was done on two different combinations of data and model. The results are shown in Table 3. The available data consist of the arc curvatures, screen orientations, and effective velocities at the two epochs of simultaneous observations, as listed in Table 1 and Table 2, where the results from the CCF analysis were used. In addition, we used all the data that had previously been used for the parameter inference in Sprenger et al. (2022), which consists of Effelsberg observations from 2019 to 2022. To differentiate, we refer to the two simultaneous epochs as the “new data”, while the combination of both sources is labeled “all data”. The new data were incorporated into the updated likelihood function as independent Gaussian random variables.

First, a two-screen model was fit to all the data. As described in Section 4, the screen orientation evolved over time. Hence, we introduce the two new parameters α_x,2023 and α_x,2024 for the screen orientation at the two simultaneous epochs in addition to the parameters defined in Sprenger et al. (2022). The projected screen velocity, V_x,∥, was derived from the absolute value and orientation that were used as fit parameters. The obtained distance, D_x, of the closer screen is much larger than the previous constraint of 127.1 ± 1.5 pc. This deviation is partially offset in the multidimensional space of parameters with the χ²/ndf increasing from 0.95 to 1.1. Nevertheless, this finding points toward the assumption of a constant screen distance possibly being wrong. A shift of ~5 pc over three years where the line of sight only travels ~0.0001 pc along the screen is rather unlikely, which makes a nonzero thickness of the screen a more physical alternative than an inclined thin screen. The screen’s velocity is still in agreement with a constant value. The orientation, α_x, of the screen is evolving in a consistent direction and remains incompatible with being constant in a two-screen model.

After including the two simultaneous observations, the inferred parameters for the second screen are the distance from the pulsar of D_y,s = 88 ± 47 pc, whose posterior does not follow a Gaussian (see Fig. 10), the orientation of α_y = −41.7 ± 3.0°, and the projected velocity of V_y,∥ = −29 ± 66 km/s. These values do not represent tighter constraints than those previously known. The second screen is also phenomenologically absent in the cross-correlation, which is compatible with a single 1D screen which is the only model predicting nearly exact identity of the two dynamic spectra beyond the shift. However, the image modulations described by Sprenger et al. (2022) and corresponding horizontal features on the arcs are equally visible in these observations and clearly point to the continuous presence of the second screen.

For comparison, we investigated a model using the new data only and setting the screen parameters orientation, distance, and velocity to be independent at each observation. Since the second screen was not constrained in this case because its shift and arc curvature were not measured, a one-screen model was employed. The same prior information on the pulsar parameters and the screen velocities was used as in the two-screen model. The results are shown in the third column of Table 3.

In both models, the pulsar parameters cannot be further constrained. Sampled posteriors for the screen angles and distances in the two-screen model and the one-screen model are presented in Fig. 10 and Fig. 11, respectively. The four contours go up to the 2σ region in equal steps as is the default setting by corner (Foreman-Mackey 2016). In the one-screen model, the screen distance is compatible with being constant and with the previously inferred one, although the uncertainties are too large to exclude an evolution on the scale hinted at by the two-screen results.

Fig. 10 Posterior distribution from 676 000 samples for the two-screen model. Left: Screen orientations at the two epochs of simultaneous observations. The red lines show the 1σ range of the screen orientation (${\alpha }_x^{SS}$) from Sprenger et al. (2022) and the black line marks identical angles. Right: Distance, Dx, of the first screen from Earth against the distance, Dys, of the second screen from the pulsar.

6 Conclusions

We performed two simultaneous observations of PSR B1508+55 at Effelsberg and with FAST. For both pairs of dynamic spectra, we formed intensity cross-secondary spectra, imaged the screen, determined the curvature parameter of scintillation arcs, and performed a cross-correlation analysis. At both epochs, the pulsar was scattered by an ISM screen into a straight line on the sky with only minor deviations. The effects of a known second screen were visible, but with a negligible impact on the aims and conclusions of this study.

We have successfully demonstrated that the scattered image of the pulsar can be imaged with 0.1 mas resolution, which corresponds to ~0.01 AU in the ISM, without using visibilities and employing only two telescopes. This is possible because of the inherently interferometric information present in scintillation as well as the evolving projected baseline due to the rotation of the Earth and the high declination of the source. Inverted arclets are not required for this method but a distinct scintillation arc is necessary because the region along the arc needs to be identified. As a result, this method works equally for thin arcs and arcs with arclets, but will struggle with fuzzy secondary spectra and with arclet apexes too far offset from the parabolic arc. Thus, it can also be extended to cases of multiple arcs, each treated independently in a one-screen approximation. The greatly reduced logistical requirements compared to classical VLBI imaging make this a useful method for future studies of scattering screens.

The screen orientation at both epochs could be measured with a very high precision that exceeds that obtained from several years of observations of annual arc curvature variations with Effelsberg alone. This precision enabled by the long baseline and the presented cross-correlation method proves that the screen orientation is evolving between observations. The new measurements have been combined with data from Sprenger et al. (2022) to allow for a likelihood analysis, leading to improved constrains on the orientation, distance, and velocity of the screen.

These measurements yielded hints to an evolving distance or thickness of the closer screen on a scale of a few parsecs, while remaining in agreement with a model of two strongly scintillating, thin, 1D screens. The scattered image of the pulsar on the closer screen is about 2 AU×0.05 AU in terms of physical size. At a distance of around 125 pc and a frequency of 1.4 GHz, this indicates (in numbers) that under the stationary phase condition, the local electron column density contains gradients corresponding to a change of 0.01 pc/cm³ in dispersion measure per astronomical unit in one direction, while at no point exceeding 0.0003 pc/cm³ per astronomical unit in the other direction. Studies such as the one presented here can inform models of the small-scale structure in the ISM, answering the question of whether structures corresponding to 1D screens are common or, in fact, exceptional.

Fig. 11 Posterior distribution from 240 000 samples for the one-screen model. The red lines show the 1σ range of the screen orientation (${\alpha }_x^{SS}$) from Sprenger et al. (2022) and the black line marks equality of the parameters. Left: screen orientations at the two epochs of simultaneous observations. Right: screen distances at the same epochs.

Data availability

Movies associated to Figs 6 and 7 are available at https://www.aanda.org

Acknowledgements

Based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. This work made use of the data from FAST (Five-hundred-meter Aperture Spherical radio Telescope). FAST is a Chinese national mega-science facility, operated by National Astronomical Observatories, Chinese Academy of Sciences. We thank GuangXing Li and JinChen Jiang for their help with data transfer. XS thanks Jumei Yao and Weiwei Zhu for their collaboration on the earlier FAST observations of B1508+55, which laid the foundation for the simultaneous observations. XS is supported by NSFC grant No. 12373025. TS thanks Nataliya Porayko for her useful comments being the internal referee at the MPIfR.

References

All Tables

All Figures

	Fig. 1 Telescope positions as viewed from the pulsar throughout the observations. Color represents the time since the start of the observation. Position of the land masses at the beginning (blue) and end (red) of the observations and that of the baseline between the telescopes are also shown.
In the text

	Fig. 2 Dynamic spectra of the two simultaneous observations. The data have been normalized to a mean of unity in every time bin and the band has been cut to the optimal subband. The time axis is synchronous and measured from the beginning of simultaneous observations.
In the text

	Fig. 3 Definitions relevant to this study. The directions of the arrows are examples and do not define the positive direction of signed quantities. Alternative notations in the one-screen case are given in parentheses.
In the text

Fig. 4 Cross-secondary spectra (top row) and cross-correlations (bottom row) from December 4–5, 2023. The observation has been split into four parts of equal length in time, shown in chronological order. The colors denote the complex phase of the cross-secondary spectra, while their brightness denotes the power. The cross-correlation functions are obtained as the Fourier transform of the cross-secondary spectra. They illustrate the temporal shift corresponding to the linear phase evolution in the cross-spectra.

In the text

	Fig. 5 Cross-secondary spectra (top row) and cross-correlations (bottom row) from 20 May 2024. The setup is identical to Fig. 4.
In the text

	Fig. 6 Image of the scattering screen at December 4–5, 2023. The central point cannot be imaged and is instead shown as fully saturated to highlight the position of the pulsar. The color scale is linear and normalized to the highest value. A resolution of about 0.1 mas is achieved. The evolution of the image over the course of the observation can be found as an online movie.
In the text

	Fig. 7 Image of the scattering screen at 21 May 2024. See Fig. 6 for more details. An online movie is available.
In the text

Fig. 8 Arc curvature measurement. Bottom row: secondary spectra with equal color scale gauged by the 99.95 percentile. Mid row: sum of eigenvalues for all subsets in frequency as a function of subsets in time and trial ζ. To row: corresponding mean over all subsets in time. The best-fit value for ζ or the scintillation arc is shown in red. In the top and bottom plots its 1σ region is shown in red and the range of possible values identified by eye is shown in yellow as a cross-check.

In the text

	Fig. 9 Evolution of the cross-correlation at December 4–5, 2023 (top) and from May 20, 2024 (bottom). Assuming a regular grid, the corresponding time shift between the telescopes is shown. Positive values mean that the scintles arrive at FAST first. The horizontal white line marks vanishing shift between the scintillation.
In the text

Fig. 10 Posterior distribution from 676 000 samples for the two-screen model. Left: Screen orientations at the two epochs of simultaneous observations. The red lines show the 1σ range of the screen orientation (${\alpha }_x^{SS}$) from Sprenger et al. (2022) and the black line marks identical angles. Right: Distance, Dx, of the first screen from Earth against the distance, Dys, of the second screen from the pulsar.

In the text

	Fig. 11 Posterior distribution from 240 000 samples for the one-screen model. The red lines show the 1σ range of the screen orientation (${\alpha }_x^{SS}$) from Sprenger et al. (2022) and the black line marks equality of the parameters. Left: screen orientations at the two epochs of simultaneous observations. Right: screen distances at the same epochs.
In the text