SIZE OF THE VELA PULSAR'S EMISSION REGION AT 18 cm WAVELENGTH

The data were recorded in single sidebands; thus, the spectra contained 8192 channels, each with bandwidth 1.95 kHz. The cross-power spectra are complex. The phase includes instrumental effects, primarily observational and instrumental delays and rates (Thompson et al. 1986), and effects of scintillation (Desai et al. 1992).

We found that the time period from 19:10 to 21:13 UT on the Mopra–Tidbinbilla baseline contained data with uniform high gain and low noise, absence of interference or gaps in correlation, and little change in length and orientation of the baseline. The coordinates projected perpendicular to the source direction vary over the range of (u, v, w) from (422.8, 1416.9, −392.67) to (733.5, 1193.7, −614.5) μs. The projected baseline length is approximately 428 km.

We edited the data to remove times and channels with interference or corrupt recording. We identified channels that showed evidence of interference such as high amplitude or high noise. This amounted to 27 channels for IF1 and 12 channels for IF2. We also identified time records with excessively low or high amplitudes, or with low correlation amplitudes, and removed those.

We corrected for average delay and rate by fringe fitting (see Thompson et al. 1986). Delay represents a phase slope with frequency, and rate a phase slope with time. We fringed the central 7168 channels of gate 2, leaving "guard zones" 512 channels wide on each end, for the regions where passband gain rolled off and instrumental phase varied most rapidly. We formed a two-dimensional discrete Fourier transform to find the fringe rate for each eight-sample (16 s) time interval, using the traditional "fringe" algorithm (Thompson et al. 1986).

We then used the fringe rate and delay from gate 2 to remove the corresponding phase slopes from the other "slaved" pulsar gates. In tests with pairs of gates that contained strong signal, we verified that delay and rate were the same for all gates, to the accuracy permitted by signal-to-noise ratio. However, we found that the residual interferometer phase depended on pulsar gate. We therefore used the phase estimated from each gate to correct that gate. For this paper, the primary purpose of fringe fitting was to remove instrumental phase, leaving only the effects of scintillation and those of statistical noise in the data.

The calibrated data take the form of complex cross-power spectra, sampled as a function of time, in six pulse gates and two IFs. These data were gathered for several baselines, as discussed in Gwinn et al. (2012). In this paper, we focus on the relatively short Mopra–Tidbinbilla baseline. As an example, Figure 2 shows the real part of the visibility, for a short span of time and frequency in three gates. The relative variation in amplitude between the gates has been removed by calibration, using the average amplitude over the spectral range for the time span of all the data, as shown in Figure 1. The ratio of calibration factors was 24:36:13 for gates 1, 2, and 3. The observed spectra differ because of the effects of noise and because of variations in the relative amplitudes of individual pulses at the different pulse phases (Krishnamohan & Downs 1983; Johnston et al. 2001; Kramer et al. 2002; Gwinn et al. 2012; Johnson & Gwinn 2012). The scattering medium is expected not to change between gates, which after all integrate over the same time interval; thus, the scattering pattern should be the same, allowing for variations in noise and pulse-to-pulse variability. An interesting question is whether differences might additionally reflect changes in the structure of the pulsar's emission region between gates. Because the effects of noise and amplitude variations are random, this issue can only be treated statistically, using the correct descriptions of noise and amplitude variation. Complicating this comparison is the fact that effects of source size are largest when scintillation leads to small flux densities (Gwinn et al. 1998). The remainder of this paper makes such a statistical comparison.

Zoom In Zoom Out Reset image size

We reduce the observational data to a histogram of measurements of the real part ${\mathcal P}_N$ and a histogram weighted by the mean square imaginary part ${\mathcal Q}_N$. The subscript "N" denotes that these distributions reflect measured histograms rather than model distributions. Mathematically, these histograms correspond to sums over the observed points V(ν, t), restricted to bins of width w about the bin centers, at real part X_k:

$$\begin{eqnarray} {\mathcal P_N}(X_k) &=& \sum _{\nu,t} 1 \quad{\rm for\ } X_k - w/2 \le {\rm Re}[V(\nu,t)] < X_k +w/2 \nonumber\\ {\mathcal Q_N}(X_k) &=& \sum _{\nu,t} {\rm Im}[V(\nu,t)]^2 \nonumber\\ && \quad{\rm for\ } X_k - w/2 \le {\rm Re}[V(\nu,t)] < X_k +w/2. \end{eqnarray} \tag{ 1 }$$

We sought to make the histograms from sufficiently narrow spectral ranges so that the amplitude did not vary greatly across the spectrum, while including enough points for robust statistics. Most of the spectral variation arises from the shape of the pulse and pulse dispersion, as Figure 1 and Figure 1 of Gwinn et al. (2012) suggest. We adopted spectral ranges of 1024 channels within each 8192-channel spectrum from each gate. We dropped the first and the last 1024 channels, to avoid effects of gain rolloff and phase variations near the edges of the observed band. Consequently, our data are indexed by IF number (1 or 2), by gate (1 through 5), and by channel range in increments of 1024, beginning at 1024 to 6144. Figure 1 shows the centroids of the 1024-channel ranges for gate 2.

Figure 3 shows examples of the measured histograms, for IF1, gate 1, and channels 4096–5120. We use these distributions, and others like them for other spectral and gate ranges, to fit for the parameters of our theoretical models, as discussed in Section 3. We used intervals of w = 0.0002 for the histograms in this paper. Wider bins would average over the smooth, rapid variations of the distribution among bins.

Zoom In Zoom Out Reset image size

Projection has a number of advantages. One-dimensional distributions are much easier to visualize and fit than two-dimensional distributions. Projection increases the number of samples per cell, reducing Poisson noise. For the short Mopra–Tidbinbilla baseline, amplitude variations from scintillation affect primarily the real part, broadening it along the real axis. Phase variation from scintillation affects the imaginary part, and noise affects both real and imaginary parts. Consequently, the second moment of the imaginary part in ${\mathcal Q}$ provides a useful constraint on noise and the effects of finite baseline length. For a short baseline, visibility is concentrated near the real axis, so most of the information is contained in these two distributions.

In the absence of noise or amplitude variations, the interferometric visibility for a scintillating point source is the product of complex Gaussian random variables. The resulting distribution of visibility peaks sharply at the origin; indeed, the distribution is singular at that point (Gwinn 2001). Furthermore, the distribution has strong exponential wings that extend to large values. This large dynamic range complicates evaluation.

For an extended source, also in the absence of noise or amplitude variations, the peak of the distribution softens and shifts toward greater real part, relative to that for a pointlike source. However, the variance decreases in both the real and imaginary directions. Thus, estimation of source size requires accurate calculation of the structure near the origin on small scales and of the behavior of the wings at large values. Figure 4 shows distributions for sources with zero size and with small size, on a short baseline, calculated using the results of this section; the parameters are comparable to those we find for the Vela pulsar in Section 4 below. We discuss effects of noise in Section 3.3 and effects of amplitude variability in Section 3.4.

Zoom In Zoom Out Reset image size

Explicitly, the normalized distribution of visibility for a pointlike, scintillating source, without noise and with constant amplitude, is given by (Gwinn 2001)

$$\begin{eqnarray} P(V) &=& \frac{2}{\pi \kappa _0^2}\frac{1}{(1-\rho ^2)}K_0\left(\frac{2}{(1-\rho ^2)}\frac{|V|}{\kappa _0}\right) \nonumber\\ && \times {\rm exp}\left(\frac{2\rho }{(1-\rho ^2)}\frac{\mathrm{Re}[V]}{\kappa _0}\right). \end{eqnarray} \tag{ 2 }$$

Here, ρ is the normalized covariance (that is, the normalized interferometric visibility of the source), and κ₀ is the scale. This distribution is that of a product, xy*, of correlated, complex Gaussian random variables. It describes the distribution of visibility, after averaging over an ensemble of electric-field measurements; it does not include effects of background noise, self-noise, or intrinsic source variability.

For a scintillating source viewed through an isotropic scattering screen (Gwinn 2001, Equations (2), (4), and (6))

$$\begin{equation} \rho = \exp \left\lbrace -{\frac{1}{2}} {\frac{(2 \pi)^2}{8 \ln 2}} {\frac{\theta _{\rm H}^2 b^2}{\lambda ^2}} \right\rbrace. \end{equation} \tag{ 3 }$$

Here, θ_H is the FWHM of the scattering disk, b is the baseline length, and λ is the observing wavelength. The mean correlated flux density of a pointlike source is, therefore,

$$\begin{equation} \langle V\rangle = \kappa _0 \rho = \kappa _0 \exp \left\lbrace -{\frac{1}{2}} (k \theta b)^2 \right\rbrace. \end{equation} \tag{ 4 }$$

Here, the wavenumber is k = 2π/λ, and the angular broadening is $\theta = \theta _{\rm H}/\sqrt{8\ln 2}$. Angular brackets 〈...〉 denote an average over many scintillations.

For short baselines, ρ ≈ 1, and the distribution of visibility is concentrated near the positive real axis; the imaginary part for a given value of Re[V] has variance that scales proportionately with Re[V]. For longer baselines, ρ → 0, and the distribution is circularly symmetric about the origin. Figure 2 of Gwinn (2001) shows examples.

At zero baseline, ρ ≡ 1, visibility V becomes intensity I, and P(V) becomes the well-known exponential distribution of intensity for a scintillating point source:

$$\begin{equation} P(I) = {\frac{1}{I_0}} \exp \lbrace - I/I_0\rbrace. \end{equation} \tag{ 5 }$$

To connect the two distributions, note that a "zero-baseline interferometer" can, in principle, measure complex visibility and will have complex noise; however, the visibility of the source is confined to the real axis, even with scintillation. Thus, Equation (5) is the projection of Equation (2) onto the positive real axis, in the limit ρ → 1.

The distribution of visibility for an extended scintillating source, without noise and with constant amplitude, is the convolution of a number of copies of Equation (2) (Gwinn 2001). For a small, scintillating source, the distribution is the convolution of three such copies, two of them with scales κ related to the size of the source along two orthogonal axes ξ, η:

$$\begin{eqnarray} \kappa _{\xi } &=& \kappa _0 (k M \theta \sigma _{\xi })^2 \nonumber\\ \kappa _{\eta } &=& \kappa _0 (k M \theta \sigma _{\eta })^2. \end{eqnarray} \tag{ 6 }$$

Here, M = D/R is the effective magnification of the scattering disk, where D is the distance from observer to scatterer, and R is the distance from scatterer to source. The dimensions of the source are σ_ξ, σ_η. We assume that the source is small in the sense that (kMθσ_ξ), (kMθσ_η) ≪ 1. If this condition does not hold, then additional terms are important; the convolution involves more distributions. The covariances, corresponding to ρ, change for the subsidiary distributions as well:

$$\begin{eqnarray} \rho _{\xi } &=& (1-(b_{\xi } k \theta)^2) \exp \left\lbrace -{\frac{1}{2}} (k \theta |b|)^2 \right\rbrace \nonumber\\ \rho _{\eta } &=& (1-(b_{\eta } k \theta)^2) \exp \left\lbrace -{\frac{1}{2}} (k \theta |b|)^2 \right\rbrace. \end{eqnarray} \tag{ 7 }$$

These are nearly equal to ρ for short baselines.

We have not found a simple analytic expression for the result of the convolution for ρ < 1. Moreover, the convolution is challenging to reproduce numerically, because P(V) has both a sharp peak at the origin and high skirts that extend to large |V|. Our strategy is therefore to proceed as far as possible via analytic calculations and then evaluate the remaining integrals numerically.

To achieve our first numerical reduction, note that we can easily convolve visibilities that are drawn from identical distributions. The average of N such visibilities is distributed according to

$$\begin{eqnarray} P_N(V) &=& \frac{2}{\pi } \frac{N^{N+1}}{(N-1)!} \frac{1}{(1-\rho ^2)\kappa ^2} \left(\frac{|V|}{\kappa } \right)^{N-1} K_{N-1} \nonumber\\ && \times \left(\frac{2 N}{1-\rho ^2} \frac{|V|}{\kappa } \right) \exp \left(\frac{2 N \rho }{1-\rho ^2} \frac{\mathrm{Re}[V]}{\kappa } \right). \end{eqnarray} \tag{ 8 }$$

We can use this result to calculate the convolution of visibilities drawn from different distributions if we first use Feynman parameters (Srednicki 2007) in Fourier space to symmetrize the corresponding conjugate product of functions. Because the visibility is complex, convolution of N visibilities requires a 2(N − 1)-dimensional integral. However, to symmetrize the convolution requires a single Feynman parameter for each visibility. The Feynman parameters also have an overall δ-function constraint, so the convolution is reduced to an (N − 1)-dimensional integral.

For a small, circular source, the observed distribution is a convolution of three visibilities, two of which are drawn from identical distributions. These two distributions are parameterized by Equations (6) and (7) with σ_ξ = σ_η and b_ξkθ ≪ 1, b_ηkθ ≪ 1. This symmetry allows elimination of an additional degree of freedom, so the original four-dimensional visibility integral is reduced to a one-dimensional integral. Explicitly, we assume that κ₀, ρ₀, κ₁, and ρ₁ are arbitrary, but κ₂ = κ₁ and ρ₂ = ρ₁. For convenience, we also introduce parameters a_i ≡ κ_i(1 − ρ_i)/2 and p_i ≡ 2[πκ²(1 − ρ²)]⁻¹. Then, the distribution of visibility is given by

$$\begin{eqnarray} P\left(V\right) &\equiv & \left(P_0\ast P_1 \ast P_2\right)\left(V\right) \nonumber\\ &=& \pi ^{2} p_0 p_1^2 a_0^2 a_1^4 |V|^2 \int _0^1 ds \ (1-s)f_1(s)K_{2} \nonumber\\ && \times (f_2(s)|V|)e^{f_3(s) {\rm Re}[V]}, \end{eqnarray} \tag{ 9 }$$

where we have defined, for convenience, the three functions

$$\begin{eqnarray} f_1(s) &\equiv & \big[\big(a_0^2 s + a_1^2(1-s)\big)\nonumber\\ && \times \big(a_0^2 s + (1-s)\big(a_1^2 - s [a_1\rho _0-a_0\rho _1]^2\big)\big)\big]^{-1} \nonumber\\ f_2(s) &\equiv & \frac{\sqrt{a_0^2 s + (1-s)\left(a_1^2 - s \left[a_1\rho _0-a_0\rho _1\right]^2\right)}}{a_0^2 s + a_1^2(1-s)} \nonumber\\ f_3(s) &\equiv & \frac{a_0 \rho _0 s + a_1 \rho _1 (1-s)}{a_0^2s+a_1^2(1-s)}. \end{eqnarray} \tag{ 10 }$$

This one-dimensional integral is suitable for efficient numerical evaluation.

We project the model distributions, including effects of noise and amplitude variations, onto the real axis. We calculate the projected probability density and the summed squared imaginary part:

$$\begin{eqnarray} \mathcal {P}(X_k) &=& \int _{X_k-w/2}^{X_k+w/2} d{\rm Re}[V]\; \int _{-\infty }^\infty d{\rm Im}[V]\; P (V) \nonumber\\ \mathcal {Q}(X_k) &=& \int _{X_k-w/2}^{X_k+w/2} d{\rm Re}[V]\; \int _{-\infty }^\infty d {\rm Im}[V]\; {\rm Im}[V]^2\; P (V). \quad\quad \end{eqnarray} \tag{ 13 }$$

These are scaled models for the projections ${\mathcal P_N}$, ${\mathcal Q_N}$ of the data as described in Section 2.3.4 and plotted in Figure 3. Note that $\mathcal {Q}(X_k)$, the mean square imaginary part in each bin, is a second moment of Im[V] and so is expected to be noisier than $\mathcal {P}(X_k)$, the zeroth moment of Im[V]. The higher noise for ${\mathcal Q}$ in Figure 3 reflects this.

Effects of noise and amplitude variations must be calculated numerically for the full two-dimensional distribution and then projected. Noise changes the value of each measurement and thus spreads the corresponding probability distribution function over a surrounding region. It smoothes a spike, corresponding to a single deterministic value, into an elliptical Gaussian distribution centered at that point, with variances given by the noise polynomial, Equation (12), evaluated at that point. The effect of noise on the distribution of visibility thus resembles a convolution of the probability distribution of the scintillating source, with the Gaussian distribution of noise. It is not a convolution, because the noise depends on the signal: this operation is sometimes called a "convolution with varying kernel." The Gaussian noise, with dependence on signal strength, is the kernel in this case. Because the kernel varies, we must project after convolution. Because the distribution is relatively concentrated along the real axis, the notion of convolution after projection has intuitive value, but we do not use this construction in our analysis for ρ < 1; we calculate the full variation of the noise kernel.

The effect of amplitude variations on timescales longer than the integration time is to superpose distributions with different κ₀ but the same noise polynomial. We implement this effect numerically, by averaging a number of distributions with varying κ₀. Fortunately, integration over 2 s reduces the intrinsic amplitude variations, as discussed below.

A sample calculation illustrates our method and shows the origin of the W-shaped signature of finite source size in the difference of projected distributions $\mathcal {P}$. Figure 5 shows how the best-fitting zero-size and finite-size models differ, for the span of data shown in Figure 4. The upper panel of the figure shows the projected distribution $\mathcal {P}$ without noise or amplitude variations. For zero size, the projected distribution shows a sharp cusp at Re[V] = 0, with an exponential decline along the positive real axis of V and a more sharply declining exponential along the negative real axis. The negative-side exponential declines so sharply because ρ = 0.986 ≈ 1 in this example. For a model with finite size, with the same normalization and mean, the bulk of the projected distribution is shifted toward the positive real axis, with a rounded peak. For the same mean amplitude, the peak has a narrower spread.

Zoom In Zoom Out Reset image size

Note that the point-source model has greater probability at large and small Re[V] than does the extended-source model, but the extended model has greater probability in between. This behavior arises because we require that the two distributions have the same mean and normalization. The shift of the maximum toward greater Re[V] for the finite-size model then requires a reduction of κ₀ (see Gwinn 2001, Equation (31)). The largest exponential scale, κ₀, dominates the behavior of the distribution away from the origin, so that a finite-size source concentrates the distribution of visibility.

Noise broadens the distributions, as the middle panel in Figure 5 illustrates. The best-fitting noise parameters are {b₀, b₁, b₂} = {0.000062, 0.0037, 0.078} for the zero-size model and {0.000058, 0.0045, 0.181} for finite-size. The noise at V = 0 has standard deviation $\sqrt{b_0}$, or approximately 0.008. This value is comparable to the width of the distribution, so that the details apparent in the top panel of the figure are blurred. Moreover, the noise increases away from V = 0, as given by the higher-order coefficients.

Intrinsic amplitude modulation also changes the distributions. The mean square of the intrinsic amplitude modulation is m²_s = 0.009 for the best-fitting zero-size model and m²_s = 0.007 for the best-fitting finite-size model, for our example of spectral and pulse-gate range. In both cases the degree of modulation is small as expected after integrating 22 or 23 pulses.

The lower panel in Figure 5 shows the difference of the zero-size and finite-size models. Even after including effects of noise and variability, the underlying features of the two distributions persist: the finite-size distribution has smaller probability density for Re[V] < 0, because of its lesser density near Re[V] = 0 in the top panel, and at large Re[V], because of its more rapid decline there. It has greater density near the mean amplitude. Thus, we expect the signature of a finite-size source to be a W-shaped difference of ${\mathcal P}$ for the best finite-size model from the best zero-size model, after effects of noise and intrinsic amplitude variations have been included.

The situation for the distribution of mean square imaginary part $\mathcal Q$ is somewhat similar to that for $\mathcal P$, but the fractional difference between models is less. Indeed, nearly all of the spread in the imaginary part results from noise, for values of ρ of interest for this baseline, so that the distribution $\mathcal Q$ functions more as a constraint on the noise model than as a carrier of information about source size. We present plots for $\mathcal Q$ below, in Section 4.1.4.

Our calculation of the model distribution involves four nested iterative loops. We calculate $\mathcal {P}$ and $\mathcal {Q}$ in parallel, on a grid of points in the real part of visibility, X_k. Each histogram bin has width w = 2 × 10⁻⁴. This size is narrow enough to track the behavior of the distribution accurately, but wide enough to contain enough visibilities so that the Poisson noise does not obscure the model differences.

At the lowest level, we integrate over s to find the probability density P(V) at a point V in the complex plane, using Equation (9). We integrate this expression iteratively using Simpson's rule (Press et al. 2007) with 2 × 3^N segments of equal width on the Nth iteration. The integration terminates when successive iterations have a fractional difference of less than 0.1%, or N > 10. Because the error bound for this integration scheme goes as N⁻⁴, this criterion yields an expected accuracy of ∼0.001%—a small fraction of the magnitude of measurable source size effects.

We broaden the probability density at each grid point by its corresponding elliptical Gaussian distribution of noise. The polynomial coefficients {b₀, b₁, b₂} (Equations (12)) parameterize the noise, and the complex visibility at the grid point determines the scale and orientation of the ellipse. We use the analytic expressions for the projection of a Gaussian distribution and the projected second moment of a Gaussian distribution to evaluate the contribution of each point, with noise, to the projections $\mathcal {P}$ and $\mathcal {Q}$ on our grid of points.

At the second level, we integrate the contributions to $\mathcal {P}$ and $\mathcal {Q}$ over the imaginary part of V. We integrate from the real axis out to the 4σ standard deviation for the distribution in the absence of size effects, as calculated from the second moment of Equation (2). Again, we use Simpson's rule iteratively with 2 × 3^N segments of equal width but now require a fractional change of less than 10⁻⁸ in $\mathcal {P}(X_k)$ between two successive iterations. We monitor $\mathcal {P}(X_k)$ because it is more sensitive than $\mathcal {Q}(X_k)$ to behavior near the real axis, where the distribution varies most rapidly, so that $\mathcal {P}(X_k)$ converges more slowly.

At the third level, we integrate the values of $\mathcal {P}$ and $\mathcal {Q}$ over sub-bins within each histogram bin along the real axis. This step accounts for effects of finite histogram bin width in discrete representations of a probability distribution function. It is most important when there is large-amplitude nonlinear structure, as exhibited by the rapid rise of P at the origin and the cusp of the noise-free distribution near the origin. In the limit of small bin width, w, the histogram representation of a function f(x) will systematically overestimate the function at x₀ by w²f''(x₀)/24. We set the number of integrated sub-bins to 3; comparison with computations using five and seven sub-bins for the fits that obtained smaller sizes yielded insignificant differences.

At the fourth and highest level, we account for the intrinsic amplitude variability of our 2 s integrations. We calculate each visibility distribution five times, with different source amplitudes, as given by a Gaussian distribution centered at 1, scaled by an amplitude-variation parameter. We then superpose these five distributions to produce a distribution including effects of variations of amplitude of the source on the final distribution. The procedure is that used by Gwinn et al. (2000) to describe the distribution of scintillation in the presence of amplitude variations. We found no difference when using more finely divided distributions of source amplitude. Effects of amplitude variations would be significant if the amplitude approached zero, but because Vela does not null and variations of individual pulses are only weakly correlated, amplitude variations after 2 s integrations are small.

For zero baseline, ρ = 1, the distribution of visibility is a sum of exponentials on the positive real axis and zero elsewhere (Gwinn et al. 1998). This form obviates the lowest two levels of the fourfold nested loops above, so the model is faster to calculate and thus is useful for exploration of parameter space. For this exponential model, we assumed a circular source, as we do for our ρ < 1 calculation. Tests showed that the two calculations were equivalent for ρ → 1: a useful check of our numerical integration.

Seven parameters describe the model. Three of these are the coefficients of the noise polynomial, Equation (12): {b₀, b₁, b₂}. The fourth and fifth parameters describe the distribution in the absence of noise and amplitude variations: the amplitude scale κ₀ and the normalized mean interferometric visibility ρ. The scale κ₀ is proportional to flux density. The mean interferometric visibility ρ describes the baseline length and angular size of the scattering disk (Equation (3)). We do not provide a parameter for overall normalization: we constrain the normalization of the model distribution to equal that of the observed distribution. The sixth parameter describes intrinsic amplitude variations (m²_s) = (δI²/〈I〉²) on timescales longer than an integration time; because of such variations, the observed distribution involves a superposition of different values of κ₀. The parameter (m²_s) gives the variance of this distribution of amplitude fluctuations, normalized by the mean amplitude. The seventh parameter gives the size of the source, usually expressed as κ_ratio = κ₁/κ₀ = (kMθσ)². In principle, additional parameters give correlations μ, ν for the subsidiary distributions due to source structure and the intrinsic elongation of the source (Gwinn 2001, Equation (28)); however, these are expected to be nearly equal to ρ for a short baseline.

Among the effects we did not include in our model are elongation of the scattering disk and elongation of the source. For a constant orientation of the baseline, elongation of the scattering disk has no effect on the distribution of visibility, for the appropriate parameters ρ, κ₀, and κ_ratio (Gwinn 2001). The rotation of our baseline was small during the test interval considered here. Effects of elongation of the source are more subtle but appear in our simple model only when the source is nearly resolved by the scattering disk, or the baseline is relatively long.

Note that ρ is a property of scattering and should be constant over the entire pulse. Finite source size will change it slightly (Gwinn 2001, Equation (18)). A source with spatial coherence, over scales comparable to those detectable via scintillation, can also change ρ by illuminating only part of the scattering disk (Gwinn et al. 1998).

Our fits minimize the weighted mean square difference between the data and the model. Individual bins in the histogram ${\mathcal P}$ measure a number of counts N and are expected to display Poisson statistics. For ${\mathcal Q}$ the values are squares drawn from a nearly Gaussian distribution and should display statistics analogous to that of the Dicke equation, Equation (11). In both cases, the distribution of noise is nearly Gaussian, with variance proportional to the bin value. We adopt this weighting above a threshold of N = 100 counts, with constant weighting below that value so as to reduce the influence of bins with zero or few counts.

The two histograms ${\mathcal P}$ and ${\mathcal Q}$ have different dimensions and different vertical scales. We seek to combine the residuals for both into one figure of merit. A reasonable conversion factor is simply the quotient of the integrated areas under the distributions, which has the correct dimensions. However, we expect the second moment to be noisier, as Figure 3 shows. To quantify this noise in ${\mathcal Q}$, we fit the distributions with simple models to find smoothed average values and then difference adjacent points to determine the noise. We find that the standard deviation of noise in ${\mathcal Q}$ determined in this way is typically three times that in ${\mathcal P}$. We therefore weight the residuals in ${\mathcal Q}$ by 1/9 of the ratio of the areas under the two histograms. This has the effect of making the mean square residual roughly equal for the two, for our models.

We find that the mean square residual is close to 1 for ${\mathcal P}$ with this weighting, for the fits described in Section 4 below. This indicates that the statistics are indeed Poisson: number of counts limits the accuracy of the histogram values in our narrow bins. As expected for our scaled weighting, the mean square residual is close to 1 for ${\mathcal Q}$, as well.

As a further test, we also fit to ${\mathcal P}$ and ${\mathcal Q}$ independently. We find that fits to only ${\mathcal P}$ found minima close to those found using both ${\mathcal P}$ and ${\mathcal Q}$, with equal or sometimes larger intrinsic size of the source. Fits to only ${\mathcal Q}$ typically do not converge; apparently this distribution does not contain enough information to determine model parameters independently. Fits with uniform rather than Poisson weighting yield similar results for size but tend to converge less quickly.

In order to examine parameter space over large scales and to provide initial parameters for our fits for arbitrary ρ, we performed a grid search using an exponential model, as described in Section 3.6.2. Because such a model can be calculated quickly, a grid search can survey large regions of parameter space efficiently. In our grid searches, we fit for three parameters: two noise coefficients b₀, b₁ and the amplitude κ₀. We searched a grid of parameters in the noise coefficient b₂ and the size parameter (kMθσ)². We demanded that b₂ > 0.030, as dictated by the number of samples correlated (Gwinn et al. 2012), and searched 0.030 < b₂ < 0.4. We ignored effects of longer-term intrinsic amplitude fluctuations between integrations: m_s ≡ 0.

The grid search indicated that the summed, weighted squared residuals vary smoothly over the parameter space defined by the five varying parameters. We found only one minimum in all IF, spectral, and gate ranges. The best-fitting size of the source σ defined by this minimum was usually significantly different from 0, with a reduction in the weighted residuals from the best-fitting zero-size model comparable to that found from the more sophisticated fits discussed below. The size parameter (kMθσ)² usually agreed to within 30% for the two kinds of fits. The amplitude-modulation parameter m_s is responsible for much of the difference in size. The minimum found by the grid search provided a useful starting point for the much-slower fits with ρ < 1 and m_s > 0, discussed in Section 4 below.

We use calculations of ${\mathcal P}$ and ${\mathcal Q}$ as described in Section 3.5 to fit for model parameters with finite size and ρ < 1, using the Levenberg–Marquardt algorithm (Press et al. 2007). We fit for the six parameters described in Section 3.7.1, using weighting described in Section 3.7.2. We initialize these parameters using the results of the grid search with an exponential model.

We expect ρ = 0.986, based on the parameters reported by Gwinn et al. (1997): angular broadening of (3.3 × 2.2) mas (FWHM intensity), with the major axis at position angle 92°, in observations at wavelength λ = 13 cm. We scale these to our observing frequency by θ∝ν² and use the length and orientation of our baseline as described in Section 2.3.1 to find ρ, using Equation (28) of Gwinn (2001). In tests, we found that the mean square residual for this data set was not particularly sensitive to ρ, for 0.85 < ρ < 1, so that our observations do not provide a good way to fit for ρ: the baseline is too short to provide much information. Our results for the fitted size, in particular, were not sensitive to ρ within this range: using ρ = 0.92 rather than ρ = 0.986 changed the inferred size parameter (kMθσ) by less than 10%. In several frequency ranges where the pulse was strong (channels 1024 through 4096 of gate 2 of both IFs), fits including ρ strongly favored ρ > 0.9, with ρ = 0.92 sometimes favored within that range. We expect ρ to be almost the same for all of the data, in all channels and gates. Differences from changes in frequency were too small to detect. Effects of source size are expected to be small for our short baseline. Indeed, we expect from the moments of the distribution that ρ and size parameter (kMθσ)² should have little covariance (M. D. Johnson & C. R. Gwinn 2012, in preparation). We adopt ρ = 0.986, noting that revisions may cause small changes in the inferred size.

The lower panels of Figure 3 show the residuals to a fit with zero source size, in histogram form. Structure remains in the residual histogram, apparent as variations near zero visibility. An independent fit, to the same data but with a finite source size, removes most of these variations. To illustrate this, we display the difference of the finite-size model from the zero-size model, as smooth curves. The difference curve for ${\mathcal P}$ is the same as that shown in the lower panel of Figure 5. Like the residuals, this curve shows the W-shaped signature of finite source size.

The mean squared residuals in ${\mathcal P}$, after fitting for finite source size, are approximately those expected for Poisson-noise-dominated errors. The data contain 437 bins with more than 100 counts. The sum of the squared Poisson-weighted residuals, after fitting, is 434. Thus, the mean square errors are approximately as expected. The fit to a model for a source with zero size leads to a sum of the squared Poisson-weighted residuals of 621, significantly greater.

For ${\mathcal Q}$, the sum of the squared residuals, after correction by the relative areas of ${\mathcal P}$ and ${\mathcal Q}$, and the factor of three from comparison of differences is 705 for a finite-size model and 785 for a zero-size model. The residual errors in ${\mathcal Q}$ appear noiselike in the plot, suggesting that the factor of three estimated from differences in other spectral and gate samples may be low for this sample of data. Nevertheless, the finite-size model improves the fit for ${\mathcal Q}$.

The reduction in summed mean-square-weighted residuals for ${\mathcal P}$ and ${\mathcal Q}$ together is 19%. This difference is highly statistically significant at more than the 40σ level, for our sample size and number of parameters, according to the F-ratio test (Bevington & Robinson 2003). At this high level, finite sampling limitations are unlikely to dominate errors in our estimates of emission size. However, this figure demonstrates that the signature of finite size appears in the data, as inspection of the figures suggests. The best-fitting size parameter is (kMθσ)² = 0.0423, or scaled size (kMθσ) = 0.21. This corresponds to a source size of approximately 180 km (standard deviation of the Gaussian distribution), or to approximately 420 km for the FWHM, as we discuss in Section 4.1.3 below.

Standard errors are small, because of the high significance of the fits as expressed by the F-ratio test. We present the best-fitting parameters for this IF, gate, and channel range in Table 1. These values are typical for our fits. As we discuss in Sections 4.1.2 and 4.2 below, errors in the estimated size are probably dominated by systematic errors.

We fit our model to spectral ranges of 1024 channels in each of the five gates, for both IFs. Our fits indicate that the pulsar has a rather large size at the start of the pulse, decreases in size over the first half of the pulse, and then increases in size again. Figure 6 shows the results of fits, giving the scaled size (kMθσ) and mean amplitude (κ₀ + 2κ₁) as a function of pulse phase. Note that linear diameter is proportional to the square root of the fitted size parameter (kMθσ)². All of the fits are independent; each point represents an independent sample of the pulsar's emission, fit completely independently using a priori parameters from independent grid searches (Section 3.7.3). The two IFs are shown by different symbols (crosses for IF1 and circles for IF2), which represent distinct frequency ranges and thus completely different scintillation patterns. At some phases, points from different gates overlap; at these points, a higher-frequency range in one gate coincides with a lower-frequency range in the previous gate. Quite often the noise parameters for these different samples are quite different; nevertheless, size and amplitude track one another with pulse phase, despite the disparate origin of the samples. Some of the fits do sample identical samples of the diffraction pattern from interstellar scattering, as Figure 2 suggests: however, these often lead to completely different results, indicating that the results arise from the pulsar, rather than from scattering. For example, data that include the three panels shown in Figure 1 yield scaled sizes of (kMθσ) = 0.08, 0.09, and 0.40. (Of course, the gray scale in Figure 2 emphasizes peaks of visibility, whereas the most critical information on source size is found near zero visibility, as Figure 5 shows.) The dependence of fitted size on pulse phase, rather than on the frequency or details of the scintillation pattern, suggests that the effect arises from the pulsar rather than the scintillation pattern.

Zoom In Zoom Out Reset image size

Our sample comprises 60 independent fits, over the two IFs, five gates, and six spectral ranges. A fit with finite size yielded significantly smaller residual for most of these: for the 43 fits with fitted size parameters of (kMθσ)² > 0.017, the reduction of the residuals was greater than 3%. This is significant at better than the 5σ level according to the F-ratio test. Unsurprisingly, gate and spectral ranges with small fitted sizes and regions with small amplitudes at the beginning and end of the pulse show the least significant results. We do not display statistical errors from the fits on this figure; they are smaller than the plotting symbols.

The scatter of independent measurements at a given pulse phase provides a measure of the precision of our results. The residuals are dominated by variations among nearby bins, consistent with Poisson noise, whereas the models predict rather more slowly varying differences across the histograms for ${\mathcal P}$ and ${\mathcal Q}$, as Figure 5 would suggest. We made the bins narrow so as to follow the rapid variation of the histograms with Re[V], as shown in Figure 3; we cannot broaden the bins without destroying this important information by averaging. We discuss averaging of the post-fit residuals over bins to visualize effects of the model in Section 4.1.4 below.

The angular width of the scattering disk θ is important for the conversion of the fitted size parameter (kMθσ)² to σ in km. The angular width appears directly in this expression and also affects the inferred magnification M as we discuss below, so that the linear size of the emission region σ depends strongly on θ. Unfortunately, our data do not provide a good determination of ρ, as discussed in Section 3.7.5 above, and so do not allow accurate determination of θ. Consequently, we adopt a simple conversion based on previous work, while noting that more accurate measurements of θ and M may revise the inferred conversion of (kMθσ) to σ, in kilometers.

Gwinn et al. (1997) found angular broadening of (3.3 × 2.2) mas (FWHM intensity), with the major axis at position angle 92°, in observations at wavelength λ = 13 cm. We adopt this value here. For simplicity in converting the observed (kMθσ) to size σ, we assume a circular scattering disk with FWHM of θ = 5.2 mas; this is the mean square of the major and minor axes. For the scintillation bandwidth, we adopt the value of Δν_s = 15 kHz from Section 4.2.4 below. With a distance to the pulsar of D = 287 pc, as measured from the parallax of the pulsar (Dodson et al. 2003), the combination of angular broadening and scintillation bandwidth yields a magnification of M = D/R = 3.1 (Gwinn et al. 1993). Consequently, the standard deviation of the Gaussian model for the source is σ_θ = (kMθσ) × 859 km. In this expression, the subscript θ serves to indicate the dependence of the inferred source diameter on the angular size of the scattering disk and the magnification. The FWHM diameter of the fitted model is $\sqrt{8 \ln 2} \sigma _{\theta }$. The right-hand vertical axis on the lower panel in Figure 6 reflects these conversion factors. The size is as large as hundreds of km, expressed as standard deviation of a circular Gaussian distribution.

Changes in the assumed parameters for angular broadening will modify the inferred emission size. The angular scale of the scattering disk θ sets the scale of the scintillation pattern at the source and thus the resolution of the lens, and θ sets the location of the scattering material along the line of sight and so the magnification factor M. If the value assumed for θ changes, the inferred size σ changes, even though the fitted parameter (kMθσ) remains unchanged. The net resulting effect is approximately proportional. Thus, halving the assumed angular broadening reduces the inferred size expressed in kilometers by a factor of two, and so on. The previous measurement of θ by Gwinn et al. (1997) covered an incomplete arc in the (u, v) plane and depended on simpler approximations to the distribution of visibility for a scintillating source than the models used here; the measurement should be repeated. Refractive scintillation can change the angular broadening with time, but the variation is expected to be only about 85 μas for this line of sight (Romani et al. 1986; Narayan 1992). If scattering material is distributed along the line of sight rather than in a thin screen and some of it is located close to the source, it can increase the magnification M without a change in the measured angular broadening. This would require scattering close to the pulsar, within the Vela supernova remnant, where density is expected to be low.

We visualize the significance of our fits, by binning our residuals. This reduces effects of small-scale noise while allowing us to present the signatures of source size in the data. Philosophically, re-binning the residuals is similar to subtracting a model from the data. For example, we subtract a model for source position, baseline length, and Earth orientation from the phases in the correlator, so that the "sky" fringe rate can be reduced enough for time integration to recover the fringe with high signal-to-noise ratio. Fringing reveals the inaccuracies of the subtracted model. Here, analogously, we remove a model for the distribution of visibility from the histograms of data, reducing variations between neighboring bins from noise. We then average the residuals of nearby bins so that the effects of source structure are more easily apparent. This averaging reveals the degree to which our model fits match the observations.

Figure 7 shows the result of binning residuals, again for the data in IF1, gate 1, channels 4096–5120. After fitting zero-size and finite-size models, the residuals were binned and averaged in groups of 20, for resulting widths of 20 w = 0.0040 in Re[V]. Averaging (rather than summing) keeps the normalization the same. The dotted histogram shows the residuals from the fit with zero size for the source. We also display the differences between the model for finite source size and for zero size as solid lines. The curves are the same as those shown in Figures 3 and 5. A zero-size source would show a flat, dotted histogram. A model that perfectly corrected the deficiencies of the zero-size model would track the dotted histogram perfectly, allowing for the finite widths of the bins. Clearly, adding a size parameter explains most of the slowly varying residuals. Quantitatively, after binning the mean square residual is reduced by a factor of 13, by including a size parameter. This is much larger than the reduction observed before binning. As Figure 3 suggests, Poisson noise dominates the residuals without binning and is greatly reduced by binning.

Zoom In Zoom Out Reset image size

The model fit is imperfect. Some systematic variations remain, more so in some ranges of pulse phase than others. Figure 8 shows residuals for the zero-size model, and differences between finite-size and zero-size models, for representative spectral ranges in three gates. Fits to the same spectral range and same gate, but different IF, resemble one another strongly, as expected because they lead to nearly the same fitted size. Noise and amplitude parameters vary between IFs, and even more between gates and channel ranges, so that the binned residuals are not identical.

Zoom In Zoom Out Reset image size

The upper pair of panels of Figure 8 show IF2, gate 1, channels 4096–5120: these data are equivalent to those in Figure 7 in pulse phase and spectral range but are for the other IF. The data are thus completely independent. Introduction of a size parameter reduces the mean square residual of $\mathcal P$ by 92% and $\mathcal Q$ by 12%. The best-fitting size parameter is κ₁/κ₀ = (kMθσ)² = 0.043, or scaled size (kMθσ) = 0.21. The IFs agree closely in residuals and in fitted size. This range lies on the rising part of the pulse profile, before the peak.

The middle pair of panels show IF1, gate 2, channels 2048–3072. This range lies near the peak of the pulse. The best-fitting size parameter is κ₁/κ₀ = (kMθσ)² = 0.006, among the smallest results we obtain. Before binning, introduction of the size parameter reduced the Poisson-weighted residuals by less than 1%; this reduction is significant at the 2σ level according to the F-test. After binning, the reduction in mean square residual for ${\mathcal P}$ is 32%; however, zero-size and finite-size models both fit $\mathcal Q$ almost equally well. Some systematic effects appear, particularly the sharp downward spike in $\mathcal P$ near Re[V] = 0.

The lower pair of panels show IF2, gate 3, channels 4096–5120. This range lies on the trailing side of the pulse, near where it briefly flattens into a plateau. At shorter observing wavelength, this plateau becomes a second component that arises close to this pulse phase. Introduction of a size parameter decreases the residuals, by 10% in $\mathcal P$ and by 4% in $\mathcal Q$ before binning, and by 54% and 24% after binning. Again the influence of $\mathcal Q$ on the overall fit is small, because zero-size and finite-size models both fit $\mathcal Q$ well. The best-fitting size parameter is (kMθσ)² = 0.039 or (kMθσ) = 0.20. Although the fit to $\mathcal P$ is good, some systematic variation appears as a shifting of the histogram relative to the best-fitting model; our circular Gaussian model does not accommodate such a shift. We are investigating a variety of simple models for the systematic variations, as we discuss briefly in Section 4.3 below.

Further significant changes in the magnification factor from changes in pulsar distance seem unlikely. The distance from parallax is far more reliable than that from dispersion measure. In principle, a model that included significant scattering in the immediate neighborhood of the pulsar, as well as scattering in the surrounding Vela supernova remnant, could increase the effective magnification by bringing the lens closer to the pulsar. Temporal broadening at low frequencies suggests that scattering material is concentrated into a thin screen (M. D. Johnson & C. R. Gwinn 2012, in preparation). Previous measurements of the angular broadening θ should be revisited, as described in Section 4.1.3 above. Further studies of the angular broadening, using long baselines with different orientations, can help to improve on the earlier results. This may change the inferred linear size σ.

Instrumental effects are unlikely to reproduce the signature of source size on the distribution of visibility, as shown in Figure 7, because they tend to vary with pulse gate and frequency within the passband. For example, for both IFs, the noise parameter b₀ shows a similar pattern for each pulse gate, reflecting variations in the gain and noise in the passband, with an overall offset for each gate that reflects changes in quantization noise (see Gwinn et al. 2012, Section 5.1). Effects of quantization, in particular, can be expressed in the spectral domain as a gain and a change in noise and so are absorbed into the amplitude parameter and b₀. Effects of variation of noise parameters would be expected to appear as differences between overlapping gates.

The direction of linear polarization varies smoothly across the pulse. We observe only left-circular polarization, the fraction of which varies little. Scintillation in the interstellar medium is nearly polarization independent (Minter & Spangler 1996; Spangler 2001). Defects in separation of polarizations at the antennas could produce artifacts that vary with pulse phase, but these can be represented accurately as pulse-phase-dependent gains. They would not produce the distortion of the distribution of visibilities that we observe as the signature of source size.

Saturation effects can be important at large antennas such as Tidbinbilla. For our observations, these might affect particularly strong pulses, such as the giant pulses reported by Johnston et al. (2001). However, such pulses carry an insignificant fraction of the total flux density of the pulsar. Saturation effects can also reduce visibility for the very strongest scintillations, at least for baselines involving two large antennas (Gwinn et al. 2000, Section 4.5). These effects are expected to distort the distribution of visibility at the largest values but are not likely to reproduce the W-shaped signature of finite size we observe.

We perform relatively little calibration on the data. We do not change amplitude, particularly within a scan or with frequency. We do change phase, in fringing. Misfringing is possible for observations of a scintillating pulsar in the speckle limit, because the phase and amplitude vary on the frequency and timescales of scintillation (Desai et al. 1992), so that the model used for fringing can be challenging to fit. Misfringing would tend to increase phase variations, while keeping amplitude the same. Such effects are reduced for a short baseline, such as the Mopra–Tidbinbilla baseline, because the phase variation is less. If the phase errors were very large, it could potentially alter the distribution of visibilities in our projections, by altering the relative sizes of the real and imaginary part. We reduced errors from misfringing by fitting our fringe model over approximately 2 decorrelation times and 800 scintillation bandwidths. The wide bandwidth reduces effects of scintillation-based phase variations. The short time was chosen to match any possible ionosphere- or atmosphere-induced phase variations.

Effects of errors in the fringe model on size would be expected to reproduce among gates, since the fringe model from gate 2 was applied to all gates (with the exception of a single phase across the band, fitted for each gate); however, size varies with pulse phase, rather than by channel range within gates, and is near zero for some spectral ranges near the pulse peak. For channel ranges with sufficient signal-to-noise ratio, such as the later channels of gate 1 and most channels of gate 3, we found that fringing to the data within that gate, rather than using the model from gate 2, produced identical results. Similarly, applying the model from other strong gates to gate 2 yielded nearly the same results.

Effects of the finite sampling of the diffraction pattern are small and are expected to be constant across the pulse. Sampling over finite intervals of time and frequency blurs together correlated samples of the diffraction pattern and so is indistinguishable from effects of source size (Johnson & Gwinn 2012). For given scintillation timescale t_ISS and scintillation bandwidth Δν, effects of sampling can be calculated analytically (Gwinn et al. 2000, Section 2.2).

We estimate the scintillation timescale and bandwidth from our observations on the Mopra–Tidbinbilla baseline, using data near the pulse peak, IF2, gate 1, channels 2048–3072. We find correlation functions for the observed real part of interferometric visibility and fit expected functional forms to them (see Gwinn et al. 2000, Equations (5) and (7)). Figure 9 displays the data and results. We have normalized the correlation functions using the amplitude and offset of the best-fitting functional forms. The data point at zero lag is omitted, because it is elevated by effects of noise and variations of the amplitude with time. Both correlation functions extend to large lag, off the right of the displayed panels; these are included in the fits but not displayed here.

Zoom In Zoom Out Reset image size

A fit of the expected Gaussian function to the temporal correlation function yields a scintillation timescale of t_ISS = 8.99 ± 0.17 s. Our sampling time of 2 s is then expected to decrease the mean square flux density by a factor of f_t = 0.994 (Gwinn et al. 2000, Equation (8)), while leaving the average flux density unchanged. The variation of f_t from the standard error of the fit is less than 0.001.

A fit of the expected Lorentzian function to the frequency correlation function yields a scintillation bandwidth of Δν = 15.2 ± 0.5 kHz. Our channel bandwidth of 1.95 kHz is then expected to decrease the mean square flux density by a factor of f_f = 0.998 (Gwinn et al. 2000, Equation (6)), leaving the mean flux density unchanged. The observed correlation function shows systematic departures from the best-fitting Lorentzian function, with a sharper peak and elevated correlation at intermediate lag, with a peak near a lag of 70 kHz. A Kolmogorov rather than a square-law structure function and a departure of the scattering disk from isotropy can produce such effects (Gwinn 2001; Johnson & Gwinn 2012). Temporal variation of the emission on short timescales will also introduce correlations in frequency (Gwinn & Johnson 2011). Models including such effects will fit better. However, as the figure shows, the correlation falls to half its peak value at a lag of approximately 14 kHz, providing a clear characteristic scale. The variation of f_f from the standard error of t_ISS is less than 0.001. Even with Δν = 12 kHz we still find f_f = 0.998.

Together, the effects of finite sampling in time and in frequency will increase the fitted size. Under the assumption that the source is pointlike, so that (kMθσ) = 0 and we misinterpret the effects of averaging in time and frequency as source size, our sampling parameters lead to a fictitious size parameter of (kMθσ) = 0.06, where we have related effects of finite sampling to size using the expressions for modulation index from averaging (Gwinn et al. 2000, Equation (9)) and for source size (Gwinn 2001, Equations (34) and (35)). This value is indeed nearly the minimum that we observe, near the pulse peak, as Figure 6 shows.

The effects of finite sampling in time and frequency are constant with pulse phase, for a pointlike source. This is clearly inconsistent with our results discussed in Section 4 above. Variation of amplitude with time can affect estimates of the scales of scintillation: variability on long timescales can affect the correlation with time, and on short timescales it can affect the correlation in frequency (Gwinn & Johnson 2011; Gwinn et al. 2011). Erroneous estimates will not produce variations of source size with pulse phase. Direct effects of source variability on estimates of source size are discussed elsewhere in this paper (Sections 3.4 and 3.5.2).