ULTRA-HIGH-RESOLUTION INTENSITY STATISTICS OF A SCINTILLATING SOURCE

Scintillation of compact astrophysical radio sources has proven to be a remarkable probe of both the turbulent interstellar plasma and the extreme conditions of the emission regions of these sources. The present work is primarily motivated by observations of pulsars at meter and decimeter wavelengths, for which the AU-scale scattering resembles a stochastic "lens" with a potential resolving power of less than a nanoarcsecond. However, our results apply to any compact object that exhibits diffractive scintillation.

The emission and scattering physics involve the interaction and mixing of many random variables. At the heart is the amplitude-modulated noise emitted by the source (Rickett 1975). Indeed, even the amplitude modulation is stochastic, with individual pulses showing complex broadband profile structures, from bursty microstructure to strong and often non-Gaussian pulse-to-pulse variations (Rickett et al. 1975; Cairns et al. 2001; Kramer et al. 2002). This emission then propagates through the dilute, turbulent plasma of the interstellar medium (ISM), where it undergoes dispersion and Faraday rotation, in addition to being strongly scattered. Finally, an observer samples the electric field in the presence of additive background noise, which often dominates the signal. Furthermore, the statistics of this random process evolve as the Earth, pulsar, and ISM move relative to each other.

Since the size of the scattering disk (Σ) scales approximately as the square of the observing wavelength (λ), the angular resolution afforded by the scintillation is θ_diff ∼ λ/Σ∝λ⁻¹. Thus, longer wavelengths permit superior resolving power. However, at longer wavelength each scintillation element contains far fewer samples because of the relatively smaller timescale (Δt_d∝λ⁻¹) and bandwidth (Δν_d∝λ⁻⁴) of the scintillation pattern. It therefore becomes subtle to decouple the scintillation from the source and background noise.

To address these difficulties, we present analytical expressions for the probability density function (PDF) of intensity after arbitrary temporal averaging of the spectra. We describe the case when single-pulse spectra are formed and then are averaged in groups of N adjacent pulses; we exclusively use N to denote such temporal averaging. We account for the decorrelation of the scintillation pattern during this averaging and derive the resulting modification to the intensity PDF. We also account for spatially extended source emission and show that the leading order PDF correction is identical, up to scale, to that of the decorrelation of the scintillation pattern within the averaging time. In addition, we demonstrate that spatially extended emission affects the measured PDF of intensity, even in the limit of N = 1 averaged spectra. This feature presents the remarkable opportunity of measuring the emitting region sizes for individual pulses, or of subclasses of pulses.

The sources and variety of noise are diverse, stemming from the background noise and the amplitude-modulated noise emitted by the source. Measured quantities such as intensity then correspond to random variables, and their stochastic variation introduces additional noise. Applied to the source, this noise is the familiar self-noise (Dicke 1946; Kulkarni 1989; Zmuidzinas 2003; Gwinn & Johnson 2011), but the background and scintillation parameters have similar limitations. Such noise is generally unbiased and scales with intensity. Our derivations fully account for all these types of noise. We also account for artifacts arising from instrumental limitations, such as quantization of the analog signal and quadrature downconversion. These types of effects introduce additional noise as well as bias and distortion.

To supplement our calculated PDFs, we also analyze expressions based on the moments of the distributions. In particular, we extend the modulation index, a traditional estimator of source emission geometry, to include the effects of self-noise and pulsar amplitude variations, and we analyze a technique that estimates self-noise. Our expressions refine these simple tools for intrinsic emission inference.

An observation well suited for comparison with our models has an observing bandwidth B ≫ Δν_d, an observation duration t_obs ≫ Δt_d, and a spectral accumulation time much longer than the pulse-broadening timescale. For each pulse period, the on-pulse and off-pulse regions yield estimates of the signal and noise intensities, respectively. These estimated parameters fully define a model of the intensity PDF (see Section 3), which is readily compared with a histogram of the measured data.

The residual structure in this comparison contains information about the spatial extent of the emission region and the evolution of the scintillation pattern, as well as artifacts from instrumental limitations and errors in the estimated on-pulse and off-pulse intensities. The on-pulse and off-pulse averages once again parameterize a model for this residual structure (see Section 4), which is quantified by a single fitted parameter for amplitude. This fitted parameter yields the transverse size of the emission region, if the scattering geometry is specified.

Most past work has focused on the regime of many samples of the spectrum averaged together: N → ∞. For such data, the effects of noise can be approximated as Gaussian. Gwinn et al. (1998) gave the expected distribution of intensity for a strongly scintillating source, including the effects of a spatially extended emission region. Gwinn et al. (2000) extended these forms to account for Gaussian self-noise at small intensity, and estimated the contributions of averaging in time and frequency by deriving their effects on the modulation index. In an unpublished manuscript, Cordes (2000) presented expressions for the PDF of scintillation gain in terms of a Karhunen–Loève expansion with the coefficients defined by a Fredholm equation, and suggested approximate forms determined by χ² distributions with the appropriate number of degrees of freedom. Gwinn (2001) gave the distribution of interferometric visibility for a pointlike or extended source in strong scintillation. Gwinn et al. (2011) presented expressions for self-noise, suitable for large N.

The present work focuses on the complementary regime of single spectral realizations, and averages of N such spectra. We connect with these previous works through the asymptotic forms of our distributions. Similarly, van Straten (2009) analyzed statistics of Stokes parameters when N = 1 and showed, for instance, that the degree of polarization is undefined in this case. His results describe the high signal-to-noise regime with no effects of propagation or source emission structure, as appropriate for the high-frequency observations of giant pulses from the Crab pulsar that motivated his work (Hankins et al. 2003).

In Section 2, we briefly review the theoretical descriptions of pulsar emission and interstellar scintillation, and establish our assumptions about each. Then, in Section 3, we present an expression for the exact intensity PDF of a scintillating point source, as well as useful approximations. We next apply these approximations in Section 4 to estimate modifications to the PDF arising from the decorrelation of the scintillation pattern within the averaging time, effects of finite source emission size, and observational constraints such as finite scintillation averages; we then use these expressions to derive the limiting resolution for emission size inference. In Section 5, we derive the modulation index and analyze a technique that identifies self-noise. Finally, in Section 6, we summarize our results and outline the prospects for pulsar observations.

Pulsar radio emission is thought to arise from coherent curvature radiation from a relativistic electron–positron plasma, flowing outward along open field lines from the stellar surface to the light cylinder. Charge acceleration occurs in a "gap" of depleted charge density where the force-free state cannot be maintained. Such sites initiate a pair-production cascade, which forms a secondary plasma (Sturrock 1971; Ruderman & Sutherland 1975). This secondary plasma generates the radio emission, which is then beamed, ducted, and refracted as it propagates through the upper magnetosphere (Barnard & Arons 1986; Arons & Barnard 1986; Lyutikov & Parikh 2000). The gap may be located close to the polar cap region (the "polar gap"; Arons & Scharlemann 1979), or along the boundary of the open and closed field lines (the "slot gap"; Arons 1983). An "outer gap" (Cheng et al. 1986; Romani 1996), situated close to the light cylinder radius, is the favored site for much high-energy emission and may also contribute to the radio emission. Differing locations of the gap may be responsible for some of the emission variety between pulsars.

The superposition of many independent radiators produces white Gaussian noise. Thus, although the radio emission process is fundamentally coherent, by groups of many particles, the bulk emission may still be effectively spatially and temporally incoherent. Intrinsic physical processes (and the pulsar's rotation) then modulate the envelope of this random field, motivating an "amplitude-modulated noise" description of the emission (Rickett 1975). Our results assume this type of emission, confined to a region that is a small fraction of the diffractive scale r_d ≡ λD/(2πΣ) ∼ 10⁶ m, where D is the characteristic Earth–scatterer distance.

Density inhomogeneities in the dilute, turbulent plasma of the ISM scatter the radio emission, leading to multipath propagation. For most pulsars observed at decimeter wavelengths, the phases between pairs of paths differ by many turns and are therefore uncorrelated; this is the regime of "strong" scattering. Equivalently, the diffractive scale is much smaller then the Fresnel scale, $r_{\mathrm{F}} \equiv \sqrt{\lambda D} \sim 10^9 \rm {\ m}$.

In general, the timescale for the evolution of the diffraction pattern at the observer ranges from seconds to hours, and so is effectively static over the scattering timescale. In this case, the effect of scattering is to convolve the pulsar signal with a "propagation kernel." The characteristic width of the average propagation kernel is the pulse-broadening timescale t₀. Dispersion contributes an additional frequency-dependent delay that is readily reversed (Hankins 1971). The motivation and limitations of approximating scattering as a convolution are given in much greater detail by Gwinn & Johnson (2011).

Although this picture of strong scattering of a point source is quite general, a description of the effects of a spatially extended emission region requires specification of the scattering geometry. We present explicit results for the case of thin-screen scattering, as is observationally motivated (Williamson 1972; Komesaroff et al. 1972; Hill et al. 2005). We give additional reductions by assuming a square-law phase structure function; the generalization to Kolmogorov or otherwise is straightforward and changes only the scale of the PDF modifications that we derive.

In practice, spectra must be formed over a finite accumulation time t_acc; the convolution representation of scattering requires t_acc ≫ t₀ (see Gwinn & Johnson 2011 for additional details). The unique nature of pulsars and their generally short duty cycles allow a particularly elegant limiting option: the formation of single-pulse spectra that contain all pulsed power. Our descriptions are exact in this limit, and appropriate whenever t_acc ≫ t₀.

We assume that such spectra are formed and then averaged in groups of N adjacent pulses. We also assume that the averaging timescale is much shorter than the scintillation timescale—the "snapshot image" of Narayan & Goodman (1989)—which subsumes the much weaker condition that the scattering material is approximately static during the accumulation of a single spectrum.

Frequency averaging is analogous but inherits additional information from non-stationary signals. The intrinsic pulse shape incurs spectral mixing via its associated frequency-domain convolution; the average is then of correlated intensities (see Appendix A.2.2). Thus, a description of frequency averaging requires the specification of individual pulse profiles, whereas a description of temporal averaging requires only the phase-averaged source intensity for each pulse.

Lastly, we assume that the data explore a large representation of the full ensemble of diffractive scintillation: (B/Δν_d) × (t_obs/Δt_d) ≫ 1. Comparison with observations is most powerful if the stronger condition B ≫ Δν_d is satisfied so that the intrinsic amplitudes of individual pulses can be estimated.

Despite the number and variety of random variables involved in the emission and scattering processes, much statistical homogeneity exists; our notation highlights this foundational uniformity. We use z_x to denote a circular complex Gaussian random variable with unit variance, indexed by x. Some common examples of such variates are a single realization of the scalar electric field or the complex gain from propagation. The PDF of z_x is exponential in its squared norm, $P(z_x) = \pi ^{-1} e^{-|z_x|^2}\!$, with respect to the standard complex metric dRe[z_x]dIm[z_x]. This norm is implicit for all PDFs presented with respect to complex random variables.

We also encounter many exponential random variables, usually in connection with a scintillation "gain," which are thus denoted G or G_i, with P(G) = e^−G for G > 0. Note that we use P() to generically denote a PDF with respect to the given variables and parameters, as in the previous examples.

The pulsar emits white Gaussian noise _i, modulated by a time-varying envelope $\sqrt{I_{\mathrm{s}}} f_i$. Here, I_s is a constant characteristic scale of the source intensity, _i is normalized to have unit variance, and the modulation f_i is power preserving. To account for pulse-to-pulse variations in the phase-averaged flux density, we also include a dimensionless pulse amplitude factor A_j, where j indexes the spectrum (i.e., the pulse number). We require no assumptions about the timescale of variability of the envelope f_i or the nature of the pulse-to-pulse variations A_j. Our treatment therefore accommodates arbitrary variability of the pulsar, such as the possibility of correlated pulse-to-pulse variations, log-normal amplitude statistics, or nanosecond-scale bursts, which merely modulates the Gaussian self-noise.

The propagation then acts to convolve this intrinsic emission with a scattering kernel, as discussed in Section 2.2. The kernel similarly consists of a power-preserving envelope g_i that modulates white Gaussian noise η_i of unit variance. The propagated signal is then superimposed with the background noise $\sqrt{I_{\mathrm{n}}}\beta _i$, where I_n is the amplitude of the noise and β_i is white Gaussian noise of unit variance. We assume that the background noise level is effectively constant over the averaging time, although this restriction is easily relaxed.

The observed scalar electric-field time series x_i and its Fourier-conjugate spectrum are thus given by

$$\begin{eqnarray} x_i &= \sqrt{A_j I_{\mathrm{s}}} \left[ \left(f \epsilon \right) \ast \left(g \eta \right) \right]_i + \sqrt{I_{\mathrm{n}}}\beta _i \nonumber\\ \Rightarrow \tilde{x}_i &= \sqrt{A_j I_{\mathrm{s}}} (\tilde{f} \ast \tilde{\epsilon })_i (\tilde{g} \ast \tilde{\eta } )_i + \sqrt{I_{\mathrm{n}}} \tilde{\beta }_i. \end{eqnarray} \tag{ 1 }$$

Here, we denote Fourier-conjugate variables with a tilde. In short, this equation describes the amplitude-modulated noise of the pulsar ($\sqrt{I_{\mathrm{s}}} f_i\epsilon _i$) convolved with the stochastic propagation kernel (g_iη_i) and added to background noise ($\sqrt{I_{\mathrm{n}}} \beta _i$). Conventionally, these terms are labeled by their noise-free limits: the pulse profile |f_i|², pulse-broadening function |g_i|², gated signal I_s, and background I_n.

Of course, $\tilde{\epsilon }_i$, $\tilde{\eta }_i$, and $\tilde{\beta }_i$ are mutually independent (circular complex Gaussian) white noise. We see that a single spectral sample $\tilde{x}_i$ is of the form $\sqrt{A_j I_{\mathrm{s}}} z_{\mathrm{f}} z_{\mathrm{g}} + \sqrt{I_{\mathrm{n}}} z_{\mathrm{b}}$, where z_f, z_g, and z_b are each circular complex Gaussian random variables with unit variance. If the scintillation is held fixed (e.g., z_g = const.), then the intensity $|\tilde{x}_i|^2$ is drawn from an exponential distribution with scale $\bar{I}_j \equiv A_j I_{\mathrm{s}} |z_{\mathrm{g}}|^2 + I_{\mathrm{n}}$. The PDF of the average of N such intensities is then given by Equation (A1):

$$\begin{eqnarray} P(I;N | z_{\mathrm{g}}) &= N \sum _{j=1}^N \left(\frac{\bar{I}_j^{N-2} }{\prod \limits _{{{\scriptstyle {\begin{array}{ll}\ssty \ell = 1\nonumber\\ \ssty \ell \ne j\end{array}}}}}^N (\bar{I}_j - \bar{I}_\ell )} \right) e^{-N I/\bar{I}_j}\\ \nonumber &= \frac{N}{ (I_{\mathrm{s}} |z_{\mathrm{g}}|^2 )^{N-1} } \\ & \quad \times \sum _{j=1}^N \frac{(A_j I_{\mathrm{s}} |z_{\mathrm{g}}|^2 + I_{\mathrm{n}} )^{N-2} }{\prod \limits _{{{\scriptstyle {\begin{array}{l}\ssty \ell = 1\\ \ssty \ell \ne j\end{array}}}}}^N (A_j - A_\ell )} e^{-\frac{N I}{(A_j I_{\mathrm{s}} |z_{\mathrm{g}}|^2 + I_{\mathrm{n}} )}}.\quad\ \end{eqnarray} \tag{ 2 }$$

Observe that when we designate N as a distribution parameter, we implicitly include the background I_n and source amplitudes {A_jI_s}.

All that remains is to incorporate the distribution of the scintillation random variable z_g. We assume that z_g explores a representation of the full diffractive ensemble. We introduce a "scintillation gain" parameter G ≡ |z_g|², so P(G) = e^−G. Integrating the scintillation ensemble through G then provides the PDF of intensity:

$$\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}{l}\ssty \ell = 1\\ \ssty \ell \ne j\end{array}}}}}^N (A_j - A_\ell )} e^{-\frac{N I}{\left(A_j I_{\mathrm{s}} G + I_{\mathrm{n}} \right)}}. \end{eqnarray} \tag{ 3 }$$

This equation gives the expected PDF of intensity for the average of N pulses with amplitudes {I_sA_i} and fixed scintillation, in the presence of white background noise with mean intensity I_n, if the spectra explore a representative ensemble of diffractive scintillation. Examples of this distribution for various N are shown in Figure 1.

Zoom In Zoom Out Reset image size

Numerical evaluation of the PDF can be challenging because of the singularity at G = 0, although the integrand is not divergent. To address this difficulty, we employ variable-precision arithmetic libraries; the additional precision is essential for N ≳ 10.

Although we have derived the exact expression for the intensity PDF (Equation (3)), we now derive several useful approximations. These approximations resolve the singularity at G = 0, uncover the N → ∞ behavior, and present a simplified framework for deriving and interpreting perturbations of the PDF arising from effects such as a finite size of the emission region, decorrelation of the scintillation pattern within the averaging time, and errors in model parameters.

3.2.1. The i.i.d. Approximation

We derive our first approximation by replacing each set of N pulse amplitude factors A_j by its average A ≡ 〈A_j〉. This replacement simplifies the statistics because the intensities that are averaged under such assumptions are i.i.d.; thus, we refer to this approximation scheme as the i.i.d. approximation.

We emphasize that the i.i.d. approximation does not ignore amplitude variability of the pulsar, but does reduce the variability to a moving average. The difference is critical because, for example, the i.i.d. approximation is exact when N = 1.

After substituting the local amplitude average A, the PDF of intensity during a fixed scintillation gain G (i.e., Equation (2)) is given by an Erlang distribution (Equation (A2)):

$$\begin{eqnarray} &&P(I;N|G) = \frac{N^N}{(N-1)!} \frac{I^{N-1}}{\left(A I_{\mathrm{s}} G + I_{\mathrm{n}}\right)^N} e^{-\frac{N I}{A I_{\mathrm{s}} G + I_{\mathrm{n}}}}. \end{eqnarray} \tag{ 4 }$$

As for Equation (3), the scintillation ensemble can then be included via its exponential density: P(I; N) = ∫^∞₀dG e^−GP(I; N|G).

3.2.2. The Gamma Approximation

We now extend the i.i.d. approximation to partially account for the effects of variation in the pulse amplitudes A_j within sets of N pulses. To achieve this extension, we replace N in Equation (4) by an effective number of degrees of freedom N_eff, chosen such that the mean and variance match those of the exact PDF. Because N_eff is not necessarily an integer, the factorial in the Erlang distribution must be replaced by a gamma function; this modified form is the gamma distribution. The mean and variance of this gamma distribution are

$$\begin{eqnarray} \mu _{\rm gamma} &= A I_{\mathrm{s}} G + I_{\mathrm{n}},\nonumber\\ \sigma _{\rm gamma}^2 &= \frac{(A I_{\mathrm{s}} G + I_{\mathrm{n}})^2}{N_{\mathrm{eff}}}. \end{eqnarray} \tag{ 5 }$$

For the exact distribution, given by Equation (2), the mean is simply the average of the marginal means, and the variance is the average of the marginal variances, divided by N:

$$\begin{eqnarray} \mu _{\rm exact} &= \frac{1}{N}\sum _{j=1}^N (A_j I_{\mathrm{s}} G + I_{\mathrm{n}}), \nonumber\\ \sigma _{\rm exact}^2 &= \frac{1}{N^2}\sum _{j=1}^N (A_j I_{\mathrm{s}} G + I_{\mathrm{n}})^2. \end{eqnarray} \tag{ 6 }$$

Matching Equations (5) and 6 gives the correspondence appropriate for the gamma distribution:

$$\begin{eqnarray} A &= \frac{1}{N} \sum _j A_j, \nonumber\\ N_{\mathrm{eff}} &= \frac{\big(\sum _{j=1}^N 1 + A_j I_{\mathrm{s}} G/I_{\mathrm{n}} \big)^2}{ \sum _{j=1}^N (1 + A_j I_{\mathrm{s}} G/I_{\mathrm{n}} )^2} \le N. \end{eqnarray} \tag{ 7 }$$

The scintillation ensemble is again included via integration over G with its exponential weight.

We refer to this extension of the i.i.d. approximation as the gamma approximation. The gamma approximation partially accounts for the pulse-to-pulse variability, and substantially improves the i.i.d. approximation for N ≫ 1.

3.2.3. The Gaussian Approximation

These three approximations are particularly useful because they reduce the N + 1 degrees of freedom in a single averaged spectra to three: an overall scale partitioned into source and background contributions, and a measure of the averaging. In fact, the i.i.d. approximation has only two degrees of freedom because the degree of averaging is N; this reduction comes at the expense of ignoring pulse amplitude variance.

Figure 2 demonstrates the relative quality of these three approximations in the presence of substantial amplitude variations. The Gaussian approximation is inferior to both the i.i.d. and gamma approximations, particularly for small N; indeed, the Gaussian approximation includes non-zero probability for negative intensities. The i.i.d. and gamma approximations, on the other hand, give excellent representations for arbitrary N and are much easier to implement numerically than the exact representation of Equation (3), particularly for large N.

Zoom In Zoom Out Reset image size

Our principal use of these approximations is to derive analytical estimates for PDF modifications arising from effects such as finite emission size and decorrelation of the scintillation pattern within the averaging time. We favor the relatively simple i.i.d. approximation for this purpose. While the small errors in the approximations might be significant in some measured PDFs, we expect the similarly small errors on the expected ≲ 1% PDF modifications to be negligible. Furthermore, the most intriguing degree of averaging is N = 1, which allows single-pulse studies and complete decoupling of emission size from decorrelation; both the i.i.d. and gamma approximations are, of course, exact for this case.

We therefore recommend the exact representation for PDF comparisons (or the gamma approximation if its accuracy suffices), and now derive PDF modifications using the i.i.d. approximation.

The scintillation pattern slowly evolves in response to the relative motions of the Earth, pulsar, and scattering medium, thereby affecting statistics of averaged intensities. We now derive the effects of this evolution on the observed distribution of intensity, in terms of the temporal autocorrelation function of the observed scintillation pattern.

Although the scintillation pattern evolves somewhat similarly with frequency, the effects of frequency averaging fundamentally differ from those of temporal averaging. Namely, for temporal averaging, the intensities within each average are independent, with correlated scales (i.e., the respective scintillation gains); whereas for frequency averaging, both the intensities within each average and their scales can be correlated because of the intrinsic amplitude modulation. As noted in Section 2.3, spectral statistics after frequency averaging depend on individual pulse profiles, whereas spectral statistics after temporal averaging only depend on the phase-averaged flux densities.

Procedurally, we must index the scintillation factors G by pulse number. The resulting set {G_i} is then a collection of correlated exponential random variables, with covariance determined by the autocorrelation function of the measured dynamic spectrum:

$$\begin{eqnarray} \Gamma _{ij} &\equiv \langle G_i G_j \rangle -1\nonumber\\ &= \frac{\langle \left[ I(t,\nu) - I_{\mathrm{n}} \right][ I(t+\Delta t_{ij},\nu) - I_{\mathrm{n}} ] \rangle }{\left\langle I(t,\nu) - I_{\mathrm{n}} \right\rangle ^2} - 1. \end{eqnarray} \tag{ 8 }$$

Details of this type of multivariate exponential distribution are given in Appendix A.2.

In the "snapshot image" regime, the averaging timescale is much shorter than the decorrelation timescale of the scintillation pattern, so Γ_ij ≈ 1. The leading order decorrelation effects in (1 − Γ_ij) then provide an excellent approximation. Using Equation (A5), we obtain

$$\begin{eqnarray} &&P(I;N,\lbrace \Gamma _{ij}\rbrace) - P(I;N)\nonumber\\ &&\quad \approx - \sum _{i<j} (1 - \Gamma _{ij} ) \int _0^\infty dG\ G e^{-G} \left.\frac{\partial ^2 P(I;N|\lbrace G_i\rbrace)}{\partial G_i \partial G_j}\right\rfloor _{G}\!.\qquad \end{eqnarray} \tag{ 9 }$$

Note that the partial derivatives are easily evaluated using the characteristic function of P(I; N|{G_i}): $\varphi (k) = \prod _\ell (1-i \bar{I}_\ell k/N)^{-1}$, where $\bar{I}_\ell \equiv A_j I_{\mathrm{s}} G_j + I_{\mathrm{n}}$. For example, applying the i.i.d. approximation (A_j ≡ A) yields the following relationship:

$$\begin{eqnarray} &&\left.\frac{\partial ^2 P(I;N|\lbrace G_i\rbrace)}{\partial G_i \partial G_j}\right\rfloor _{G} = \frac{1}{N(N+1)} \frac{\partial ^2 P(I;N|G)}{\partial G^2}. \end{eqnarray} \tag{ 10 }$$

We can also reduce the sum over Γ_ij in Equation (9) by approximating the evolution of the scintillation pattern as relative motion of the pulsar and observer with respect to a "frozen" random screen. The decorrelation in time is then simply expressed in terms of the phase structure function of the scattering medium (Goodman 1985; Cornwell et al. 1989). For example, a square-law phase structure function has Γ_ij = exp [ − (Δt_ij/Δt_d)²], even in the presence of axial anisotropy (Rickett 1977; M. D. Johnson et al. 2012b, in preparation). If the averaging timescale is much less than the decorrelation timescale, then ∑_{i < j}(1 − Γ_ij) ≈ N²(N² − 1)/(12Δτ²_d), where Δτ_d is the decorrelation timescale in pulse periods. The quartic dependence on N leads to a rapid onset of decorrelation artifacts in the observed PDF after only modest averaging.

Combining this approximation for the decorrelation behavior with the general form for the PDF modification (Equation (9)), and applying the reduction afforded by the i.i.d. approximation (Equation (10)) then provides the following estimate:

$$\begin{eqnarray} &&P(I;N,\lbrace \Gamma _{ij}\rbrace) - P(I;N) \nonumber\\ && \quad \approx - \frac{N(N-1)}{12 \Delta \tau _{\mathrm{d}}^2} \int _0^\infty dG\ G e^{-G} \frac{\partial ^2 P(I;N|G)}{\partial G^2}. \end{eqnarray} \tag{ 11 }$$

A spatially extended emission region smooths the fluctuations from scintillation by superimposing many slightly offset copies of the diffraction pattern at the observer. In fact, extended source emission is merely one example in a broad class of physical effects that alter the distribution of scintillation gain. Other possibilities arise from modified scattering assumptions, such as weak scintillation or Lévy-type scattering statistics (Boldyrev & Gwinn 2003). To account for these cases, we introduce a generalized scintillation gain random variable $\mathcal {G}$, which is no longer constrained to follow exponential statistics. Specification of $P(\mathcal {G})$, combined with $P(I;N|\mathcal {G})$ given by Equation (3), then provides the intensity PDF P(I; N).

Gwinn et al. (1998) examined $P(\mathcal {G})$ for thin-screen scattering of a Gaussian distribution of source intensity. They demonstrated that, if the source emission region is small relative to the magnified diffractive scale, $\mathcal {G}$ is well approximated by the convolution of three independent exponential random variables. The scales of these random variables reflect the size of the emission region. Gwinn (2001) derived an analogous result for interferometric visibility.

By extending the description of Section 3.1, we now demonstrate that an equivalent form arises for any spatially incoherent, compact emission region in strong scattering. We parameterize the transverse coordinates of the source by $\textbf {s}$. Each emitting location has independent source noise, and the emission envelope at each location may also vary. However, the propagation kernel is correlated over the emission region. A single spectral sample thus takes the form

$$\begin{eqnarray} &&\tilde{x}_i = \left\lbrace \int d^2 \textbf {s} \ \sqrt{A(\textbf {s}) I_{\mathrm{s}}} z_{\mathrm{f}}(\textbf {s}) z_{\mathrm{g}}(\textbf {s}) \right\rbrace + \sqrt{I_{\mathrm{n}}} z_{\mathrm{b}}. \end{eqnarray} \tag{ 12 }$$

We specify our source coordinates, $\mathbf {s}$, to be centered on the spatial mean of source intensity: $\int d^2\textbf {s}\ \textbf {s} A(\mathbf {s}) = \mathbf {0}$.

If the emission is confined to a small fraction of the magnified diffractive scale, then the scintillation random variable, $z_{\mathrm{g}}(\textbf {s})$, will be strongly correlated over the region of source intensity. We therefore expand to linear order: $z_{\mathrm{g}}(\textbf {s}) \approx z_{\mathrm{g}}(\textbf {0}) \,{+}\, (\textbf {s} \,{\cdot}\, \nabla) z_{\mathrm{g}} \rfloor _{s=0}$. The source term in Equation (12) is then a convolution of three complex Gaussian random variables. Because of our choice of origin for $\textbf {s}$, these three random variables are mutually uncorrelated during a fixed scintillation pattern. The marginal variances therefore add to give $\mathcal {G}$ (up to overall normalization). In addition, the scintillation random variable, $z_{\rm g}(\textbf {s})$, is uncorrelated with its spatial derivatives, so the scales of the three variances are mutually independent. Thus, $\mathcal {G}$ is the convolution of three independent exponential random variables:

$$\begin{eqnarray} &&\mathcal {G} \equiv \frac{1}{1 + \gamma _{\mathrm{s},1} + \gamma _{\mathrm{s},2}} (G + \gamma _{\mathrm{s},1} G_1 + \gamma _{\mathrm{s},2} G_2 ). \end{eqnarray} \tag{ 13 }$$

The dimensionless subsidiary scales γ_{s, i} ≪ 1 contain information about the transverse extent of the source emission. The scaling prefactor is merely chosen so that $\langle \mathcal {G} \rangle = 1$ and could equally well be absorbed into the definition of source intensity.

Using the results of Gwinn et al. (1998), we now explicitly relate the parameters γ_{s, i} to the emission geometry by assuming a Gaussian distribution of source intensity and thin-screen scattering. In this case, the γ_{s, i} give the squared size of the source in orthogonal directions $\hat{x}_i$, in units of the magnified diffractive scale:

$$\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{ 14 }$$

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 scattering. If the source emission arises from a circular Gaussian intensity profile, then the FWHM of the emission region is $2\sqrt{\ln 4} \sigma \approx 2.35\sigma$, where σ ≡ σ_i.

Efforts to infer emission size often attempt to isolate the distribution $P(\mathcal {G})$ through sufficient spectral averaging; the modulation index m then quantifies the effects of emission size on this distribution (Salpeter 1967; Cohen et al. 1967; Gwinn et al. 1998):

$$\begin{eqnarray} m^2 &= \frac{\langle \mathcal {G}^2 \rangle }{\langle \mathcal {G} \rangle ^2} - 1\nonumber\\ &= \frac{1 + \gamma _{\mathrm{s},1}^2 + \gamma _{\mathrm{s},2}^2}{(1 + \gamma _{\mathrm{s},1} + \gamma _{\mathrm{s},2} )^2} \approx 1 - 2(\gamma _{\mathrm{s},1} + \gamma _{\mathrm{s},2} ). \end{eqnarray} \tag{ 15 }$$

In Section 5.1, we discuss the modulation index in depth and extend Equation (15) to account for the effects of self-noise, temporal decorrelation within the averaging time, and pulsar amplitude variability.

To determine the PDF modification resulting from a finite emission size, we first expand the PDF for a fixed scintillation pattern to linear order in γ_{s, i}:

$$\begin{eqnarray} &&P(I;N,\gamma _{\mathrm{s}}|\mathcal {G}) - P(I;N,\gamma _{\mathrm{s}}\,{=}\,0|G)\nonumber\\ &&\quad \approx [ G_1 \gamma _{\mathrm{s},1} + G_2 \gamma _{\mathrm{s},2} - G(\gamma _{\mathrm{s},1} + \gamma _{\mathrm{s},2} ) ] \frac{\partial P(I;N,\gamma _{\mathrm{s}}\,{=}\,0|G)}{\partial G}.\nonumber\\ \end{eqnarray} \tag{ 16 }$$

Of course, the scintillation ensemble must now be included. Integrating each G_i with its respective weight $e^{-G_i}$ gives

$$\begin{eqnarray} &&P(I;N,\gamma _{\mathrm{s}}|G) - P(I;N,\gamma _{\mathrm{s}}\,{=}\,0|G) \nonumber \\ && \quad \approx 2\gamma _{\mathrm{s}} \left(1 - G\right) \frac{\partial P(I;N,\gamma _{\mathrm{s}}\,{=}\,0|G)}{\partial G}. \end{eqnarray} \tag{ 17 }$$

Here, we have eliminated the degeneracy that appears at linear order between γ_{s, 1} and γ_{s, 2} by setting them both equal to a single characteristic size γ_s. If the effects of source emission size are substantial, then additional emission information, such as elongation or core/cone geometry, can be inferred from the higher order terms, although higher terms in the original gain expansion must also be retained. To facilitate a comparison with the decorrelation effects derived in Section 4.1, we integrate by parts after combining with the final scintillation weight:

$$\begin{eqnarray} &P(I;N,\gamma _{\mathrm{s}}) - P(I;N,\gamma _{\mathrm{s}} \,{=}\, 0) \nonumber\\ & \quad \approx - 2\gamma _{\mathrm{s}} \int _0^\infty dG \ G e^{-G} \frac{\partial ^2 P(I;N,\gamma _{\mathrm{s}}\,{=}\,0|G)}{\partial G^2}. \end{eqnarray} \tag{ 18 }$$

See Figure 3 for a comparison of these effects of emission size for various N.

Zoom In Zoom Out Reset image size

The approximate form of Equation (18) must be applied with some care; it fails in the limit N → ∞, for example, because of the resulting divergent derivatives of the delta function that characterizes the measured intensity. However, the large-N limit is still easily accessible numerically by applying the full distribution of scintillation gain, as given by Equation (30) in Gwinn et al. (1998). Substituting γ_s = γ_{s, 1} = γ_{s, 2} in their expression and scaling to normalize the mean gives

$$\begin{eqnarray} &&P(\mathcal {G}) = \frac{1+2\gamma _{\mathrm{s}}}{\left(1-\gamma _{\mathrm{s}}\right)^2}\nonumber\\ &&\quad \times\bigg [ e^{-(1 + 2\gamma _{\mathrm{s}})\mathcal {G}}- \left(\frac{\mathcal {G}}{\gamma _{\mathrm{s}}} + (1+\mathcal {G}) - 2\mathcal {G}\gamma _{\mathrm{s}} \right) e^{-\left(2 + \frac{1}{\gamma _{\mathrm{s}}} \right)\mathcal {G}} \bigg ].\quad\quad \end{eqnarray} \tag{ 19 }$$

When combined with the conditional PDF $P(I;N|\mathcal {G})$ (Equation (2)), this equation produces the scintillation-averaged PDF of intensity. The disadvantage of this method is that the PDF must be calculated separately for every value of γ_s. Fits for emission size effects are then nonlinear and computationally expensive, in contrast with the simple linear fit enabled by Equation (18). Fortunately, the linear expansion of Equation (18) is easily adequate for the full expected range of N, even exceeding 10⁶, if γ_s ≪ 1. For γ_s ≳ 0.1, the linear approximation begins to differ substantially from the exact expression, but this region requires the appropriate quadratic term added to the original gain expansion as well.

Comparison of Equation (18) with Equation (11) shows that, at linear order and to the accuracy of the i.i.d. approximation for the pulse amplitudes, the effects of an extended emission region on the PDF of intensity are identical to the effects of temporal decorrelation, up to overall scale; this relationship is further analyzed in Section 4.3. Moreover, the modification arising from extended emission increases with the number of samples averaged, as was the case for temporal decorrelation (see Figure 3). These effects also depend on the gated signal-to-noise S ≡ 〈A〉I_s/I_n of the observation; see Figure 4. Indeed, the dependence on S, when combined with the Poisson noise, determines the resolution limit of an observation. We derive this limit explicitly in Section 4.6.

Zoom In Zoom Out Reset image size

4.2.1. Effects of a Non-static Emission Region

As we have already observed, comparison of Equation (11) with Equation (18) demonstrates that the linear corrections to the intensity PDF for effects of finite emission size and temporal decorrelation are perturbations of identical shape, and differing weights, at least within the excellent i.i.d. approximation for the averaged intensities. The relative strength of these two perturbation weights is

$$\begin{eqnarray} &&\frac{ \mathrm{size}}{\mathrm{decorrelation}} \approx 24 \frac{\Delta \tau _{\mathrm{d}}^2}{N(N-1)} \gamma _{\mathrm{s}}, \end{eqnarray} \tag{ 20 }$$

where we have substituted the square-law autocorrelation given in Section 4.1. Equation (20) allows immediate assessment of the relative importance of effects. For example, if the source extends over 5% of the diffractive scale, then effects arising from the finite emission size will dominate those from decorrelation of the scintillation pattern as long as the averaging is over less than a quarter of the diffractive timescale.

In practice, the estimates of parameters such as the source and background intensity are random variates, and are therefore subject to both noise and bias; we now derive the consequences of such errors. These effects, in contrast with effects in the previous sections, do not alter the theoretical PDF of intensity. However, they lead to similar changes in the model PDF, which mimic the effects of an extended emission region. They can thus lead to errors in the inferred parameters (γ_s and Δτ_d) that are distinct from the fundamental limits posed by sampling.

We have seen that the PDF modifications that arise from both an extended emission region and decorrelation of the scintillation within the averaging time are naturally represented in terms of partial derivatives of P(I; N|G) with respect to G, prior to integration with P(G). To express all the results of this section in the same way, we again analyze departures using the i.i.d. approximation for P(I; N|G). We also give the expected bias to γ_s that arises from each type of parameter error.

4.4.1. Effects of Biased Background or Source Amplitude

Many effects can bias the estimates of the background or source amplitudes. For example, leaked pulse power from quadrature downconversion is frequency-reversed and can be dispersed into the off-pulse region (Demorest 2007). Moreover, the analog-to-digital conversion of the signal requires quantization of the signal; this process introduces an offset of the background noise that varies with signal intensity (Jenet & Anderson 1998; Gwinn 2006).

A bias δI_n in the estimated background noise will result in a bias of the estimated pulse amplitudes δA_j = −δI_n/I_s, leading to a bias of the estimated mean intensity within a scintillation element $\delta (\bar{I} \equiv A I_{\mathrm{s}} G + I_{\mathrm{n}}) = \left(1 - G \right) \delta I_{\mathrm{n}}$. To obtain the incurred PDF modification, we apply the i.i.d. approximation and use that P(I; N|G) then only depends on the scale $\bar{I}$ and the degree of averaging N. We first expand to leading order in $\delta \bar{I}$, then re-express the derivative using $A I_{\mathrm{s}} \partial /\partial \bar{I} = \partial /\partial G$, and finally integrate G with its exponential weight to obtain

$$\begin{eqnarray} && P(I;N,\delta I_{\mathrm{n}} \,{=}\, 0) - P(I;N,\delta I_{\mathrm{n}})\nonumber\\ && \quad \approx - \frac{\delta I_{\mathrm{n}}}{A I_{\mathrm{s}}} \int _0^\infty dG\ \left(1 - G \right)e^{-G} \frac{\partial P(I;N,\delta I_{\mathrm{n}}\,{=}\,0|G)}{\partial G}.\nonumber\\ \end{eqnarray} \tag{ 21 }$$

Comparison with Equation (17) shows that this incurred modification is identical to a change in emission size γ_s → γ_s − δI_n/(2AI_s). A bias in pulse amplitudes is equivalent, with the substitution δI_n = −I_sδA_j. Thus, at leading order, emission size effects are indistinguishable from a biased estimate of the source or background intensities, regardless of averaging N. However, because the source emission size and the noise bias divided by signal-to-noise likely vary differently with pulse amplitude, segregating the pulses by strength provides a means to resolve this degeneracy.

4.4.2. Effects of Background or Source Parameter Noise

We similarly analyze the consequence of unbiased noise in the estimated source and background model parameters. This type of parameter error introduces spurious modulation in the corresponding model PDF, which must be compensated by artificial size effects. We therefore expect a positive bias on the inferred size.

Such parameter noise is frequently dominated by the source power estimate because of the limited sampling of the scintillation ensemble in each spectrum. This noise δA will be effectively independent of N, because we have assumed that the scintillation pattern remains approximately constant during the averaging, so each of the N consecutive pulses will be similarly biased. The mean intensity within each scintillation element is therefore biased: $\delta \bar{I} = \left(G - 1 \right) I_{\mathrm{s}} \delta A$. However, because this parameter noise is generally unbiased, the leading order correction after an ensemble average is quadratic in $\delta \bar{I}$. Again translating to derivatives with respect to G and averaging over an ensemble of pulse amplitudes, we obtain

$$\begin{eqnarray} && P(I;N,\delta A \,{=}\, 0) - P(I;N,\delta A)\nonumber\\ &&\quad \approx - \frac{1}{2} \left\langle \left(\frac{\delta A}{A} \right)^2 \right\rangle \int _0^\infty dG\ G^2 e^{-G} \frac{\partial ^2 P(I;N,\delta A\,{=}\,0|G)}{\partial G^2}.\nonumber\\ \end{eqnarray} \tag{ 22 }$$

This departure is similar to, though not strictly degenerate with, the effects of a finite emission size. This noise δA increases the inferred dimensionless emission size by γ_s → γ_s  +  〈(δA/A)²〉/4. The normalized variance is given by the reciprocal of the number of averaged scintillation elements, N_scint ∼ B/Δν_d, so γ_s experiences a positive bias ∼1/(4N_scint).

The analogous emission size bias from noise in the estimated background intensity is ∼1/(4S²N_b), where S is the gated signal-to-noise, and N_b is the number of samples averaged to estimate the background intensity. Observe that parameter noise for the background intensity, unlike parameter noise arising from scintillation, is independent for each spectrum and can therefore be mitigated by averaging over series of pulses.

Similarly with effects of bias in the model parameters, the near-degeneracy with effects of finite emission size is broken by comparing multiple analyses. For example, varying the bandwidth used to estimate the model parameters changes their noise; the evolution of measured residual structure with analyzed bandwidth thus precisely quantifies the influence of parameter noise.

Many instrumental distortions depend on the signal amplitude; they can therefore be modeled by imposing appropriate distortions on the distribution of scintillation gain P(G). The influence on the PDF can be modeled by combining the modified P(G) with the conditional distributions P(I; N|G) derived above. Alternatively, the modulation index m (Equation (15)) provides a simple measure to quantify changes in P(G); we therefore use m to characterize effects in this section.

Instrumental distortions can arise, for example, from saturation of the observing system. Because strong pulsars can be highly variable and dominate backgrounds when they are on, optimal quantization levels are often difficult to determine, and even systems with many bits may saturate. A simplified model of saturation imposes a voltage cutoff on the observed scalar electric field in the time domain:

$$\begin{eqnarray} \hat{x}_i = \left\lbrace \begin{array}{@{}l l@{}}x_i & \quad {\rm if\ |x_i| < \sigma _{\mathrm{max}}},\\ \sigma _{\mathrm{max}}\mathrm{sign}(x_i) & \quad {\rm else.}\\ \end{array} \right. \end{eqnarray} \tag{ 23 }$$

The effect of this cutoff is twofold: power from the strongest samples is spread across all channels in the spectral domain, and the observed signal intensity is underestimated. At leading order, the second effect predominates, so that the estimated off-pulse noise overestimates the effective on-pulse noise. In terms of the discussion in Section 4.4.1, saturation leads to a bias δI_n < 0 and thus, an inferred emission size that is upward biased.

Likewise, if the pulse is heavily dispersed, then the saturation predominantly affects samples with the highest scintillation gain G. As a simplified model of this effect, we replace the distribution of scintillation gain P(G) = e^−G by a truncated version: $\hat{P}(G) \equiv e^{G}\theta \left(G_{\mathrm{max}} - G \right) + e^{-G_{\mathrm{max}}} \delta \left(G-G_{\mathrm{max}} \right)$, where θ(x) is the Heaviside step function. The corresponding modulation index is

$$\begin{eqnarray} &&\hat{m}^2 = \frac{\langle G^2 \rangle _{\hat{P}}}{\left\langle G \right\rangle _{\hat{P}}^2} - 1 = \frac{\sinh (G_{\mathrm{max}}) - G_{\mathrm{max}}}{\cosh (G_{\mathrm{max}}) - 1} \le 1. \end{eqnarray} \tag{ 24 }$$

Comparison with Equation (15) shows that even a relatively large cutoff G_max = 5 gives suppression equivalent to modest emission size effects, γ_s ≈ 0.014. Furthermore, the exponential character of the incurred suppression leads to an extremely rapid onset of saturation artifacts once the largest intensities clear the threshold for saturation; the effects on the strongest pulses may therefore be quite pronounced relative to those on the weaker pulses.

Poisson noise presents a fundamental limit on the resolving power afforded by scintillation, as a fraction of the magnified diffractive scale. This limit depends on both the gated signal-to-noise ratio S ≡ 〈A〉I_s/I_n and the observational parameters B and T_obs. Figure 4 illustrates the variation of PDF modification with S, for N = 1. The modification from extended emission size is maximal at zero intensity, with fractional strength

$$\begin{eqnarray} &&\left|\frac{P(I\,{=}\,0;N\,{=}\,1,\gamma _{\mathrm{s}}) - P(I\,{=}\,0;N\,{=}\,1,\gamma _{\mathrm{s}}\,{=}\,0)}{P(I\,{=}\,0;N\,{=}\,1,\gamma _{\mathrm{s}}\,{=}\,0)} \right| \nonumber\\ &&\quad = -2 \gamma _{\mathrm{s}} \left[ 2 + \frac{1}{S} + \left(1+S\right) \frac{e^{-1/S}}{\mathrm{Ei}\left(-1/S \right)} \right], \end{eqnarray} \tag{ 25 }$$

where Ei(x) is the imaginary error function. Equating Equation (25) with the level of Poisson noise η in the observed PDF then provides an excellent assessment of resolving power.

We apply an approximation for Equation (25) to obtain the following estimate for the achievable resolution, $\sqrt{\gamma _{\mathrm{s},\mathrm{min}}}$, in units of the magnified diffractive scale:

$$\begin{eqnarray} &&\sqrt{\gamma _{\mathrm{s},\mathrm{min}}} \approx \sqrt{\frac{\eta }{S}} \exp [0.2 + 0.09 \left(\ln S\right)^2 + 0.005 \left(\ln S \right)^3 ].\qquad \end{eqnarray} \tag{ 26 }$$

This approximation is good to within 10% for 10⁻³ < S < 2. Simulations indicate that substituting $\eta \equiv ({1.3}/{\sqrt{N_{\mathrm{tot}}}})$ effectively estimates the standard errors from the full distribution fits, where N_tot is the total number of sampled points used to estimate the PDF of intensity.

We now apply this result to estimate the resolving power for typical observational parameters. A $(1\mathrm{\ {\rm hr}})\,{\times}\,(4\mathrm{\ MHz})$ observation of a pulsar with 5% duty cycle gives ∼10⁹ on-pulse samples of the flux density. A gated SNR S = 1 then gives a 3σ detection limit of about 1% of the magnified diffractive scale, S = 0.1 gives about 10% of the magnified diffractive scale, and S = 0.01 gives about 50% of the magnified diffractive scale. Since the detection limit is much more sensitive to S than to the number of sampled points, an analysis of a subset of strong pulses may provide optimal resolution.

In Section 4.2, we demonstrated that the size of the emission region affects the modulation index m. However, self-noise, temporal decorrelation, and background noise also affect the modulation index. These effects are readily evaluated, either by analysis of the random variables that make up each sample of intensity (Equation (1)), or by appropriate derivatives of the corresponding characteristic function.

Here, we analyze only the case for which the pulse amplitudes are i.i.d.: 〈A_iA_j〉 = 〈A〉² + δ_ij〈δA²〉. This condition is different than the assumptions for the i.i.d. approximation derived above, which assumes the stronger condition that the spectral statistics are i.i.d. Although many pulsars appear to have i.i.d. pulse amplitudes, others display striking departures. For example, phenomena such as nulling and mode-changing break statistical isotropy, and the brightest pulses from the Vela pulsar tend to occur in groups (Palfreyman et al. 2011).

To explicitly account for the temporal decorrelation of the scintillation pattern, we again assign the autocorrelation appropriate for square-law scattering, quantified by the decorrelation timescale (in pulse periods) Δτ_d. We incorporate effects of an extended emission region in terms of the dimensionless size parameter γ_s, defined by Equation (14), with the assumption that γ_s ≪ 1. Modification of the scattering assumptions will introduce minor scaling alterations for these contributions.

To leading order in the effects of emission size and temporal decorrelation, the modulation index is then given by

$$\begin{eqnarray} m^2 &\equiv& \frac{ \langle (I - I_{\mathrm{n}})^2 \rangle }{ \left\langle I - I_{\mathrm{n}} \right\rangle ^2} - 1 \nonumber \\ \nonumber&\approx & \ 1 - \overbrace{4 \gamma _{\mathrm{s}} \left(1 + \frac{1}{N} \right)}^{\mbox{\hspace{-2.45274pt}{Emission Size\hspace{-3.45274pt}}}} \ +\ {\overbrace{\frac{2}{N}}^{\mbox{\hspace{-3.45274pt}{Self-Noise}\hspace{-3.45274pt}}}} \ +\ \overbrace{\frac{1}{N S}\left(2 + \frac{1}{S} \right)}^{\mbox{\hspace{-12.45274pt}{Signal-to-Noise}\hspace{-12.45274pt}}}\\ && + \underbrace{4 \left(1 - 2\gamma _{\mathrm{s}} \right) \frac{\langle \delta A^2 \rangle }{N \left\langle A \right\rangle ^2}}_{\mbox{\hspace{-28.45274pt}{Intrinsic Variability}\hspace{-28.45274pt}}} \ -\ \underbrace{\left(1 - 4 \gamma _{\mathrm{s}} \right) \frac{N^2-1}{6 \Delta \tau _{\mathrm{d}}^2}}_{\mbox{\hspace{-28.45274pt}{Decorrelation}\hspace{-28.45274pt}}}.\qquad \end{eqnarray} \tag{ 27 }$$

Here, we have grouped the contributions to the modulation index, to explicitly illustrate the relative influence of the various effects. We also have used γ_{s, 1} = γ_{s, 2} ≡ γ_s because of the degeneracy at linear order. We see that the relative contributions of an extended emission region and decorrelation of the scintillation pattern give the previously derived relationship between their respective PDF modifications (Section 4.3).

The modulation index corresponds to the zero time and frequency lags in the autocorrelation function of intensity; analysis of non-zero lags eliminates many of the terms in Equation (27). For example, non-zero spectral lags eliminate the "Signal-to-Noise" contribution and mitigate the "Self-Noise" contribution (although intrinsic modulation induces correlations in the spectral self-noise that preserve the contribution). Non-zero temporal lags fully eliminate all terms other than "Emission Size" and "Decorrelation," which is enhanced.

The self-noise of a signal is a useful diagnostic, which can be used to identify intrinsic variability on timescales shorter than the spectral accumulation time. Intrinsic modulation induces correlations among frequencies, without modifying the mean spectrum (Gwinn & Johnson 2011). If spectra are formed according to the criteria in Section 2.3, then the power in each spectral channel is an exponential random variable. Because such variables are completely characterized by their mean, the single-channel statistics are therefore immune to intrinsic modulation. However, the intrinsic modulation does introduce correlations in the noise of nearby channels. We now derive the expected self-noise in the spectral domain, thereby enabling detection of such correlations.

One method for estimating self-noise is to compare pairs of nearby samples that are assumed to be within the same scintillation element. We will assume that the pair consists of samples with uncorrelated self-noise (e.g., pairs from different pulses, or from the same pulse with negligible intrinsic modulation). Because self-noise is heteroscedastic (see Section 3.2.3), we calculate the noise δI² ≡ (I₁ − I₂)²/2 in pairs of samples {I₁, I₂} as a function of their mean I = (I₁ + I₂)/2; we denote the resulting noise estimate δI²(I). The ensemble average is given by

$$\begin{eqnarray} &&\delta I^2(I) = \frac{\int _0^{2I} 2(I - I^{\prime })^2 P(I^{\prime }, 2 I - I^{\prime };N) dI^{\prime }}{\int _0^{2I} P(I^{\prime }, 2 I - I^{\prime };N) dI^{\prime }}, \end{eqnarray} \tag{ 28 }$$

where P(I₁, I₂; N) represents the probability of sampling the pair {I₁, I₂}, including an ensemble average over scintillation. Because we assume that the samples are independent but are drawn from the same scintillation element, P(I₁, I₂; N) = ∫^∞₀dG P(I₁; N|G)P(I₂; N|G).

As usual, we apply the i.i.d. approximation. Although the scintillation integral cannot be evaluated in closed form, Equation (A2) readily yields the identity

$$\begin{eqnarray} && \int _0^{2I} 2(I - I^{\prime })^2 P(I^{\prime }, 2 I - I^{\prime };N) dI^{\prime }\nonumber \\ &&\quad = \frac{I^2}{N + \frac{1}{2}} \int _0^{2I} P(I^{\prime }, 2 I - I^{\prime };N) dI^{\prime }. \end{eqnarray} \tag{ 29 }$$

This relationship immediately unveils the simple relationship between the estimated mean signal and noise:

$$\begin{eqnarray} &&\delta I^2(I) = \frac{I^2}{N + \frac{1}{2}}. \end{eqnarray} \tag{ 30 }$$

In fact, because the identity given by Equation (29) holds regardless of the integration over G, this property of the estimated self-noise is independent of the degree of sampling and the properties of scintillation. We see that this method for estimating the self-noise gives the same form as the exact expression (see Section 3.2.3), except that N → N + 1/2.

If the averaged intensities are from different pulses, then pulse-to-pulse variations contribute additional noise. If the averaged intensities are from the same pulses, then intrinsic variations on timescales shorter than t_acc induce correlations in self-noise, and thereby decrease the measured noise; if N > 1 then pulse-to-pulse variations within the averaging will increase the noise. Each of these effects can be estimated by applying the appropriate substitution N → N_eff. Similar tests were applied by Gwinn et al. (2011) to infer short-timescale variability (<300 μs) of PSR B0834+06.

Consider an exponential random variable I: $P(I;\bar{I}) = \bar{I}^{-1} e^{-I/\bar{I}}$ for I > 0, where $\bar{I} \equiv \left\langle I \right\rangle$. The corresponding characteristic function is $\varphi (k) = (1 - i \bar{I} k)^{-1}$. The distribution $P(I;\lbrace \bar{I}_j\rbrace )$ of the average of N such variables from independent distributions with different scales $\lbrace \bar{I}_j\rbrace$ is then easily derived by inverting the partial fraction decomposition of the product of characteristic functions:

$$\begin{eqnarray} &&P(I;\lbrace \bar{I}_j\rbrace ) = N \sum _{j=1}^N \left(\frac{\bar{I}_j^{N-2} }{\prod \limits _{{{\scriptstyle {\begin{array}{@{}l@{}}\ssty \ell = 1\\ \ssty \ell \ne j\end{array}}}}}^N (\bar{I}_j - \bar{I}_\ell )} \right) e^{-N I/\bar{I}_j},\quad I>0.\nonumber\\ \end{eqnarray} \tag{ A1 }$$

In other words, the distribution of the average of exponential random variables is a weighted sum of the marginal distributions. This sum is the hypoexponential (or generalized Erlang) distribution.

The representation of Equation (A1) is singular if any of the scales are equal. In particular, if the random variables are i.i.d., each with scale $\bar{I}$, then their average is drawn from an Erlang distribution, which is simply a gamma distribution with integer shape parameter:

$$\begin{eqnarray} &&P(I;\bar{I},N) = \frac{N^N}{(N-1)!} \frac{I^{N-1}}{\bar{I}^N} e^{-\frac{N I}{\bar{I}}},\quad I>0. \end{eqnarray} \tag{ A2 }$$

Goodman (1985) discusses both these distributions in the similar context of polarized thermal light.

We now approximate the distribution of N correlated exponential random variables {I_i}. The exact distribution follows by first constructing a complex multivariate normal distribution with the desired covariance matrix; the exponential random variables are the squared norms of each circular complex Gaussian random variable. We restrict ourselves to the case in which each random variable has unit mean. The characteristic function of P({I_i}) can then be written (Mallik 2003)

$$\begin{eqnarray} &&\varphi \left(\lbrace k_i \rbrace ; \boldsymbol{\Sigma } \right) = \frac{1}{\mathrm{det}[ \textbf {I} - i \mathrm{diag}\left\lbrace k_i \right\rbrace \boldsymbol{\Sigma } ]}. \end{eqnarray} \tag{ A3 }$$

Here, $[\boldsymbol{\Sigma }]_{ij} = \sqrt{ \langle I_i I_j \rangle - 1 } \equiv \sqrt{\Gamma _{ij}}$.

We now expand in $\boldsymbol{\epsilon } \equiv \textbf {J} - \boldsymbol{\Sigma }$, where $\textbf {J}$ is the N × N unit matrix: $[\textbf {J}]_{ij} = 1\ \forall \ i,j$. Since $[\boldsymbol{\epsilon }]_{ij} \approx (1 - \Gamma _{ij} )/2$, we can approximate the characteristic function to leading order as

$$\begin{eqnarray} \varphi \left(\lbrace k_i \rbrace ; \boldsymbol{\Sigma } \right) &=& \frac{1}{1 - i \sum _{j=1}^N k_j} + \frac{\sum _{i<j} k_i k_j (1 - \Gamma _{ij})}{\big(1 - i \sum _{j=1}^N k_j\big)^2}\nonumber\\ && + \mathcal {O}([ 1 - \Gamma _{ij} ]^2 ). \end{eqnarray} \tag{ A4 }$$

The inverse Fourier transform of $\varphi \left(\lbrace k_i \rbrace ; \boldsymbol{\Sigma } \right)$ then gives our approximation to P({I_i}). Although this inverse can be written in terms of delta functions and derivatives, the most natural form expresses the action as a weight function:

$$\begin{eqnarray} &&\int d^N I_i \ f(\lbrace I_i \rbrace) P(\lbrace I_i \rbrace)\nonumber\\ &&\quad = \int _0^\infty dI \ \left[ f(\lbrace I\rbrace) e^{-I} - \sum _{i<j} (1 - \Gamma _{ij} ) I e^{-I} \left.\frac{\partial ^2 f}{\partial I_i \partial I_j}\right\rfloor _I\right.\nonumber\\ &&\qquad \left.+\, \mathcal {O}([ 1- \Gamma _{ij} ]^2 ) \right]. \end{eqnarray} \tag{ A5 }$$

Here f({I}) denotes f({I_i = I ∀ i}).

A.2.1. Example: Bivariate PDF of Correlated Exponential Random Variables

The joint distribution of a pair of correlated exponential random variables is particularly simple. In this case, the characteristic function (A3) is easily inverted:

$$\begin{eqnarray} P(I_1,I_2;\Gamma) &=& \frac{1}{1-\Gamma } I_0\left(\frac{2\sqrt{\Gamma }}{(1-\Gamma)} \sqrt{I_1 I_2} \right) \nonumber\\ &&\times \exp \left[ -\frac{1}{(1-\Gamma)} \left(I_1 + I_2 \right) \right],\quad I_1,I_2>0.\nonumber\\ \end{eqnarray} \tag{ A6 }$$

Here, Γ = 〈I₁I₂〉 − 1, and I₀(x) is the modified Bessel function. Using Equation (A5), we can approximate Equation (A6) as

$$\begin{eqnarray} P(I_1,I_2;\Gamma) &=& e^{-\hat{I}} \left\lbrace \delta (\Delta I) + \left(1 - \Gamma \right) \left[ \frac{2 - \hat{I}}{4} \delta (\Delta I) + \hat{I} \delta ^{\prime \prime }(\Delta I) \right] \right.\nonumber\\ &&\left.+\, \mathcal {O}[(1 - \Gamma )^2 ]\right\rbrace ,\quad \hat{I}>0, \end{eqnarray} \tag{ A7 }$$

where $\hat{I} \equiv (I_1 + I_2)/2$ and ΔI ≡ I₁ − I₂.

A.2.2. Convolution of Correlated Exponential Random Variables

The gamma distribution is an excellent approximation for the sum of independent exponential random variables that are not necessarily isotropic (see Section 3.2) and can also be used to approximate the convolution of correlated exponential random variables. This approximation is most effective for correlations that are particularly strong or weak. For example, consider the average of two exponential random variables, each with unit mean, and with covariance 〈I₁I₂〉 = 1 + Γ. The gamma distribution that matches the mean and variance of the exact average has a corresponding effective number of degrees of freedom N_eff = 2/(1 + Γ). These distributions agree exactly for the limits Γ → 0, 1, but suffer sizeable discrepancies for intermediate correlations. For example, the third central moment of the average of two correlated exponential random variables is (1 + 3Γ)/2 but is 2/N²_eff = (1 + Γ)²/2 for the corresponding gamma distribution—a difference of ∼10% for intermediate correlation.

More generally, for N averaged exponential random variables with correlation matrix Γ_ij, as defined above, and unit mean, the appropriate number of degrees of freedom based on variance matching is N_eff = N(1 + 2N⁻¹∑_{i < j}Γ_ij)⁻¹. Even in the case of large N, the disparity in skewness between the exact average and the appropriate gamma distribution remains. For instance, if Γ_ij ≡ Γ for i ≠ j then,

$$\begin{eqnarray} &&\mu _3 = \left\lbrace \begin{array}{@{}ll@{}}2\Gamma ^{3/2} + \dsty\frac{6\Gamma }{N} (1 - \sqrt{\Gamma } )&\\ \quad + \dsty\frac{2}{N^2} (1 - 3 \Gamma + 2 \Gamma ^{3/2} ) & {\rm Exact\; Average,} \\\ms 2\Gamma ^2 + \dsty\frac{4\Gamma }{N} \left(1 - \Gamma \right) + \frac{2}{N^2}\left(1-\Gamma \right)^2 &\!\!\!\!{\rm Approximation\; with }\\ &N_{\mathrm{eff}} = \frac{N}{1 + (N-1)\Gamma }{,} \end{array}\right.\nonumber\\ \end{eqnarray} \tag{ A8 }$$

where μ₃ denotes the third central moment. Again, we see that these expressions agree for Γ → 0, 1 as expected, but differ substantially for intermediate correlation, even as N → ∞.