Introduction

The propagation of pulses from radio pulsars is affected by temporal and spatial variations in integrated electron column number density and by small-scale inhomogeneities in the ionized plasma of the interstellar medium (ISM; e.g., Armstrong et al. 1995). All ISM effects scale rapidly with observing frequency as ν -γ (with γ = 2 and γ ∼ 4 for dispersion and scattering/scintillation, respectively; e.g., Cordes et al. 1985;Goldreich & Sridhar 2006), and radio astronomy of pulsars has allowed for the number density of the ionized plasma of the ISM to be directly probed and turbulence to be inferred along different lines of sight (e.g., Armstrong et al. 1981Armstrong et al. , 1995;;Levin et al. 2016;Reardon et al. 2020;Alam et al. 2021;Turner et al. 2021).

Time-varying inhomogeneities in the line-of-sight electron column give rise to two important effects: scattering and scintillation. Multipath scattering causes temporal broadening of the pulses and angular broadening of the image, analogous to a convolution of the intrinsic pulse shape with an asymmetric broadening function (e.g., described in detail in Williamson 1972Williamson , 1973Williamson , 1974; and see recent work by Geyer &Karastergiou 2016 andMcKee et al. 2018), which is canonically assumed to be exponential, although a Kolmogorov medium with a small inner scale will depart from this (see, e.g., Lee & Jokipii 1975;Lambert & Rickett 2000).

A phenomenon related to scattering is scintillation; since the pulsar has an angular size that is small enough for it to be considered a point source, the varying inhomogeneities along the line of sight cause the scattered rays to interfere with each other constructively and destructively, giving rise to an interference pattern. As the Earth moves across this interference pattern, the observed flux of the pulsar is seen to vary as a function of time and observing frequency. Studies of the scintillation properties of pulsars have proved to be a very sensitive tool for probing the ISM, studying binary pulsars (Lyne 1984; and see recent work by Reardon et al. 2020), and identifying structures within the ISM. Stinebring et al. (2001) were the first to demonstrate that a two-dimensional Fourier transform of a dynamic spectrum produces scintillation arcs in the "secondary spectrum", which are interpreted as arising from discrete scattering screens along the line of sight to the pulsar. The curvaturesfoot_0 of these arcs are highly dependent on the velocity as well as the geometry of scattering screens, and studying their variation with time has allowed very precise measurements of scattering regions (e.g., Simard 2019;Reardon et al. 2020).

The origin of these scattering screens has been variously attributed to known sources, including circumstellar structures in the vicinity of hot stars (Walker et al. 2017), H II regions along the line of sight to the pulsar (Stinebring et al. 2000), and expanding supernova remnants (e.g., Lambert & Rickett 2000). However, it has not always been possible to associate pulsars with any known structures making the origin of some scattering screens mysterious. In these cases, the inhomogeneities have instead been explained by structures and boundaries within the Local Bubble (Bhat et al. 1998), corrugated current sheets, coincidentally angled close to the line of sight to the pulsar (Pen & Levin 2014), and thin but highly elongated lenses which lead to an apparent overdensity due to their projection along the line of sight (Romani et al. 1987). Such explanations are particularly necessary in the cases where multiple scattering screens are identified within a small region between the Earth and the pulsar (see Putney & Stinebring 2006). Observations of several pulsars have shown evidence of highly anisotropic scattering, where the scattered images of the pulsars are highly elongated along particular axes (Walker et al. 2004;Brisken et al. 2010;Stinebring et al. 2019;Rickett et al. 2021). Such anisotropies in scattering may result from anisotropic plasma structure in the ISM, e.g., corrugated plasma sheets viewed from a grazing angle (Pen & Levin 2014), or "noodles" of plasma which are stabilized by magnetic reconnection (Gwinn 2019).

A notable pulsar where multiple scintillation arcs have been detected is PSR B1133 + 16 (known as J1136+1551 in the J2000 equinox; see Cordes et al. 2006;Putney & Stinebring 2006). We summarize some of this pulsar's important astrometric parameters in Table 1. This is a nearby pulsar (372 ± 3 pc; 1σ uncertainties are used throughout), with a very high transverse velocity of -+ 659.7 4.5 4.2 km s -1 , where both quantities have been measured using very-long-baseline interferometry (VLBI; Deller et al. 2019). 5 This combination of the pulsar's close distance and high speed leads to rapid variations in its observed scintillation, which has been studied extensively (e.g., Trang & Rickett 2007). PSR B1133 + 16 is one of the brightest known pulsars (20 ± 10 mJy at 1382 MHz; Jankowski et al. 2018), which makes it an excellent target for studying the low signal-to-noise ratio (S/N) secondary-spectra parabolas. PSR B1133 + 16 is one of a few pulsars for which the distances to the scattering screens have been estimated (Putney & Stinebring 2006). We note, however, that the distances to the PSR B1133 + 16 screens reported by Putney & Stinebring (2006) assumes a screen orientation such that the major axis of the anisotropic image is aligned with the decl. axis (see Section 4.3) and therefore represent a lower distance limit. Precise estimation of the screen distance requires this angle to be known, either through direct VLBI measurements, as in Brisken et al. (2010), or by modeling the contribution from the Earth's orbit, which we present in detail in Section 4.3.

The remainder of this paper has the following structure. In Section 2, we outline the relevant theory for describing scintillation variations and their applications in isolating scattering-screen properties. In Section 3, we describe our observations and data-reduction techniques. We present the results of our analysis in Section 4 and discuss the implications of our findings in Section 5. Finally, we make concluding remarks in Section 6.

Scintillation Theory

We follow the theoretical treatment of thin-screen scattering as presented in, e.g., Walker et al. (2004) and Cordes et al. (2006). The pulsed flux from pulsars is affected by interstellar scintillation such that the observed intensity, I, varies as a function of time, t, and observing frequency, ν, i.e., I(ν, t), known as the dynamic spectrum. A two-dimensional Fourier transform (represented by a tilde) recasts variations in observed intensity in terms of time delay, τ, and Doppler frequency, f D , and we sum over pairs of interfering images (represented by j, k) with magnification, μ, relative to the unlensed image:



The above equation defines a secondary spectrum, the twodimensional power spectrum of the dynamic spectrum. We characterize scattering screens by their effective distance and velocity vector, D eff and V eff , respectively. The effective distance of the scattering screen, D eff , is defined in terms of its fractional distance, s, from the pulsar to the observer (i.e., a small value of s indicates that a scattering screen is close to the pulsar),

where

, 3 screen psr i.e., 0 < s < 1. The effective velocity of the interference pattern in the plane of the observer is a combination of the pulsar, Earth, and scattering-screen velocities (Cordes & Rickett 1998;Gupta et al. 1994):

where all of the velocities above refer to the components perpendicular to the line of sight (as opposed to the full threedimensional velocities). When the scattering screen is highly anisotropic, i.e., when it is composed of a group of parallel and highly elongated plasma lenses, the resulting image of the pulsar is also highly elongated, with its major axis perpendicular to the lenses on the scattering screen. In this scenario, only the effective velocity parallel to the image can be extracted (see Equation (6)).

The effective distance and velocity allow us to define the axes of the secondary spectrum in terms of the temporal and spectral delays, given by Note. All parameters are measured using VLBI, and are presented in Deller et al. (2019). Note that this reference lists an uncertainty of zero for the pulsar distance, and we have instead calculated the distance from the more precise parallax values listed on the PSRPI website. and

respectively, where λ is the observing wavelength, and θ describes the observed angle of a scattered image. In highly anisotropic scattering, the interference between the scattered images and the line-of-sight image of the pulsar leads to parabolic arcs in the secondary spectrum, described by t h = f D 2 (Stinebring et al. 2001), with curvatures given by

, and α is the angle between the effective velocity and the major axis of the image. The parabola curvature is seen to vary over time due to the proper motion of the pulsar, the Earth's orbital motion, and the movement of the scattering screen (and the pulsar orbit, when in a binary system).

We note that although the component of the pulsar velocity parallel to the line of sight does not contribute to the scintillation time it does, of course, change the distance to the pulsar. In the case of PSR B1133 + 16, this change is negligible: a parallel velocity of ∼100 km s -1 leads to a change in pulsar distance of ∼0.0035 pc (∼700 au) over our 34 yr data set, and we use the assumption that the screen properties have not changed significantly in that time. However, in the case where the pulsar is very fast moving, the length of the data set is very long, and a screen is very close to the pulsar, the parallel component of the velocity becomes significant. A pulsar moving relative to a fixed screen experiences a change in screen placement over time of

for a pulsar at an initial distance d , psr 0 with a parallel velocity component of V ∥ . In the case where the distance between the screen and the moving pulsar remains constant (e.g., the bow shock nebula associated with PSR J0437-4715, where the screen moves with the pulsar; Bell et al. 1993), the change in screen placement over time becomes . This assumes that other properties of the screen remain constant, which may not be the case. For example, significant differences in the distance of the shock contour in the Guitar Nebula are observed, due to changes in the ambient density (Chatterjee & Cordes 2004;Ocker et al. 2021), which would further complicate modeling of the screen.

Observations and Data Reduction

Our data set is composed of observations taken over a 34 yr period using the Arecibo telescope. Between 1980 and 2015, observations were taken sporadically, with gaps of many years in the data set (the longest being ∼10 yr). Observations were taken using the Mock spectrometer backendfoot_3 , at Nyquist sampling rates of typically 6.4-38.4 μs, and with data recorded using 10 s subintegrations. Using the approach outlined in Hill et al. (2003), spectra were determined by recording the flux densities, S, of the on-pulse and off-pulse phase regions in each frequency channel, ν, and in each subintegration, t i , and weighting by the off-pulse average from the entire observation, i.e.,

In addition to the historical data, our data set includes observations taken as part of a dense campaign, carried out between 2015 January and July. During this time, 21 multifrequency observations were taken with an approximately weekly cadence. Observations were typically 45 minutes in length, primarily at a center frequency of 432 MHz with 26.7 MHz of bandwidth. Immediately after most of these observations, the pulsar was observed again, using a different frequency: either 327 MHz with 53.3 MHz bandwidth, or 1450 MHz (160 MHz bandwidth). The number of spectral channels used in the 327, 430, and 1450 MHz observations were, respectively, 4096, 2048, and 2048, resulting in a Nyquist sampling rate of 38.4 μs for the two lower frequencies, and 6.4 μs for 1450 MHz. Our observational configurations are summarized in Table 2.

The high sensitivity of the Arecibo telescope, together with the fine spectral resolution of our data set, has enabled highquality dynamic spectra to be generated from all of our observations. An example of a dynamic spectrum from an observation at 1450 MHz is shown in the upper panel of Figure 1. All of our dynamic spectra were converted to secondary spectra as the squared modulus of the twodimensional Fourier transform (described in Section 2).

Results and Analysis

We observe scintillation arcs in the secondary spectra from all of our observations and see that many observations show multiple simultaneous arcs. These arcs are mostly characterized as very sharp and well-defined parabolas, which indicates that they arise from highly anisotropic thin screens (Figure 1). In a small number of observations, especially those below 700 MHz, we see very clear inverted arclets in the secondary spectra (Figure 2), further confirming that the scattering in this pulsar is highly anisotropic.

Arc Curvature Measurement

The high-quality secondary spectra from our observations allow for precise curvature measurements for the scintillation arcs. Our analysis of the secondary spectra arc curvatures made use of the SCINTOOLSfoot_4 Python module (described in Reardon 2018Reardon , 2020)). This tool set has recently been used in measuring the orbital parameters of PSR J0437-4715 via the variations in arc curvatures modulated by the pulsar orbit and the Earth orbit (Reardon et al. 2020), providing an independent analysis to those of pulsar timing measurements (Perera et al. 2019). It has recently been suggested that particular scatteringscreen morphologies will cause the apex of the parabolas in the secondary spectrum to be offset from the origin (Shi 2021). We do not see any evidence for this in our data set, and we assume that the apex of all parabolas is at the origin.

The dynamic spectra were resampled in terms of wavenumbers rather than frequencies (i.e., k = 2π/λ; see Fallows et al. 2014) before secondary spectra were computed, to allow arc curvatures from our multifrequency observations to be directly compared in a consistent manner. This resampling also eliminates the wavelength dependence in the scintillation arc curvature, which makes the arcs sharper over a finite bandwidth. The curvature for a scintillation arc in the resampled spectra therefore takes the form

where η is the curvature of the scintillation arc without resampling, defined in Equation (7). The scintillation arc Figure 1. An example of a PSR B1133 + 16 dynamic spectrum (above, linear grayscale), taken at a center frequency of 1450 MHz. The observation was made on MJD 57098 (day number 119 in Figure 5), with darker colors corresponding to greater flux density. One subintegration containing radio frequency interference has been removed. The associated secondary spectrum (below, using a logarithmic grayscale, and with the center fringe frequency removed), the squared modulus of a two-dimensional Fourier transform of the dynamic spectrum, displays at least three separate arcs (labeled C, D, and E, in order of increasing curvature), for which we have measured curvatures of 0.0041(5), 0.0054(3), and 0.013(2) s 3 (scaled to a frequency of 1400 MHz). The arcs correspond to persistent discrete scattering screens, which we observe throughout the length of our data set. The sharply defined arcs indicate the presence of weak scintillation in this observation. Note that there is a hint of a lower-curvature arc, which is consistent with the Arc B curvature, but the S/N is too low for the curvature to be precisely measured.

curvatures were then measured by an application of a generalized Hough transform (see Bhat et al. 2016). This was done by evaluating a path integral along different trial curvatures, producing a distribution of the integrated power as a function of trial curvature. In the resulting distribution the optimum arc curvature and its uncertainty were determined by fitting parabolas to the peaks in integrated power (see the example given, shown in Appendix A, Figures A1 and A2). Since many of the PSR B1133 + 16 secondary spectra display multiple scintillation arcs (e.g., Figure 1), the above process was repeated over several ranges of η to avoid identifying only the brightest arc at the global maximum of the integrated power. Final η values were checked by overplotting their resulting parabolas on the secondary spectra, and η values with large uncertainties, e.g., due to low S/N, were excluded. In some cases, anomalously small uncertainties in the curvature measurements were recorded, which appeared to be significantly underestimated. To reduce the impact of these measurements on the model fitting (see Section 4.3), we have scaled the uncertainties in quadrature for each arc separately. After a first iteration of model fitting, a factor E quad was found that produced a reduced χ 2 = 1 when scaling the uncertainties as

scaled initial 2 quad 2 1 2 . We note that several observations appear to have faint scintillation arcs for which we have not been able to precisely measure curvatures and so were not included in our analysis. As the ranges of curvatures that we have measured are consistent for decades (Figure 3, and see below), this implies that all of the arcs are likely to be present in all observations, albeit at a S/N that is too low for us to reliably detect. Taking this to be true, it is not unreasonable to expect that there may be additional arcs that we are not sensitive to, either due to low S/ N or being at curvatures that are too high or low to easily distinguish in our observations. A targeted campaign with highly sensitive, next-generation telescopes would potentially uncover many more screens than the six we have presented in this work.

Arc Curvatures Associated with Specific Screens

When plotting our measured arc curvatures as a function of time we see that several distinct populations emerge. By noting which arcs appear in secondary spectra with multiple arcs (e.g., Figure 1), we are able to unambiguously assign curvatures to particular populations. We show this in Appendix B (Figure B3), and find that a minimum of six populations are required, which we refer to as Arcs A-F in ascending order of curvature (Figure 3). We note that, with the exception of Arcs A and E, all pairs of arcs have been detected in the same observation (see Table B3). Remaining curvatures were assigned to populations based on these ranges, and the annual variation of the curvatures (discussed in detail in Section 4.3). Arc C is almost always the brightest arc, and in instances where only one arc is detectable this is always Arc C. There is a large gap between detections of the highest curvatures, but we assume the simplest case where they are all part of the same population (Arc F).

As we outlined in Section 2, each range of curvatures arises from an individual scattering screen; we have therefore demonstrated there are at least six distinct screens along the 372 ± 3 pc line of sight to PSR B1133 + 16.

Annual Variation of Arc Curvatures and Determination of Screen Astrometry

As mentioned in Section 2, the orbital motion of the Earth, together with the proper motion of the pulsar and the scattering screen, leads to annual variations in observed arc curvatures (under the assumption that the scattering geometry remains static throughout a 1 yr period). As the proper motion and distance of PSR B1133 + 16 are known to high precision, this allows for the astrometric parameters of the scattering screens to be extracted when using the contribution of the Earth's orbit to disentangle the velocities, if the annual signal is well sampled with curvature measurements of sufficiently high precision.

When folding all curvature measurements into a single yearlong period, the annual signals for Arcs C-F immediately become apparent. For each scintillation arc, we fit the curvature measurements with three parameters: the fractional distance, s, the transverse velocity, V screen , and the angle between the major axis of the anisotropic image and the decl. axis on the associated scattering screen, α s . As we noted in Section 2, only the effective velocity parallel to the major axis can be probed in highly anisotropic scattering. Hence, instead of fitting a vector velocity for the scattering screen, we fit a single scalar along the major axis of the image. Practically, we fit the curvature measurements using the following functional form:

where

is the known part of the effective velocity, t is the time of the year, and α v (t) is the angle between A(t) and the decl. axis (so α s -α v (t) is the angle between A(t) and the major axis of the image), which varies due to changing Earth motion. The screen velocity is canonically expected to be comparable to the speed of sound of interstellar plasma, which is ∼10 km s -1 (Goldreich & Sridhar 1995). However, this is only true in the absence of external factors driving ISM turbulence. Some examples of this are evaporated material shed from objects in the vicinity of hot stars (Walker et al. 2017), expanding supernova remnants (Lambert & Rickett 2000), stellar winds, and high-velocity (i.e., supersonic) neutron stars moving in the local ISM and producing shocks. Additionally, material at high galactic latitudes would have an acceleration contribution from their orbital motion in the galactic potential. This means that although screen velocities consistent with the speed of sound might be expected, there are many scenarios that can explain departures from it.

Applying Equation (12), we extract the unknown screen parameters: s, α s , and V screen . We calculate the pulsar velocity in R.A. (α) and decl. (δ) using the coordinates, proper motion, and distance obtained from VLBI measurements (presented in Deller et al. 2019), as =a -+ V 129.404 0.018 0.050 km s -1 and = d -+ V 642.89 0.10 0.13 km s -1 . We made use of the ASTROPYfoot_5 package (Astropy Collaboration et al. 2013, 2018) to calculate the observer motion in the pulsar reference frame, and relative to the solar system barycenter using the coordinates of the Arecibo Observatory.

The screen parameters in Equation (12) are highly covariant, particularly when the amplitude of the periodic curvature variations is smaller than the measurement uncertainties (as is the case for Arcs A and B). We therefore explore the parameter space of the models using a Markov Chain Monte Carlo (MCMC) technique. As mentioned above, there are values of α s every 180°that produce solutions to Equation (12). We chose priors for the model parameters such that only one of the α s solutions was being tested, which we arbitrarily chose to be the one closest to 0°. The resulting corner plots from this MCMC analysis are shown in Figure 4, and we display the models together with the annual curvatures in Figure 5. The parameters that we have extracted from the fits to the model (Equation ( 12)) are listed in Table 3.

We find that our model provides a good description of the observed curvature variations in Arcs C-F, due to the sufficiently high-precision detections of arcs associated with these screens, compared to the amplitude of their variations in curvature. As mentioned above, we do not attempt to fit a model to Arc A, owing to the few detections we have made of curvatures associated with its screen. Although we have modeled the periodic variations in Arc B, we were not able to find models with a consistent phase, indicating that the amplitude of the annual variation is of comparable or smaller size to our measurement uncertainties.

Since Arc B is well sampled throughout the Earth's orbit, it is visible in observations separated by decades (Figure 3), and its range of measured curvatures is relatively low (∼0.0014-0.0035 s 3 , at a reference frequency of 1400 MHz), we can be confident that the curvatures we have assigned to this population are indeed associated with a single scattering screen. The annual variations are clearly detectable in Arcs D and E, but the magnitude of the measurement uncertainties prevents precise parameters from being extracted. At our measurement precision, we are only able to very vaguely place Screens B, D, and E at a fractional distance of ∼0.5.

In general, screens closer to the Earth give rise to higher amplitudes in the annual variations of arc curvatures, which implies that Screen B is much closer to PSR B1133 + 16 than the upper limit of our distance uncertainty. This also suggests that, at our measurement precision, precise astrometry of Screens A and B (and any that are undiscovered) would require much higher S/N and greater delay and fringe frequency Figure 3. Variations in arc curvature, η, as a function of time (scaled to a reference frequency of 1400 MHz), for our entire data set. Multiple arcs are resolvable, corresponding to specific ranges of η which we attribute to specific scattering screens, and are color-coded. The detection of these arcs over a long baseline demonstrates that the underlying screens persist for decades, over the length of our entire data set. resolutions (i.e., longer observations with higher frequency resolution) than is available with our data set, or coherent phase retrieval (Walker et al. 2008;Baker et al. 2021). Observations of PSR B1133 + 16 with the new generation of high-gain telescopes will certainly improve the measurement precision of Arcs B, D, and E to a level where precise astrometric values of the underlying screens can be extracted.

We find precise solutions for the parameters of the screens associated with Arcs C and F (Figure 4). We see from Figure 5 that the models offer an excellent description of our measured curvatures for these two arcs. Although we have relatively few detections of Arc F, the large amplitude of its annual variation, due to its close proximity to the Earth, allows for more precise parameters to be extracted from fits to the data. We note that the amplitude of the annual variation predicted by the best-fit models varies by approximately one order of magnitude (Figure 5), indicating that the precision we have achieved here could be drastically improved, given further measurements throughout the Earth's orbit.

Corner plots showing the results of Markov Chain Monte Carlo exploration of the parameter space for our screen models. In each case, 300,000 models were generated using flat priors, in a range of α s corresponding to one specific solution (see text). Dashed lines represent the median, 0.16, and 0.84 quantiles, while solid lines indicate the distribution mean. While precise screen parameters have been determined for Screens C and F, the Screen B, D, and E parameters are not well constrained.

Discussion

Screen Properties

We have precisely measured the parameters of Screens C and F (Table 3). The distances we measure are in general agreement with those presented in Putney & Stinebring (2006), although the values in their work used single observations and assumed a screen angle α s = 0°, i.e., they did not use the knowledge of the orbit to break degeneracies between parameters. Therefore, their reported values represent lower limits on s. Since Putney & Stinebring (2006) did not measure temporal variations in arc curvatures, they do not present measurements of the screen angles or velocities.

Although the solar system is moving at ∼26 km s -1 with respect to the local ISM (e.g., Wood et al. 2015), and so velocities in the range of approximately ± 30 km s -1 would be consistent with the ISM speed of sound ( ∼ 10 km s -1 ), there appears to be a departure from this in the case of Screen C. While the precision of the screen velocity measurement is quite low, there is a strong preference for a velocity far from that expected from thermal motion of the ISM, which would indicate that the scattering medium has been energized by an external source.

The Local Interstellar Medium

The close distance of PSR B1133 + 16 combined with the persistence of several discrete scattering screens over decades indicate that the local ISM contains several distinct regions, which we illustrate in Figure 6. It perhaps seems surprising that there are at least six distinct ISM regions within only 372 ± 3 pc. We do not find evidence for any known structures in the region around PSR B1133 + 16, using the Aladin Sky Atlasfoot_6 (Boch & Fernique 2014), such as the H II region that has been identified as the dominant source of the observed scattering in PSR B1642-03 (Gupta et al. 1994) and PSR J1643-1224 (Mall et al. 2022), nor do we find any reported nearby hot stars.

The lack of associated objects along the line of sight lends weight to the predictions by Romani et al. (1987) and Pen & Levin (2014), that thin elongated lenses and corrugated plasma sheets, respectively, formed within the ISM give rise to scattering. However, as mentioned above, we measure a velocity for Screen C that is significantly different from the speed of sound in the ISM, which may exclude this explanation. PSR B1133 + 16 has traveled ∼0.0035 pc throughout the duration of our data set, which places a lower limit on the size of the scattering regions we probe.

Remaining explanations are ones where density fluctuations within the local ISM are responsible. The Earth is located in the Local Bubble: an underdense region of interstellar space driven by the aggregated effects of nearby supernovae over the last ∼10-20 Myr (e.g., Frisch 2007). Maps of the Local Bubble produced by measurements of interstellar absorption Figure 5. Annual variations in arc curvature, η, for the PSR B1133 + 16 arcs, labeled as Arcs A-F in order of ascending curvature. The curvatures have been scaled to a reference frequency of 1400 MHz. Black lines show 24 randomly selected models from the results of our Markov Chain Monte Carlo parameter estimation (Figure 4), with the exception of Arc A, due to the small number of detections of this arc. Note that the Arc F plot uses a logarithmic scale due to the high amplitude of the annual variation, while the other arc plots use a linear scale. Our models provide a good description of the periodic variations in Arcs C-F. Note. Due to the small number of detections, a model is not fit for Screen A.

There are solutions at integer multiples of 180°from the listed α s , for which the velocity alternates between positive and negative values. Here we arbitrarily list the solution with the angle closest to zero. We use the astronomer's convention of the angle being measured north through east. (Lallement et al. 2014) show that the line of sight toward PSR B1133 + 16 cuts through several density boundaries.

Although the distance we measure to Screen C is relatively imprecise, it is consistent with the approximate location of an overdensity within a cavity of the GSH238 + 00 + 09 super bubble (Heiles 1998). Screens B, D, and E could possibly be attributed to this region of the Local Bubble, as well. While it is possible that Screen F also arises from structure within the Local Bubble, it is extremely nearby, at distances below the resolution of current Local Bubble maps, but unlikely to be described by the same explanation as Screens B-E.

Future Prospects

The properties of Screens C and F we have measured here are precise and demonstrate the utility of pulsar scintillation as a means of probing the local ISM. We briefly discuss some ways that this work can be built upon to uncover further information about the local ISM.

There are a number of new facilities that have commenced operations in recent years that will improve the prospects for identifying more pulsars with multiple well-resolved scattering screens, and for making high-precision measurements of arc curvatures. The Canadian Hydrogen Intensity Mapping Experiment (CHIME; detailed in CHIME/Pulsar Collaboration et al. 2021) observes almost every known pulsar in the Northern sky with an approximately daily cadence, in the 400-800 MHz band, making it an ideal instrument for measuring temporal variations in ISM properties. The Five-Hundred-Meter Aperture Spherical Radio Telescope (FAST; detailed in Lu et al. 2019) has recently been used to measure scattering screens in PSR B1929 + 10 and PSR B1842 + 14, which was not possible with the previous generation of radio telescopes (Yao et al. 2020). MeerKAT (Bailes et al. 2020), the precursor to the Square Kilometre Array (SKA), has recently been utilized in making scattering measurements of a small number of pulsars in the Southern sky (Oswald et al. 2021) and will be further improved with the completion of the full SKA (Keane 2018). Independent determination of the screen angle through direct VLBI measurements (in the same way as Brisken et al. 2010) could potentially offer better constraints on the allowed parameter space of the curvature variation models. This would allow a more precise distance to the screens to be determined.

In this work, we do not detect all six arcs in every one of our observations, and in the cases of Arcs A and F, we have very few detections. The results of our analysis tell us that all six arcs are probably present in the secondary spectra. There are likely three reasons for these nondetections:

1. Stochastic variations in the scattering properties of a screen (e.g., as a result of intermittent turbulence; see, e.g., Stinebring et al. 2001;Walker et al. 2004;Cordes et al. 2006) may also result in nondetections at a particular epoch. This is expected in the Pen & Levin (2014) corrugated-sheet model. 2. Some of the arcs are simply too faint to be detected in every observation (we note that many different configurations of receivers and backends have been used throughout our data set, and the presence of radio frequency interference can quickly diminsh the S/N of secondary-spectra arcs). 3. In the cases of Arcs A and F, curvatures at certain parts of the Earth's orbit are too great or small to be identified in our secondary spectra (in Figure 5 we see that the predicted amplitude of the Arc F variation is up to an order of magnitude greater than our highest measured curvature).

Therefore observations with the new generation of highly sensitive telescopes with broadband recording will allow for both higher precision curvatures to be measured and more consistent detections of the fainter arcs. PSR B1133 + 16 is special in that six scintillation arcs are now known in its secondary spectra, more than the other pulsars with known multiple scintillation arcs (Putney & Stinebring 2006). However, there is seemingly nothing unique about the line of sight to PSR B1133 + 16 to which these multiple scattering screens could be attributed. At the time of writing, this pulsar is the 29th brightest (20 ± 10 mJy at 1400 MHz) listed in the Australia Telescope National Facility (ATNF) pulsar catalog 10 (Manchester et al. 2005), which makes it likely that it is the intrinsic brightness of the pulsar (and therefore, the scintillation arcs), and the high sensitivity of the Arecibo Telescope, rather than its line of sight, which gives rise to a large number of detectable scintillation arcs. It then follows that highersensitivity observations of PSR B1133 + 16 will uncover previously unknown scintillation arcs and that observations of other pulsars will uncover multiple arcs. Such a census of this pulsar and others will allow for the distinct regions of the local ISM to be mapped. As well as providing a robust test of the predictions for the origin of scattering regions by Romani et al. (1987) and Pen & Levin (2014), this will provide a new picture of the way the local ISM is organized, on a distance scale far below that of current local ISM maps.

Conclusions

We have used 34 yr of high-sensitivity observations with the Arecibo Telescope to study the scintillation properties of PSR B1133 + 16. By measuring the arc curvatures in the secondary spectrum generated from our observations, we have identified six scintillation arcs that are consistent for decades and which we have attributed to discrete scattering screens. Using the assumption that the properties of the scattering screens have not changed significantly throughout our data set, we have modeled the variations in arc curvatures, finding the locations of five of these six scattering screens, and have achieved subparsec precision in one case. We speculate that structures within the Local Bubble are responsible for the underlying scattering screens, and suggest that some of the screens may be associated with the boundary of a known region of the Local Bubble, although this is unlikely for the screen closest to the Earth.

Table B1 Number of Simultaneous Detections of Screens A-F in our Data Set Screen A B C D E F A L 2 1 2 0 1 B 2 L 27 9 15 6 C 1 2 7 L 9 2 8 1 2 D 2 9 9 L 4 3 E 0 1 5 2 8 4 L 6 F 1 6 1 2 3 6 L