CONSTRAINING THE VELA PULSAR'S RADIO EMISSION REGION USING NYQUIST-LIMITED SCINTILLATION STATISTICS

Owing to their extraordinary emission, pulsars are invaluable probes of a vast array of extreme physics, from the supranuclear density interiors and the ∼10¹² G magnetospheres of neutron stars, to the turbulent plasma of the interstellar medium (ISM). The regularity of their time-averaged emission has been applied to sensitive tests of general relativity (Taylor & Weisberg 1989; Kramer et al. 2006), and the achieved timing precision holds promise for direct observation of gravitational waves (Detweiler 1979; Jenet et al. 2006; Hobbs et al. 2010). Yet, remarkably, no consensus exists for the precise origin of this emission; its geometry, stability, and physical mechanism remain enigmatic, even after nearly a half century of study (Melrose 2004; Verbiest et al. 2009). The incredible compactness of the pulsar and its surrounding plasma prohibit traditional imaging of the emission process, but a number of methods can achieve the requisite resolution, using the interstellar plasma as an effective "lens" (Lovelace 1970; Backer 1975; Cordes et al. 1983; Cornwell & Narayan 1993; Gwinn et al. 1998; Shishov 2010).

Variations of this stochastic lens give rise to the observed scintillation, and statistics of these variations convey information about the spatial structure of the emission. Longer observing wavelengths enhance the signature of this structure because the stronger scattering effectively enlarges the aperture. However, the increased spatial resolution comes at the expense of spectral resolution, and it becomes difficult to decouple the scintillation from source noise and variability. Indeed, the pulse-broadening timescale can easily exceed that of intrinsic pulsar variation.

We can circumvent these difficulties by constructing spectra over accumulation times that include all pulsed power. Intrinsic sub-pulse variability then affects the correlation of spectral noise but not the mean spectrum or single-channel statistics (Gwinn & Johnson 2011). Furthermore, we can distinguish the effects of extended source geometry from the scattering evolution and pulsar variability by varying the degree of temporal averaging; with no averaging, only the source emission structure contributes. Additionally, we can identify the emission sizes of individual pulses. We have presented a detailed mathematical treatment separately (Johnson & Gwinn 2012, hereafter JG12).

We constructed spectra in this way and then compared histograms formed from the spectral samples with appropriate models for a strongly scintillating source. Such models benefit from minimal assumptions about the nature and distribution of both the scattering material and the underlying emission; this generality is achieved at the expense of resolution in pulse phase. The inferred emission size is therefore a characteristic size of the emission region, quantified via a dimensionless parameter γ_s, which depends on the spatial standard deviation of source intensity, weighted by integrated flux density. Translating γ_s to physical units at the pulsar additionally requires specification of the scattering geometry.

Several previous investigations have also quantified scintillation statistics in order to estimate the size of the Vela pulsar's emission region. In general, for a strongly scattered point source, the scintillation acts upon intensity as a multiplicative gain $\mathcal {G}$ that is the squared modulus of a complex random walk. The probability density function (PDF) of $\mathcal {G}$ is therefore exponential, $P(\mathcal {G}) = e^{-\mathcal {G}}$, with an associated modulation index of unity: $m^2 \equiv \langle \mathcal {G}^2 \rangle / \langle \mathcal {G} \rangle ^2 - 1 = 1$. A spatially extended emission region smoothes this variation, and so the modified $P(\mathcal {G})$ or reduced modulation index can be used to estimate the emission size (Salpeter 1967; Cohen et al. 1967; Gwinn et al. 1998).

The principal difficulty in detecting these modifications is the additional noise that arises from the noiselike nature of the source, known as self-noise, and from the background. Pulsars are particularly challenging targets because their dramatic intrapulse and interpulse variations also contribute modulation. All of these sources of noise can easily imitate or mask the signature of an extended emission region. Previous efforts have attempted to minimize these sources of noise by spectral averaging in both frequency and time. In contrast, the present technique fully accounts for the effects of self-noise and intrinsic pulsar variability.

For example, Gwinn et al. (1997) sought to detect the modification of $P(\mathcal {G})$ through the PDF of the correlation function amplitude on the Tidbinbilla-Parkes baseline at 2.3 GHz. Because this baseline was much shorter than the scale of the diffraction pattern, the observational setup was effectively a zero-baseline interferometer. Their best-fitting model corresponded to an emission region with an FWHM of 460 ± 110 km. Gwinn et al. (2000) then extended this analysis to three gates across the pulse, also at 2.3 GHz. They improved the model PDF and quantified the effects of their averaging in frequency and time. They estimated a declining size of the emission region from 440 ± 90 km to less than 200 km across the pulse.

Rather than analyzing distribution functions, Macquart et al. (2000) used the modulation index to quantify the effects of emission size for single-dish data. They measured m² ≈ 0.9 at 660 MHz and thereby inferred the extent of source emission to be no more than 50 km (FWHM = 120 km). Including the contribution of self-noise leads to an inferred FWHM of the emission region of ∼200 km (JG12).

More recently, Gwinn et al. (2012) analyzed visibilities on the Tidbinbilla-Mopra baseline at 1.65 GHz and fit models to the projections of the real part and imaginary variance. They found the size of the emission region to be large both early and late in the pulse (∼400 km and ∼800 km, respectively) but near zero in the central portion, where the pulse is strong.

An analysis at lower frequencies has the advantage of increasing the signature of size but the disadvantage of incurring more mixing between the scintillation and the pulsar noise. We demonstrate that this difficulty is not prohibitive, and that a careful statistical description of observed power spectra at low frequencies can yield fresh insight into the physics of both interstellar scattering and pulsar emission.

In Section 2, we briefly review the theoretical descriptions of pulsar emission and interstellar scintillation. Then, in Section 3, we describe our observations and the subsequent construction of both the observed and model PDFs of intensity; we also present the expected residual structure arising from an extended emission region and outline the model fitting procedure. Next, in Section 4, we give our analysis for various pulse selection cuts and derive our estimates for the size of the emission region. In Section 5, we quantify possible biases and systematic errors in our inferences of the size of the emission region. Then, in Section 6, we interpret our measurements in terms of beamed, dipolar emission and derive the corresponding emission altitudes. Finally, in Section 7, we summarize our results, their impact on our understanding of pulsar magnetospheres, and the prospects for future work.

A full understanding of pulsar emission is elusive, but proposed scenarios typically invoke a relativistic electron–positron plasma that flows outward from the pulsar polar caps to the light cylinder (Goldreich & Julian 1969). This outflow then generates radio emission via coherent curvature radiation (Sturrock 1971; Ruderman & Sutherland 1975; Mitra et al. 2009). The emission is beamed, ducted, and refracted as it propagates through the upper magnetosphere (Barnard & Arons 1986; Arons & Barnard 1986; Lyutikov & Parikh 2000).

The precise altitude at which the emission is generated is uncertain, and various analyses suggest emission heights from a few stellar radii (Arons & Scharlemann 1979) to hundreds of stellar radii (Karastergiou & Johnston 2007) for Vela-like pulsars. In the classification and interpretation of Rankin (1990), the Vela pulsar is a core-single ($\mathrm{\textbf {S}}_\mathrm{t}$) pulsar, with emission arising effectively at the stellar surface.

2.2.1. Scintillation Physics

Density inhomogeneities in the dilute, turbulent plasma of the ISM scatter the radio emission, leading to multipath propagation. At 760 MHz, the Vela pulsar is deep within the regime of "strong" scattering, for which the Fresnel scale is much larger than the diffractive scale. Furthermore, the timescale for evolution of the diffraction pattern at the observer is several seconds, so each scintillation element is effectively static for the duration of a single pulse (∼5 ms). The scattering then acts to convolve the pulsar signal with a stochastic propagation kernel, which reflects the scintillation and temporal broadening (Hankins 1971; Williamson 1972; Gwinn & Johnson 2011).

The spatial extent of the emission region influences the scintillation. A pair of dimensionless parameters, {γ_{s, 1}, γ_{s, 2}}, quantify the emission region and give its squared characteristic size in an orthogonal, transverse basis $\lbrace \hat{\textbf {x}}_1,\hat{\textbf {x}}_2 \rbrace$, in units of the magnified diffractive scale (Gwinn et al. 1998):

$$\begin{eqnarray} &&\gamma _{\mathrm{s},i} = \left(\frac{D}{R} \frac{\sigma _i}{\frac{1}{2\pi }\frac{\lambda }{\theta _i} } \right)^2\!. \end{eqnarray} \tag{ 1 }$$

Here, D is the characteristic observer–scatterer distance, R is the characteristic source–scatterer distance, λ is the observing wavelength, θ_i is the angular size of the scattering disk along $\hat{x}_i$, and σ_i is the standard deviation of the distribution of source intensity along $\hat{x}_i$. This representation accommodates anisotropy of both the source emission and the distribution of scattering material.

Nevertheless, the effects of anisotropy in the emission and scattering only manifest at quadratic order in γ_s; thus, we use γ_s ≡ (γ_{s, 1} + γ_{s, 2})/2 to avoid parameter degeneracy. Likewise, we state our results in terms of the standard deviation σ_c ≡ σ₁ = σ₂ of a circular Gaussian intensity profile. The FWHM of such a profile is $2\sqrt{\ln 4} \sigma _{\mathrm{c}} \approx 2.35\sigma _{\mathrm{c}}$.

2.2.2. Interpretation of the Inferred Emission Size

2.2.3. Scattering Geometry of the Vela Pulsar

We observed the Vela pulsar for one hour on 2010 October 22 and one hour on 2010 October 25 using the Robert C. Byrd Green Bank Telescope (GBT) and the recently built Green Bank Ultimate Pulsar Processing Instrument (GUPPI). For flux calibration, we also observed 3C 190 (0758+143) for approximately 20 minutes prior to each pulsar observation.

We used GUPPI in 8 bit, 200 MHz mode, and our observation spanned 723.125 − 923.125 MHz. Rather than utilizing the standard real-time dedispersion of GUPPI, we recorded the channelized complex baseband voltages directly to disk for offline processing. The collected data constitute approximately 8 TB.

We processed the baseband data using the digital signal processing software library dspsr (van Straten & Bailes 2011), which first performed coherent dedispersion and then formed dynamic spectra via a phase-locked filterbank. This technique constructs each spectrum from a fixed-length time series beginning at the corresponding pulse phase, rather than windowing the pulse into a fixed number of temporal phase bins. Our spectra had 262,144 channels per 25 MHz, corresponding to 95.4 Hz resolution. The ∼10.5 ms spectral accumulation time easily contained all of the power in one pulse.

We exported the coherently dedispersed spectra from dspsr into psrfits format and then performed flux calibration for each polarization using the software suite psrchive (Hotan et al. 2004). Subsequent processing utilized the psrchive libraries, which we integrated into our custom software.

Although reversing the Faraday rotation is straightforward, we chose to analyze the measured linear polarization streams without modification in order to avoid introducing possible artifacts from calibration errors and mixing. Consistency of our results between the two polarizations is, therefore, a meaningful and encouraging indication that our inferences reflect properties of the pulsar rather than instrumental limitations.

3.3.1. RFI Excision

3.3.2. Background Noise and Bandpass Shape Estimation

3.3.3. Rejection of Pulses with Quantization Saturation

3.4.1. Histogram Formation

3.4.2. Effects of Finite Histogram Resolution

In order to compare a theoretical PDF P(I) with an observed PDF, we must account for the finite bin width of the histogram. Specifically, the theoretical (normalized) histogram value H(I; ΔI) for a bin of width ΔI centered on I is given by

$$\begin{eqnarray} &&H(I;\Delta I) = \frac{1}{\Delta I}\int _{I - \Delta I/2}^{I + \Delta I/2} P(I^{\prime }) dI^{\prime }. \end{eqnarray} \tag{ 2 }$$

Expanding P(I') about I and replacing the second derivative with its discrete representation gives

$$\begin{eqnarray} H(I;\Delta I) &\approx P(I) + \frac{P(I + \Delta I) - 2 P(I) + P(I - \Delta I)}{24}.\qquad \end{eqnarray} \tag{ 3 }$$

This approximation corrects for the finite bin width to quadratic order in ΔI, while only requiring calculation of the theoretical PDF at the histogram bin centers. For our data, the next order correction, ΔI⁴P⁽⁴⁾(I)/1920, was well below our Poisson noise. The effects of the finite bin width were small when N = 1 (typically ≲ 0.005% of the PDF amplitude), but they became appreciable for larger N.

We independently constructed observed PDFs for each scalar electric field (i.e., each linear polarization) and with various degrees of temporal averaging N, bandwidth, and pulse selection cuts. In all of these cases, the mean source and background intensities of the pulses parameterized the point-source model PDF of intensity; this parameterization accounted for the intrinsic pulse-to-pulse variability as well as the quasistatic evolution of receiver gain and noise levels. Because this treatment incorporated the observed set of pulse amplitudes rather than appealing to assumptions about their distribution, our models were robust to arbitrary pulse-to-pulse variations, including strongly correlated behavior such as nulling or heavy-tailed distributions of pulse amplitude, as are observed for some pulsars (Lundgren et al. 1995; Cairns et al. 2001).

3.5.1. Model PDF Parameter Estimation

3.5.2. Construction of the Model PDF

We now describe the basis and formation of our model PDFs of intensity. First, consider N consecutive pulses with gated source intensities A₁I_s, ..., A_NI_s and mean background intensity I_n. After averaging the N spectra formed from these pulses, the PDF of intensity is (JG12, Equation (3))

$$\begin{eqnarray} P(I;N) &=&N \int _0^\infty dG \ \frac{e^{-G}}{ \left(I_{\mathrm{s}} G\right)^{N-1} } \nonumber \\ &&\times \sum _{j=1}^N \frac{(A_j I_{\mathrm{s}} G + I_{\mathrm{n}})^{N-2} }{\prod \limits _{{{\scriptstyle {\begin{array}{c}\ell = 1\\ \ell \ne j\end{array}}}}}^N (A_j - A_\ell)} e^{-\frac{N I}{(A_j I_{\mathrm{s}} G + I_{\mathrm{n}})}}. \end{eqnarray} \tag{ 4 }$$

The integral over G arises because we assume that each spectrum effectively explores the full scintillation ensemble.

We calculated this theoretical PDF for every sequence of N pulses that contributed to an observed PDF, using the estimated model parameters {A_jI_s, I_n}. We then averaged these distributions to yield the "model PDF" of intensity. Note that the model PDF required no fitted parameters and was independently constructed for each observed PDF.

The PDF residual $\mathcal {R}(I)$ (i.e., the observed minus model PDF) contains signatures of extended emission size, decorrelation of the scintillation within the averaging time, and errors in the model parameters {A_jI_s, I_n}. These perturbations have nearly identical shape, so a single fitted parameter quantifies their combined effect. We can break this degeneracy by varying the reduction parameters such as N and the analyzed bandwidth. For example, the temporal decorrelation of the scintillation pattern does not contribute when N = 1.

We approximate the expected cumulative residual arising from the emission size and the temporal decorrelation as (JG12, Equations (9), (10), and (18))

$$\begin{eqnarray} \mathcal {R}(I) &\approx& - \left[2\gamma _{\mathrm{s}} + \frac{\sum _{i<j}(1 - \Gamma _{ij})}{N(N+1)} \right]\nonumber \\ && \times \int _0^\infty dG \ G e^{-G} \frac{\partial ^2 P(I;N|G)}{\partial G^2}. \end{eqnarray} \tag{ 5 }$$

This general form expresses the influence of the temporal autocorrelation function of the scintillation pattern:

$$\begin{eqnarray} &&\Gamma _{ij} \equiv \frac{\langle \left[I(t_i, f) - I_{\mathrm{n}}(t_i, f)\right][ I(t_j,f) - I_{\mathrm{n}}(t_j,f)] \rangle }{\langle I(t,f) - I_{\mathrm{n}}(t,f) \rangle ^2} - 1.\quad \end{eqnarray} \tag{ 6 }$$

We will quantify alternative sources of residual structure, such as parameter errors, in terms of their incurred bias on γ_s. As was required for the model PDF, this model residual was calculated for every sequence of N pulses that contributed to the corresponding observed PDF and subsequently averaged.

If the N averaged intensities are drawn from consecutive pulses, then the scale $\mathcal {R}_0$ of the residual structure can be approximated as

$$\begin{eqnarray} &&\mathcal {R}_0 \equiv \gamma _{\mathrm{s}} + \frac{\sum _{i<j}(1 - \Gamma _{ij})}{4 N(N+1)} \approx \gamma _{\mathrm{s}} + \frac{1}{24} \frac{N(N-1)}{\Delta \tau _{\mathrm{d}}^2 }, \end{eqnarray} \tag{ 7 }$$

where Δτ_d ≫ N is the decorrelation timescale in pulse periods, and we have assumed the standard square-law autocorrelation $\Gamma _{ij} = e^{-\left[ (i-j) / \Delta \tau _{\mathrm{d}} \right]^2}$. Thus, the relative influence of temporal decorrelation to emission size is approximately quadratic in N.

For each observed PDF, we calculated the corresponding model PDF and residual shape. We then performed a weighted least-squares fit to the measured residual after correcting for the finite histogram resolution; we assigned weights by assuming that the histogram errors reflected Poisson noise. The only fitted parameter was the amplitude scale $\mathcal {R}_0$ of the residual. In some cases, we performed multiple reductions on the same set of pulses by varying the averaging (N) or the lag between the averaged pulses (Δτ); we then fit all these reductions simultaneously to the pair of parameters {γ_s, Δt_d} by assuming a square-law autocorrelation structure.

To utilize the full statistical content of our observations, we first analyzed all pulses concurrently, in groups only determined by frequency and polarization. We constructed each observed PDF from a (1 hr) × (3.2 MHz) on-pulse block of data: ∼10⁹ spectral samples and ∼40, 000 pulses. The scintillation pattern evolved substantially over both the span of time and frequency; each observed PDF contained samples from ∼10⁶ independent scintillation elements.

Figure 1 demonstrates one example of the comparison of our model and observed PDFs of intensity; the plotted data correspond to the stronger linear polarization for the 3.2 MHz sub-band centered on 757.5 MHz. The agreement extends over the nearly eight decades spanned. The tail of the distribution arises from the combined effects of intrinsic variability, noise, and scintillation; the success of the model indicates that it correctly accounts for each effect.

Zoom In Zoom Out Reset image size

Our models also correctly characterize the effects of temporal averaging, N. The decorrelation timescale in pulse periods is Δτ_d ∼ 50, so the correction for temporal decorrelation is quite small for the plotted reductions: N = 1, 3, and 5 (see Equation (7)).

For N = 1, the residual structure in the stronger polarization has an amplitude of ≲ 10⁻⁴, which is comparable to the Poisson noise in the histogram. A weighted linear least-squares fit to the lowest 200 samples (I < 10) in the observed PDF estimates γ_s = (2.0 ± 1.0) × 10⁻⁴ with reduced chi-squared χ²_red = 1.13; the standard deviation of χ²_red is $\sqrt{2/200}=0.1$. Furthermore, because intrinsic modulation correlates noise in adjacent spectral channels, Poisson weights underestimate the histogram errors, so χ² ≳ 1 is expected. Thus, a conservative upper bound on the effects of size from this single sub-band is γ_s ≲ 5 × 10⁻⁴, which is equivalent to σ_c ≈ 8 km at the pulsar with the conversion of Section 2.2.3.

For the adjacent sub-bands (centered on 751.25 MHz and 763.75 MHz), the fits give γ_s = (− 2.2 ± 1.1) × 10⁻⁴ and γ_s = (− 1.2 ± 1.0) × 10⁻⁴. Combining these results refines the 3σ size upper bound to approximately 4 km at the pulsar. Narrower sub-bands give similar results, with larger errors because of the increased Poisson noise and a larger bias from the finite bandwidth (see Section 5.1). In all cases, our results suggest a characteristic size of ≲ 10 km for the emission region.

We next examined effects dominated by temporal decorrelation of the scintillation pattern. We constructed PDFs with N = 2 that averaged pairs of pulses separated by Δτ pulse periods; the amplitude of the residual structure (Equation (5)) then depends roughly quadratically on Δτ. We approximated the autocorrelation function as Gaussian, $\Gamma (\Delta \tau) = e^{-\left(\Delta \tau /\Delta \tau _{\mathrm{d}} \right)^2}\!$, to simultaneously fit all lags to the pair of parameters, {γ_s, Δτ_d}.

The residuals match their expectations quite well; Figure 2 shows one example, again from the stronger linear polarization for the 3.2 MHz sub-band centered on 757.5 MHz. This excellent agreement demonstrates the quality of our analytical representation of decorrelation of the scintillation pattern within the averaging time. The fits estimate γ_s = (9.3 ± 5.4) × 10⁻⁵, giving a 3σ upper bound on the size of the emission region as ∼7 km.

Zoom In Zoom Out Reset image size

Individual pulse profiles for the Vela pulsar show marked dependence on their strength (Krishnamohan & Downs 1983). To examine covariance between pulse strength and emission size, we segregated the pulses into quantiles by their mean flux density and then analyzed each subset independently.

This approach has several utilities. Because the resolving power improves with the signal-to-noise ratio, a subset of strong pulses can maximize the potential resolution. Furthermore, many instrumental biases vary differently with pulse strength than their effects on the inferred size do (see Section 5.2); the variation of γ_s between quantile subsets can therefore suggest the presence of such a bias.

We established deciles by pulse strength. Figure 3 shows the PDFs and their residual structure for the top and bottom deciles; the mean gated signal-to-noise 〈A_jI_s/I_n〉 in these deciles is 1.39 and 0.48, respectively. The two residuals are comparable, and both suggest an emission size that is substantially smaller than the magnified diffractive scale. The limit is most stringent for the pulses in the top decile because of their higher signal-to-noise ratios: γ_s = (1.1 ± 1.8) × 10⁻⁴, or an 11 km upper bound on the characteristic emission size at a 3σ confidence level. For the bottom decile, we obtain γ_s = (− 7.9 ± 5.5) × 10⁻⁴, which gives a 12 km upper bound at the same confidence.

Zoom In Zoom Out Reset image size

Lastly, we examined the statistics for individual pulses. For each pulse, our model requires only two parameters to describe the spectral statistics: the mean intensities of the background (I_n) and gated signal (AI_s). Poisson noise generally limits each estimate of emission size to a sizeable fraction of the magnified diffractive scale; the standard error δγ_s in the fitted γ_s is well approximated by (JG12, Equation (26))

$$\begin{eqnarray} &&\delta \gamma _{\mathrm{s}} \approx \frac{1.3}{\sqrt{N_{\mathrm{tot}}} S} \exp [0.4 + 0.18 \left(\ln S\right)^2 + 0.01 \left(\ln S \right)^3 ]. \end{eqnarray} \tag{ 8 }$$

Here, N_tot is the total number of points used to estimate the observed PDF, and S ≡ A_jI_s/I_n is the gated signal-to-noise ratio. Thus, for a strong pulse (S ∼ 2), γ_s can be determined to within ±0.006, whereas for a weak pulse (S = 0.1), γ_s can only be established to within ±0.2. Interestingly, for observations of highly variable pulsars, the resolution afforded by strong single pulses may exceed that of the entire pulse ensemble analyzed jointly.

Figure 4 shows the fitted sizes for ∼30, 000 individual pulses as a function of their gated signal-to-noise ratios; values from both linear polarizations are plotted. This sample of pulses included several "micro-giant" pulses, as well as tens of pulses with the late "bump" component identified by Johnston et al. (2001). We measured no confident detection of emission size effects; the inferred values of γ_s are consistent with the expected standard errors for observations of a source with point-like emission, as given by Equation (8).

Zoom In Zoom Out Reset image size

Noise in the model parameters introduces spurious modulation. Subsequent fits to the spectral statistics will then exhibit a positive bias in γ_s to quench this excess modulation. Fortunately, we can readily characterize the expected noise δA_j and δI_n and can therefore anticipate the influence of such a bias.

For example, the noise δA_j in the source amplitudes is dominated by the finite sampling of the scintillation pattern for a single spectrum; the corresponding bias is γ_s → γ_s + 1/(4N_scint) (JG12, Section 4.4.2), where N_scint is the number of averaged scintillation elements used to estimate A_j.

On the other hand, because the background noise is white and approximately stationary, δI_n only depends on the number of samples N_b that are averaged to estimate I_n. However, the incurred bias to γ_s also depends on the gated signal-to-noise ratio S: γ_s → γ_s + 1/(4S²N_b).

For each 3.2 MHz sub-band, we had N_scint ∼ (3.2 MHz)/Δν_d ≈ 3, 000, and N_b ∼ 10⁵. Both biases to γ_s were therefore within our quoted confidence. Furthermore, each bias scales as the reciprocal of the analyzed bandwidth. We observed no evolution in the inferred values of γ_s for 1.6 MHz and 0.8 MHz sub-bands and thus confidently established that this class of parameter error did not substantially influence our results.

Instrumental limitations can also lead to errors in the model parameters. For example, interleaved sampling results in a negatively dedispersed pulse "echo"; for highly dispersed pulsars, such as the Vela pulsar, this process can therefore leak pulsed power into the off-pulse region, leading to a positive bias in the estimated background noise: I_n → I_n + δI_n. Alternatively, quantization introduces signal-dependent noise (Jenet & Anderson 1998; Gwinn 2006); the off-pulse region then underestimates the effective on-pulse noise.

Of course, any bias in the estimated background noise incurs a corresponding bias in the estimated source amplitude, and vice versa: δI_n = −δA_jI_s. These errors then affect the inferred emission size as γ_s → γ_s − δI_n/(2AI_s). Because the dependence of emission size with pulse strength may vary differently than the noise bias, segregating pulses by their gated signal-to-noise ratios can identify the influence of this type of bias. Because we did not observe variation in γ_s between decile subsets, this type of parameter error is unlikely to have significantly corrupted our emission size inference.

We further tested for the presence of a noise bias by comparing the estimates of I_n as a function of pulse strength. Leaked power, for instance, would have introduced a bias that was positively correlated with pulse strength. We detected a slight negative trend in our estimates of I_n with pulse strength; the variation in 〈I_n〉 between the top and bottom deciles (as shown in Figure 3) is approximately 0.15%. If this variation reflects a bias in the estimated noise, then the incurred contribution to γ_s is at approximately the 3σ errors we quote. However, the covariance between the fitted values of γ_s and 〈I_n〉 is marginal, suggesting that this variation in 〈I_n〉 is not responsible for any appreciable bias in the emission size inference. The variation may instead reflect noise in I_n, which weakly anti-correlates with the gated signal-to-noise ratio of the corresponding pulse. Then, for the full pulse ensemble, the error estimates of Section 5.1 are still appropriate, although the segregation into quantiles will somewhat exaggerate this effect.

Overall, these tests suggest that our 8 bit quantization scheme and high-quality backend (GUPPI) effectively controlled any systematic errors in the model parameters. Remarkably, the determination of γ_s is then limited by the Poisson noise in the observed PDF of intensity (constructed from ∼10⁹ samples). Because the transverse size at the source is proportional to $\sqrt{\gamma _{\mathrm{s}}}$, our resolving power scales as the fourth root of the observing duration and bandwidth.

Thus far, we have assumed thin-screen scattering with a square-law phase structure function, as is observationally motivated for the Vela pulsar (Williamson 1972; Komesaroff et al. 1972; Roberts & Ables 1982; Johnston et al. 1998; M. D. Johnson et al. 2012, in preparation). However, even with relaxed scattering assumptions, the same form of Equation (1) still applies, with modified scaling.

For example, the translation to a general power-law structure function $D_\phi (\textbf {x}) \propto | \textbf {x} |^\alpha$ is given by (M. D. Johnson et al. 2012, in preparation)

$$\begin{eqnarray} &&\gamma _{\mathrm{s},\alpha } = \left[ \frac{2^{3-\alpha }}{\alpha \Gamma \left(\frac{\alpha }{2} \right)} \gamma _{\mathrm{s}} \right]^{2/\alpha }\!. \end{eqnarray} \tag{ 9 }$$

Hence, if the phase structure function is Kolmogorov (α = 5/3), then the effects of a small emission size γ_s ≪ 1 will be somewhat more pronounced than for square-law scattering. The inferred transverse sizes will therefore be smaller if we assume α = 5/3 than those if we assume α = 2. Our given limits γ_{s, α = 2} ≲ 10⁻⁴ translate to γ_{s, α = 5/3} ≲ 2.3 × 10⁻⁵, so this alternative scattering environment roughly halves our inferred limits.

The translation from γ_s to a dimensionful emission size also depends on the bulk scattering geometry (see Equation (1)). In this case, the uncertainty in the magnification D/R directly affects the inferred emission size. However, the excellent determination of both the temporal and angular broadening of the Vela pulsar leads to relatively small uncertainties for the magnification (Gwinn et al. 2000), and thus the uncertainties for the inferred emission size are within the standard errors from Poisson noise.

Because γ_s represents a characteristic size of the emission region, effects such as the proper motion of the pulsar can be important. However, the Vela pulsar has a transverse velocity of only ∼60 km s⁻¹ in its local environment (Dodson et al. 2003)—quite slow relative to other young pulsars and the typical neutron star birth velocity (Lyne & Lorimer 1994). Over our 10 ms accumulation time, the pulsar travels only 600 m. Moreover, since size effects are weighted by pulse amplitude, the relevant timescale is the mean width of individual pulses, typically ∼2 ms, contributing an effective displacement of ∼100 m—much smaller than our present detection threshold.

A highly simplified, conventional model describes the geometry of pulsar emission at centimeter and longer wavelengths. Dyks et al. (2004) provide an overview of this model, with a discussion of its assumptions and parameters. We now outline the salient features for completeness. The model assumptions include a dipolar magnetic field, radio-wave emission tangent to field lines, and a single emission altitude for identifiable features of the average-pulse profile. This model follows from theoretical predictions of an ultra-relativistic electron–positron wind, streaming along the open dipolar field lines in the co-rotating frame (Goldreich & Julian 1969; Ruderman & Sutherland 1975; Melrose 2004), so we assume that the emission is beamed tangent to the field lines into an infinitesimal solid angle. In this case, an observer receives radiation from a single point in the pulsar's magnetosphere at each pulse phase. The variations of the spatial displacement of this point with pulse phase then determine the apparent size of the emission region.

Because we merely seek a rough correspondence between the size and altitude of the emission region, we ignore effects of relativistic terms, modifications to the dipolar field structure, and magnetospheric ducting and reprocessing (Blaskiewicz et al. 1991; Phillips 1992; Kapoor & Shukre 1998; Hibschman & Arons 2001; Arons & Barnard 1986; Lai et al. 2001).

A pair of parameters define the orientation of the pulsar: the inclination of the observer relative to the rotational axis (ζ) and the impact angle of the magnetic pole with the line of sight (β). We also employ standard angular coordinates: the longitude or rotational phase of the pulsar (ϕ) and the angle between the point of emission and the magnetic pole (θ_em(ϕ)). For a fixed emission altitude r_em, the spatial displacement of the point of emission ${\bm D}_{\rm em}(\phi)$ then follows from the dipolar geometry: $|{\bm D}_{\rm em}| \approx r_{\rm em} \tan ^{-1}(\theta _{\rm em}/2)$. The cosine formula of spherical trigonometry relates θ_em to pulse longitude ϕ and the orientation parameters {β, ζ}, and the well-known formula for the swing of position angle with ϕ gives the direction of ${\bm D}$ (Komesaroff 1970). Lyne & Graham-Smith (2005, Figure 16.14) and Lorimer & Kramer (2004, Figure 3.4) provide useful figures and discussion.

The Vela pulsar benefits from several features that permit confident inferences of the orientation parameters. Fits to the X-ray emission from its pulsar wind nebula yield ζ ≈ 63° (Ng & Romani 2004, 2008). In addition, the swing of the position angle of polarization across the pulse indicates that β ≈ −6.5° (Radhakrishnan & Cooke 1969; Johnston et al. 2005). Thus, the Vela pulsar is a nearly orthogonal rotator, with the line of sight close to the plane of rotation. Note that the duty cycle of the Vela pulsar is quite small. The full angular extent of the emission is only ρ ≈ 16° (Johnston et al. 2001), and the effective angular extent, expressed as the standard deviation in pulse phase weighted by flux density, is only ρ_eff = 3.0°, from our measurements.

We now relate the effective emission size σ_c to r_em, within the context of this geometrical model. A perfectly orthogonal rotator (β = 0) provides an enlightening demonstration. In this case, the standard deviation of the transverse spatial displacement, weighted by flux density, is σ_eff ≈ ρ_effr_em/3. The factor of 1/3 reflects that the angular difference between the tangent and radial directions of the dipolar field is roughly one-third the colatitude of the tangent direction. The translation to an effective circular emission region σ_c then incurs an additional factor: $\sigma _{\rm eff} = \sqrt{2}\sigma _{\rm c}$. Applying our most stringent limits for the size of the emission region, σ_c < 4 km, we obtain $r_{\rm em} < 3 \sqrt{2} \sigma _c/\rho _{\rm eff} \approx 320\ \mathrm{km}$.

Applying the measured geometry for a nearly orthogonal rotator, appropriately weighted by the individual pulse profiles for the stronger linear polarization at 757.5 MHz, we find r_em ≈ 86σ_c < 340 km. Moreover, if we assume that the phase structure function of the scattering medium is Kolmogorov (see Section 5.3), then we obtain r_em < 170 km.

Our results are comparable with the typical conclusions derived from pulse profile analysis and polarimetry. For example, Kijak & Gil (1997) calculated pulsar emission heights by measuring pulse profile widths at 0.4 and 1.4 GHz and assuming that the emission spanned the open field-line region at a single altitude. For the Vela pulsar at 757.5 MHz, their model predicts r_em ≈ 360 km.

More generally, the emission may arise from field lines that span only a portion of the open field-line region. Researchers typically characterize this possibility with the fractional colatitude of active field lines $f_{\theta } \approx \theta _{\rm max} \sqrt{R_{\rm LC}/r_{\rm em}}$. The inferred emission altitude is quite sensitive to the assumed emission structure, scaling as ∼f⁻²_θ. However, with strict symmetry assumptions, f_θ can be inferred in pulsars with clearly delineated core and cone components (Gupta & Gangadhara 2003; Dyks et al. 2004). Such measurements suggest that f_θ generally lies within the range 0.2 − 0.8 for conal emission.

Alternatively, the aberrational shift in position angle indicates the emission height (Blaskiewicz et al. 1991; Phillips 1992; Hibschman & Arons 2001). This shift is insensitive to f_θ and thus gives a robust r_em without constraining f_θ. For the Vela pulsar, the position angle offset indicates r_em ∼ 100 km (Johnston et al. 2001).

A measurement of σ_c concurrently determines f_θ and r_em, regardless of the profile morphology or the scatter broadening. Thus, it can probe the emission structure at low frequencies, where both the size of the emitting region and the effects on the scintillation are most pronounced.

The sizes that we infer for the emission region are substantially smaller than those found by the previous scintillation analyses discussed in Section 1.2; however, our lack of resolution in pulse phase precludes a direct comparison. For example, if the emitting regions are largest when the pulsar is weak, then the effects of emission size will be diminished. Indeed, the phase-resolved analysis at 1.66 GHz suggests precisely this scenario. Furthermore, our observations span lower frequencies than previous analyses. Although the radius-to-frequency mapping (Ruderman & Sutherland 1975; Cordes 1978) suggests that the emission height decreases with increasing frequency, the Vela pulsar possibly does not conform to this picture. Alternatively, the pulse profile shows considerable intrinsic evolution at 1 GHz, and the disparity between the measured emission size may indicate the introduction of a new, offset emission component at these higher frequencies. Pending observations should resolve this intriguing discrepancy.