A Cyclic Spectroscopy Scintillation Study of PSR B1937+21
I. Demonstration of Improved Scintillometry

0000-0002-2451-7288]Jacob E. Turner

Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA

0000-0001-8885-6388]Timothy Dolch
Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA
Eureka Scientific, 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA

0000-0002-4049-1882]James M. Cordes
Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

0000-0002-4941-5333]Stella K. Ocker
Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91101, USA
The Observatories of the Carnegie Institution for Science, Pasadena, CA 91101, USA

0000-0002-1797-3277]Daniel R. Stinebring
Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA

0000-0002-2878-1502]Shami Chatterjee
Cornell Center for Astrophysics and Planetary Science and Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

0000-0001-7697-7422]Maura A. McLaughlin
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA

0000-0002-4925-8403]Victoria E. Catlett
Green Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA

0000-0002-4188-6827]Cody Jessup
Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA
Department of Physics, Montana State University, Bozeman, MT 59717, USA

0009-0009-7809-3335]Nathaniel Jones
Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA

0009-0007-2434-2776]Christopher Scheithauer
Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA

§ INTRODUCTION

Pulsar scintillation is a valuable tool for understanding the small scale structure of the ionized interstellar medium (ISM), while also acting as a significant source of noise in pulsar timing efforts. In many cases, gaining a better understanding of the former can lead to better mitigation efforts of the latter. It is also crucial to understand the frequency dependence of propagation through the ISM, as the effects seen in pulsar scintillation evolve with these frequencies.
As a consequence of interactions between pulsar emission and free electrons along the signal's propagation path, the emission will experience multipath propagation and its image will undergo angular broadening, resulting in interference between multiple ray paths at the observer and a broadening of the observed pulsar signal <cit.>. The effects of this interference can be seen in pulsar dynamic spectra, which track the intensity of the observed signal over frequency and time. The widths of bright patches in these spectra, known as scintles, provide insight into the rapidity of these intensity fluctuations and thereforeprovide a means to estimate the level of scattering and the corresponding arrival time delay a given pulse may experience. Additionally, if the broadened pulse can be deconvolved into the intrinsic pulsar signal and the pulse broadening function (PBF), which has also been referred to in the literature as the impulse response function (IRF), although the former is an ensemble average of the latter, of the ISM – the signal imparted by the ISM onto the pulse as a result of their interaction – fitting this PBF to an exponential decay function provides another way to quantify the level of scattering present in an observation. These delays, $\tau_{\rm d}$, are expected to scale as $\nu^{-4.4}$, where $\nu$ is the observing frequency, under the assumptions of thin-screen scattering andthat the ISM can be accurately characterized by Kolmogorov turbulence and that the smaller length scales of the wavenumber spectrum dominate <cit.>. Under the assumption of single-screen scattering, these delays are also inversely proportional to the characteristic scintle width in frequency of a dynamic spectrum, also known as the scintillation bandwidth, $\Delta \nu_{\textrm{d}}$ <cit.>. Consequently, $\tau_{\rm d}$ can be inferred from $\Delta \nu_{\textrm{d}}$, which is expected to scale with observing frequency as $\nu^{4.4}$ under the same assumptions mentioned above.
Additional information can be gained about the structure of the ISM along a given line of sight (LoS) by examining the the 2D Fourier transform of the dynamic spectrum, also known as the secondary spectrum, of an observation, where interference patterns in the dynamic spectrum can give rise to parabolic structures known as scintillation arcs, whose curvature, in conjunction with the pulsar's LoS velocity across the screen, $\textbf{V}_{\textrm{eff}, \perp}$, and distance, can be used to determine the location of the corresponding scattering screen along the LoS <cit.>. Structures within these arcs can also be used to gain insight into the structure of the scattering region <cit.>. Under the assumption that pulsar scattering occurs at a thin screen, the arc curvatures are expected to scale with frequency as $\nu^{-2}$ <cit.>.

Delays from the most highly scattered pulsars are the most pertinent for mitigation in pulsar timing efforts. Additionally, the LoSs through the ISM to these pulsars are often the most interesting to study. However, observations that use pulsar timing array (PTA)-style observing setups, as well as those that use traditional polyphase filterbanks and incoherent dedispersion, typically lack the frequency resolution necessary to see scintillation structures within their dynamic spectra, particularly at lower observing frequencies. This often leaves observers unable to effectively quantify the level of scattering present or resolve scintillation arcs. However, cyclic spectroscopy is a technique that takes advantage of the periodic and amplitude-modulated nature of the pulsar signal to drastically improve the resolution of baseband observations. Cyclic spectroscopy, used for years within various engineering communities, was first introduced to pulsar timing in <cit.>, where the technique made use of the phase information in the voltage data and exploited the fact that the scintillation time spans many pulse periods, allowing the PBF to be estimated. The cyclic spectrum is defined as

\begin{equation}
\label{act_cyc}
{S_E}(\nu, \alpha_k) = \langle E(\nu + \alpha_k/2)E^{*}(\nu - \alpha_k/2)\rangle,
\end{equation}

where $E(\nu)$ is the electric field of the signal and $\alpha_k = k/P$ is the cyclic frequency at the $k^{\textrm{th}}$ harmonic for pulse period $P$ <cit.>. To account for changes in the pulsar signal due to the ISM, over integration times similar to a pulsar's scintillation timescale, we can write the cyclic spectrum as

\begin{equation}
\label{prac_cyc}
{S_E}(\nu, \alpha_k) = \langle H(\nu + \alpha_k/2)H^{*}(\nu - \alpha_k/2)\rangle{S_x}(\nu, \alpha_k),
\end{equation}

where $H(\nu)$ is the transfer function of the ISM, which is the Fourier transform of the ISM PBF, and ${S_x}(\nu, \alpha_k)$ is the Fourier transform of the intrinsic pulse profile <cit.>.
As a result of the additional information made available in the pulsar signal through this technique, while maintaining a given number of pulse phase bins, one can acquire significantly higher frequency resolution than accessible through the use of traditional polyphase filterbank channels, allowing for the emergence of previously unresolved features in periodic and dynamic spectra while maintaining high pulse phase resolution. This resolution improvement has already been used to great effect in millisecond pulsar (MSP) observations at low radio frequencies <cit.> and could lead to pulsar timing procedures where scintles and scintillation arcs (assuming observing times where scintillation timescales are obtainable) are resolved in all epochs. This would allow for scintillation arcs to be studied with preexisting monthly observations of many pulsars without the need to alter the observing parameters crucial to PTA feasibility, and providing many LoSs through the ISM over which to examine arc evolution over frequency and time <cit.>. In particular, these high cadence observations are crucial for tracking arc curvature over time, allowing for tighter constraints on scattering screen distances.
Additionally, this scintle resolution would make it possible to estimate scattering delays at all epochs, provided flux density is sufficient for visible scintles in dynamic spectra. This would (perhaps dramatically) advance the feasibility of correcting for scattering delays in the pulsar timing models used by PTAs. Arguably, the main benefit of CS is its ability to resolve pulse broadening functions, providing a direct approach to mitigating scattering delays across individual scintillation timescales. However, the capability of high-resolution dynamic spectra also allows for the monitoring of scattering delays in highly scattered pulsars, potentially informing red noise modeling efforts to remove these effects in instances when we cannot recover pulse broadening functions. While PTAs monitor these delays over time <cit.>, they do not account for them in timing models. This results in a noise floor of tens to hundreds of nanoseconds at L-band, depending on the pulsar, and even higher at lower observing frequencies <cit.>.
There are attempts within PTAs to indirectly account for the effects of scattering. For example, some PTAs fit $\nu^{-4}$ delays to TOAs <cit.>, but this approach ignores that scaling indices have been observed to vary across both LoSs and epochs. Ignoring this variability could result in improper correction of these scattering effects, thus introducing unnecessary noise into PTA data <cit.>. For highly scattered pulsars, which can have many scintles across an observing band of only a few hundred MHz, the frequency resolution improvement provided by CS could allow for more accurate, frequency-dependent scaling law measurements on individual epochs. This improved accuracy could potentially inform PTA noise modeling by providing priors to red noise analyses that attempt to mitigate delays from interstellar scattering.
The ideal situation when performing frequency-dependent scintillation analyses of pulsar emission propagation through the ISM is to obtain all observations at the various frequencies simultaneously so that the same LoS is being sampled throughout (although, note the small discrepancies in effective column density that can arise due to frequency-dependent scattering volumes; see  for details.). Typically when frequency-dependent analyses are performed, they are limited by the frequency resolution and bandwidth of the observations, as well as the observing cadence between observations at different frequencies. Some of these analyses used a large frequency range, but only 2$-$4 measurements <cit.>, with these measurements sometimes taken on different days <cit.>. In general, this should not affect results provided they were taken within a pulsar's refractive timescale and the effective velocity remains similar over the course of the observations, but nearly simultaneous measurements are preferred. Other methods make simultaneous measurements across many frequencies <cit.>, although all such studies have been performed at lower frequency ranges <cit.>. Studies at L-band have taken the approach of breaking up large bandwidth observations into smaller frequency slices as well, although this appears to have been limited to 200 MHz slices across 800 MHz bands <cit.>.
In this paper we use cyclic spectroscopy to acquire fine frequency resolution across a relatively small bandwidth, allowing for eight simultaneous measurements per epoch for frequency-dependent analyses. In Section <ref> we present the data used for this work. In Section <ref> we discuss the data processing and analyses performed on these data. Section <ref> discusses the results of these analyses and their comparison to theory and previous literature. In Section <ref> we briefly discuss benefits CS can provide for PTAs.Finally, in Section <ref> we discuss our conclusions and the potential for future work that can take full advantage of the cyclic spectroscopy techniques discussed in this work.
While this paper focuses on the benefits of cyclic spectroscopy (CS) in the context of analyses of the ISM in PTA observing setups, as well as providing parameter estimations and scaling law fits of these data, accompanying papers (Dolch et al., in prep, Turner et al., in prep) will focus on a more detailed interpretation of these measurements.

§ DATA

Our observations of PSR B1937+21 (P2627 PI Stinebring) took place on MJDs (and fractional years) 56183 (2012.70), 56198 (2012.74), and 56206 (2013.01). They spanned approximately 2.5 hours each and were taken in baseband mode with the FPGA-based PUPPI spectrometer at the Arecibo Observatory using 200 MHz of bandwidth with 6.25 MHz wide filterbank channels centered at 1373.125 MHz. Our data were recorded on eight 25 MHz-wide banks, which we maintained for our subband analyses. Our raw baseband data were then processed via cyclic spectroscopy using [<http://dspsr.sourceforge.net/>] <cit.> with 1024 pulse phase bins and 1024 cyclic channels per filterbank channel, for a final channel resolution of around 6.1 kHz. We used blocks of data averaged over approximately 28 second intervals for our analyses. An example dynamic spectrum processed from these baseband data using the typical NANOGrav frequency resolution of 1.5 MHz <cit.> compared with the same dynamic spectrum processed with cyclic spectroscopy is shown in Figure <ref>.
The CS technique makes use of the phase coherent treatment of radio signals. However, it may be just as effective to simply treat the pulsar as a continuum source (i.e., non-periodic) and sacrifice time resolution to achieve the desired frequency resolution in accordance with the Gabor limit on time-frequency sampling. In Appendix <ref> we discuss why this is an inferior approach to true CS. In principle one could also achieve their desired frequency resolution with baseband data. However, this would be at the expense of time resolution, requiring two data products to achieve what cyclic spectroscopy can achieve in one, resulting in more efficient processing and more sustainable usage of available data storage.

A 100 MHz slice of the dynamic spectrum on MJD 56198 showing data processed using Fourier spectroscopy that yields 1.5 MHz resolution with the standard frequency resolution available to NANOGrav (1.5 MHz channels, left) and using cyclic spectroscopy (6.1 kHz channels, right), both using 1024 pulse phase bins. Cyclic spectroscopy greatly improves the frequency resolution of the dynamic spectrum while maintaining high time resolution, allowing for high time and frequency resolution within the same data product. The darker bands in the image on the right are the gaps between the filterbanks.

§ ANALYSES

To examine the frequency dependence of scintillation effects in our data, we divided the total band into eight 25 MHz wide subbands, which are wide enough to contain a large number of narrow scintles. This frequency slice size was chosen to allow for a high degree of robustness in our analyses while limiting the effects of both the finite scintle effect and the filterbank gaps.

§.§ Dynamic Spectra

After processing our baseband data using CS to improve frequency resolution, for each observation, we summed the two orthogonal polarizations and created a dynamic spectrum, defined as

\begin{equation}
\label{dynspec}
\end{equation}

where $S$ is the intensity of the pulsar signal at each observing frequency, $\nu$, and time, $t$, ${P}_{{\rm{bandpass}}}$ is the total power of the observation, and ${P}_{{\rm{on}}}$ and ${P}_{{\rm{off}}}$ are the power in all on- and off-pulse components, respectively, at each frequency and time (see Figure <ref>). The on-pulse region was determined by first creating a template from the average profile of the corresponding subband and then performing an overlap integral of $P(\nu,t)$ with that template.

To determine the scintillation bandwidth for a given frequency slice, we first acquired to corresponding dynamic spectrum's 2D autocorrelation function (ACF), from which we retrieved its 1D frequency and time ACFs. In virtually all studies that use an ACF analysis to determine scintillation parameters, these ACFs are then fit with a Gaussian and then the scintillation bandwidth, $\Delta \nu_{\rm d}$, is measured by determining the half width at half maximum of its fit and the scintillation timescale, $\Delta t_{\rm d}$, is measured by determining the half width and $e^{-1}$ of the time ACF fit. However, given that the ISM's PBF, a time domain function, is characterized by a one-sided decaying exponential, it makes more mathematical sense to fit the corresponding frequency domain ACF using a Lorentzian, as the two functions are Fourier pairs. In order to maintain consistency with the majority of existing scintillation studies while also being as mathematically thorough as possible, we performed both Gaussian and Lorentzian fits on the frequency ACFs when determining  scintillation bandwidths. To examine a given dynamic spectrum's characteristic scintle widths in frequency and time, we acquired its autocorrelation functions (ACFs) in both frequency and time using pypulse <cit.> to extract that spectrum's 2D ACF. We then took slices at zero time lag and zero frequency lag to acquire time and frequency ACFs, respectively. As a consequence of the high S/N in our data, noise spikes in both the frequency and time ACFs were insignificant enough that no modification to the ACFs was necessary prior to fitting. To determine the scintillation bandwidth for a given subband, we took its dynamic spectrum and retrieved its 1D frequency ACF. We then fit its frequency ACF with a Lorentzian and measured the scintillation bandwidth, $\Delta \nu_{\rm d}$, by finding the half width at half maximum of its fit. Similarly, to determine the scintillation timescale for a given subband, we took its dynamic spectrum and fit its time ACF with a decaying exponential and measured the scintillation timescale, $\Delta t_{\rm d}$, by finding the half width at $1/e$ of the fit <cit.>. Many scintillation studies of this pulsar, and many scintillation studies in general, have performed this ACF fit with a Gaussian <cit.>. However, given that the ISM's PBF, a time domain function, is characterized by a one-sided decaying exponential, it makes more mathematical sense to fit the corresponding frequency domain ACF using a Lorentzian, as the two functions are Fourier pairs. In order to maintain consistency with the majority of existing scintillation studies on this pulsar while also being as mathematically thorough as possible, we performed both Gaussian and Lorentzian fits on the frequency ACFs when determining  scintillation bandwidths. The corresponding scattering delays $\tau_{\rm d}$, were then found via

\begin{equation}
\label{band_delay}
2\pi \Delta \nu_{\text{d}} \tau_{\rm d} = C_1,
\end{equation}

where the dimensionless value ${C}_{{\rm{1}}}$ can range between $0.6-1.5$ depending on the geometry of the medium and spectrum of the electron density fluctuations therein <cit.>. In this work we chose ${C}_{{\rm{1}}}=1$, as is common in many scintillation studies <cit.>.

As with all pulsar scintillation studies, our precision in estimating scintillation bandwidths and timescales from ACFs was limited by the number of scintles visible in a given observation. As such, our uncertainties were dominated by the finite scintle effect, and are quantified by

\begin{equation}
\label{finite_scintle}
\begin{split}
\epsilon &  = \frac{\Delta \nu_{\rm d}}{2\ln(2) N_{\rm scint}^{1/2}} \\
& \approx \frac{\Delta \nu_{\rm{d}}}{2\ln(2)[(1+\eta_{\text{t}}T/\Delta t_\text{d})(1+\eta_\nu B/\Delta \nu_\text{d})]^{1/2}},
\end{split}
\end{equation}

with ${N}_{{\rm{scint}}}$ being the number of scintles observed over the course of the observation, $T$ and $B$ being the total integration time and total bandwidth, respectively, and ${\eta }_{{\rm{t}}}$ and ${\eta}_{\nu}$ being filling factors that can range from $0.1$ to $0.3$ depending on how one defines the characteristic timescale and bandwidth, which for this paper are both set to 0.2 <cit.>. Throughout all of our observations, this resulted in 90$-$160 scintles in the lowest frequency subbands and around 40$-$120 scintles in the highest frequency subbands. The equivalent finite scintle error for scintillation timescale is identical to the scintillation bandwidth error aside from the factor of $\ln(2)$. These uncertainties were then added in quadrature with the uncertainties derived from the covariance matrices of ACF fits to obtain the total uncertainty for a given scintillation parameter estimation.
A complete list of all ACF-measured scintillation parameters can be seen in Table <ref>. An example dynamic spectrum, along with its 2D ACF and 1D ACF fits, are shown in Figure <ref>, while an example set of all frequency ACFs from a single epoch is shown in Figure <ref>. A full table of all scintillation bandwidth and timescale measurements will be made available in the other paper accompanying this work (Dolch et al., in prep, Turner et al., in prep).Also of note in our observations, and visible in Figure <ref>, is the strong slant in the 2D ACF, indicating a non-zero scintillation drift rate, $\rm{d}\nu/\rm{d}t$. , which in many cases would be indicative of strong refraction in our observations. However, due to this pulsar's abnormally small transverse velocity of around 6 km s$^{-1}$ <cit.>, scattering angles due to refractive scintillation are only on the order of a microarcseconds, indicating refraction to be negligible in these observations.

Subband Scintillation Bandwidth and Timescale Measurements
MJD     Frequency     $\Delta \nu_{\textrm{d}}$     $\Delta t_{\textrm{d}}$      $N_{\textrm{scint}}$
(MHz)     (MHz)     (min)

56183     1285     0.42 $\pm$ 0.04     3.56 $\pm$ 0.24     95

56183     1310     0.47 $\pm$ 0.05     3.28 $\pm$ 0.25     100

56183     1335     0.49 $\pm$ 0.06     4.15 $\pm$ 0.29     78

56183     1360     0.42 $\pm$ 0.04     3.58 $\pm$ 0.24     102

56183     1385     0.45 $\pm$ 0.05     3.77 $\pm$ 0.30     93

56183     1410     0.51 $\pm$ 0.06     4.34 $\pm$ 0.30     72

56183     1435     0.42 $\pm$ 0.04     4.04 $\pm$ 0.30     92

56183     1460     0.66 $\pm$ 0.09     5.40 $\pm$ 0.49     48

56198     1285     0.53 $\pm$ 0.06     3.96 $\pm$ 0.24     91

56198     1310     0.69 $\pm$ 0.07     2.95 $\pm$ 0.22     93

56198     1335     0.52 $\pm$ 0.05     2.88 $\pm$ 0.16     122

56198     1360     0.72 $\pm$ 0.09     4.07 $\pm$ 0.29     67

56198     1385     0.62 $\pm$ 0.06     3.34 $\pm$ 0.21     92

56198     1410     1.04 $\pm$ 0.16     4.50 $\pm$ 0.36     45

56198     1435     1.10 $\pm$ 0.17     4.78 $\pm$ 0.47     41

56198     1460     0.81 $\pm$ 0.10     3.89 $\pm$ 0.31     63

56206     1285     0.32 $\pm$ 0.03     2.88 $\pm$ 0.19     154

56206     1310     0.28 $\pm$ 0.02     3.11 $\pm$ 0.20     164

56206     1335     0.46 $\pm$ 0.04     2.73 $\pm$ 0.21     117

56206     1360     0.41 $\pm$ 0.04     2.98 $\pm$ 0.21     120

56206     1385     0.46 $\pm$ 0.05     3.36 $\pm$ 0.23     98

56206     1410     0.54 $\pm$ 0.06     3.62 $\pm$ 0.26     79

56206     1435     0.41 $\pm$ 0.04     2.94 $\pm$ 0.20     123

56206     1460     0.44 $\pm$ 0.04     3.22 $\pm$ 0.22     106
Scintillation bandwidth and timescale measurements, as well as the calculated number of scintles, across all subbands and epochs. All errors shown are 1$\sigma$.

(a) A dynamic spectrum from MJD 56206 taken at 1360 MHz. (b) The corresponding 2D autocorrelation function. (c) The 1D frequency slice (blue) of the 2D autocorrelation function with its corresponding Lorentzian fit (green). (d) The 1D time slice of the 2D autocorrelation function (blue) with its corresponding one-sided decaying exponential fit (green). Despite the strong slant in the 2D ACF, due to the low transverse velocity of PSR B1937+21, contributions to scattering due to refraction are negligible in these observations <cit.>.

The scintillation bandwidth and scintillation timescale scaling indices, ${x_{\Delta \nu_{\rm d}}}$ and $x_{\Delta t_{\rm d}}$, at each epoch were determined by performing a weighted linear least-squares fit of the form

\begin{equation}
\label{scattering_index_fit}
\Delta \nu_{\textrm{d}}=a_{\Delta \nu_{\rm d}}\nu^{x_{\Delta \nu_{\rm d}}},
\end{equation}

\begin{equation}
\label{timescale_index_fit}
% \log_{10}\Delta \nu_{\textrm{d}}=\xi\log_{10}\nu+b
\Delta t_{\rm d} =      a_{\Delta t_{\rm d}}\nu^{x_{\Delta t_{\rm d}}},
\end{equation}

with the weights being the squared inverses of the scintillation bandwidths and timescales, respectively, across all scintillation bandwidth and timescale measurements from a given MJD.
We do not see any evidence of the filterbank gaps influencing either our ACF analyses using 25 MHz subbands. as well as the full-band time ACFs.  However, their effects are present in full-band frequency ACFs. Many peaks emerge at non-zero lag in the full-band frequency ACFs that are a significant fraction of the height of the zero-lag peak. However, efforts to remove the filterbank gaps from the dynamic spectra by setting them equal to the spectra background noise yielded no significant change in our scintillation bandwidth estimations, full-band or otherwise. As a result, we can conclude that the filterbank gaps make no significant difference in our scintillation bandwidth measurements.

Frequency ACFs for all observing frequencies on MJD 56206

Our observations are separated by 15 and eight days (see Table <ref>.) By calculating refractive timescales for each epoch at each frequency and then taking a weighted average, we can check whether our scintillation measurements are independent from each other. From <cit.>, the refractive timescale, $t_r$, can be approximated by

\begin{equation}
t_r\approx\frac{4}{\pi}\Big(\frac{\nu\Delta t_{\rm d}}{\Delta \nu_{\rm d}}\Big).
\end{equation}

Using the above equation, we find weighted average refractive timescales of 10.2 $\pm$ 0.5, 6.3 $\pm$ 0.3, and 9.2 $\pm$ 0.4 days for epochs 56183, 56198, and 56206 respectively. Based on these values, we can conclude that MJDs 56183 and 56198 do not fall within the same refractive timescale and are consequently independent from each other. MJDs 56198 and 56206 may marginally fall within the same refractive timescale, although the refractive timescale of 6.3 days for MJD 56198 renders this claim inconclusive.

§.§ Secondary Spectra

To create the secondary spectrum for each observation, we took the squared modulus of the two-dimensional Fourier transform of the corresponding dynamic spectrum (see Figure <ref>.

Secondary spectra computed from the two dynamic spectra shown Figure <ref>. The two order-of-magnitude difference in the y-axis ranges is a consequence of the different delay Nyquist limits, which are caused by the different frequency resolutions of their corresponding dynamic spectra.

At each frequency slice, to measure the curvature of each arm in a given arc, we followed the approach described in <cit.> whereby we divided each secondary spectrum along the center of its fringe frequency axis and fit the maximum arc value at each delay using $\tau=\eta f_D^2$, where $\tau$ is the differential time delay (Fourier conjugate variable to observing frequency), $\eta$ is the arc curvature, and $f_D$ is the fringe frequency (Fourier conjugate variable to time).

Following <cit.>, to measure scintillation arc scaling indices, which we refer to as $x_\eta$, we performed a weighted linear least-squares fit of the form

\begin{equation}
\label{arc_scale_fit}
%\log_{10}\eta = \alpha\log_{10}\nu+\beta,
\eta = a_\eta \nu^{x_\eta},
\end{equation}

with the weights being the squared inverses of the uncertainties on each arc curvature, using the measured curvatures of a given arm at all frequencies within a given MJD. This fit was performed separately for arc measurements on the left and right sides of the fringe frequency axis to acquire scaling indices $x_{\eta,L}$ and $x_{\eta,R}$, respectively.

§ RESULTS & DISCUSSION

§.§ Scintillation Bandwidth and Timescale Measurements & Scaling

Although scattering delays are not currently accounted for explicitly in PTA timing models, the ISM introduces time-variable delays on the order of tens to hundreds of nanoseconds at at center frequencies around 1.4 GHz in MSPs timed by various PTAs <cit.>. Observing with larger frequency bandwidths can improve precision of timing model parameters, such as correcting for the frequency-dependent effects of DM. However, wider bandwidths also mean a wider range of frequencies over which effects of the ISM can vary in influence. Additionally, pulse profiles evolve with frequency and scintillation introduces further frequency-dependent modulations in pulse intensity, adding stochastic noise to this chromatic behavior. In the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) PTA collaboration, current practice corrects for time-variable dispersion measure (DM) using a piece-wise fitting procedure <cit.>. If a similar approach were to be used to account for scattering delays at a given epoch, we might assume the medium exhibits Kolmogorov turbulence and use a $\nu^{-4.4}$ dependence for scattering delays across the band. However, it has been shown that many LoSs frequently deviate from this dependence and have scaling indices that fluctuate across epochs <cit.>. There is also a strong covariance between DM and scattering delay determined this way as a consequence of DM fits absorbing scattering delays <cit.>. As such, this approach could result in TOA misestimations of tens to hundreds of nanoseconds depending on the pulsar. For this reason, it becomes crucial to understand how scintillation bandwidth, and, consequently, scattering delay, vary across these wide observing bands when attempting to correct for scattering effects.
Thanks to our ability to fit eight simultaneous measurements across the observing band, our scintillation bandwidth scaling analyses are robust in the context of a PTA-style, i.e., with an emphasis on time resolution, usually at the expense of frequency resolution, observation at L-band, as similar efforts on data with current NANOGrav frequency resolution yield scintillation bandwidths estimations that are 2$-$3 times larger than the true scintillation bandwidth and scaling indices that are consistent with zero within error.  The results are shown in Table <ref>, with the corresponding fits shown in Figure <ref>. Taking a weighted average of each subband across all three epochs and then performing the fit results in a cross-epoch scaling index of 2.58 $\pm$ 1.09. This weighted average fit is shown in Figure <ref>.

[Scintillation bandwidth scaling index fit on MJD 56183.]
[Scintillation bandwidth scaling index fit on MJD 56198.]
[Scintillation bandwidth scaling index fit on MJD 56206.]
Scintillation bandwidth scaling indices for all three epochs, with scaling index fits indicated by the blue lines. 1$\sigma$ error bars are shown. Both axes are in log space on all plots.

Fitted Pulsar Scintillation Bandwidth Scaling Indices
MJD      $x_{\Delta \nu_{\rm d}}$
56183     1.25 $\pm$ 1.11

56198     4.26 $\pm$ 1.69

56206     3.09 $\pm$ 1.40

Fitted scintillation bandwidth scaling indices to the eight frequency ACFs across the observing band. Errors are 1$\sigma$ fit uncertainties.

Scintillation bandwidth scaling index fit (blue line) using the weighted average of all scintillation bandwidth measurements in a given subband across all three epochs, with 1$\sigma$ weighted uncertainties shown. Both axes are in log space.

Fitted Pulsar Scattering Delay Scaling Indices

MJD      Index

56183     $-$2.40 $\pm$ 0.67

56198     $-$1.75 $\pm$ 1.48

56206     $-$3.83 $\pm$ 0.99

Fitted scattering delay scaling indices to the eight frequency ACFs across the observing band. Errors are 1$\sigma$ fit uncertainties. Our measurements are all consistent with a shallower power law than that expected of a Kolmogorov medium, agreeing with previous measurements in the literature.

The indices show some noticeable variation between epochs, with the value reported on MJD 56183 being noticeably shallower than MJDs 56198 or 56206. Since these indices can fluctuate from epoch to epoch as the LoS changes, and presumably similar LoSs would occur the closer in time measurements were taken, this may at least partially explain why the indices on MJDs 56198 and 56206 are more similar to each other than MJDs 56198 and 56183, as the former two are only eight days apart, while the latter two are 15 days apart. While the final two epochs have scaling indices consistent with Kolmogorov turbulence, the first epoch does not, although it does have a large uncertainty.

Our scaling indices all agree with each other within error and are consistent with those previously reported in the literature for this pulsar, with sources generally quoting indices between 3$-$3.6 within 1$\sigma$ uncertainty both on indices averaged over many epochs as well as for individual epochs <cit.>, although <cit.> reported indices on individual days as shallow as 1.3 $\pm$ 0.3 and as steep as 4.9 $\pm$ 0.8. Given the relation between the Lorentzian function and the structure of the ISM's PBF, it is possible that the generally shallower indices measured using this function indicate the true scattering delay scaling relation of this pulsar may be shallower than previously reported, as previous scintillation bandwidth measurements of this pulsar have generally relied on Gaussian ACF fits.

The treatment of this LOS as having an infinite screen is likely not appropriate, as it has been shown that shallower scaling indices such as the ones seen for this pulsar are more characteristic of finite thin screens <cit.>. It has also been observed that higher DM pulsars in general tend to exhibit shallower scaling laws <cit.>, although this may also be a consequence of a truncated scattering region. These scaling indices and the clumpiness of the scintillation arcs in our data (see Section <ref>) seem to suggest we are dealing with an extended medium, or possibly a thick screen.
Generally, it can be challenging to determine whether such deviations from Kolmogorov turbulence in data are a consequence of assumptions in analyses, lack of S/N, or are truly astrophysical in origin. Fortunately, scintillation arcs can provide insight through which we can distinguish between these possibilities. As discussed in Section <ref>, all arcs seen in our data exhibit patchiness, indicating an inhomogeneous scattering screen or a double lensing effect, providing evidence that our measured deviations from power laws consistent with Kolmogorov turbulence are astrophysical in nature, at least for our final two epochs.

As a consequence of the longer observing times required to measure scintillation timescales for most pulsars, few studies examine the corresponding scaling behavior. The studies that do exist use measurements taken on separate days <cit.>, and simultaneously at multiple frequencies <cit.>, although the latter has primarily been done at observing frequencies considerably below 1.4 GHz, where scintillation timescales are much shorter. While a pulsar like PSR B1937+21 has a relatively short scintillation timescale (around 7 minutes at 1.4 GHz per <cit.> and around 3-5 minutes in this work), making it possible to track using PTA-length observing campaigns, the large scattering delays make resolving scintles difficult under these setups. Our scintillation timescale scaling indices are shown in Table <ref>, with the corresponding fits in Figure <ref>. Taking a weighted average of each subband across all three epochs and then performing the fit results in a cross-epoch scaling index of 1.31 $\pm$ 0.65. This weighted average fit is shown in Figure <ref>. Our measurements for the final two epochs are consistent with indices measured in the three studies mentioned above. At least some discrepancies between our results and the aforementioned studies may be due in part to the different methodologies used to estimate scintillation timescales. While we fit a one-sided decaying exponential to time ACFs, <cit.> and <cit.> fit Gaussians, while <cit.> looked at the structure functions of flux density variations to obtain scintillation timescales. Like with our scintillation bandwidth indices, the scintillation timescale indices from MJDs 56198 and 56206 agree more with each other than either does with 56183, further supporting the notion that observations taken close together in time should have similar LoSs, consequentially resulting in similar scaling indices.

Under the assumption of a Kolmogorov medium with a negligible inner scale, one would expect scintillation timescale to scale as $\Delta t_{\rm d} \propto \nu^{1.2}$ <cit.>. However, the three studies cited above find a shallower index than what would be expected if the medium exhibited Kolmogorov turbulence, agreeing with the findings of most scintillation bandwidth index measurements. Two of our three epochs have scintillation timescale scaling indices that are consistent with a Kolmogorov medium well within 1$\sigma$, with the final epoch being consistent within 1.8$\sigma$.

Fitted Pulsar Scintillation Timescale Scaling Indices
MJD      $x_{\Delta t_{\rm d}}$
56183     2.38 $\pm$ 0.86

56198     1.75 $\pm$ 1.58

56206     0.90 $\pm$ 0.70

Fitted scintillation timescale scaling indices to the eight frequency ACFs across the observing band. Errors are 1$\sigma$ fit uncertainties.

[Scintillation timescale scaling index fit on MJD 56183.]
[Scintillation timescale scaling index fit on MJD 56198.]
[Scintillation timescale scaling index fit on MJD 56206.]
Scintillation timescale scaling indices for all three epochs, with scaling index fits indicated by the blue lines. 1$\sigma$ error bars are shown. Both axes are in log space on all plots.

Scintillation timescale scaling index fit (blue line) using the weighted average of all scintillation timescale measurements in a given subband across all three epochs, with 1$\sigma$ weighted uncertainties shown. Both axes are in log space.

§.§ Scintillation Arcs

With our high frequency resolution measurements, we can obtain secondary spectra with conjugate frequencies $\tau$ (also called delay) that extend far beyond what is traditionally possible for PTA-style observations, i.e., time and frequency resolutions that are limited by the pulse width and sampling rates. This is a consequence of the Nyquist limit for delay, which is given by $\tau_{\rm Nyq} = N_{\rm chan}/2\rm{BW}$, where $N_{\rm chan}$ is the number of frequency channels and BW is the observing bandwidth. The difference in both the visible range and detail within spectra can be seen quite clearly when comparing secondary spectra created using data processed with a frequency resolution typical of NANOGrav observations and secondary spectra processed using CS, as in Figure <ref>. Of particular significance is that in standard NANOGrav observations we are only sensitive to $\tau$ values less than 0.3 $\mu$s, whereas our CS-processed data extends out to over 80 $\mu$s. However, it is important to note that we have limited the range in our plots to 20$\mu$s, as the range over which arcs are visible in our data only extend out to close to this delay. Comparing these two arcs, we can see that the secondary spectrum processed at current NANOGrav frequency resolution only contains features associated with the power near the origin of the secondary spectrum in the CS-processed data, and so neither of the arms of the arc, much less features within the arms, are visible.
Additionally, we have the ability to analyze the intra-epoch evolution of scintillation arcs across a wide range of frequencies, as we were able to resolve a scintillation arc in each of our 25 MHz slices across 200 MHz of bandwidth in all epochs. This represents a substantial improvement in arc sensitivity and indicates the use of cyclic spectroscopy in future studies of MSPs could lead to important developments in our understanding of both inter- and intra-epoch arc evolution, particularly when observed over the large bandwidths used in PTAs. Some example scintillation arcs with overlaid fits can be seen in Figure <ref>.

[Secondary spectrum on MJD 56183 at 1310 MHz, with the fits to both arms in the scintillation arc shown in orange.]
[Secondary spectrum on MJD 56198 at 1310 MHz, with the fits to both arms in the scintillation arc shown in orange.]
Example scintillation arcs with overlaid fits.

§.§.§ Arc Features & Evolution Over Frequency

Scintillation arcs for this pulsar have been previously shown once in the literature at 427 MHz also utilizing cyclic spectroscopy <cit.> and once at L-band using Fourier spectroscopy filterbank methods <cit.>. While the arc power seen in <cit.> is roughly symmetric, the arcs in both <cit.> and this work all display noticeable asymmetry. Additionally, at least two arcs can be seen in both their data and ours. However, due to the patchiness of these arcs and small size of calculated refraction angles, this asymmetry is likely not the result of refraction, but rather the result of various features in the arms moving across the arcs over time.An important distinction between the arcs in <cit.> and those in both our work and <cit.> is the significant difference in the scale of our fringe frequency axes. Using a by-eye comparison and the relation between fringe frequency, observing wavelength, and angle on the sky, $\theta \propto f_D\nu$ <cit.>, we have shown the arcs in our data probe an angular scale 6$-$10 times smaller than the arcs in <cit.>. Additionally, due to the significantly larger delay axes in <cit.>, all structures in the arcs we observe are fully contained within the central feature of their arcs.
Arclets, resolved as patchy features in our scintillation arcs, which have been attributed to a variety of astrophysical phenomena, including AU-scale inhomogeneities within the scattering screen <cit.> or a double lensing effect from material behind the scattering screen relative to Earth <cit.>, can be seen in the secondary spectra in our data, and indeed can be seen in all of the observed scintillation arcs in both arms on MJDs 56183 and 56198 and in the left arm on MJD 56206. These were also seen in the arcs shown in <cit.>, although only visible in the right arm, which is also the brighter arc in those observations, but not in <cit.>. They mention their arcs are largely featureless, although they do report some discrete structures in a single epoch. This makes our observations the first reported detections of arclets in this pulsar at L-band. As in <cit.>, we expect arcs at lower frequencies to be wider and have more diffuse features. While difficult to discern throughout the evolution in the arcs in our data, when comparing with the secondary spectrum seen in the 427 MHz observations in <cit.> the arcs displayed in that spectrum are clearly wider and more diffuse compared to the arcs in our data when examined by eye (around 6-8 mHz in width for their data compared to around 2-3 mHz in width in ours), agreeing with this expectation.
At all three epochs, we find that the arcs dim with increasing observing frequency. The positions of the arclets remain constant with frequency and dim from the edges to the center of the arcs with increasing frequency. Regardless of observing frequency, we expect the features closest to the centers of arcs to remain the brightest, as they correspond to constructive interference in the larger scintle structure in the dynamic spectrum. The first structures that disappear are the dim portions of the arc higher up the delay axis, while the only remaining structures at the highest frequencies are the apexes of the arclets, which were in fact the brightest structures in the arc aside from the concentrated power at the base of the fringe  frequency axis, which also remains visible. It is well established that pulsar flux decreases with frequency over the ranges we observe <cit.>, and we can quantitatively see this in the dynamic spectra in our data, as the mean and median flux density of the 25 MHz frequency slices of these spectra clearly get progressively smaller at higher frequencies. As a result, it may be that we simply are unable to see some dimmer features in the higher frequency arcs as a consequence. Additionally, this pulsar is known for having a particularly steep spectral index <cit.>, making this effect much easier to see. That being said, these observations only span 200 MHz, and so changes in flux density over frequency likely only partially explain what we are seeing in our data. An example of this evolution can be seen in Figure <ref>. Note that, although all plots in this figure display the same dynamic range for ease of visual comparison, these features are still present when we allow their dynamic ranges to vary.

Scintillation arc evolution on MJD 56183 in increasing 25 MHz intervals across and down beginning at 1285 MHz (top left) and ending at 1460 MHz (bottom right). Features in the arcs dim as we go to higher frequencies, with only the power near to base of the parabolas and centers of the arclets remaining visible at the highest frequency.

§.§.§ Curvature Measurements & Scaling

The results of our scintillation arc scaling analysis are shown in Table <ref>, with example fits from MJD 56198 shown in Figure <ref>. We can see that only some of our scaling indices agree with  the $\eta \propto \nu^{-2}$ theory laid out in <cit.>, indicating the curvature of these arcs scale may not perfectly scale with frequency following the same power law as the angular deflection of the scattered rays. Other arc curvature analyses of MSPs have been performed <cit.>, but only one () has attempted to fit for frequency-dependent arc curvature evolution in an MSP. Unlike our fits, for which all eight measurements were taken simultaneously within the same observation, their fits used two measurements which were taken two weeks apart. They found a scaling index of $-2.35 \pm 0.07$, steeper than predicted by <cit.>, which they attribute to a change in $\textbf{V}_{\textrm{eff}, \perp}$ between their two epochs, an issue that we do not need to be concerned with for our observations, as all of our measurements for a given scaling fit are taken on the same day, and therefore have the same $\textbf{V}_{\textrm{eff}, \perp}$. As such, the results in some of our epochs mark the first agreement with the predicted arc curvature scaling relation for an MSP. Additionally, our results laid out in this subsection, along with the single arc visible in our observations, strongly suggests that the scattering along this LoS is dominated by a single scattering screen, making this the furthest and highest DM MSP for which the scattering along the LoS has been found to be dominated by a single, localized region.

Scintillation arc curvature scaling indices in both arms on MJD 56198. The fits agree well with the frequency dependence of angular deflections as detailed in <cit.>. The offset between the two fits may partially be attributable to difficulties in fitting curvatures to patchy arcs.

Fitted Pulsar Scintillation Arc Curvature Scaling Indices
MJD     $x_{\eta,L}$     $x_{\eta,R}$
56183     $--$2.05 $\pm$ 0.17     $--$1.49 $\pm$ 0.33

56198     $--$2.38 $\pm$ 0.32     $--$1.89 $\pm$ 0.27

56206     $--$1.60 $\pm$ 0.38     $--$1.29 $\pm$ 0.28

Fitted arc curvature scaling indices for both left and right arcs. Different power laws and curvatures of the left and right arms on a given day may partially be a consequence of the difficulty fitting curvatures to discontinuous arcs comprised mostly of patchy features.

While our scaling indices are somewhat consistent with $-2$, two epochs display a noticeable asymmetry in the measured index between the left and right arms. As discussed earlier, there are strong asymmetries in the structure of all arcs as well as significant patchiness which are unlikely to be the result of refractive scintillation.Despite the high resolution of these observations, this patchiness may make it difficult to properly measure curvature in all subbands, particularly if the apexes of some arclets are offset from the main parabolic fit, as in <cit.>. As a result, measured curvatures may slightly deviate from the true curvature of the arcs, which may at least in part contribute to the power law asymmetries we have measured. Consistent offsets in the curvatures of the left and right arms are also visible across all frequencies and all epochs, with curvatures in the right arms always being smaller than those on the left, even when there is more visible structure in the right arm than the left. Interestingly, over the course of these observations, the arc with most power and structure is initially the left one, but this power transitions to the right arc over the course of our observations, whilethe arc curvatures in right arm transitions from consistently smaller than the left arm to almost identical. While this may partially be the result of lensing or anisotropic scattering, the origin of this asymmetry remains unclear.An example of these offsets can also be seen in Figure <ref>.  as a $-2$ scaling index only assumes scattering as a consequence of diffractive effects.ccc

Fitted Pulsar Scintillation Arc Curvature Scaling Indices

MJD     Scaling Index Left Arc     Scaling Index Right Arc

56183     $-$2.05 $\pm$ 0.68     $-$1.84 $\pm$ 0.35

56198     $-$1.76 $\pm$ 0.16     $-$1.83 $\pm$ 0.36

56206     $-$2.02 $\pm$ 0.59     $-$1.86 $\pm$ 0.39

Fitted arc curvature scaling indices for both left and right arcs. We can see that the majority of our fits agree with arc curvatures scaling in the same way as angular deflection of the scattered rays <cit.>. Different power laws and curvatures of the left and right arms on a given day can largely be attributed to high levels of refraction occurring during these observations.

Fitted Pulsar Scintillation Arc Curvature Scaling Indices

MJD     Scaling Index Left Arc     Scaling Index Right Arc

56183     $-$2.05 $\pm$ 0.68     $-$1.84 $\pm$ 0.35

56198     $-$1.76 $\pm$ 0.16     $-$1.83 $\pm$ 0.36

56206     $-$2.02 $\pm$ 0.59     $-$1.86 $\pm$ 0.39

Fitted arc curvature scaling indices for both left and right arcs. We can see that the majority of our fits agree with arc curvatures scaling in the same way as angular deflection of the scattered rays <cit.>. Different power laws and curvatures of the left and right arms on a given day can largely be attributed to high levels of refraction occurring during these observations.

Once curvature measurements are obtained, they can be used in conjunction with $\textbf{V}_{\textrm{eff}, \perp}$ to estimate screen distances and potentially associate scattering screens with known structures along the LOS to a pulsar. In cases where the pulsar velocity dominates, we can ignore the velocities of Earth and the ISM when estimating these distances. However, the low velocity of PSR B1937+21 (9 km s$^{-1}$, per <cit.>) requires us to take these additional velocities into account. Unfortunately, as a consequence, degeneracies in our estimated screen distance solutions (either around 0.3 or 2.1 kpc from Earth, assuming a screen orientation of zero degrees) necessitate high cadence, long-term arc measurements to properly constrain these distances <cit.>. These long term observations can properly account for the orbital motion of Earth, which causes periodic variations in the observed curvature, and can help break the degeneracy between screen distance and screen orientation angle. However, it is unlikely that there would be any obvious structures, such as an HII region or molecular cloud, along the LOS to this pulsar which we could associate with a screen in the first place. The LOS is mostly devoid of significant structure and predominantly consists of the inter-arm region between the Orion-Cygnus and Carina-Sagittarius spiral arms of the Milky Way, although a small fraction of the material along the LOS can be attributed to the Orion-Cygnus arm <cit.>. This likely makes interstellar plasma filaments the most plausible candidate for this pulsar's screen(s) (Stock et al., in prep). In general, as more scintillation arc distances are constrained, there is an increasing need to identify such diffuse structures along LoSs with scattering <cit.>.That being said, one of the degenerate solutions places a screen around 150$-$200 pc away from Earth, fairly close to the edge of the local bubble <cit.>.

§ CONCLUSIONS & FUTURE WORK

We have used cyclic spectroscopy to improve the frequency resolution of L-band baseband observations of PSR B1937+21 compared to what is feasible in PTA observing setups that rely on limits of non-cyclostationary time and frequency sampling, yielding highly detailed dynamic and secondary spectra and allowing for scintillation parameter estimations within 25 MHz subbands. These results have yielded, in some epochs, the first arc curvature scaling fits for an MSP that agree with theory, as well as scintillation bandwidth scaling indices that match previous observations and the first detection of arclets in this pulsar at L-band. We found arc curvature scaling indices for different arms within individual observations to exhibit noticeable asymmetries. , which we have attributed to refractive scintillation in these epochs., indicating screen anisotropy or a double lensing effect in this frequency range.An examination of arc feature evolution over frequency shows a dimming of prominent features as we go to higher frequencies, at least in part likely due to this pulsar's steep spectral index. , possibly originating from anisotropic scattering or lensing. The thorough analysis in this work made possible through cyclic spectroscopy demonstrates the potential for larger-scale efforts using this technique and adds to a growing case towards its eventual adoption by PTAs  due to its ability to calculate high-resolution dynamic and secondary spectra. Traditional polyphase filterbanks have limited the combined time and frequency resolution of scintillation studies with pulsars, particularly in studies of scintillation arcs because the secondary spectrum is Nyquist limited. This in turn may make more pulsars accessible for scintillation studies at lower observing frequencies. Additionally, the ability to use cyclic spectroscopy to recover a pulsar's PBF through deconvolution <cit.> or other methods <cit.> adds another element to these studies, especially given that the PBF may present a more accurate picture of the scattering present in a given observation than the traditional ACF fitting approach <cit.>.
Efforts are currently underway at the Green Bank Observatory to implement a real-time or near real-time cyclic spectroscopy backend, with development nearing completion and testing on real data to begin shortly. This will provide exciting opportunities for new science to all astronomers studying periodic or pulsed emission using this facility, as well as aid in the timing efforts of the NANOGrav PTA collaboration, for whom the 100 meter telescope at the observatory plays a significant role in the ongoing study of low-frequency gravitational waves.

§ ACKNOWLEDGEMENTS

We gratefully acknowledge support of this effort from the NSF Physics Frontiers Center grants 1430284 and 2020265 to NANOGrav. All of the data processing in this work utilized the resources of the Link computing cluster at West Virginia University. We thank the NANOGrav Collaboration for helpful discussions and suggestions, particularly the Noise Budget Working Group. We thank the AO staff for their assistance during the observations reported in this work, in particular P. Perillat, J. S. Deneva, H. Hernandez, and the telescope operators. TD, MAM, and DRS were partially supported through the National Science Foundation (NSF) PIRE program award number 0968296. TD acknowledges NSF AAG award number 2009468 as well as sabbatical and summer leave funding from Hillsdale College.

Software: scintools <cit.>, dspsr <cit.>, astropy <cit.>, pypulse <cit.>, scipy <cit.>, numpy <cit.>, and matplotlib <cit.>.

§ DISCUSSION OF CYCLIC SPECTROSCOPY IN COMPARISON WITH INCREASED FREQUENCY RESOLUTION USING TRADITIONAL TECHNIQUES

Suppose that we had a one-dimensional time series of pulses and a very high spectral
resolution. Could we still form a similar data product to the dynamic spectra we have acquired for this
paper? In this context, we are operating in what would be the $0^{\rm th}$ lag from a cyclic spectroscopy perspective and would end up with $1/T$ frequency resolution, with $T$ being the sampling time. In other words, until we start shifting by $1/P$ (or integer multiples thereof) subband spacing, we are not really doing cyclic spectroscopy.What distinguishes cyclic spectroscopy from the technique described above is that in cyclic spectroscopy-based processing, everything is focused on the periodicity of the signal, whereas the above method treats the pulsar like a continuum source. While we can indeed achieve high frequency resolution in this way, and can still generally exploit that the scintillation pattern if the resulting dynamic spectrum is sustained over multiple resolutions, a key difference is that in the CS approach the specific on-pulse/off-pulse distinction is not necessary. the power spectrum will not always be positive definite and the resulting periodic spectrum can have negative values, whereas the power spectrum would be positive definite using the process described above.

Additionally, CS is required to achieve this resolution in the context of data quality required for pulsar timing arrays. By averaging over a slice of the frequency-time plane that is smaller than the uncertainty principle dictated by the Gabor limit of signal processing, cyclic spectroscopy allows one to overcome that limit for a cyclostationary signal like that of a pulsar as seen from Earth. Exceeding this resolution limit is also what results in negative values in power spectra described in the previous paragraph. Pulsar timing arrays require high time resolution, usually on the order of microseconds, to accomplish their scientific goals. In the scenario described above, achieving such time resolution would be impossible with the desired frequency resolution. More specifically, to achieve 1 MHz frequency resolution, one would be limited to 1 $\mu$s time resolution. Since CS allows for both high time and frequency resolution, the dynamic spectra shown throughout our paper would only be possible through CS in data that simultaneously requires the high timing resolution needed for pulsar timing arrays.

Perhaps most importantly, there does appear to be a fundamental difference between a dynamic spectrum retrieved using CS compared with conventional techniques, even at the same time and frequency resolution. It appears the CS approach actually is better than attempting to form a dynamic spectrum of equivalent resolution using traditional methods. In the particular case of B1937+21, to achieve the approximately 6.1 kHz frequency resolution in the dynamic spectra for this paper using standard processing techniques, we would be limited to a time resolution of approximately 0.16 milliseconds. This would limit us to four phase bins for this pulsar, assuming we are only allowing for phase bins that are powers of two. The resulting dynamic spectrum is shown in Figure <ref> alongside the dynamic spectrum with the exact same resolution but produced using cyclic spectroscopy.

(Dynamic spectra from data with 6.1 kHz frequency channel resolution. The spectrum on the left is from data processed using traditional methods and so is limited to four phase bins, while the spectrum on the right is from data processed using cyclic spectroscopy, and therefore has no such limitations on phase resolution.

While the scintillation pattern is essentially recovered using conventional processing, the spectrum
appears to be broken up into smaller bands within each polyphase filterbank (PFB) that are not present in the CS-processed data, resulting in additional filterbank-style gaps. This effect can be seen even more clearly if we examine these dynamic spectra after they have been integrated in time, shown in Figure <ref>.

Dynamic spectra from non-CS (green) and CS (blue) processed data from Figure <ref> integrated in time. The banding effect seen in the full non-CS dynamic spectrum is even more apparent after this integration. The values in the CS-processed dynamic spectrum have been scaled down for visual aid.

This effect is caused by the periodic modulation of the on/off sections of the pulsar signal, which is not properly accounted for in the non-CS approach. Consider the time domain situation: the radio signal is multiplied by a periodic windowing function that has many periods, all of the same form. The convolution theorem tells us that the the radio spectrum will be convolved with the Fourier transform of the resulting long, periodic modulating function. This demonstrates why it is essential to do the CS-style processing, even if one is only going to use the zero-lag (positive definite) part of the signal. Additionally, if one only wants or needs the amplitude of the spectrum, then a CS-processed (i.e. pulse-phase cognizant) zero-lag spectrum is all that is needed. That being said, CS becomes particularly important if one wants to extract any of the electromagnetic phase from the signal. This was demonstrated in <cit.>, where it was shown that, except for really low S/N situations, one is always better off doing the full CS processing and using some reconstruction technique to get the PBF. This is particularly true in the case of most interest to NANOGrav, i.e., where there are few scintles across the band. In this case, an ACF-inferred scattering delay is a poor estimator compared the the CS-derived estimator.

Artifacts are also seen in the resulting secondary spectra, although they are a bit more subtle, as shown in Figure <ref>. The main distinguishing feature is the segmented bright spike that arises at 0 mHz, which, again, is a consequence of not fully treating the signal as a pulsed source.

Secondary spectra from the analogous dynamic spectra in Figure <ref>. The spectrum on the left is from data processed using traditional methods, while the spectrum on the right is from data processed using cyclic spectroscopy.