Resolving the Emission Regions of the Crab Pulsar's Giant Pulses II. Evidence for Relativistic Motion

The Crab Pulsar is the prime example of an emitter of giant pulses. These short, very bright pulses are thought to originate near the light cylinder, at $\sim\!1600{\rm\;km}$ from the pulsar. The pulsar's location inside the Crab Nebula offers an unusual opportunity to resolve the emission regions, using the nebula, which scatters radio waves, as a lens. We attempt to do this using a sample of 61998 giant pulses found in coherently combined European VLBI network observations at $18{\rm\;cm}$. These were taken at times of relatively strong scattering and hence good effective resolution, and from correlations between pulse spectra, we show that the giant pulse emission regions are indeed resolved. We infer apparent diameters of $\sim\!2000$ and $\sim\!2400{\rm\;km}$ for the main and interpulse components, respectively, and show that with these sizes the correlation amplitudes and decorrelation timescales and bandwidths can be understood quantitatively, both in our observations and in previous ones. Using pulse-spectra statistics and correlations between polarizations, we also show that the nebula resolves the nanoshots that comprise individual giant pulses. The implied diameters of $\sim\!1100{\rm\;km}$ far exceed light travel-time estimates, suggesting the emitting plasma is moving relativistically, with $\gamma\simeq10^{4}$, as inferred previously from drifting bands during the scattering tail of a giant pulse. If so, the emission happens over a region extended along the line of sight by $\sim\!10^{7}{\rm\;km}$. We conclude that relativistic motion likely is important for producing giant pulses, and may be similarly for other sources of short, bright radio emission, such as fast radio bursts.


INTRODUCTION
The Crab Pulsar (PSR B0531+21) is the remnant of supernova SN 1054 and the central star in the Crab Nebula. It powers the pulsar wind nebula (PWN) that fills the interior of the Crab Nebula. The PWN is expanding into the freely expanding supernova ejecta (Chevalier 1977), sweeping up ejecta into a dense thin shell of material. This shell is subject to Rayleigh-Taylor instabilities leading to the filamentary structure seen in optical images of the Crab Nebula (e.g., Chevalier & Gull 1975;Jun 1998) and the observed acceleration of the Crab filaments (Trimble 1968).
The pulsar itself was discovered by Staelin & Reifenstein (1968) through the detection of individual radio pulses and has since been studied extensively (for a review, see . Its mean radio profile shows seven components, the dominant components of which are the main pulse (MP) and the interpulse (IP) at frequencies < 4 GHz, which are separated by ∼ 145 • in rotation phase.
The Crab Pulsar is unusual among pulsars in that it shows "giant pulses", extremely narrow and bright pulses that occur randomly within the phase windows of its MP and IP components (Hankins et al. 2003). At lower observing frequencies (< 4 GHz), the properties of the MP and IP giant pulses are quite similar, with pulses typically lasting overall a few microseconds, but comprised of numerous nanoshots that have durations down to the time resolution with which they have been observed, some clumping together in microbursts (Sallmen et al. 1999;. At radio frequencies above ∼ 4 GHz, there is also an interpulse (Moffett & Hankins 1996). But it is offset by ∼ 10 • in phase from the IP at low frequencies (∼ 0.3−3.5 GHz) and has properties sufficiently different that it is almost certainly a different pulse component (e.g., Hankins et al. 2016).
The emission mechanism of giant pulses is unknown. However, the strong alignment of the MP and (low frequency) IP components from radio to γ observations to within ∼ 2 ms (Moffett & Hankins 1996) offers empirical evidence that both the radio and high energy emissions emanate from the same spatial region. For the MP, this is strengthened by the correlation between radio giant pulses and optical (Shearer et al. 2003(Shearer et al. , 2012) and X-ray pulses (Enoto et al. 2021). Since pair production strongly absorbs γ-rays near the polar cap, giant pulse emissions likely occurs in other magnetospheric gaps (Romani & Yadigaroglu 1995;Muslimov & Harding 2004;Harding et al. 2008) or in regions beyond the light cylinder, r LC = cP/2π ≈ 1600 km, (Lyubarskii 1996;Philippov et al. 2019) where plasma instabilities and magnetic reconnection can occur.
To help understand the emission mechanism, it would be useful to have constraints on the locations and sizes of the emission regions. For this purpose, it may be possible to use the fact that pulsars scintillate, displaying intensity variations in time and frequency because of multipath propagations of their radio emission: if one can retrace those paths, they act like an interstellar interferometer with which the emission can be studied at extreme angular resolution. Cordes et al. (1983) first used scintillation in this way to probe pulsar magnetospheres, by inferring limits on the transverse separations of the emitting regions of PSR B0525+21 and PSR B1133+16. Since then, this method has been applied to several other pulsars (Wolszczan & Cordes 1987;Smirnova et al. 1996;Gupta et al. 1999;Pen et al. 2014).
For the Crab Pulsar, observations of its scintillation show a fairly stable angular broadening, and a highly variable temporal broadening (Rankin & Counselman 1973). VLBI measurements favor two screen locations with angular broadening originating from the interstellar medium (ISM) and (most) temporal scattering originating at the filaments inside the Crab Nebula (Vandenberg 1976). For resolving the Crab's emission, it is the nebular screen that is most relevant, as it is closer to the pulsar and thus gives higher spatial resolving power. Hence, we focus on the nebular screen in this paper.
For individual Crab giant pulses, if they originated from the same physical location, one would expect that their radiation would follow the same paths through the scattering regions and thus that they would have imprinted on them the same impulse response function (IRF). Hence, one would expect power spectra of pulses close in time to correlate strongly. Indeed, Cordes et al. (2004) found that pulses close in time correlated strongly, with an average correlation coefficient of 1/3, as would be expected for pulses that had different intrinsic frequency structure -associated with the random nanoshots that they are comprised off -but that were imprinted with the same IRF. Similar levels of correlation were also found by Karuppusamy et al. (2010) for spectra of MPs and IPs in the same pulse rotation, as well as for spectra between different microbursts of individual MPs.
In contrast, Main et al. (2021, hereafter Paper I) found a surprisingly low correlation coefficient of ∼ 2% when correlating power spectra of nearby MP-MP pairs and MP-IP pairs. They suggested this could result if individual giant pulses arose in different parts of an extended emission region, one larger than the resolution of the scattering screen. If so, then the nebular screen must have had lower effective resolution during the observations of Cordes et al. (2004) and Karuppusamy et al. (2010); this is indeed consistent with those observations having had smaller scattering time when scaled to our observing frequency (see derivation in Section 2). With the large sample of MPs and IPs in Paper I, it was also possible to see differences between the scintillation patterns of the two components, hinting at a projected physical separation between the emitting regions of 50−400 km.
Another result found in Paper I, the importance of which we realized only later, was that the mean power in giant pulse spectra was approximately equal to the standard deviation. This is not consistent with multiple nanoshots going through the same screen, as in that case the standard deviation should be √ 3 times larger than the mean (making the usual assumption that the IRF has Gaussian statistics; Rickett 1990). Instead, it suggests that the emission locations of the individual nanoshots are also separated sufficiently for them to be resolved by the nebular screen, and that, therefore, the scintillation pattern imprinted on each nanoshot is different. This has the paradoxical implication that the nanoshots, which occur within a few µs in the giant pulse and thus must be causally related to each other, appear to have separations on the sky in excess of the ∼ 300 km resolution found in Paper I, i.e., far more than can naively be understood from light-travel time arguments.
In this paper, we use another dataset to get a more complete picture of the Crab's emission regions. We begin with a brief review on how scattering screens close to pulsars can be seen as magnifying lenses in Section 2, and infer the expected resolution for the Crab Pulsar. In Section 3 we de-scribe our multi-telescope data sets and in Section 4 their reduction, including the coherent combination of the data from different telescopes. We measure scattering times from our data in Section 5, and correlations between pulses and polarizations in Section 6. The latter show that in our data both giant pulses and their constituent nanoshots are imprinted with different IRFs. In Section 7, we discuss the ramifications for the sizes of the emission regions and for the screen. We address the paradox of the resolved nanoshots and conclude it is most easily resolved if the plasma emitting the pulses moves highly relativistically, as has also been suggested by Bij et al. (2021) based on drifting frequency structure in the scattering tail of giant pulses. We finish with an outlook in Section 8.

INTERSTELLAR INTERFEROMETRY
A scattering screen can be seen as a lens, which yields a certain resolution at the pulsar, over which the interference pattern changes by order unity. The resolution depends on where the scattering occurs, improving as the screen is placed closer to the pulsar. If the probability of scattering is normally distributed around the line of sight -which yields an exponential scattering tail as is often observed -the variance of the distribution is related to the exponential decay time τ by, where d s is the distance to the scattering screen and d eff = d p d s /(d p − d s ) the effective distance, with d p the distance to the pulsar. For these normally distributed scatterers, the corresponding angular resolution is λ/2πσ L , where λ is the observing wavelength. Hence, the physical resolution at the pulsar is given by (Cordes & Rickett 1998), In the Crab Nebula, as discussed in Paper I, the only conceivable location for scattering is in the optically emitting filaments, as only their densities are high enough (of order n e ≈ 10 3 cm −3 Osterbrock 1957). Lawrence et al. (1995) and Martin et al. (2021) show that these filaments reside within 0.5 − 2 pc of the pulsar and Trimble (1973) give a range of distances to the pulsar from 1.4 − 2.7 kpc. For an estimate of the resolution σ x of the lens at the pulsar, we use a nominal pulsar distance of d p = 2 kpc and distance of the screen from the pulsar d p − d s = 1 pc along with a geometric time delay τ = 1 µs at λ = 18 cm (Paper I), to infer σ x ≈ 290 km.
The above only holds if the scattering time is dominated by delays in the nebula. It is known, however, that there is a contribution from the interstellar medium, which may be important at times that the scattering in the nebula is relatively weak (Rankin & Counselman 1973;Vandenberg 1976). We can estimate its contribution by assuming that the interstellar scattering time does not vary much, so that an upper limit to its contribution is set by the lowest observed scattering times. Losovsky et al. (2019) find that at 111 MHz the scattering time varies between 10 and 115 ms, while Serafin-Nadeau et al. (in preparation) find that at 600 MHz the lowest scattering times are ∼ 10 µs. Scaling with ν −4 , both results imply an upper limit of ∼ 0.2 µs at λ = 18 cm. If the observed scattering time is close to this range, the resolution due to the nebular screen will be poorer than inferred from Equation 2.
We can verify the above result by noting that the interstellar screen dominates the angular broadening of the Crab Pulsar, and at λ = 18 cm, Rudnitskii et al. (2016) measure full width at half maximum 0.5 w 1.3 mas. The range likely reflects that the screen is not isotropic, but for an estimate we nevertheless assume an isotropic screen with a normal distribution, such that the lens size is σ L = d s (w/ √ 8 log 2). For a screen halfway to the pulsar, Equation 1 then would imply a range in scattering time of 0.2 τ 1.5 µs. In reality, since the Crab is relatively far above the Galactic plane, the interstellar screen is likely closer to us than halfway, which would make the result more consistent with the upper limit above. Assuming that the actual value of the interstellar scattering time is close to the upper limit, τ = 0.2 µs, then to get the maximum observed angular broadening, w = 1.3 mas, requires a screen distance of 0.24 kpc. Without any nebular contribution, this geometry would imply a resolution σ x 80, 000 km.

OBSERVATIONS
We analyse a total of 7.31 hr of European VLBI Network (EVN) dual-polarization data, taken at four epochs between 2015 Oct and 2017 May (see Table 1). For our analysis, we use data only from the up to 8 telescopes that had relatively clean signal in both polarizations and covered the full frequency range of 1594.49 − 1722.49 MHz. At each telescope, real-sampled data in both circular polarizations were recorded in either 2-bit MARK 5B or VDIF format, covering the frequency range in either eight contiguous 16 MHz wide bands or four contiguous 32 MHz wide bands. During each observation run, the telescopes regularly switched to calibrator sources resulting in short breaks in our data.
The use of multiple telescopes at baselines of up to 10, 000 km as an interferometer allows us to achieve a resolution of ∼ 4 mas. For our purposes here, this is useful as it resolves the radio-bright nebula of angular size ∼ 6 × 4 , the dominant source of noise. The angular resolution is not high enough, however, to resolve the nebular or interstellar scattering screens (which at our frequency have sizes of ∼ 0.005 In order to combine individual telescope data coherently and obtain sensitive measurements of giant pulses, we first need to align the voltage data in both time and phase. The largest delays come from differences in path length travelled by the signal to each telescope. We used the software program CALC10 1 (Ryan & Vandenberg 1980) through a wrapper from the Super FX Correlator (SFXC; Keimpema et al. 2015) to calculate these geometric delays. The geocentric frame was chosen as the reference frame and the geometric delays were generated at one second intervals which we interpolate with an Akima (1970) spline in our beamformer pipeline. The second largest delays come from differences between each telescope's local clock. We obtained the clock offset and rate information from the post-observation VEX file 2 which were obtained from standard VLBI clock searching techniques. With these corrections, giant pulses common between telescopes are aligned to within ∼ 10 ns (see Figure 1). Lastly, cables and electronic components, and the atmosphere also introduce time delays and phase rotations. Since the time-averaged emission of the Crab Pulsar is not very bright relative to the nebula, we determined these time delays and phase rotations using just the giant pulses themselves.
We present our beamformer pipeline in Section 4.1. In Section 4.2 and Section 4.3 we describe how we identified giant pulses from our data, and how we then use these giant pulses to find fringe solutions respectively. We present results of our coherently combined data in Section 4.4.

Beamformer Pipeline
Our pipeline closely follows that used for the Large European Array for Pulsars (LEAP), as described by Bassa et al. (2016); Smits et al. (2017) and SFXC as described by Keimpema et al. (2015).
We first bring the signal from each telescope to the geocentric frame using the pre-determined geometric and clock delays. Given a geocentric time t geo , the corresponding time at a telescope is where τ geo + clock, tel is the geometric and clock delay for the telescope of interest. We use BASEBAND (Van Kerkwijk et al. 2020) and PULSARBAT (Mahajan & Lin 2022) to read in the baseband data to the nearest integer time sample, flip lower-sidebands so that frequencies are all in increasing order, and convert the real-sampled data to complex 3 . This realto-complex conversion shifts the frequency from the edge of each sub-band to the center. Since the integer geometric and clock delay compensation was applied when reading in the data (before the real-to-complex conversion), we need to shift the integer delay compensation frequency to the center of the sub-band. This was done by applying the shift where N int geo + clock, tel is the geometric and clock delay in integer number of samples. The residual fractional delay manifests itself as a phase error across the frequency band and is removed by rotating the signal by where N frac geo + clock, tel is the fractional delay, ν is the baseband frequency and SR is the sampling rate. Due to varying radial velocities between the different telescopes, the signals are Doppler shifted relative to each other; we correct for this frequency shift by a phase rotation in the time domain, whereτ geo + clock, tel is the time derivative of the geometric and clock delays (which encodes the Doppler shift), and ν sky is the central observing frequency of the sub-band. Lastly, the above time shifts in the down-converted baseband signal does not correct for the phase rotation of the radio waves at each telescope at the central observing frequency. We corrected for this by rotating the phase of the signal by, φ = 2πτ geo + clock, tel ν sky , a process called fringe stopping.
In the corrections outlined above, only single values of the delay and its rate are used. This is only acceptable at short time scales where the delay is approximately constant. Thus, we apply our corrections in data chunks of no more than 128 samples (8 us), which ensures drifts in delay of less than 0.000192 sample (0.012 ns) even for the maximum, ∼ 1.5 µs/s delay rate on terrestrial baselines.
After applying the geometric and clock delays, we removed radio frequency interference (RFI) from our data in chunks of 2 24 samples (1.048576 s). We start by removing obvious RFI and those mentioned in each telescope's log files. To detect residual RFI, we first channelized our data into 8192 frequency channels per sub-band. We then normalized these by the square root of the time-averaged power spectrum, thus correcting for the fact that the passband is not perfectly flat, with roll-offs at the edges of the band and some other structures. RFI spikes and highly variable channels were flagged by comparing with a 128-channel median-filter (0.25 MHz), removing all signals above a 5σ cut-off. We also normalized time-variability resulting from instrumentation by dividing by the square root of the frequency averaged power spectrum smoothed over 8.192 ms.
Next, we coherently de-dispersed the data. For the dispersion measure (DM), we started with initial guesses from Crab monitoring data 4 (Lyne et al. 1993), and then adjusted the value to ensure the profiles of bright giant pulses were aligned in frequency, leaving us with final values listed in Table 1.
At this point, we incoherently summed the data from each telescope and searched for giant pulses (see Section 4.2) as the telescope data are aligned to within a sample in time.
Once we have modelled the delays and phase rotations resulting from instruments and atmospheric variations using giant pulses (using Effelsberg as a reference; see Section 4.3), we applied the fringe delay and phase solutions to the baseband data. After all the delay and phase corrections are applied we weighed the data from each telescope by the amplitude of their gains to ensure maximum S/N in the coherently added baseband data (given that the nebula dominates the system temperature, the different telescopes have very similar gain; see below). A final normalization of the coherently summed voltage data was performed such that the intensity of the noise level is unity in each sub-band and polarization. Finally, we applied a parallactic angle correction to the coherent beam, as appropriate for the Crab at our reference telescope (Effelsberg). The rotation measure (RM) towards the Crab Pulsar is ∼ 45 rad/m 2 (Sobey et al. 2019). Thus, at our observing wavelength of 18 cm the expected phase rotation across our total bandwidth is ∼ 0.1 rad. Since this is small and we do not detect RM or cable delay in our sub-bands, we do not correct for it.

Giant Pulse Search
We searched for giant pulses in two separate passes, the first using incoherently summed data to determine fringe solutions and the second using coherently added data. In both passes we summed the power over both polarizations and all eight sub-bands. Peaks above 8σ in a 16 µs wide running average of the intensity time stream (with the window size roughly matched to the typical width of a giant pulse) were flagged as potential giant pulses. We define the start time of a giant pulse as the start of the first window where we detect it.
We created a folded pulse profile from our list of potential giant pulses using polyco files generated with TEMPO2 (Hobbs & Edwards 2012). The polyco file contains a polynomial model of the pulsar phase as a function of time at the geocenter for our central observing frequency and observation window. From the folded pulse profile, we determined the MP and IP phase windows and any potential giant pulses that did not fall within these phase windows were discarded. As we had removed most of the RFI from our data, there were only about ∼ 10 false detections per observation; visual inspection showed these were clearly RFI. Given the narrowness of the pulse windows, our sample should thus contain no spurious pulses.
In Table 1, we list the resulting number of MP and IP detections for each observation in our coherent data, as well as their rate of occurrence. One sees that the rates vary significantly between the observations, something which has been seen before (Bera & Chengalur 2019), but is not understood.

Fringe Solution
To determine the delay and phase models for coherent combination, we used giant pulses with a signal-to-noise S/N > 50 (as measured on the incoherently summed data). We chose Effelsberg to be the reference location and time standard because of its relatively clean signal.
In the first processing stage, after correcting for the geometric delays and clock offsets (see 4.1), there will be further delays due to instrumental effects and ionospheric variations. To measure these, we correlated giant pulses observed at each individual telescope with the reference telescope in voltage, A third-degree polynomial fit of the amplitude is represented by the solid red line. Note that the amplitude is very close to unity, as expected given that the dominant source of noise for each telescope comes from the Crab Nebula. Middle: Relative fringe phases determined from performing an eigenvalue decomposition on matrices of visibilities of the individual giant pulses. The red line shows our fit. Bottom: Residuals of our phase fit, with the root-mean-square scatter indicated by the pink-shaded range. and fit the resulting offsets using a polynomial across each observation, weighting the results for each pulse by its S/N.
Since each polarization and sub-band may be affected by instrumentation and the atmosphere differently, we modelled them separately. We achieved fringe delay solutions for each baseline to Effelsberg with average root-mean-square deviation of 5 ns, i.e., better than 8% of a 62.5 ns time sample. An example is shown in Figure 1.
In the second processing stage, after correcting for the fringe delays, we determined and modelled the fringe amplitudes and phases. To solve for the time-varying complex telescope gains in each polarization and sub-band, we first crosscorrelate the voltage series of giant pulses between pairs of telescopes, integrating over time, to form complex visibilities. A matrix of complex visibilities was then created for each giant pulse, where V i,j denotes the visibility of telescope pairs and the telescopes are numbered from 1 to n. Note that autocorrelation data were not included (the diagonals are set to 0) and a minimum of 4 telescopes were used to create the visibility matrix in order to satisfy both closure phase and closure amplitude (Thompson et al. 2017), ensuring that any residual telescope based errors cancelled out. We then performed an eigenvalue decomposition on these visibility matrices. The eigenvectors corresponding to the dominant eigenmode give the relative complex gain between telescopes up to a complex constant. The gain amplitudes were modelled by fitting a polynomial across the whole observation and the complex gains were modeled by a sum of sinusoids which we then convert to phases. We show an example of our gain calibration in Figure 2. The average root-mean-square deviation for our phase model is 0.11 rad (implying S/N 9 per sub-band, polarization and telescope, consistent with what is expected given S/N 50 for the incoherently summed signal over 8 telescopes, 8 sub-bands and 2 polarizations, taking into account that this is measured over 16 µs, while most signal is within the scattering time of ∼ 5 µs).

Combined Data
Our pipeline achieves the expected coherence, as can be seen from the example giant pulse shown in Figure 3: the voltage series from multiple telescopes align well after delay and phase corrections, and the coherently summed giant pulse profile has S/N higher by the number of telescopes than the profiles of the single dishes, much better than the incoherently summed giant pulse profile where the S/N increase only by the square root of the number of telescopes. In the single-dish data, one sees that the intensities are very similar in units of the off-pulse noise, independent of telescope aperture. This reflects that the total system noise is roughly the same, as it is dominated by the radio emission from the Crab Nebula itself, which has a flux density of S CN ≈ 833 Jy at our observing frequency (Bietenholz et al. 1997), while the nominal system equivalent flux (SEFD) for a telescope ranges between S tel ≈ 19 and 450 Jy 5 , with an average S tel 300 Jy. Thus, the SEFD for our tied-array beam, which largely resolves out the nebula, can be estimated as (S CN + S tel )/N ≈ 140−160 Jy (depending on the number of telescopes N used).  Bottom: Pulse profile of this pulse (which happens to be a MP) as detected at each telescope (coloured lines), and after summing incoherently (dotted black line) and coherently (solid black line). For these profiles, the intensities in all sub-bands and both polarizations were summed, and were divided by the background, so that the profile is in S/N units. The time resolution is 62.5 ns. EK036 A 5.08 ± 0.15 0.76 ± 0.05 3.75 ± 0.14 0.43 ± 0.07 EK036 B 8.90 ± 0.19 0.63 ± 0.08 5.8 ± 0.3 0.34 ± 0.05 EK036 C 6.17 ± 0.16 1.00 ± 0.10 5.1 ± 0.2 0.62 ± 0.15 EK036 D 8.45 ± 0.14 2.31 ± 0.06 6.8 ± 0.2 1.94 ± 0.18 NOTE-The scattering timescales, τ , and intrinsic widths, p, of stacked MP and IP giant pulse profiles are presented here. The uncertainties from bootstrapping are likely underestimated, mostly because of the presence of echoes in the scattering tail (see Section 5).
In Figure 4, we highlight an IP and some other interesting giant pulses that we found. Two pulses with distinct multiburst structure are shown in Figure 5; we will return to the properties of pulses with multiple bursts in Section 6.3.

SCATTERING TIMESCALE
Emission from the Crab Pulsar passes through two scattering screens, one originating at the optical filaments in the  Crab Nebula and one in the ISM. The temporal broadening is usually dominated by the nebular screen (Vandenberg 1976), and one sees its effect in Figures 3-5 in the roughly exponential tail shared by all pulses, with a timescale of ∼ 5 µs. One also sees that the pulses have internal structure on a similar timescale.
To measure the timescale more precisely -and thus help determine the resolution at the pulsar (Section 2) -we construct average pulse profiles for the MP and IP components separately, by aligning and stacking giant pulses with S/N above 100 (removing profiles with obvious multiple bursts like those shown in Figure 5). For the alignment, we fit the pulse intensities with an exponentially modified Gaussian of the form, where A, t 0 , and p are the amplitude, centroid, and standard deviation of the Gaussian, respectively, τ is the scattering timescale, and C is the background intensity. We then create our stacked pulse profiles by adding MPs and IPs aligned using the individual t 0 , and fit the same exponentially modified Gaussian profile above to the stacks. Our pulse stacks and fits are shown in Figure 6 and relevant fit parameters are given in Table 2. One sees that the fits are not formally acceptable: they do not capture the relatively slow rise in intensity, suggesting that a Gaussian is not a good model for the average giant pulse profile, and the scattering tail shows bumps inconsistent with smooth exponential decay. Hence, the formal errors of the fits have little meaning, and we instead use bootstrapping, determining the standard deviations in the fit parameters from 8192 sets of pulse profiles created using random selection with replacement from our input profiles.
Note that the bad fits are not surprising: the individual pulses are made up of nanoshots and their average distribu- 10 20 Time (µs) Figure 6. Stacked MP (blue line) and IP (red line) components normalized to the peak of the exponentially modified Gaussian fit (dark blue and dark red lines, respectively) and centred so that t0 is at 0 µs. Fit values are given in Table 2. We also show the residuals between the best-fit model and the data for each of our observations. A: Stacked pulses from EK036 A. B: Stacked pulses from EK036 B. It is clear that the same echo features, which are indicated by arrows, (around ∼ 6, 12 and 18 µs) appear in both MP and IP components. C: Stacked pulses from EK036 C. D: Stacked pulses from EK036 D showing an echo feature around ∼ 8 µs as well as possibly one near a few µs that causes the profile to be substantially broader than is the case in the other epochs. tion does not have to resemble a Gaussian, and the scattering tail can have structure arising from individual scattering filaments. Indeed, the scattering tail for EK036 B shows clear evidence for echoes (see Figure 6B), with both the MP and IP profile showing similar structure in the scattering tail. These echo features appear on all high S/N giant pulse profiles in EK036 B (at lower S/N, it is difficult to make out echoes by eye). Echoes have previously been seen at lower frequencies, from 327 to 1400 MHz (Backer et al. 2000;Lyne et al. 2001;Crossley et al. 2007;Driessen et al. 2019), so the present study extends it to 1660 MHz. Echoes were not seen at 4.9 GHz during a period when echoes were present at 1.4 GHz (Crossley et al. 2007), but this may just reflect that echoes are frequency dependent, as expected if they arise in refraction in nebular structures.
In EK036 A and EK036 C, the modified exponential fits the scattering tail quite well with no apparent echoes. In EK036 D, however, there are deviations from an exponential scattering tail. Furthermore, the inferred intrinsic profile is much broader, with a larger σ. Since it seems unlikely the actual intrinsic average giant pulse profile changed, this may reflects an echo as well, one that is at a delay of only a few microseconds.
For all our observations, we see that the sharp rise of the model does not fit the intrinsic pulse structure of giant pulses well, especially for the MP profiles. Comparing MP and IP profiles, one sees that the former are systematically broader, with larger fitted σ and τ . The larger σ suggests MPs have intrinsically longer durations over which nanoshots are emitted than IPs. The longer scattering times likely also reflect some intrinsic difference in emission, e.g., that the nanoshots in MPs have a more skewed intensity distribution, falling off more slowly at the tail end, and that this skewed distribution leads to fits with τ biased high.
From the above, we conclude that among our measurements of the scattering time, those from the IP are probably more reliable, although even for those the errors are likely underestimated. For the purpose of a rough estimate, however, taking the scattering time in all epochs to be τ 5 µs should be good to about 50%.

CORRELATIONS
For regular pulsars, auto-correlations of dynamic spectra can be used to infer the scintillation bandwidth and timescale, and to learn about possible spatial offsets between pulse emission regions. For the Crab Pulsar's randomly occurring giant pulses, we follow Cordes et al. (2004) and Paper I and construct the normalized time and frequency correlation by first correlating pairs of giant pulse power spectra and then binning the correlations by the time separation between the giant pulses.
The correlation coefficient ρ(P 1 , P 2 ) between two power spectra P 1 and P 2 sampled at k frequencies can be estimated with, where m and s are estimates of the average and standard deviation of the power spectra, respectively. The second term accounts for noise biases (see Paper I and Appendix A), using that the mean and standard deviation of the noise power spectra are 1 as we have normalized our data by the background noise.
NOTE-The scintillation timescale tscint and frequency decorrelation νdecorr are defined as the values where the correlation function drops to 1/e and 1/2, respectively. For EK036 C, there were too few MP-IP pairs to derive a meaningful correlation. For Cordes et al. (2004), the values are approximate, as they were interpolated between their 1.48 and 2.33 GHz observations (assuming νdecorr ∝ ν 4 and tscint ∝ ν).
To construct our correlations, we create power spectra of giant pulses, discarding the 15% of each sub-band near each edge where little signal is detected because the passband rolls off. We then correlate spectra in each sub-band separately over frequency, taking into account that for larger frequency offsets fewer points contribute, giving individual estimates of the correlation as a function of frequency offset ∆ν and time offset ∆t. Next, we construct the normalized 2-D correlations by binning all pulse pairs by time separation, using a bin width of 1 s. We average pairs across sub-bands and polarizations using optimal weights, appropriate for our case where the correlations are low (see Appendix A, Equation A10). Below, in Section 6.1, we discuss the correlations between MP-MP and MP-IP pairs of pulses and in Section 6.2 we take a look at differences between the correlations for the first and second halves of pulses. Next, in Section 6.3, we focus on the correlations between multiple bursts within a given pulse, and in Section 6.4, we look at correlations between the left and right circular polarizations. We make qualitative comparisons with previous work and with expectations based on the scintillation screen's resolution, but leave the interpretation in terms of physical properties of the emission regions to Section 7. The binned correlations for all epochs are fitted with bivariate Gaussians, with as parameters the amplitude (A), the time (σ t ) and frequency (σ ν ) widths, the correlation between time and frequency (ρ f,t ), and, for the MP-IP correlations, possible time (∆t 0 ) and frequency (∆ν 0 ) offsets (which are relative to the MP, i.e., a positive sign of the time offset indicates that the IP trails the MP),

Correlations of Pulse Pairs in Time and Frequency
The results of the fits for all four epochs are collected in Table 3, where we also list values from Cordes et al. (2004) and Paper I. One sees that the biggest differences between our data and the ones from the literature are seen for the correlation amplitudes, which are very low for all our observations, The correlations are fitted with a two-dimensional Gaussian (as described in Section 6.1), with the white contour representing 1σ away from the peak (fitted values are listed in Table 3). The attached panels show the average correlation within the corresponding regions marked by dashed white lines and the red line shows the fit, also averaged over the dashed white lines. The red bars indicate the best-fit time and frequency offsets in the MP-IP correlation. Here, IP precedes MP, though the difference from zero is not significant. The complete figure set (4 images) is available in the online journal. at A 0.6% on average, which is about three times lower than found in Paper I, and a factor ∼ 50 lower than what was seen by Cordes et al. (2004). In contrast, the scintillation timescale of about t scint ≈ 10 s on average is similar to that measured in Paper I, but substantially shorter than the ∼ 27 s inferred by interpolating in the measurements of Cordes et al. (2004), while the typical decorrelation bandwidth we find, of ν decorr 0.4 MHz on average, is about half what was found in Paper I, yet similar to what one infers from Cordes et al. (2004). We will argue in Section 7 that these apparently contradictory trends can be understood from considering differences in effective resolution of the scintillation screens between the different observations. Like in Paper I, we find that MP-IP correlations are weaker and broader than the MP-MP correlations. For EK036 B, we also find a marginal time offset, of ∆t 0 = −0.4 ± 0.4 s, as well as a slightly more significant frequency offset, of ∆ν 0 = +0.05 ± 0.02 MHz. Both have opposite sign of those found in Paper I. The measurements for epochs A and D are less reliable, but the small frequency offsets are again inconsistent with those of Paper I.

Differences in Correlation between Pulse Halves
The high detection rate and large number of strong giant pulses in EK036 B allows us to look for differences in correlation with time of the pulse, and thus use signals that travelled to the scattering screen at increasing distance from the line of sight. To do this, we divide our giant pulses in half and create spectra from the first 8 µs and second 8 µs, resulting in 125 kHz channels. We then correlate spectra of the first half (corresponding to the inner regions of the scattering screen) and second half (corresponding to the outer regions of the scattering screen) separately and bin in time.
The resulting averaged correlations are shown in Figure 8 and the fit parameters listed in Table 3. From those, most striking is the change in frequency decorrelation width, which nearly halves. This is not unexpected, as the interference pattern from the outer scattering regions probed at later time will have finer structure. The decrease in amplitude by about 10% between first and second half also seems qualitatively consistent with an increase in resolution, as it means two separate pulses are less likely to be in the same resolution element. Clearly, however, even in the first half, the correlation amplitude of ∼ 1% is still much less than the 1/3 expected if giant pulses pass through the same part of the scattering screen. In contrast to the frequency and amplitude, the time decorrelation remains the same. That the time scale does not decrease suggests again that it is determined not by the size of the effective resolution element produced by the scintillation screen but rather by the size of the region in which giant pulses occur.
In an attempt to better constrain the MP-IP correlation offset in time, we also correlated the first halves of MP-IP pairs. We found that while these yielded a slightly higher amplitude than what we inferred from the correlation of spectra from the full, 16 mus windows on the giant pulses, the uncertainties on the offsets were not smaller (see Table 3).
We also attempted to correlate IP-IP pairs for EK036 B, but as there are few giant pulse pairs with sufficiently high S/N, we do not trust the fit parameters and we do not list them in Table 3.  Figure 9. Correlation of microburst pairs within each pulse component, with MP pairs (387 pairs) in blue and IP pairs (8 pairs) in orange. We omitted any correlations for which the uncertainty exceeded 0.15. The opacity of the individual points scales with the uncertainty of the correlation. The level of correlation is as low as seen for nearby pairs in different rotations, as expected given that the pairs likely not causally related but rather due to chance coincidence.

Correlation of Nearby Microbursts
Sometimes, multiple microbursts occur within one rotation. So far, we have excluded these from our analysis as we wished to avoid the risk of contamination of the power spectra, but we can check whether they correlate differently. We focus on EK036 B, where we found 1302 MP and 163 IP cases where we detected multiple giant pulses from the same rotation (with S/N greater than 8 and separated from each other by at least 16 µs): 1185 MP and 151 IP with 2 microbursts, 105 MP and 12 IP with 3 microbursts, 11 MP with 4 microbursts and 1 MP with 5 microbursts.
The expectation for the correlation depends on whether the multiple bursts in a rotation are causally related, i.e., consequences of the same physical events, or whether they are due to multiple observable giant pulses happening to be emitted in a single rotation. To estimate the fraction of multiples expected by chance, we can use that for EK036 B we detect a MP or IP giant pulse roughly every 12 and 48 rotations, respectively. Chance coincidences of N pulses should happen roughly at those rates raised to the power N . Given that we observed for about 177,404 periods, we would thus expect to see about 1299 MP and 78 IP pairs, about 111 MP and 2 IP triples, about 10 MP and 0.03 IP quadruples and about 1 MP and 0.0007 IP quintuples. For the MP, the observed numbers are roughly consistent with the expected ones, but for the IP they are significantly higher, suggesting some causal connection between multiple IP bursts. If most pairs are not causally related, we expect the correlation between them to be similar to what is seen between giant pulses that occurred nearby in time. As can be seen in Figure 9, this is indeed the case for the MP: the average weighted correlation amplitude of microbursts within the MP component is 1.8 ± 0.4%, a bit higher than but still consistent with the MP-MP correlations in Section 6.1. The average weighted correlation amplitude of microbursts within the IP component is nominally higher, at 6 ± 3%, perhaps again suggesting some causal connection between IP microbursts.
For both MP and IP, the correlation amplitude is much less than what was found by Karuppusamy et al. (2010), consistent with the suggestion that they observed at a time when the scattering time was shorter, the screen's resolution poorer, and the scintillation pattern for different bursts thus more similar.

Frequency Correlation Between Polarizations
At high frequencies, the MP giant pulses are observed to consist generally of multiple nanoshots, which have different center frequencies, bandwidths, and polarization ). If the same held at lower frequencies and if all nanoshots passed through the same scattering screen, one would expect the power spectra of giant pulses to have a variance three times larger than the mean squared. As mentioned in Section 1, this is not the case for the observations of Paper I and it is not the case for ours either: the variance of the power spectra equals their mean squared (see . This suggests that if a burst is indeed composed of multiple nanoshots, the shots are imprinted with different IRFs.
Correlating the two polarizations provides a way to verify this conclusion independently of the number of nanoshots, only relying on the fact that scintillation should not depend on the polarization. This is because the number of nanoshots influences both the degree of polarization of the giant pulse as a whole and the expected strength of the correlation, with the latter depending on whether the nanoshots are imprinted with the same or with different IRFs.
One can see that this would be the case by first considering a single nanoshot: since individual nanoshots are usually very highly polarized, one expects a high degree of polarization. Furthermore, this single signal is unlikely to be resolved by the screen, so its two polarizations should correlate perfectly. As one adds more nanoshots, all highly polarized but in roughly random directions, the polarization will start to average out, down to zero for a large number of nanoshots. The correlation between polarizations will also decrease, but differently depending on whether the nanoshots are all imparted with the same IRF, in which case the limiting correlation coefficient would be 1/3, or whether they have different IRF, in which case the correlation would tend to zero.
In order to look for this, we calculated Stokes parameters and correlation coefficients between left and right polarization for all pulses with S/N > 100. We follow the PSR/IEEE convention ) and the Stokes parameters are defined by where E − and E + are the left (−) and right (+) polarized electric field, {. . . } and {. . . } indicate the real part of the complex values and the angular brackets · · · denote averages over a pulse. In terms of the Stokes parameters, the total, linear, and circular degree of polarization of a pulse are given by, respectively, where we subtract 2 from the intensities to correct for the contribution from background noise (as appropriate given that we normalized our data by the intensity of the noise level in each sub-band and polarization; see Section 4.1).
In Figure 10, we show the inferred polarization properties. One sees that the degree of polarization varies over the full range, but with most pulses being only modestly polarized, with a root-mean-squared value of d 2 1/2 0.45. If individual nanoshots are fully but randomly polarized, one would expect that d 2 1/n s , where n s is the number of nanoshots, so we infer that, typically, n s 5 nanoshots contribute (or a few more taking into account that the nanoshots will likely have a range in brightness; see Appendix B). This is of the same order seen by Hankins et al. (2016) at higher frequencies, where nanoshots can be resolved.
One also sees that most pulses have little circular polarization, much less than the d V /d = 1/2 expected for (an average over) randomly polarized impulses. This is especially true at high degree of polarization, where a pulse is more likely to be dominated by just a few nanoshots. It is consistent with what is observed at high frequency, where individual nanoshots have been observed to be strongly linearly polarized (Jessner et al. 2010;Hankins et al. 2016).
Turning now to the correlation between left and right circular polarization, we first compare in Figure 11 pulse profiles and power spectra in left and right polarization for one highly polarized and one weakly polarized pulse. One sees that, as expected, for the highly polarized pulse, left and right are very similar, while for the weakly polarized one the power spectra in particular are quite different.
For the case that each nanoshot is imprinted with a different IRF, it turns out that it is possible to write the correlation coefficient directly in terms of Stokes parameters (see Ap- where we once again subtract 2 from the intensities to correct for the contribution from background noise. Hence, the expected dependence is on something close to the degree of lin-  Figure 13. Average correlation between left and right polarization of MP (blue) and IP (orange) giant pulses as a function of frequency, with fits using the sum of two Gaussians overdrawn. The narrow peak is due to signals from individual nanoshots shared between polarizations, and its width is as expected from the scattering time. The wide peak reflects correlation between the nanoshots; the small amplitude and large width show that these are resolved on the sky by the scattering screen. ear polarization, not the degree of polarization itself (though for our case, where the giant pulses are mostly linearly polarized, V 2 is small and one expects r(P − , P + ) d 2 ). For the case that all nanoshots are imprinted with the same IRF, no such direct relation is possible, but, as noted above, the expectation is that towards zero polarization, the average correlation will be 1/3 (see Appendix B.1).
We show the observed correlation coefficients as a function of (Q 2 + U 2 )/(I 2 − V 2 ) in Figure 12. We find that the correlation scales nearly perfectly as Equation 21, approaching zero for pulses with small polarization. This confirms that individual nanoshots are imprinted with different IRFs. As for pulse pairs, we construct weighted averages of the correlations between left and right polarization as a function of frequency, both for MP and IP. As one can see in Figure 13, these show a narrow, strong component and a wide, weak one, which can be reasonably well fit by a combination of two Gaussians (see Table 4). The narrow component is due to the part of the signal from the nanoshots that is shared between polarizations (i.e. the non-circularly polarized parts of individual nanoshots), with amplitude A narrow 15% and width ν narrow 0.037 MHz. This component's width is consistent with the scattering time (1/2πν narrow 4 µs). The wide component reflects the correlation between different nanoshots: its small amplitude A wide 2% and large width ν wide 0.32 MHz show again that the emission region is resolved. The fact that the amplitude is larger and the width smaller than what is found for pulse pairs, however, suggests that the region in which the nanoshots of a given pulse arise is less resolved and thus smaller than that spanned by pulse pairs.

EMISSION REGIONS
Our analysis suggests that the Crab Pulsar's emissions regions can be resolved by the scintillation screen. Since the amount and geometry of the scattering material changes over time scales of the order of months (Vandenberg 1976), the extent to which the emission regions are resolved will change as well, and thus different observations can give complementary information.
In this section, we aim to sketch out what the emission regions may look like as projected on the sky, trying to follow what a consistent picture would look like starting from the measurements of Cordes et al. (2004), and then including those of Paper I and ourselves. Below, we focus on the two emission regions from which the giant pulses that comprise the MP and IP originate. In addition, we can constrain the properties of the parts of the regions from which the nanoshots originate that make up a given giant pulse, since we found in Section 6.4 that these parts are also resolved in our observations. We show our sketch of the constraints on the emission regions in Figure 14, beginning with the results of Cordes et al. (2004) in panel A, and ending with what we infer from all observations in panel D; a side view of the latter is shown in Figure 15. In Figure 16 we show how the MP-MP correlation parameters allow us to place constraints on the emission region size from which MPs originate.
For our numerical estimates, we calculate resolutions from the scattering time τ using Equation 2, assuming d p − d s 1 pc (which is much smaller than d p 2 kpc, and thus the term d p /d s 1). We will ignore dependencies of the prefactor in Equation 2 on the degree of anisotropy of the scattering screen. We will also use the following relations, derived in Appendix C, that link the scattering time τ , resolution σ x , the size of the emission region σ s , the typical number n s of nanoshots per giant pulse, and the relative velocity between the pulsar and screen v rel with our observables, the amplitude A, decorrelation bandwidth ν decorr and scintillation timescale t scint : Here, all sizes are defined as standard deviations of normal distributions. The intermediate quantities σ c and n g are the effective spatial correlation scale and effective number of resolution elements, respectively; for unresolved emission, σ c σ x and n g 1, while for highly resolved emission, σ c 2σ s and n g (2σ s /σ x ) 2 . We write σ c, as a reminder that the scintillation time depends on the sizes parallel to the direction of motion.
For the pulsar velocity, we use the value inferred by Kaplan et al. (2008) from the proper motion (12.5 mas yr −1 at an angle of 290 • east of north): v p = 120(d p /2 kpc) km s −1 . For the screen velocity, we will assume that it has negligible average velocity with respect to the pulsar, but that it can vary, with velocities ranging up to ∼ 70 km s −1 like for the optical filaments (Trimble 1968).

Constraints from Unresolved Observations
From the average correlation amplitude of ∼ 1/3 of MP spectra at small time delays, Cordes et al. (2004) suggests that the MP emission region is one where giant pulses are comprised of multiple nanoshots which pass through the same scintillation element, i.e., the emission region is not resolved.
In the limit of unresolved emission, the scintillation timescale is determined by the speed with which the pulsar crosses a scintillation resolution element, t scint (σ x, / √ 2)/v p (see Equation 26), where σ x, is the resolution along the direction of motion. Cordes et al. (2004) measure t scint = 25 ± 5 and 35 ± 5 s at 1.48 and 2.33 GHz, respectively, which yields σ x, = 4240 and 5940 km, respectively. Since σ x scales linearly with wavelength (see Equation 2), the two results are consistent with each other.
Another estimate of the resolution can be made from the scattering time τ (Section 2). Cordes et al. (2004) do not measure τ directly, but for an unresolved source it can be inferred from the decorrelation bandwidth. At 1.48 GHz, they can only set a limit, ν decorr 0.8 MHz, which implies τ 1480 1/2πν decorr 0.2 µs, but at 2.33 GHz, they mea-sure ν decorr 2.3 MHz, which implies τ 2330 0.07 µs. 6 If this scattering were dominated by the nebula, one would infer, from Equation 2, σ x 730 km and σ x 790 km at 1.48 and 2.33 GHz, respectively (again consistent with each other given the dependence on wavelength). This is much smaller than inferred above from the scintillation time. This might be taken to suggest a very anisotropic screen, but more likely is the assumption that the nebula dominates scattering does not hold. As estimated in Section 2, the interstellar screen should have τ 0.2 µs at 1.66 GHz. Hence, the interstellar screen likely contributed to the scattering (though it cannot have dominated it completely, since that would imply a scintillation time of several minutes).
We conclude from the Cordes et al. (2004) observations that the emission region is small enough to not contribute to the scintillation time. From Equation 22, one sees that for given size σ s of the region in which giant pulses are emitted, the relevant ratio that determines whether it contributes is 2σ s /σ x (see Appendix C and Gwinn et al. (1998) for derivations). Taken this ratio to be smaller than unity and using the smaller of the two constraints on σ x, inferred from the scintillation timescale, we conclude that the emission region has a size σ s, 4240/2 = 2120 km.

The Emission Regions Resolved
The observations presented in Paper I and here were both taken at times that the scattering time was substantially longer than those of Cordes et al. (2004), and therefore dominated by the nebular screen. The scattering times of τ 1660 1 and 5 µs, respectively, imply resolutions of σ x = 290 and 130 km. In both observations, the maximum correlation amplitudes are greatly reduced compared to the 1/3 found by Cordes et al. (2004), to ∼ 1.8 and 0.6%, respectively, and the decorrelation bandwidths are substantially larger, at 1.1 and 0.4 MHz respectively, than expected from the scattering time (1/2πτ = 0.16 and 0.03 MHz, respectively).
Both these indicate the emission regions are resolved. If two giant pulses, even close in time, arise from a region much larger than the scintillation elements, they are imprinted with typically different IRFs and hence the correlation amplitudes will be low. Given that they will correlate best at the start  Cordes et al. (2004) suggests that at the time of their observations the resolution element covered the full MP emitting region. For simplicity, we assume that the emission region is circular. Since no information about the IP emitting region size/location is known but there is an IP emission regions somewhere, we show this region in grey. B: The clues from Paper I suggest that the MP and IP emission regions are resolved and that the size of the MP and IP region along the direction of the pulsar motion is larger than the physical resolution at the pulsar, ∆x. In addition, the slight positive time offset, ∆t0, indicates that the MP emission region precedes the IP emission region but that the regions likely overlap. C: In all our EVN observations, the physical resolution at the pulsar is smaller than in the previous observations. We also find that the size of the IP regions is larger than the MP region along the pulsar motion. D: In order to have sign changes in ∆t0 and ∆ν0, the emission regions cannot be oriented along the direction of the pulsar velocity. Rotating them, we arrive at our final picture of the emission regions. We use smaller circles within the larger MP and IP emission regions to indicate the regions where nanoshots of a giant pulses may arise. ∼ ∼ ∼ Figure 15. Sketch of the giant pulse emission regions as seen from the side (not drawn to scale). The black circle on the left indicates the pulsar. Blobs of highly relativistic material are ejected from around the light cylinder, at radius, RLC , and emit nanoshots at some specific locations along their trajectory, as indicated with the orange circles. Projected on the sky, these occur within a smaller part, of ∼ 1100 km, of the full ∼ 2400 km emission region. Along the line of sight, giant pulse emission occurs over ∼ 10 7 km. of the scattering tail, with the paths still close to the line of sight and the effective resolution poor, and poorly correlated at later times when the resolution is higher, only the coarser frequency structure in the power spectra will correlate and thus the decorrelation bandwidth will be relatively large (as derived for resolved sources by Gwinn et al. 1998). Following the same logic, the even lower amplitude of the MP-IP correlation suggests that the IP emission region is even more resolved and thus somewhat larger.
We also found that the regions from which nanoshots of given giant pulses arise are resolved. This is surprising, since the nanoshots are causally related and occur within a few microsecond. We return to this in Section 7.5, but here recall that if they were not resolved, i.e., if all nanoshots were imprinted with the same IRF, the variance of the power spectra should be three times larger than the mean squared. In contrast, in both our observations and those of Paper I, the variance is equal to the mean squared as expected if nanoshots are imprinted with different IRFs (see Appendix A). Furthemore, the correlations between the left and right circular polarization approach zero for pulses with low polarization (see Section 6.4), which can only be understood if the nanoshots are imprinted with different responses, i.e., if their physical separations as projected on the sky are resolved (see Appendix B).
We can use Equation 24 to estimate the number of scintillation resolution elements n g from the observed amplitudes and the number of nanoshots per giant pulse, n s 5 (see Section 6.4). We find that the emission region covers approximately 55 and 165 resolution elements in the data from Paper I and presented here, respectively.
As a consistency check, we can use the decorrelation bandwidth, which for low correlation amplitude depends on the amplitude and scattering time as ν decorr (1/ √ A)/2πτ (see Equation 25). Inserting numbers, we find ν decorr = 1.2 and 0.4 MHz, respectively, reasonably consistent with what we observe (1.1 and 0.4 MHz, respectively; see Table 3).
Following the same logic for the correlations between the two polarizations of individual pulses, for which we found an amplitude of ∼ 2% for the wide component for epoch B,  Tables 2  and 3). Overdrawn with a solid line is what is expected (Eqs 24-26 and Appendix C) for an emission region size σMP = 835 km, using a pulsar distance dp = 2 kpc and distance between the screen and the pulsar of dp−ds = 1 pc; the shading shows the effect of varying these within 1.4 ≤ dp ≤ 2.7 kpc and 0.5 ≤ dp − ds ≤ 2 pc. For comparison, the dotted and dashed blue lines show the expectations for emission sizes of σMP = 10 2.5 and 10 3.5 km. At long τ , the size of the resolution element σx ∝ 1/ √ τ is smaller than the emission region size and hence the scintillation time is constant, while the amplitude and decorrelation bandwidth decrease as 1/τ and 1/ √ τ , respectively. Conversely, at short τ , the resolution is poorer and the scintillation time is proportional to 1/ √ τ , while the amplitude approaches 1/3 and the decorrelation bandwidth approaches 1/2πτ . The model's assumption of negligible contribution from interstellar scattering breaks down at short τ : the solid gray region indicates where interstellar scattering likely dominates and the hatched grey one where it contributes significantly. A significant contribution from interstellar scattering means the model cannot reproduce the measurements of Cordes et al. (2004). we infer that the regions in which nanoshots from individual pulses arise cover about 50 resolution elements. The implied decorrelation bandwidth is ν decorr = 0.23 MHz, reasonably consistent with the observed width (∼ 0.32 MHz; see Table 4), although not as close as we found for the pulse pairs.

Emission Region Sizes
If the emission regions are larger than the resolution elements, the time decorrelation scale should reflect the time required for an emission region to move by its own size, independent of resolution. Indeed, we find that the time decorrelation scale is similar between the observations of Paper I and ours, at t scint, MP-MP 10 s. Thus, along the direction of motion, the size of the MP emitting region should be σ MP, v p t scint, MP-MP / √ 2 ≈ 850 km (see Equation 26). Note that here we ignore any contribution of the screen motion, which, at up to 70 km s −1 , introduces an uncertainty of about 30%. Indeed, variations in screen velocity likely are responsible for the fact that the individual scintillation times are not quite the same within the measurement errors (see Table 3, in particular epoch C).
Above, we found that the resolution elements had size σ x 290 and 130 km in the observations of Paper I and here, respectively, and that the MP emitting region spans about n g 55 to 165 resolution elements, respectively. Since in the limit of well-resolved emission, n g (2σ MP /σ x ) 2 (Equation 23), one infers σ MP 1 2 σ x √ n g = 1075 and 835 km, respectively. These numbers are pleasingly close to what we inferred from the scintillation time above, with the small differences perhaps due to changes in screen distance within the nebula, or deviations from isotropy. The low amplitude and increased time and frequency widths of the MP-IP correlations in both our datasets and that of Paper I suggests that the MP emission size is slightly smaller than the IP emission size. We can estimate the IP emission size along the direction of the pulsar velocity from the scintillation times for MP-MP and MP-IP, by σ IP, = v p (t 2 scint, MP-IP − 1 2 t 2 scint, MP-MP ) 1/2 1040 km. Unfortunately there are insufficient pairs of IP pulses to construct an IP-IP correlation meaningful enough for a complementary measurement, though for completeness we note that the very noisy result is certainly consistent with what is implied by the above size: n g (2σ IP, /σ x ) 2 255 and thus A IP-IP 1/n g 0.4% and ν decorr,IP-IP √ n g /2πτ 0.5 MHz As noted, the correlations between polarizations imply that the screen also resolves the region in which the nanoshots that comprise individual giant pulses originate. Using n g 50, we infer σ nano 460 km.
Overall, we infer that the MP and IP emission regions, if described by normal distributions, have sizes σ MP 835 km and σ IP 1040 km, or in terms of full width at half maximum (FWHM), about 2000 and 2400 km, respectively. The constituent nanoshots arise in somewhat smaller regions, σ nano 460 km (or FWHM of ∼ 1100 km). In Figures 14D  and 15 we show how we envisage the nanoshot emission regions to fit within the larger ones.

Orientation and Separation between Emission Regions
Time and frequency offsets in the MP-IP correlations can only occur if the respective emissions regions are offset as projected on the sky. Time offsets are easiest to interpret, as they reflect mostly how one region will move through the same resolution element at a slightly different time than the other. In both our observations and those of Paper I, marginally significant time offsets up to ∼ 1 s are detected. Treating those as upper limits suggests that the MP and IP emission regions are coincident along the pulsar motion to within a few 100 km. If we instead take them as measurements, we need to consider that the sign is different. This is possible only if the two regions are not aligned along the pulsar motion and the scattering screen is anisotropic on both occasions, but with a different orientation (see Figure 14D). In that case, for a given angle ψ between the short axis of the resolution element (long axis of the scattering screen) and the proper motion, the projected separation would be ∼ v p ∆t 0 / cos ψ, i.e., ∼ 170 km for a typical angle ψ = 45 • .
Independent evidence for anisotropy comes from the measurement of significant time-frequency covariance ρ f,t in the correlations. Furthermore, frequency offsets can only happen for an anisotropic screen. For those, one requires knowledge of the screen to infer a spatial offset, but a lower limit is given by (w/2)∆ν 0 /ν decorr (where w is the full width at half maximum of the resolution element), which implies offsets of 80 and 20 km for Paper I and epoch B, respectively.
Overall, we conclude that the MP and IP emission regions overlap substantially, but are likely offset by a few 100 km, i.e., about 10% of their diameter, in a direction different from the direction of the proper motion.

Evidence for Superluminal Velocity
From our observations, it seems clear the individual nanoshots that make up a giant pulse are resolved by the scattering screen, arising in a region on the sky with width σ nano 460 km. Yet, the nanoshots in a given giant pulse are clearly causally related, and arrive within σ t 0.6 µs ( Table 2), which at face value suggests a physical separation of 0.18 km.
As mentioned in Section 1, a similar problem was encountered by Bij et al. (2021), who found that in order to reproduce the drifting bands seen during the scattering tail of one particular giant pulse, the screen had to resolve the constituent nanoshots, which therefore had to be separated by ∼ 60 km. Bij et al. suggest that the discrepancy could be resolved if the nanoshots are emitted by blobs of material moving highly relativistically, with γ 10 4 . Highly relativistic motion would also naturally explain why the nanoshots are resolved in our observation; they require γ 3 × 10 3 .
One implication of relativistic motion is that the nanoshots are emitted over a region that is extended along the line of sight by another factor γ, i.e., several Gm, or of order 10 3 light cylinder radii. Given that, our measured sizes would be upper limits to the sizes of the regions where the plasma causing the emission originates, and the differences in apparent size between the interpulse and main pulse that we find may reflect differences in γ rather than true size.
A better measure of the true size may be the observed pulse phase widths of ∼ 7.6 • and ∼ 9.4 • for MP and IP, respectively. These are much larger than the beaming angle implied by the relativistic motion, of ∼ 1/γ 1 , suggesting that the blobs are emitted in a small range of directions, and thus also from a range of positions. Assuming the source is near the light cylinder radius, the implied widths perpendicular to the spin axis (and also roughly perpendicular to the direction of the pulsar velocity) are ∼ 210 and 260 km, respectively.

CONCLUSIONS AND RAMIFICATIONS
We find that in our observations, when the scattering was relatively strong and dominated by the nebular screen, the physical resolution at the pulsar was ∼ 130 km. From correlations between spectra of our large numbers of giant pulses, it is clear that this resolves the giant pulse emission regions.
We infer apparent diameters of ∼ 2000 and 2400 km for the emission regions of the MP and IP components. This strongly favors emission arising beyond the light cylinder radius. Slight time and frequency offsets in the MP-IP correlations suggest that MP and IP emission regions overlap significantly but not completely, with changes in sign suggesting that they are not aligned along the direction of the pulsar motion.
The largest surprise is that we also resolve the parts of the emission region from which the ∼ 5 nanoshots that comprise a given giant pulse arise: this is clear both from the statistics of the giant pulse power spectra as well as from the dependence of the correlation between polarizations on polarization properties. From the frequency dependence of the polarization correlations, we infer that the nanoshots arise in a region with diameter of ∼ 1100 km, smaller than but of the same order as the size of the region in which giant pulses occur. Since nanoshots are causally related, the simplest solution seems to be that the plasma emitting them move at highly relativistic speeds with γ 10 4 , generalizing to all pulses what was found for a single pulse from drifting drifting bands in its scattering tail Bij et al. (2021). It thus provides additional support for models that require highly relativistic motion (e.g., Lyutikov 2021).

Implications for FRB-substructure
Several lines of evidence suggest that at least some repeating fast radio bursts (FRBs) are generated by young magnetars (see e.g. Petroff et al. 2022, for a recent review). Observational support for such a scenario is given by, e.g., the FRB-like signals from the Galactic magnetar SGR 1935+2154 (CHIME/FRB Collaboration et al. 2020;Bochenek et al. 2020) and the extreme magneto-ionic environment of FRB 20121102A (Michilli et al. 2018;Spitler et al. 2016).
Additionally, substructure on microsecond and nanosecond scales has been reported for FRB 20180916B (Nimmo et al. 2021) and FRB 20200120E (Nimmo et al. 2022), respectively. Such short timescale variation translates to light travel times of order 1 km and below, i.e., similar in extent to what one would expect given the time separation between nanoshots making up an individual giant pulse of the Crab Pulsar (Table 2, Section 7.5). The evidence for superluminal motion discussed above may also play a role in FRB emission, suggesting that relativistic effects cannot be neglected when modeling FRBs. In fact, relativistic beaming would help explain the high brightness temperatures of FRBs, and may even be directly responsible for the emission, if the plasma blobs can act as a magnetic mirror (Yalinewich & Pen 2022).

Future work
To improve our understanding of the orientation and sizes of the giant pulse emission regions, further observations may help. In particular, probing the emission regions at different scattering timescales will further constrain the size of the emission regions and place bounds on the transverse velocities of the optical filaments. Similarly, observations at higher frequencies will result in a lower resolution at the pulsar which can also help to quantify the size of the emitting regions. With additional MP-IP correlations, we can more confidently determine whether there is a spatial offset between the emission regions and also map changes in the shapes of the resolution elements.
Our beamformed data shows the greatly improved sensitivity gained using multiple telescopes. With even more telescopes, we can perhaps detect more IPs and obtain more meaningful IP-IP correlations, which will constrain the size of the IP region directly. As already mentioned in Paper I, some of the uncertainty in our determination for the resolution element size comes from not having well constrained distances to the Crab Pulsar and the optical filaments; a good parallax distance would resolve this.
Facility: The European VLBI Network (EVN) is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project codes: EK036 A-D.

APPENDIX A. OPTIMAL ESTIMATES OF INTRINSIC CORRELATIONS FROM POWER SPECTRA
For two observed power spectra P 1 (ν) and P 2 (ν) sampled at k frequencies ν i , the sample correlation coefficient is given by, where m and s are the sample mean and standard deviation 7 of the power spectra, given by m = 1 k i P i and s 2 = 1 k−1 i (P i − m) 2 , respectively.
We wish to determine the intrinsic degree of correlation, ρ, between pairs of power spectra from multiple pairs of measurements, correcting the individual estimates for biases due to noise, and combining them using optimal weights, i.e., for a set of pairs ij to determine an optimalr where in general one would expect that the weights w ij and correction factors c ij would depend primarily on the signal-to-noise ratios of the input power spectra. Paper I already addressed the correction factors, but as we need those for the optimal weights as well, we briefly repeat the logic here.
The power spectra we are considering are Fourier transforms of the measured electric field of a giant pulse, given by where S(ν) and N (ν) are the pulse signal convolved with the IRF and the measurement noise, respectively. Since we are working with complex data, the measured power is, where | . . . | indicates absolute values and {. . . } indicates the real part of the complex values. We wish to determine the correction to the correlation coefficient to make it an estimate of the intrinsic correlation ρ(S 1 , S 2 ). As noted in Paper I, the presence of (uncorrelated) noise does not bias the sample covariance, which is the numerator in Equation A1, but it does bias the standard deviations in the denominator. Hence, ρ(S 1 , S 2 ) will be related to the expectation value ρ(P 1 , P 2 ) for the correlation coefficient by, Thus, to calculate the bias we need to relate the variances σ 2 P and σ 2 S . The expectation values µ P and σ 2 P of the sample mean and variance m P and s 2 P are given by, where we assumed that S(ν) and N (ν) are independent and used µ S , µ N , σ 2 S , and σ 2 N to represent their means and variances. Note that while the cross-term between signal and noise that originates from the squaring of the voltages does not contribute to the mean, it does contribute to the variance.
With the above, one sees that one could make an estimate of σ 2 S by s 2 P − s 2 N − 2(m P − m N )m N , but as noted in Paper I, this will lead to difficulties since s P is a noisy estimate of σ P (the estimate of m P is better, and those of µ N and σ N much better, since these are based on more data). A better estimate is possible using information on the distribution of S(ν) and N (ν). In particular, if both are normally distributed, with zero mean but different variances, their powers will be distributed like scaled χ 2 distributions with 2 degrees of freedom (real and imaginary parts). Hence, one will have µ 2 S = σ 2 S and µ 2 N = σ 2 N , which implies is not χ 2 distributed, as in the case where each giant pulse is comprised of many nanoshots with the same IRF, then σ 2 S = 3µ 2 S (see Section B). We show that S(ν) is indeed χ 2 distributed in Figure 17. This allows one to estimate c 12 more simply, ρ(S 1 , S 2 ) = ρ(P 1 , P 2 )c 12 = ρ(P 1 , P 2 ) µ P,1 µ P,1 − µ N,1 µ P,2 µ P,2 − µ N,2 .
To determine the weights w ij for an optimal estimate ofr opt , we need to know the uncertainties in the sample correlation coefficients. The general case depends on both ρ and the nature of the distributions, but for the case of no intrinsic correlation (i.e., ρ = 0), the variance of ρ is simply Var(ρ) = 1/(k − 1), i.e., independent of the signal-to-noise ratio. In our case, we expect this to hold roughly as well, since we have low ρ. Hence, for our estimate of the intrinsic correlation, the expected variance is, Given this variance, if we choose weights equal to its inverse, we will have an optimal estimater opt (even for higher ρ, this likely is reasonably optimal). Note that this implies that for high signal-to-noise pulses, with µ P µ N , the weight approaches a constant: all such pulses bring equal information, independent of exactly how high their signal-to-noise ratio is. But for low signal-to-noise pulses, with µ P 2µ N (i.e., signal-to-noise ratio per channel dropping below 1), the weight decreases as they bring little information.

B. CORRELATION BETWEEN CIRCULAR POLARIZATIONS AND STOKES PARAMETERS
If individual nanoshots are fully polarized, but with a random direction, then the ensemble should be less polarized. With a random walk like picture, one expects that while the intensity will be just be the sum of the intensities of the individual pulses, I = i I i , the polarized intensity will be I p ( i I 2 i ) 1/2 and thus the degree of polarization d = I p /I ( i (I i /I) 2 ) 1/2 . For n s equal-intensity nanoshots, one would thus expect d 1/ √ n s .
For multiple nanoshots, the voltage spectra are given by, where the sum is over the number of shots n s . The corresponding Stokes parameters (see Equations 14-17) following the PSR/IEEE convention  are where we defined I ±,i ≡ A 2 ±,i , and assumed everywhere that differences in angles φ i (ν) causes the cross terms between nanoshots to cancel. For the polarized flux, one finds, where the approximate equality uses the assumption that the polarizations of the individual shots are random, so that the crossterms cancel. Hence, we reproduce the expectation above, that the expected degree of polarization is d 2 i (I 2 i /I 2 ) 1/n s . In order to see the effect of scattering of the pulses in a screen, we should include an IRF g i (ν) in Equation B12 for the voltages. We will assume that all g i (ν) have the same total power, i.e., the same µ g = |g(ν)| 2 , and that the real and imaginary parts are normally distributed, so that |g(ν)| 2 is distributed as a χ 2 distribution with 2 degrees of freedom and hence σ 2 g = (|g(ν)| 2 − µ g ) 2 = µ 2 g . For the Stokes parameters, which are frequency averages of linear additions of the nanoshots, this will introduce a factor µ g in all equations, which we will absorb in the definition of the intensities, i.e., I ±,i = µ g A 2 ±,i . Below we discuss two cases, one where all nanoshots have the same IRF (Appendix B.1) and one were they all have different IRF (Appendix B.2).

B.1. Same Impulse Response Function
If all nanoshots are affected by the same IRF, then the left and right pulse power spectra are where we have abbreviated ∆φ ±,ij (ν) = (φ i (ν) − φ j (ν)) ± (ψ i − ψ j ). The means and variances of these power spectra are where in the approximate equality we used that I ±,i I ± /n s . For a large n s , one thus has σ 2 ± = 3µ 2 ± . The covariance between the two polarizations is given by . The polarization of each nanoshot is drawn randomly from the full Poincaré sphere, and its amplitude drawn from a power law with index α = −2.8, mimicking what we see in our data and others (Bera & Chengalur 2019); the frequency bandwidth, number of channels and scattering timescale are chosen to be as in our EK036 B dataset. We show 500 randomly selected giant pulses for ns = 2, 8, and 32 (as indicated). The red lines mark expectations for large samples of pulses, with the bold markers indicating where ns = 2, 8 and 32. One sees that the simulations confirm the expectations. The scatter is much larger for the case that the response functions are the same, because the expectation is only for the average over many pulses, in which polarization properties average out. In contrast, the expectation for different response functions is for individual giant pulses, only relying on averaging out random wave angles.
where we simplified using cos(∆φ − ) cos(∆φ + ) = 1 2 cos(2∆ψ), I ± = 1 2 (I ±V ), and µ 2 g A 2 −,i A 2 +,i = I −,i I +,i = 1 4 (I 2 i −V 2 i ) = 1 4 (Q 2 i + U 2 i ). For the approximate equality we used that for random polarizations, Q 2 + U 2 i (Q 2 i + U 2 i ). Combined with the product of the variances, σ − σ + (3 − 2/n s )I − I + = 1 4 (3 − 2/n s )(I 2 − V 2 ), the expected correlation coefficient is, For large n s , ρ 1/3, as noted by Cordes et al. (2004). Of course, we do not know n s a priori, but we can use the degree of polarization as an estimate. In Figure 18, we compare simulations with the prediction for a combination akin to the degree of linear polarization,d = (Q 2 + U 2 )/(I 2 − V 2 ), which we find below is expected to be a good approximation for the correlation coefficient for the case that nanoshots are imprinted with different IRF. For our case, given that Q 2 U 2 V 2 I 2 /3n s , one expectsd 2/(3n s − 1) and thus ρ (2 +d)/3(2 −d). One sees that the simulation confirms the expectation, albeit with fairly large scatter for individual simulated pulses. This is because the correlation coefficient for a given pulse depends on the polarizations of its constituent nanoshots; the approximations only hold for the average over a large number of pulses.

B.2. Different Impulse Response Functions
Turning now to the case that each nanoshot is imprinted with a different IRF, g i (ν), the left and right pulse power spectra are ±,i |g i (ν)| 2 + 2 i<j A ±,i A ±,j |g i (ν)||g j (ν)| cos(∆φ ±,ij (ν) + ∆φ g,ij (ν)), with means and variances Here, the variances are smaller because the cross-terms between the different nanoshots now cancel. Similar cancellation happens for the covariance between left and right polarizations, leading to, where used the same simplifications as in Equation B21 as well as that cos(∆φ − + ∆φ g ) cos(∆φ + + ∆φ g ) = 1 2 cos(2∆ψ). The product of the left and right standard deviations σ − σ + = 1 4 (I 2 − V 2 ), and thus the expected correlation coefficient is, where for the approximate equality we used that Q 2 U 2 V 2 I 2 /3n s . Note that unlike for the same-IRF case, for a given pulse with measured Stokes parameters, one does not have to rely on the last approximation -the result is known independent of any approximation for the individual properties of the nanoshots. Hence, comparing with simulations in Figure 18, one sees that even for individual pulses, the predictions are good; there is only some scatter due to the random wave angles, which causes imperfect cancellation of cross-terms.

C. EFFECTS OF RESOLVING THE EMISSION REGION
The correlation between giant pulses is low, much lower even than the ρ = 1/3 expected for giant pulses consisting of large numbers of nanoshots. Furthermore, the ratio of mean to variance of pulse spectra as well as the correlation between polarizations suggests that the different nanoshots generally are imprinted with different IRFs. Here, we derive expected correlation coefficients assuming each giant pulse consists of n s nanoshots, which arise randomly in some larger emission region, which is resolved by the screen. We will assume that the nanoshots are true impulses, without spectral structure, which randomly sample the whole emission region without any spatial correlation within a pulse.
For a first, simple derivation, we assume the emission region is resolved by the screen into n g different patches, each of which have completely independent IRFs. We write the electric field as in Equation B12, as a sum of impulses with amplitudes A i seen through possibly different IRFs g i (ν). Inspecting the results above for the power spectra for the case of nanoshots that all share the same IRF or all have different ones (Equations B18-B20 and B23-B25, respectively), one sees that the mean, µ P , is always the same, but the variance, σ 2 P , will depend on the number of pairs n corr that share their IRF (and thus correlate fully with each other), as, σ 2 P = µ 2 P 1 + 4n corr n 2 s = µ 2 P 1 + 2 (1 − 1/n s ) n g .
where we assumed that all nanoshots have the same amplitude A i = µ P /n s (for different amplitudes, one could rewrite using some effective number of pulses and pairs), and where in the second equality we used that for any given pair of nanoshots, the probability that they share the same IRF is simply the inverse of the number of patches, so that n corr = n s (n s − 1)/2n g . Similarly, for the covariance between two giant pulse spectra, one finds Cov(P 1 , P 2 ) = µ 1 µ 2 n corr n s,1 n s,2 = µ 1 µ 2 1 n g , where in the second equality we again used that the probability for a given pair to share the same IRF is 1/n g and thus that n corr = n s,1 n s,2 /n g .
Combining the above, one finds that the correlation coefficient for two separate pulses is expected to be, ρ(P 1 , P 2 ) = 1 n g + 2 (1 − 1/n s ) .
Thus, for two single-shot pulses, one has the expected ρ(P 1 , P 2 ) = 1/n g and for pulses with large n s , ρ(P 1 , P 2 ) 1/(n g + 2). The above derivation is a simplification in that it assumes that the screen resolves the region in a fixed number of patches, while in reality the resolution of the screen varies throughout the scattering tail, from very poor early on, when the paths along which radiation travels are all close to the line of sight, to much better later on. Indeed, in Section 6.1, we found that spectra taken in the first half of pulses correlated more strongly with each other than those taken in the second half.
In order to take this into account, we first start by rewriting the above in terms of a fraction of the IRF that is correlated. Next, to understand not just the amplitude of the correlations at zero lag, but also the correlation at other ∆ν and thus the decorrelation width ν decorr , we consider the behaviour of the IRF in the time domain, g(t).
To take into account that part of the IRF of nanoshot i is shared with that of other nanoshots, we write the IRF in terms of a shared (correlated) part g c (ν) that contributes a fraction f c of its power, and a different (uncorrelated) part g u,i (ν) that contributes the rest, i.e., g i (ν) = g c (ν) f c (ν) + g u,i (ν) 1 − f c (ν).
Note that we implicitly assume all IRFs are normalized the same way, i.e., all have the same ν |g(ν)| 2 dν. With this, where we usef c to indicate an appropriate average over frequency of f c (ν) (note that for relatively narrow-band observations like ours, the dependence on frequency can be safely ignored). Thus, identifying n g = 1/f 2 c and n s = µ 2 P / i µ 2 i one recovers Equation C28. Similarly, for the correlation coefficient (assuming both pulses have the same effective n s ), ρ (P 1 , P 2 ) = i,j A 2 1,i A 2 2,j Cov |g 1,i (ν)| 2 , |g 2,j (ν)| 2 σ 1 σ 2 =f 2 c 1 + 2f 2 c (1 − 1/n s ) .

(C33)
Turning now to the time domain, without much loss of generality we can take g(t) to be amplitude-modulated noise, i.e., g(t) = a(t)CN (0, 1), where CN (0, 1) is the complex normal distribution with mean zero and variance one, and a(t) is normalized to give unit total power, i.e., p(t)dt = 1 with p(t) = a 2 (t). For two impulses arising from different locations, the fraction f c (t) of the response that is shared (and thus correlates) will decrease with increasing time, i.e., the part phased coherently will have reduced power p c (t) = p(t)f c (t).
We can relate this to the frequency domain by noting that generalizing the cross-correlation coefficient in Equation C33 to a function of ∆ν, the covariance term becomes a cross-covariance, which requires the auto-correlation of |g c (ν)| 2 . Using the cross-correlation theorem, we can write the latter as, |g c (ν)| 2 |g c (ν)| 2 = F F −1 |g c (ν)| 2 2 = F |g c (t) g c (t)| 2 , where F indicates a Fourier transform and F −1 its inverse. Using that the expectation for the auto-correlation power between the part of impulses that is shared in the time domain is given by the autocorrelation of their power envelopes, the expectation value for the auto-correlation is, |g c (ν)| 2 |g c (ν)| 2 = F {p c (t) p c (t)} = |p c (∆ν)| 2 .
The precise shape of the autocorrelation of |g c (ν)| 2 will depend on p(t) and f c (t), but one can gain insight using a simple assumption, that both p(t) and f c (t) are exponentials. For p(t), this is reasonable, as the pulse profiles can be modelled fairly well as the convolution of a Gaussian arising from nanoshots with an exponential scattering tail (see Section 5). It is also the expected profile for an isotropic scattering screen in which the scattering points are distributed normally. Since the profile has to be normalized, one has p(t) = exp(−t/τ )/τ (for t ≥ 0) where τ is the scattering timescale.
For the shared fraction f c (t) averaged over many pulses, an exponential distribution may also be plausible: if the locations of the individual pulses are distributed roughly isotropically with a normal distribution in each direction, then squared distances between pulse pairs will follow an exponential distribution. Hence, while for close pairs f c (t) will extend further in time than for distant pairs, on average one may expect f c (t) exp(−t/τ c ), with τ c corresponding to the time delay for which the IRF of two typical pulses will start to differ significantly. Thus, the correlated power envelope will be, p c (t) = 1 τ exp(−t/τ ) exp(−t/τ c ) = 1 τ exp(−t/τ ), with τ = τ τ c /(τ + τ c ). The autocorrelation in time is (p c p c )(∆t) = (τ /2τ 2 ) exp(−|∆t|/τ ), and hence in the frequency domain, one expects a Lorentzian, |p c (∆ν)| 2 = τ τ 2 1 1 + (2π∆ντ ) 2 . (C38) Comparing with Equation C33, we can identifyf c = τ /τ . The half-width at half maximum is given by ν decorr = 1/2πτ , which is larger than the 1/2πτ one would infer from the measured scattering time τ by a factor 1/f c . For well-resolved emission, i.e., smallf c , the dependence on number of pulses drops out, and one expects a scaling with amplitude A of the decorrelation width as ν decorr (1/ √ A)/2πτ . As noted above, an exponential for f c (t) is expected if the nanoshots are distributed roughly isotropically on the sky, with a normal distribution. Assuming they have a variance σ 2 s in each direction, differences in position between pairs will have a variance 2σ 2 s in each direction. Hence, the distance squared r 2 between pairs will follow an exponential distribution exp(−r 2 /4σ 2 s ) and the typical distance is 2σ s . Thus, like Gwinn et al. (1998), we find that 2σ s /σ x is the relevant measure of the degree to which the emission region is resolved. We confirmed using simulations that the average amplitude for a pair of pulses equals (τ /τ ) 2 = (1 + (2σ s /σ x ) 2 ) −1 .
The amplitude is inversely related to the effective number of resolution elements (see Equation C30), which suggests defining an effective spatial correlation scale σ c = (σ 2 x +(2σ s ) 2 ) 1/2 , such that one has n g = (σ c /σ x ) 2 (and thus τ /τ =f c = σ x /σ c = n −1/2 g and ν decorr = √ n g /2πτ ). The correlation scale σ c can be related to the spatial separation scint between giant pulses at which the correlation will decrease by 1/e, which is useful for comparing with the scintillation time t scint = scint /v rel . Since amplitudes vary with separation as exp(− 2 /2σ 2 c ), the correlation, which is a fourth-order product, will fall off as exp(−2 2 /σ 2 c ). Hence, one infers scint = σ c / √ 2 (which we also confirmed with simulations).