# Disentangling interstellar plasma screens with pulsar VLBI: Combining auto- and cross-correlations

###### Abstract

Pulsar scintillation allows a glimpse into small-scale plasma structures in the interstellar medium, if we can infer their properties from the observed scintillation pattern. With Very Long Baseline Interferometry, and working in delay-delay rate space (after a Fourier transform of the dynamic spectra) where the contributions of pairs of images to the interference pattern become localized, the scattering geometry and distribution of scattered images on the sky can be determined if a single, highly-anisotropic scattering screen is responsible for the scintillation. However, many pulsars are subject to much more complex scattering environments where this method cannot be used. We present a novel technique to reconstruct the scattered flux of the pulsar and solve for the scattering geometry in these complex cases by combining interferometric visibilities with cross-correlations of single-station intensities. This takes advantage of the fact that, considering a single image pair in delay-delay rate space, the visibilities are sensitive to the sum of the image angular displacements, while the cross-correlated intensities are sensitive to the difference, so that their combination can be used to localize both images of the pair. We show that this technique is able to reconstruct the previously published scattering geometry of PSR B0834+06, then apply it to simulations of more complicated scattering systems, where we find that it can distinguish features from different scattering screens even when the presence of multiple screens is not obvious in the Fourier transform of the visibilities. This technique will allow us to both better understand the distribution of scattering within our local interstellar medium and to apply current scintillometry techniques, such as modelling scintillation and constraining the location of pulsar emission, to sources for which a current lack of understanding of the scattering environment precludes the use of these techniques.

###### keywords:

pulsars:general – ISM:general – ISM:structure – techniques:interferometric^{†}

^{†}pubyear: 2018

^{†}

^{†}pagerange: Disentangling interstellar plasma screens with pulsar VLBI: Combining auto- and cross-correlations–Disentangling interstellar plasma screens with pulsar VLBI: Combining auto- and cross-correlations

## 1 Introduction

Pulsar scintillation, variation in the observed flux of a pulsar with time and frequency due to refractive and diffractive effects of the intervening plasma in the interstellar medium (ISM), can be a nuisance for high-precision pulsar timing experiments. The extra delays imparted by the ISM reduce the precision of pulsar timing, especially at low frequencies where pulsars are brightest but interstellar scattering is strongest. On the other hand, pulsar scintillation can be advantageous to studies of the ISM itself. Pulsars, as unresolved coherent radio sources, act as probes of the ionized plasma and can provide insight into the structure of the plasma on very small scales. For decades, pulsars have been used as a probe of the electron density power spectrum in the ISM (e.g. Lee & Jokipii, 1976; Armstrong et al., 1995; Xu & Zhang, 2017) and the distribution of free electrons within our galaxy (e.g. Cordes & Lazio, 2002, 2003). These models have been used to place constraints on the distances to dispersed radio sources such as RRATS and Fast Radio Bursts (e.g. Keane et al., 2011).

In the early 2000’s, Stinebring et al. (2001) noticed that the secondary spectra (the 2-dimensional power spectra of the dynamic spectra, which is the frequency spectrum as a function of time), of some pulsars contain remarkably organized parabolic structure. Since then, it has been found that the secondary spectra of many pulsars, especially those bright and nearby, show this parabolic structure (e.g. Stinebring et al., 2003; Putney & Stinebring, 2006; Stinebring, 2007), which can be explained by highly-anisotropic scattering at a thin screen localized between us and the pulsar (Walker et al., 2004; Cordes et al., 2006).

In some cases (e.g. PSR B0834+06 (Hill et al., 2005; Brisken et al., 2010) and PSR B1737+13 (Stinebring, 2007)), distinguishable inverted arclets with the same curvature as the main parabola and apexes along the main parabola are also visible in the secondary spectrum. In an anisotropic, thin-screen interpretation of scintillation, each of these arclets corresponds to the interference of a single image of the pulsar on the scattering screen with the remaining images of the pulsar, allowing one to map and track the individual images over time and frequency through wide-band, multi-epoch observations. Hill et al. (2005) followed the individual arclets in the secondary spectrum of PSR B0834+06 over 26 days, and found that the observed motion of the arclets through the secondary spectrum can be accounted for solely by the motion of the pulsar behind the screen. They suggest that the stationarity of the images in time is indicative of compact lensing regions, regions where the electron column density varies rapidly, embedded within the scattering screen. Hill et al. (2003) and Brisken et al. (2010) find that the dependence of image location on frequency is very weak, another hint that the individual lensing regions are compact. Stinebring (2007) found that the parabola curvature does not vary significantly over 20 years for some pulsars, suggesting that the screen location and orientation are constant over decades, while Hill et al. (2005) find that individual arclets persist over a 26-day observation, evidence that the electron density perturbations on the scattering screen are unchanged over this timescale.

The stationarity of the compact lensing regions within the scattering screen allows us to move beyond a statistical treatment of scintillation to modeling and predicting scintillation patterns for individual pulsars. Along this line, Simard & Pen (2018) develop a model of pulsar scintillation from corrugated refractive plasma sheets within the interstellar medium that makes clear predictions for both the flux variations and the motion of the images through the sheet with time and frequency, while Gwinn (in prep.) develop a model based on 1-D plasma structures in the ISM. By comparing the evolution of images with time and frequency to these models, pulsar scintillation can be used to investigate the characteristics of the plasma structures in the ISM, in the same way that quasar extreme scattering events are used to distinguish between various plasma lens models (e.g. Clegg et al., 1998; Bannister et al., 2016; Tuntsov et al., 2016; Dong et al., 2018; Er & Rogers, 2018). This stationarity may also allow the screens themselves to be used to achieve unprecedented spatial precision at the pulsar. In many pulsars, the scintillation pattern is different for different parts of the pulse profile, evidence that the interstellar screen resolves the pulsar emission region (e.g. Backer, 1975; Smirnova et al., 1996; Gupta et al., 1999; Gwinn et al., 2000, 2012; Johnson et al., 2012; Pen et al., 2014; Main et al., 2017). In order to place physical sizes and separations on pulse components using these measurements, one must have an understanding of the scattering geometry, particularly the distances to and the resolving powers of the scattering screens.

When parabolic arcs are present in the secondary spectrum, the curvature of these arcs, which depends on the distances to the screens and pulsar as well as the velocities of the pulsar, screens and observer in the direction of scattering (see Section 2.1), can inform the scattering geometry. However, in order to break the degeneracy between the distances, velocities, and angle of scattering, additional information is needed. Putney & Stinebring (2006) assume that the velocity of the pulsar dominates over all velocities in the system and is aligned with the direction of scattering, allowing them to determine the fractional distances, , of the scattering screens of many pulsars from the curvature of the parabolic arc in the secondary spectrum. Others (e.g. Smirnova et al., 2014; Popov et al., 2016; Shishov et al., 2017; Fadeev et al., 2018) assume that the scattering is isotropic and use the correlated flux of the scattered pulsar on global and space baselines to determine the angular size of the scattered image. Combining this with the scattering timescale, one can estimate the fractional distance to the scattering screen. Brisken et al. (2010) demonstrated that by using Very Long Baseline Interferometry (VLBI), one can map out the scattered images using the phases of arclets in the secondary cross-spectra (related to the Fourier transform of the interferometric visibility; see Section 2.2) along multiple baselines. Using this technique, they were able to reconstruct the locations of the scattered images of PSR B0834+06 in two angular dimensions on the sky and measure the angular size and orientation of the scattering screen. This can be combined with the locations of the arclets corresponding to the scattered images to measure the fractional distance to the screen without making assumptions of the scattering geometry.

Many scattering systems are much more complex than a single scattering screen. In some, (e.g. B1133+16 and B0329+54 (Putney & Stinebring, 2006)), multiple parabolas are visible, evidence of multiple scattering screens along the line of sight to the pulsar, while many pulsars show much less organized structure in the secondary spectrum. It is still unknown if this disorganized power is due to many scattering screens at different distances and orientations, or due to a distributed, isotropic component of scattering.

In this work, we present a technique of combining interferometric cross-correlations along global baselines with auto-correlations from the same observations to disentangle individual scattering screens in complex scattering systems. This method not only allows one to measure the distance to and velocity of individual scattering screens, but also allows one to map out the distribution of scattered images on each screen. This will help to inform our understanding of the distribution of scattering material in the local ISM and allow us to extend techniques of constraining separations of pulse components to systems with complex scattering.

In Section 2.1, we recall the theory of scattering by a thin, anisotropic screen. In Sections 2.2 and 2.3, we review the theories of using VLBI visibilities and intensities from simultaneous observations at multiple stations to measure the locations of images on the sky and the distance to the scattering screen. In Section 3, we introduce the technique of combining auto- and cross-correlations and explain how this can be used to investigate scattering in multi-screen systems. In Section 4 we verify that this technique can reproduce the results obtained by Brisken et al. (2010) through VLBI observations of PSR B0834+06, and in Section 5 we simulate multi-screen systems in order to test this technique further. Finally, we present closing remarks and ongoing applications of our technique in Section 6.

## 2 Theory

Term | Symbol | Description |

Quantities derived from the intensities | ||

Dynamic spectrum | The observed intensity of the pulsar at a frequency and time , which varies due to scintillation. | |

Conjugate spectrum | The Fourier transform of the dynamic spectrum. | |

Intensity cross secondary spectrum | The cross-correlation of the dynamic spectra at stations and , in the Fourier domain. Sensitive to the angular separations between pairs of images. | |

Quantities derived from the visibilities | ||

Visibility dynamic cross-spectrum | The visibility of the pulsar between two stations at a frequency and time . | |

Visibility conjugate spectrum | The Fourier transform of the visibility. | |

Visibility secondary cross-spectrum | Sensitive to the sum of the angular separations of images from the pulsar. | |

Quantities derived from the intensity and visibility cross spectra | ||

Quaternary spectra | Sensitive to , the larger angle in each pair of interfering images | |

Sensitive to , the smaller angle in each pair of interfering images |

### 2.1 Scattering by a thin screen

The anisotropic thin screen model of pulsar scintillation (Walker et al., 2004; Cordes et al., 2006) explains the observed parabolic structure in the secondary spectra of pulsars by considering the interference of many images along a line on the sky and at a single distance between us and the pulsar. In the dynamic spectrum, , the interference of any two images and with angular separations and relative to the dominant core of the scattered flux of the pulsar (which our calibration centers at the origin of the secondary spectrum, and which we will refer to as the ‘core image’) create a fringe pattern, which under the Fourier transform to the conjugate spectrum, , leads to power at

(1) |

and

(2) |

where and are the Fourier conjugate variables of time and frequency respectively, is the speed of light, and is the wavelength of observation. Note that this will also contribute power at (,) due to the symmetry about switching the indices and , so that amplitude of the conjugate spectrum is point symmetric about the diagonal. Since and are the relative Doppler shift and geometric delay between the images, we will refer to them as the Doppler shift, , and delay, , respectively. Note that is also sometimes called the delay rate, and is the time derivate of . When the pulsar is moving towards an image, the geometric delay is decreasing, and , so that features due to an individual image appear at low and move towards higher over time.

(3) |

and

(4) |

are the effective distance and velocity respectively. , , and are the velocities of the pulsar, observer, and screen respectively, and and are the distance from the observer to the pulsar and screen respectively. The scintillation pattern moves across Earth with the speed , where the subscript indicates that this is the projection onto the direction of , but in the opposite direction of the projection of onto . When the scattered images are aligned and from a single screen, both and are constant in the system, so that the interference of images with the unlensed () image of the pulsar leads to power along a parabola with curvature

(5) |

We can also consider the case where multiple screens are present along the line-of-sight to the pulsar. In this case, we will see images due to paths that experience high bending angles at more than one screen. Paths that go through two screens, which we label and , have a geometric delay, , and Doppler shift, , relative to the the central image of the pulsar of:

(6) |

and

(7) |

where the subscript indicates that the light is scattered by both screens, and are the angular separations of the image from the core image after being lensed by screens (closer to the pulsar) and (closer to the observer) respectively. , and are the distances from the observer to the pulsar, screen and screen respectively, while , and are respectively the velocities of the pulsar, screen and screen .

### 2.2 VLBI measurements of pulsar scattering

As shown by Brisken et al. (2010), if one can identify points in the secondary spectrum where , one can reconstruct the scattered image of the pulsar and the lensing geometry from the interferometric visibility dynamic spectrum, . The phase of a point in the conjugate visibility spectrum, , where indicates the 2-dimensional Fourier transform (see Table 1 for definitions of all quantities derived from the dynamic spectra), due to the interference of two images, and , contains a term due to the phase imparted by the lens as well as a term due to geometric delay, which depends on the projected baseline :

(8) |

Thus in the visibility secondary cross-spectrum, , this same point has the phase

(9) |

(We will use the subscripts and to indicate quantities derived from the visibilities and intensities respectively.) When , , and, with multiple baselines, one can measure . In order to use this technique, one must be able to distinguish regions in the secondary spectrum where , ideally from the apexes of inverted arclet, or, if these are not present, by choosing points that lie closest the main parabola. Thus, this relies on a system in which the secondary spectrum is well-organized. Once is known for these points where is assumed to be , one can use the Doppler shifts and geometrical delays of these points, equations (1) and (2), to determine and respectively. If points are from a single, highly-anisotropic screen, they can be combined statistically to improve these measurements. Once the distance to and orientation of the screen are known, the astrometry of the images on the screen can be improved by using the locations of the arclets in the secondary spectrum, allowing 100 microarcsecond precision (Brisken et al., 2010). Using this technique, Brisken et al. (2010) measured the angular locations of bright scattered images of PSR B0834+06 to 100 as precision along with the distance to and orientation of the the scattering screen.

### 2.3 Simultaneous single-dish measurements of pulsar scattering

In a single screen system, the dynamic spectra, from simultaneous observations at multiple stations can be used to measure the distance to the screen without correlating baseband voltages. Galt & Lyne (1972) measure the scintillation speed and orientation of the scattered images for PSR B0329+54 by monitoring the offset in the scintillation pattern between Jodrell Bank and Penticton, British Columbia over an entire day as the Earth rotates. At the time, the phenomenon of scintillation arcs was not yet known, but now the measured scintillation pattern speed and orientation can be combined with the curvature of the parabola, equation (5), to measure .

Instead of using the cross-correlation of the dynamic spectra at two stations to measure the scintillation speed, we construct the cross secondary spectrum, , where the subscript indicates that this is constructed from intensities as opposed to visibilities, the superscript indicates taking the complex conjugate, and and label the two stations separated by the projected baseline . In this case, the phase of is simply that due to the delay between the images at the two stations:

(10) |

Note the similarity to the phase of the visibility secondary cross-spectrum, , given in equation (9). While the visibility secondary cross-spectrum is sensitive to the sum of the scattering angles, the intensity cross secondary spectrum is sensitive to the difference.

For a single, anisotropic screen,

(11) |

where is the angle of the projected baseline in the - plane and is the orientation of the line of scattered images in the - plane. Note the similarity of this relation to that for points with in the interferometric analysis. The only difference is that when using the intensity cross secondary spectrum, this relation is generalized to all points in the secondary spectrum, including those for which . By constructing and multiple baselines, one can measure both and . In order to measure , one must have a measurement of the curvature of the parabola in the secondary spectrum, (equation (5)) as well.

This technique of measuring phase gradients in the intensity cross secondary spectrum is much more sensitive to small time delays than real-space cross-correlations of the dynamic spectra, and allows one to include only regions of the secondary spectrum that one is confident belong to the scattering screen of interest. It also has advantages over using the interferometric secondary cross-spectrum: With this technique one is able to use all points with high signal-to-noise in the secondary spectrum, rather than just those along the main arclet. In addition, the measured effective velocity is independent of the measured curvature of the parabola in the secondary spectrum, and thus this technique is much better suited to systems that don’t show clear inverted arclets.

While correlating voltages requires aligning the clocks at two stations to nanosecond precision, in this case, where typical delays between the scintillation patterns at the two stations are on the order of seconds, millisecond precision is sufficient. Finally, using this technique does not require recording raw voltages, as all that is needed is the gated dynamic spectra at each station. As a result, one can take advantage of pulsar backends at many stations, which can record over a much wider bandwidth and with higher bit depth than most VLBI systems. For the brightest pulsars, this bit depth is needed to ensure that the brightest events are not degraded due to saturation, and also allows for better characterization and removal of RFI, while the large bandwidth both increases the number of scintles sampled and therefore the signal-to-noise in the secondary spectrum, and allows us to examine changes in the scattering with frequency.

## 3 Combining interferometric measurements with simultaneous auto-correlations

In some scattering systems, the combination of scattering from many different screens results in regions in the secondary spectra where parabolic arcs cannot be distinguished. Here, we present a novel technique which combines the two techniques above in order to reconstruct the scattering geometry and scattered flux distribution in these cases.

By multiplying the visibility secondary cross-spectrum by the intensity cross secondary spectrum, we construct a quaternary spectrum, , where the phase due to the interference of two images at and is:

(12) |

Note that the phase is no longer dependent on . Meanwhile, multiplying the visibility secondary cross-spectrum by the conjugate of the intensity cross-secondary spectrum constructs with phase

(13) |

which is independent of . Thus, under the assumption that the power in a given pixel in the secondary spectrum is dominated by the interference of a single pair of images, we can measure the angular separations of both images from the pulsar, projected along the projected baseline . With multiple baselines, we can use this to determine the positions of the scattered images in coordinates.

By choosing only points where , we can create a sample of the scattered images along with their Doppler shifts and geometric delays relative to the central image. For each point in this sample, one is able to determine ( to ) and from equations (1) and (2) (or equations (7) and (6) if one suspects these images are doubly-scattered). (This can be extended to all points in the secondary spectrum, but the equations for the Doppler frequency shift and the geometric delay become more complex when considering interfering images that may be on screens at different distances from the observer.) One can then build a 3-D reconstructed ’image’ of the scattered pulsar, where the third dimension is the distance to the screen scattering each image.

Standard software correlators like DiFX (Deller et al., 2007) and SFXC (Keimpema et al., 2015) allow the user to retain the auto-correlations as well as the interferometric visibilities, making this method straightforward to implement in VLBI campaigns. This method comes with significant advantages over using solely the interferometric visibilities: Not only does this technique allow one to separate multiple anisotropic screens apparent in the secondary spectrum, but as it does not make assumptions about the anisotropy of the scattering, it may also be used to separate scattering screens by their distances or velocities when parabolic arcs are not visible in the secondary spectra.

## 4 Psr B0834+06

In order to ensure that the three methods discussed in this work are consistent with one another, we apply them to the scattering system of PSR B0834+06 observed with global VLBI in 2005 by Brisken et al. (2010). The pulsar was observed with a bandwidth of 32 MHz centered at 316.5 MHz using Arecibo (AR), the Green Bank Telescope (GB), Jodrell Bank (JB), and tied-array Westerbork (WB). We received the data after correlation with the DiFX software correlator to visibilities and auto-correlations with 244 Hz resolution and 1.25 s gated integrations. For more details on the observation, we refer the reader to Brisken et al. (2010). The baselines from Green Bank or Arecibo to Westerbork and Jodrell Bank are very similar, but visibilities on baselines to Westerbork have lower signal-to-noise than those on baselines to Jodrell Bank. Thus, we have used only Arecibo, Green Bank and Jodrell Bank in this analysis. The geometry of the baselines is shown in Fig. 1.

Brisken et al. (2010) present a very detailed analysis of the scattering screen, using not only the phase of the secondary cross-spectrum along the main parabola, but also making use of the discrete arclets visible in the secondary cross-spectrum, a method that is more robust to images offset from the main arc. Brisken et al. (2010) also examine the properties of the scattering screen, including its evolution with frequency, in great detail. Here, we only aim to use the secondary cross-spectra, cross secondary spectra, and their combination to reproduce the measurement of the screen distance, orientation, and velocity made by Brisken et al. (2010).

To construct the visibility dynamic cross-spectra, we averaged the on-pulse channelized visibilities to 5-second time integrations. We phase calibrated the data by performing a singular value decomposition (SVD) on each sub-band of each visibility dynamic cross-spectrum. We then constructed a model for each IF/baseline pair by keeping only the first two modes of the SVD. A single mode captures only variations that can be accounted for by a multiplication of slow changes in both time and frequency, and is insufficient to calibrate the data, while using more modes risks including the phases of the scintles (the bright and dark patches in the intereference pattern) in our model. The phase of each model was then removed from the corresponding visibilities.

To construct the dynamic spectra from the autocorrelations, we began by dividing the on-pulse dynamic spectrum by the off-pulse dynamic spectrum, in order to remove any time-dependent sensitivities and the bandpass. We then divided by the average intensity for each integration, in order to remove pulse-to-pulse variations. This can also decrease the signal-to-noise of the scintillation pattern, but since the coherence bandwidth of the scintillation pattern is much smaller than the observing bandwidth, such that many scintles are visible across the observing band, this loss is minimal. As pulse-to-pulse variations will correlate between the two stations on very short timescales, they can bias the measurement of the scintillation pattern delay to lower values and it is vital that they be removed. This is especially true for PSR B0834+06, which shows extreme amplitude modulation between pulses (Rankin & Wright, 2007; Gwinn et al., 2011). We normalized each spectrum by subtracting the overall median value and dividing by the root-mean-squared value of each spectrum. The JB spectrum contains RFI which we flagged by eye and removed. We then averaged the resulting dynamic spectra to 5-second integrations.

We calculated the conjugate spectra, for each visibility and intensity dynamic spectrum using a 2-D Fast Fourier Transform (FFT). In order to prevent artifacts due to the continuity assumption of FFTs, we padded the dynamic spectra to twice their size in both time and frequency using the mean amplitude of each dynamic spectrum. As this is equivalent to interpolating every second pixel in the conjugate spectrum, we then averaged adjacent pixels in both time and frequency in the conjugate spectrum. The resulting secondary spectra constructed from the intensities, , are shown in Fig. 2.

For the visibility analysis, we will work with the visibility secondary cross-spectra: . The phases of the secondary cross spectra for all three baselines, GB-AR, JB-AR, and GB-JB, are shown in Figure 3. Recall that in this case the phase at a given point in the secondary cross-spectrum dominated by interference between images and is . Along any given arclet, the phase is furthest from zero at low , where the interference is due to two images that are very nearby on the sky, and closest to zero at large , where the interference is due to two images that are roughly equal distances from the pulsar, but on opposite sides. We once again see, especially from the JB-AR and GB-JB spectra that the 1-ms feature does not follow the same phase trend as the main parabolic structure, indicating that those images are due to a scattering structure with either a different distance, different velocity, or different orientation (or some combination of all three) than the screen causing the main parabolic arc.

For the intensity analysis, we will need the intensity cross secondary spectra, , where the subscripts and indicate the two different stations. The phases of the cross secondary spectra for the three baselines are shown in Fig. 4. Recall that in this case, the phase at any point in the intensity cross secondary spectrum is , where and are the locations of the two images whose interference is leading to power at that point. Thus, for constant delay on a parabolic arc, the phase is lowest near low , where the images that are interfering are very close together, and highest at large , where the images interfering are far apart. If all features are from the same anisotropic scattering screen, we expect a linear phase gradient with . We note that the feature at negative and a delay of 1 ms, which Brisken et al. (2010) suggest is from a different screen than the main parabolic structure in the secondary spectrum, does not follow the linear gradient of the main parabola in the JB-AR and GB-JB spectra.

### 4.1 Separate Analysis of the Interferometric Visibilities and Autocorrelations

Visibility meas. | Intensity meas. | |
---|---|---|

(seconds) | (seconds) | |

-7.33 0.06 | -7.88 0.05 | |

-5.71 0.11 | -5.19 0.07 | |

-1.79 0.10 | -1.51 0.11 |

(pc) | (km/s) | (deg) | ||||

Visibilities | Intensities | Visibilities | Intensities | Visibilities | Intensities | |

AR-GB and AR-JB | 1220 20 | 1150 20 | 319 3 | 298 2 | -32.8 0.5 | -29.4 0.3 |

AR-GB and JB-GB | 1220 20 | 1140 20 | 319 3 | 296 2 | -32.1 0.4 | -33.6 0.4 |

AR-JB and JB-GB | 1170 60 | 1570 70 | 319 8 | 349 7 | -65.8 0.3 | -66.2 0.4 |

Brisken et al. (2010) | 1170 20 | 305 3 | -27 2 |

In our analysis of the interferometric secondary cross-spectra, we begin by measuring the curvature of the parabolic arc. Similar to the Hough transform used for this purpose by Bhat et al. (2016), we calculate the average power in pixels that lie along parabolas of different curvatures and fit a Gaussian to the resulting power against curvature curve using a least-squares fitting routine. We note that this method relies on the assumption that the brightnesses of the highly scattered images are all much less than the dominant core of the image.

We take the centroid of the best-fit Gaussian as the curvature and construct a parabolic mask centered on this main parabola and with a width of 0.29 mHz. We calculate the average complex value within all unmasked points that lie in each 5 consecutive bins, and take the phase of that average value. In order to place uncertainties on the measured phases, we take the uncertainty on the real and imaginary parts to be the standard deviation of the mean. The resulting phase against trends are shown for all three baselines in Fig. 5, where we see a linear trend indicative of a dominant single scattering screen. At large , we expect this to dissolve into random phases as random noise begins to dominate the measurement. The analysis of the simultaneous single-dish autocorrelations is similar; however as the single parabola and steady gradient in phase with suggest that a single screen is dominating the secondary spectra (apart from the 1-ms feature), we do not use a mask to select points to include in our average. Instead, we calculate the delay-averaged phase from all points with positive delays, again over 4 consecutive bins. The resulting trend, shown in Fig. 6, is well-described by a line for all baseline pairs.

For each baseline, we use a least-squares routine to fit a line to the phase against trend where the phases appear to be dominated by signal. In the visibility case, we use points interior to 20, 20 and 10 mHz for the GB-AR, JB-AR, and GB-JB baselines respectively while in the intensity case we use points interior to mHz for all baselines. Note that in all cases, these limits exclude the 1-ms feature from our analysis. From these fits, we calculate the delay in the scintillation pattern along all three baselines, of ; the delays are given in Table 2. The uncertainties provided are only the statistical uncertainties and do not take into account systematic effects, such as those due to calibration, and (in the case of the visibility analysis) errors in the curvature fit due to bright images offset from the dominant core of the image. Measurements of also depend on the frequency at which the curvature is measured. We take this to be the central frequency of the observation, but if the pulsar spectrum is steep (as is the case for many pulsars), lower frequencies may have a more substantial weight in the curvature. This can be corrected by Fourier transforming each frequency channel independently and scaling the Fourier frequencies for each channel by , where is the reference frequency at which further calculations will be done. As the fractional bandwidth in this observation is only %, we don’t make this correction here.

The delays measured from the intensities and visibilities agree to within 15% for all baselines, but those measured from the intensities are smaller than those measured from the visibilities in all cases. This is due to the fact that the phases where there is no signal are not randomly distributed over 2 in the intensity cross-spectra, as would be expected due to pure noise. This is likely due to residual pulse-to-pulse variations that have not been removed, in addition to artifacts from the pulsar binning that correlate between spectra. The effects of these can be mitigated with a more thorough reduction and calibration of the data, and can be avoided by not included regions that are dominated by noise in the calculation of the phase. The interferometric observations are less affected by this, as we are choosing to use only points with high signal-to-noise in the analysis. However, the visibilities are subject to their own phase uncertainties due to self-calibration of the phases, which can have similar effects.

We calculate (the orientation of the scattering screen in the plane) and the for each baseline pair using the method described in Sections 2.2 and 2.3. By combining this with the measured curvature, we calculate . The results are shown in Table 3. Again, these quoted uncertainties do not take into account systematic errors. As we can see from Table 3, the values for , , and measured using the auto-correlations taken simultaneously at various stations are consistent with those measured from the interferometric visibilities from the same baseline pair to within 10% for the AR-GB/AR-JB and AR-GB/JB-GB baseline pairs, but only consistent to within 35% for the AR-JB/GB-JB pair. We also note that the AR-JB/GB-JB pair gives a much larger measurement of . This is not unexpected - these two baselines are separated by less than , and so the speed of the scintillation pattern perpendicular to these two baselines is not well constrained. As a result, they favour an orientation of the scattering screen that is more aligned with the two baselines. The measurements of from the perpendicular baseline pairs (AR-GB/AR-JB and AR-GB/JB-GB) are consistent with the original measurement by Brisken et al. (2010) to within 25%, while the measurements of and from these same baselines pairs is consistent with the original measurement to within 5%.

While the outcomes of these two methods are similar, these methods differ significantly in their ease of implementation. Without the need for baseband voltages, cross-correlating intensity dynamic spectra allows us to take advantage of the wide-bandwidth receivers and high bit depth pulsar recorders available at many stations, and doesn’t require the nanosecond-precision temporal alignment and significant processing power necessary for correlation. This makes using the single-station intensities a much more practical method of measuring the distances to single-screen systems.

### 4.2 Combined Analysis

In addition to using the interferometric visibilities and auto-correlations separately to measure the distance to scattering screens, the visibilities and auto-correlations can be combined. In Section 5, we show with simulations how this can recover the distances to multiple screens, even when those screens are difficult to distinguish in the secondary spectrum. In this section, verify that this analysis is consistent with the separate analyses described in Section 4.1 by applying it to the same 2005 Brisken et al. (2010) observations of PSR B0834+06.

We begin with the visibility secondary cross-spectra, , and the intensity cross secondary-spectra, , shown in Figs. 5 and 6. Then, using the phases of the quartenary spectra, as described in Section 3, and assuming that each pixel in the secondary spectrum is dominated either by noise or by power due to the interference of only one pair of images, we determine and , the angular locations of the two images interfering to produce power at each pixel. This can be done for each pair of baselines; the results for the AR-GB/AR-JB pair are shown in Fig. 7. Once the angles have been measured, they can be combined with and to measure (or ) and at every point, where is now measured parallel to the separation vector between the two images. and measured using the AR-GB/AR-JB baseline pair are shown in the right-hand panels of Figs. 8 and 9 respectively. Note that where the phase is dominated by noise, the phase angles will be randomly distributed between 0 and 2 radians, resulting in large values inferred for and , which in turn place the estimated distance for that point at small , when the screen is close to the pulsar.

Using only points with mas, we then create histograms of the and and distributions. These are shown in the left-hand panels of Figs. 8 and 9 for and respectively, for the AR-GB/AR-JB baseline pair. We fit a Gaussian to these histograms in order to determine the best-fit values for and , indicated by the orange vertical lines on Figs. 8 and 9 and given for all baseline pairs in Table 4. (For comparison with the separate analysis of the visibilities and intensities, we have calculated from our best-fit values of .) To calculate the orientation of the scattering screen, , we take a subset of these points for which mas and the signal-to-noise in the absolute value of the JB-GB visibility secondary cross-spectrum (the noisiest spectrum) is greater than 10.3 (this value is empirically chosen to exclude points that are dominated by noise but for which mas). We plot these points in the , plane, and fit a line through them. The resulting values of for all baseline pairs are given in Table 4.

Next, we take points for which 1 mas, and combine with the power at each pixel in order to reconstruct the scattered image of the pulsar, shown in Fig. 10. The brightness of the points here represents the unnormalized flux in log-scale, assuming that our limit of values enforces that all points in the secondary spectrum are due to interference of the image with the core image of the pulsar. From this image, it is apparent that our geometry is reflected about the diagonal compared to that of Brisken et al. (2010). However, our mapping causes the pulsar to be moving towards images at negative , as expected since for these images the delay is decreasing so the delay rate should be negative. This is also seen in multi-epoch observations of PSR B0834+06 - arclets arise at negative and move towards positive (Simard et al. in prep.). As such, we are confident in our orientation of the screen. We also note the contribution of the 1-ms feature, which does not fall along the line traced by the other images. In plotting it here, we have ignored the fact that the phases measured for the 1-ms have likely wrapped around radians. Brisken et al. (2010) use the fact that the image locations are expected to vary very little with frequency to break this degeneracy - this is another reason that wide-bandwidth observations are vital to scintillation studies.

(pc) | (km/s) | (deg) | |
---|---|---|---|

AR-GB and AR-JB | 1200 100 | 300 20 | -25.0 0.1 |

AR-GB and JB-GB | 1200 100 | 300 20 | -25.3 0.2 |

AR-JB and JB-GB | 1200 300 | 300 30 | -25.4 0.2 |

Brisken et al. (2010) | 1170 20 | 305 3 | -27 2 |

We can also consider the 1-ms feature, and use it to constrain the distance to the second screen. First, we must take into account the fact that there may be phase wrapping at these high delays in the secondary cross-spectra. Brisken et al. (2010) look at a number of characteristics of these images, and find that phase wrapping most likely has occurred on the baselines to JB. We don’t repeat those diagnostics here, but will continue by assuming that phase wrapping has occurred on those baselines. Then, we choose points in the 1-ms feature for which mas. This gives us 32 points centered around (, ) mas in , coordinates. As these images are closely clustered on the sky, we use the central location in the following analysis. We also calculate the average and for the 1-ms feature, and find mHz and s.

Once we know the location of the 1-ms feature as well as the properties of the scattering screen causing the main parabola, we can proceed to determine the full geometry of the PSR B0834+06 system using a geometrical analysis, as shown by Liu et al. (2016). Since a parabola cannot be drawn through the apexes of the arclets in the 1-ms feature and through the origin, one can assume that this feature is due to radiation that is scattered by multiple screens. (Images that are not aligned with the images that comprise the main parabolic arc could also cause this.) We will assume that that the second screen is in front of the screen that is responsible for the main parabolic arc. (We could also assume, as Liu et al. (2016) do, that the screen is behind the screen responsible for the main parabolic arc. In this case, one would have to assume that the main screen is perturbed only in 1-dimension, so that it is only able to bend light in a direction parallel to the line of images that make up the main parabola.) To determine the distance to the second screen, screen , we require the distance to screen (the screen closer to the pulsar), as well as the scattering angles at both screens. In the case of PSR B0834+06, the pulsar’s proper motion is large and the screens are mid-way between us and the pulsar, so we expect that the pulsar’s velocity is the only one that must be considered in the system. (However, the assumption that the screen closer to the observer, screen , has no significant velocity would be sufficient to solve this system.) Under this assumption, equation (7) simplifies to

(14) |

From this, we calculate mas, and we assume that it is oriented along the same line of scattered images as the majority of the images in the main parabola. We then use equation (6) to determine the distance to the second screen, and find pc.

We’ve shown in this section that an analysis including both the intensity and visibility cross-spectra from a VLBI observation not only has the power to replicate the measurement of the scattering geometry made using the visibilities separately, but also allows one to easily map the points in the secondary spectrum that are due to interference of images with the central bright image of the pulsar, which can be used to reconstruct the scattered image of the pulsar. Tracking these images and the evolution of their locations and fluxes with frequency and time is crucial to revealing the astrophysical origin of the compact structures in the interstellar medium responsible for scintillation arcs, making this combined analysis a powerful technique that comes at little extra cost.

## 5 Simulations

All Sims. | ||

(MHz) | 326 | |

BW (MHz) | 32 | |

2048 | ||

duration (s) | 6000 | |

1024 | ||

(km) | 3000 | |

(pc) | 500 | |

(pc) | 15 | |

(pc) | 100 | |

(km s) | 400 | |

Sim. A | Sim. B | |

(km s) | -279 | -279 |

(km s) | -160 | -720 |

In Section 4, we applied the technique of combining visibility and intensity cross-spectra to measure the distance to the screens in a simple scattering system, but the true power of this technique comes when applied to systems that are more complex. In this section, we use simulations of scattering systems with two screens to show that combining the visibility and intensity cross-spectra allows one to separate features in the secondary spectrum based on the distance to the screen, and ultimately to determine the distances to multiple screens when many are present in the scattering system.

We use a 1-D simulation in which the pulsar’s proper motion, the scattering direction, and the baseline are all aligned for simplicity, although with multiple baselines this can be generalized to the 2-D case. We begin by choosing the size and channelization of our dynamic spectrum and defining the geometry of the scattering system, including the distances to both screens. Then, we randomly select 100 angular locations for the images on each screen, choosing a maximum angular separation from the pulsar so that the resulting scintillation pattern will be well-resolved with our chosen integration time and channelization. Each image is assigned a magnification by choosing a random magnification between 0 and 1 and multiplying it by a decaying exponential in angular separation from the pulsar, with a separation of a third of the maximum separation. For simplicity, we assume that the image position and flux do not vary over the bandwidth of our observation. From both theoretical predictions (Simard & Pen, 2018) and observations (Hill et al., 2005; Brisken et al., 2010), we expect these to change very slowly with frequency. We assign the pulsar an electric field magnitude of in arbitrary units, and calculate the electric field received at the two stations for each path through one or more images (from the resulting geometric delays, Doppler frequency shifts at the center frequency, and magnifications). We then construct the secondary spectra for the field observed at both stations by adding the field amplitudes and phases to the appropriate pixel in the secondary spectrum for each path, essentially creating the secondary spectrum of the electric field due to the sum of all images of the pulsar. We construct the dynamic spectra of the electric fields, , at the two stations through a 2-D inverse FFT of each secondary spectrum.

Up until this point, we have neglected the frequency dependence of . In order to include the dependence of the , we inverse Fourier transform the dynamic spectra along the time axis and scale the Fourier frequencies by , where is the frequency of each channel and is a reference frequency, chosen to be the frequency of the lowest channel. We then perform an inverse FFT to return to the dynamic spectra. At this point, we add Gaussian noise to both the real and imaginary parts of the dynamic spectra, with a mean of 0 and a standard deviation 1/10 of the standard deviation of the amplitude of the scintillation pattern. We the construct the intensities, and , and the visibilities, . From these, we then calculate the visibility secondary cross-spectrum, , and the intensity cross secondary spectrum, .

Here we show two different simulations of scattering systems with two screens: In Simulation A, the two screens have different curvatures, so that they can be distinguished by eye in the secondary spectrum, while in Simulation B, the two screens have similar curvatures. In both cases, the screen closer to the observer does not resolve the screen closer to the pulsar - all lines of sight through the screen closer to the observer are assumed to originate from the pulsar without being scattered by the intervening screen. The parameters chosen for both simulations are given in Table 5, and the secondary spectra are shown for Simulations A and B in Figs. 11 and 12 respectively.

We then apply the combined analysis of the visibilities and intensities described in Section 3 and applied to PSR B0834+06 in Section 4.2. The results for calculated at pixels with signal-to-noise greater than 13.7 in the intensity cross secondary spectrum are shown for all three simulations in Fig. 8. We use the subset of points where mas and with signal-to-noise greater than 13.7 in the amplitudes of the secondary spectrum to construct histograms of values, shown in Fig. 8. By fitting a Gaussian curve to the histograms, we can attempt to recover the input values of for the screens. For Simulation A, we recover (for a fit to the histogram between and ) and (for a fit to the histogram between and ). The input values were and . For Simulation B, although there are two noticeable clusters of in the secondary spectrum, only the screen closer to the pulsar is apparent in the histogram. Fitting a Gaussian curve between and , we find . However, with the obvious presence of two screens in the secondary spectrum, one could pick out features belonging to each screen in order to improve the measurements of both distances.

## 6 Conclusions

In this work, we present a novel technique for disentangling features in the secondary spectra of pulsars from different interstellar screens by combining the secondary cross-spectra, constructed from the visibilities, with the cross secondary-spectra, constructed by correlating the dynamic spectra from the intensities at each station. This technique allows one to determine the location on the sky of the two images that are interfering to produce power at each point in the secondary spectrum. When combined with the Doppler shift and delay of these points, one can determine the distances and velocities of the structures in the interstellar medium that are lensing the pulsar radiation, even if multiple structures are contributing and without making the assumption of anisotropic scattering screens. One particularly nice feature of this analysis is that it allows one to easily isolate features that are due to interference of a lensed image with the unlensed image of the pulsar, as at these points one of the angles is very close to zero. This makes it very simple to reconstruct the scattered image of the pulsar on the sky. This technique is simple to implement in VLBI campaigns, as many software correlators, including DiFX and SFXC, allow to user to choose to save the auto-correlations in addition to the cross-correlations, and powerful for reconciling the parabolic arcs observed in the secondary spectra of nearby pulsars with the less organized structure seen in the secondary spectra of others. In this work, we showed the application of this technique to PSR B0834+06, a system with two distinct scattering screens, but the true strength of this technique comes when screens are difficult to distinguish in the interstellar medium. One example is PSR B0329+54, where multiple screens are seen at high frequencies (Putney & Stinebring, 2006), but the emission becomes much messier at lower frequencies (Gwinn et al., 2016; Popov et al., 2017). By applying this novel analysis technique over a large frequency range, and comparing with simulations of scattering by multiple screens, we can start to tease apart the contributions to scattering at different frequencies.

In addition, we expand upon the approach of Galt & Lyne (1972), who measured the delay in the scintillation pattern of B0329+54 between two stations using the intensities recorded at the two stations over 24 hours, and used this to determine the velocity of the scintillation pattern and orientation of the scattering. We show that by working in - space, rather than - space, the delay can be measured to much greater precision, allowing a more precise measurement of the speed of the scintillation pattern, , and the orientation of the scattering screen, , especially for systems with large . Amaral et al. (in prep.) apply this technique to combine Algonquin Radio Observatory (ARO) and Dominion Radio Astrophysical Observatory (DRAO) observations of PSR B1133+16 and measure the distance to one of the intervening scattering screens, using Earth rotation to achieve multiple baselines, while Syed et al. (in prep.) apply this technique to a VLBI observation that was not correlated in order to confirm the distance to the scattering screen, and combine this with variations in the phase of the scintillation pattern across the pulse profile in order to place a physical separation on the two main pulse components of B1133+16. This technique allows measurements of the distances to scattering screens in systems that are dominated by a single screen without requiring the recording, storing, and correlating of raw voltages, reducing the computational power and disk space required for these types of analyses. As a result, it allows studies of the geometry of pulsar scattering in parabolic-arc systems to make use of the wide-bandwidth receivers being installed or recently installed at many radio observatories as well as the pulsar backends at these observatories. With wider-bandwidth and higher bit observations, we will be better able to look at the frequency evolution of structures in the secondary spectra, and to consider the evolution of bright images, which may be degraded by the 2-bit sampling used in most VLBI backends. This frequency and temporal evolution is vital to testing predictions of pulsar scintillation based on models of the scattering plasma.

We anticipate that by incorporating these two techniques, to disentangle multiple scattering screens and to measure distance to screens in single-screen systems using only auto-correlations, the study of pulsar scattering systems and their geometries will expand from the small-number of sources studied today to a large set that engenders a better understanding of pulsar scattering.

## Acknowledgements

We thank Marten van Kerkwijk, Robert Main and Franz Kirsten for helpful discussions. We are grateful for the very useful suggestions Robert Main, Franz Kirsten, J-P Macquart, Marten van Kerkwijk and Dongzi Li provided on early drafts. Research in Canada is funded by NSERC. The Dunlap Institute for Astronomy and Astrophysics is funded through an endowment established by the David Dunlap family and the University of Toronto. The Arecibo Observatory is operated by SRI International under a cooperative agreement with the National Science Foundation (AST-1100968), and in alliance with Ana G. Méndez-Universidad Metropolitana, and the Universities Space Research Association. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

## References

- Armstrong et al. (1995) Armstrong J. W., Rickett B. J., Spangler S. R., 1995, ApJ, 443, 209
- Backer (1975) Backer D. C., 1975, A&A, 43, 395
- Bannister et al. (2016) Bannister K. W., Stevens J., Tuntsov A. V., Walker M. A., Johnston S., Reynolds C., Bignall H., 2016, Science, 351, 354
- Bhat et al. (2016) Bhat N. D. R., Ord S. M., Tremblay S. E., McSweeney S. J., Tingay S. J., 2016, ApJ, 818, 86
- Brisken et al. (2010) Brisken W. F., Macquart J.-P., Gao J. J., Rickett B. J., Coles W. A., Deller A. T., Tingay S. J., West C. J., 2010, ApJ, 708, 232
- Clegg et al. (1998) Clegg A. W., Fey A. L., Lazio T. J. W., 1998, ApJ, 496, 253
- Cordes & Lazio (2002) Cordes J. M., Lazio T. J. W., 2002, preprint (astro-ph/0207156)
- Cordes & Lazio (2003) Cordes J. M., Lazio T. J. W., 2003, preprint (astro-ph/0301598)
- Cordes et al. (2006) Cordes J. M., Rickett B. J., Stinebring D. R., Coles W. A., 2006, ApJ, 637, 346
- Deller et al. (2007) Deller A. T., Tingay S. J., Bailes M., West C., 2007, PASP, 119, 318
- Dong et al. (2018) Dong L., Petropoulou M., Giannios D., 2018, MNRAS, 481, 2685
- Er & Rogers (2018) Er X., Rogers A., 2018, MNRAS, 475, 867
- Fadeev et al. (2018) Fadeev E. N., et al., 2018, preprint (arXiv/1801.06099)
- Galt & Lyne (1972) Galt J., Lyne A. G., 1972, MNRAS, 158, 281
- Gupta et al. (1999) Gupta Y., Bhat N. D. R., Rao A. P., 1999, ApJ, 520, 173
- Gwinn et al. (2000) Gwinn C. R., et al., 2000, ApJ, 531, 902
- Gwinn et al. (2011) Gwinn C. R., Johnson M. D., Smirnova T. V., Stinebring D. R., 2011, ApJ, 733, 52
- Gwinn et al. (2012) Gwinn C. R., et al., 2012, ApJ, 758, 7
- Gwinn et al. (2016) Gwinn C. R., et al., 2016, ApJ, 822, 96
- Hill et al. (2003) Hill A. S., Stinebring D. R., Barnor H. A., Berwick D. E., Webber A. B., 2003, ApJ, 599, 457
- Hill et al. (2005) Hill A. S., Stinebring D. R., Asplund C. T., Berwick D. E., Everett W. B., Hinkel N. R., 2005, ApJ, 619, L171
- Johnson et al. (2012) Johnson M. D., Gwinn C. R., Demorest P., 2012, ApJ, 758, 8
- Keane et al. (2011) Keane E. F., Kramer M., Lyne A. G., Stappers B. W., McLaughlin M. A., 2011, MNRAS, 415, 3065
- Keimpema et al. (2015) Keimpema A., et al., 2015, Experimental Astronomy, 39, 259
- Lee & Jokipii (1976) Lee L. C., Jokipii J. R., 1976, ApJ, 206, 735
- Liu et al. (2016) Liu S., Pen U.-L., Macquart J.-P., Brisken W., Deller A., 2016, MNRAS, 458, 1289
- Main et al. (2017) Main R., van Kerkwijk M. H., Pen U.-L., Rudnitskii A. G., Popov M. V., Soglasnov V. A., Lyutikov M., 2017, preprint (arXiv/1709.09179),
- Pen et al. (2014) Pen U.-L., Macquart J.-P., Deller A. T., Brisken W., 2014, MNRAS, 440, L36
- Popov et al. (2016) Popov M. V., et al., 2016, Astronomy Reports, 60, 792
- Popov et al. (2017) Popov M. V., et al., 2017, MNRAS, 465, 978
- Putney & Stinebring (2006) Putney M. L., Stinebring D. R., 2006, Chinese Journal of Astronomy and Astrophysics Supplement, 6, 233
- Rankin & Wright (2007) Rankin J. M., Wright G. A. E., 2007, MNRAS, 379, 507
- Shishov et al. (2017) Shishov V. I., Smirnova T. V., Gwinn C. R., Andrianov A. S., Popov M. V., Rudnitskiy A. G., Soglasnov V. A., 2017, MNRAS, 468, 3709
- Simard & Pen (2018) Simard D., Pen U.-L., 2018, MNRAS, 478, 983
- Smirnova et al. (1996) Smirnova T. V., Shishov V. I., Malofeev V. M., 1996, ApJ, 462, 289
- Smirnova et al. (2014) Smirnova T. V., et al., 2014, ApJ, 786, 115
- Stinebring (2007) Stinebring D., 2007, Astronomical and Astrophysical Transactions, 26, 517
- Stinebring et al. (2001) Stinebring D. R., McLaughlin M. A., Cordes J. M., Becker K. M., Goodman J. E. E., Kramer M. A., Sheckard J. L., Smith C. T., 2001, ApJ, 549, L97
- Stinebring et al. (2003) Stinebring D. R., Hill A. S., McLaughlin M. A., Becker K. M., Cordes J. M., Kramer M., 2003, in Bailes M., Nice D. J., Thorsett S. E., eds, Astronomical Society of the Pacific Conference Series Vol. 302, Radio Pulsars. p. 263
- Tuntsov et al. (2016) Tuntsov A. V., Walker M. A., Koopmans L. V. E., Bannister K. W., Stevens J., Johnston S., Reynolds C., Bignall H. E., 2016, ApJ, 817, 176
- Walker et al. (2004) Walker M. A., Melrose D. B., Stinebring D. R., Zhang C. M., 2004, MNRAS, 354, 43
- Xu & Zhang (2017) Xu S., Zhang B., 2017, ApJ, 835, 2