On Single-pulse Energies of Some Bright Pulsars Observed at 1.7 GHz

We present Arecibo observations of the bright pulsars B0301+19, B0525+21, B0540+23, B0611+22, and B0823+26 at 1.7 GHz with 100 MHz bandwidth. No giant pulses were found, except for B0823+26, where we recorded a giant interpulse with 230 times the average peak intensity. The postcursor in B0823+26 shows a symmetric double-peaked structure, indicating that it is frequency dependent. In all pulsars, for a given single-pulse peak intensity there is a range of equivalent widths up to a maximum, which becomes smaller the stronger the pulses are, thereby apparently limiting the energy output. Forming average profiles from pulses with certain equivalent widths leads to profiles with changing component characteristics and could allow exploring the magnetosphere at different heights, assuming a dipolar field geometry. We found that the normalized lag-1 autocorrelation coefficient for single-pulse energies can be over 0.5, indicating high correlations. From the first peak of the energy autocorrelation function a so-far-unobserved 15-period modulation is found for B0540+23, as well as a possible 10-period modulation for B0611+22. We also show that a fit of the Weibull distribution to the cumulative probability for the energies yields a better fit than the usual lognormal distribution. The cumulative probability distributions permit an estimate of the nulling fraction, which ranges from 0.6% for B0611+22 to 24% for B0525+21.


Introduction
Each single pulse from a pulsar looks different in amplitude, width, and exact phase, and yet synchronous adding of a sufficient number of single pulses leads to a temporally stable profile (e.g., Rathnasree & Rankin 1995). The early years of pulsar research saw many studies on periodicities in amplitude fluctuations, which ultimately were termed P 3 and had typical ranges between 2 and 15 rotation periods P. Examples include a periodicity of about 4 P in the discovery pulsar B1919+21 (Conklin et al. 1968), 5.5 P in B0823+26 (Taylor et al. 1969), 10.0 P in B1929+10 (Lang 1969), and 11.2 P in B1133+16 (Hesse 1970). Amplitude fluctuations have also been reported for some millisecond pulsars (MSPs; e.g., Jenet et al. 1998;Edwards & Stappers 2003) and hence seem to be a universal phenomenon. Several large studies of single-pulse amplitude fluctuations at different frequencies have been published since then (e.g., Rankin 1986;Weltevrede et al. 2006and Basu et al. 2016. It has been proposed that the cause for these amplitude fluctuations are subpulses that appear from pulse to pulse at a slightly phase-shifted position and hence are drifting across the average profile with a timescale that has been termed P 2 (Drake & Craft 1968;Backer 1970a), representing the mean separation between the drifting subpulses. Widely accepted today is a model of spark discharges in a gap over the polar cap in the pulsar magnetosphere (Ruderman & Sutherland 1975), where the sparks have fluctuating intensity and are circulating around the magnetic pole, thereby giving rise to the phenomenon of drifting subpulses. However, drifting is not observed in all pulsars and has been tied to the shape of the average profile, which itself is viewed as a cut through the emission cone, which is often described with a core-cone model (Rankin 1986(Rankin , 1990(Rankin , 1993a(Rankin , 1993b. According to this model, radio emission can originate from a central region (the core) and a surrounding annular region (the cone), which can also comprise nested cones. The core emission is centered on the magnetic axis, and pulsars with a simple onecomponent profile are considered core-single S t -class pulsars. Depending on the line of sight, a single-component profile may also be a conal-single S d -class pulsar, where an emission cone is cut tangentially. If both sides of the cone appear along the line of sight, then it belongs to the conal-double D class; in the case of a central line-of-sight traverse, both the cone and the core may be seen, leading to the triple T class. Combining the circulating spark gap model with the conal emission structure, one is led to the conclusion that the drifting phenomenon should be strongly present in S d -class pulsars, only weakly present in D pulsars owing to the steep angle of the line of sight, and absent in pure S t pulsars.
Many observations indeed favor this scenario, but the situation is as yet not fully understood. Examples are B0823 +26, which has clear periodic amplitude fluctuations but no prominent subpulse drifting is seen (Taylor & Huguenin 1971;Basu & Mitra 2019), and B0540+23, which was found to have signs of sporadic drifting but without amplitude periodicities (Nowakowski 1991). Using phase-resolved fluctuation spectra, drifting and related intensity fluctuations have been further divided into phase and amplitude modulation subcategories (Weltevrede et al. 2006(Weltevrede et al. , 2007Basu et al. 2016).
Occasionally, there are no pulses observed at all, which has been dubbed nulling (Backer 1970b) and has since been observed in most pulsars. The nulling fraction (NF) can be anywhere between less than 0.1% and 60% or more (e.g., Ritchings 1976;Rankin 1986;Vivekanand 1995), with S t pulsars typically nulling the least (Biggs 1992). Some pulsars undergo mode switching, where the average radio emission manifests in two or more distinct metastable profiles as originally observed in B1237+25 (Backer 1970c). B0823+26 has been discovered to switch between a radio-bright (B) and a radio-quiet (Q) mode (Sobey et al. 2015), with the B-mode dominated by core emission and the Q-mode dominated by conal emission (Rankin et al. 2020).
The intensity of emission is modulated by at least two other timescales, one at microseconds to milliseconds and the other at minutes to hours. The former is known as microstructure and originates from emission fluctuations on top of single pulses that are on the order of milliseconds in normal pulsars and has been contemplated to originate from either radial (temporal) flux modulation or longitudinal sweeping beamlets (e.g., Cordes et al. 1990;Lange et al. 1998;Kramer et al. 2002;Mitra et al. 2015). Pulsar signals are traveling through the turbulent interstellar medium, and the flux can be additionally modulated (termed scintillation) by inhomogeneities that disturb the phase of the radio waves (e.g., Rickett 1970 and references therein). One distinguishes between the weak and the strong scintillation regime, with the latter comprising a shorter diffractive and a longer refractive timescale, both of which are frequency and distance dependent (Lorimer & Kramer 2005 and references therein).
Some pulsars emit occasionally single pulses that by far exceed the mean energy, the most famous of which are the Crab pulsar B0531+21 and the MSP B1937+21, where many detailed studies exist (e.g., Jenet et al. 2001;Majid et al. 2011;Popov & Stappers 2007;Zhuravlev et al. 2013;Bera & Chengalur 2019;McKee et al. 2019). Originally called "jumbo pulses" (Sutton et al. 1971), they are now referred to as giant pulses (GP) and are usually defined as single pulses with an energy over 10 times the mean single-pulse energy (e.g., Karuppusamy et al. 2012;McKee et al. 2019).
While the distribution of regular pulses can be reasonably well fit with a lognormal distribution, the GPs follow a power-law distribution (e.g., Cairns et al. 2001;Johnston & Romani 2002;Zhuravlev et al. 2013;McKee et al. 2019). However, another study found that regular pulses in many pulsars show deviations from a lognormal distribution (Mickaliger et al. 2018). GPs on average have a tendency to become narrower the stronger they get (Popov & Stappers 2007;Majid et al. 2011;Bera & Chengalur 2019). The same tendency for stronger pulses to become narrower has been also observed in a few MSPs for regular pulses (Jenet et al. 1998;Vivekanand 2000;Liu et al. 2016). If stronger pulses are narrower, then this would imply that the area under the pulse, which is proportional to the pulse energy, is limited by a certain maximum value for a given pulsar. There are some theoretical arguments for GPs (Petrova 2004) and some observational arguments (Burke-Spolaor et al. 2012) that energy-regulating mechanisms could be present in pulsars.
Hence, pulsar radio emission is fluctuating on several different timescales and intensity levels. Due to the lack of a unified emission theory, one way forward is trying to set constraints on timescales, flux density levels, energies, and widths of single pulses. Here we present observations at 1690 and 1730 MHz of five known bright pulsars with the Arecibo Observatory's legacy 305 m radio telescope. A narrow 100 MHz bandwidth allows any frequency-specific effects to remain preserved as much as possible. The pulsar B0540+23 has been observed twice, in 2016 and in 2019, and will serve as a test of the reproducibility of results and to potentially detect intrinsic changes in the pulsar magnetosphere. Section 2 describes observations and data reduction methods, Section 3 presents analyses and results, and Section 4 discusses the findings and provides conclusions.

Observations and Methods
All observations were undertaken in 2016 and 2019 with the Arecibo Observatory's legacy 305 m radio telescope that collapsed on 2020 December 1. Table 1 lists some basic pulsar properties for the five bright pulsars that were observed, the observation dates, and frontend and backend configurations of the telescope; all observations used the L-band wide receiver.
The choice of pulsars and frequency bands was dictated by the desire to use only raw data without appreciable radio frequency interference (RFI), which turned out to be always the highest available center frequency around 1.7 GHz. Pulsar B0540+23 was observed during both projects, with the observation from 2016 labeled (A) and the one from 2019 labeled (B). The reason to include the second B0540+23 observation is to verify the reproducibility of all processing, analysis, and fitting steps and possibly to detect small emission changes that might go unnoticed otherwise. The remainder of our observed pulsars and the lower frequency bands have some time-domain RFI or radar spikes, baseline fluctuations, and contaminated frequency channels, or were simply too weak for a single-pulse analysis. The selected raw data presented here are clean to such an extent that neither clipping nor masking or channel excision was necessary, and hence any potential data altering is avoided. Figure 1 shows representative examples of time series with 50 consecutive periods, confirming the already well-known fact that D-class pulsars typically have a much larger NF than the S t -class pulsars (e.g., Ritchings 1976;Rankin 1986;Biggs 1992;Vivekanand 1995).
In Table 2 we list data reduction parameters and some calculated average profile characteristics. The initial data reduction was performed with the software package PRESTO [ascl.net/1107.017]. To reduce the number of data samples, we downsampled the data by the factors indicated in Column (3), arriving at a number of samples per pulse (i.e., phase bins) between about 3000 and 6000. Downsampling simply adds a predefined number of consecutive data samples as opposed to a resampling and avoids any potential data altering, as has been previously reported (Karastergiou et al. 2001).
Since there is usually no integer number of data samples per exact pulse period and we utilized the search mode, a correction has to be applied before folding the time series. By maximizing the signal-to-noise ratio (S/N) of the folded profile, a correction factor was determined for adding or removing individual data samples. This is a standard procedure, and for our observations this amounts to adding or removing one data sample every few pulse periods (typically between 2 and 10). The last step involves converting machine counts into a Jansky scale, which was done with the radiometer equation. We have set the digitization loss factor equal to 1 and determined a conversion factor for each individual pulse with the respective standard deviation obtained from 1000 off-pulse data samples. The system parameters for the L-band wide receiver in 2016 were T sys = 30 K, G = 10 K Jy −1 , and in 2019 they were T sys = 25 K, G = 8 K Jy −1 . In between, on 2017   September 20, Hurricane Maria had made landfall in Puerto Rico as a Category 4 storm and broke off the 96-foot line feed. The PUPPI backend suffers from a slight phase mismatch of the time-interleaved ADC samplers, leading to weak, reflected "ghost images" of the pulsar signal (Lam et al. 2019). The amplitude of the "ghost image" becomes weaker toward higher frequencies, and we found it to be unnoticeable above about 1.4 GHz, which was another reason to choose a GPU node at the high-frequency end. The Mock backend, on the other hand, does not have this problem. Figure 2 shows the average profiles of the flux-calibrated pulsars with some basic characteristics indicated. The red horizontal line is drawn at the 3% level of the profile maximum and serves to derive the width, which is labeled W3 in the figure and Table 2. We confirm the core-single S t and conaldouble D general profile shape as established by the Rankin classification (Column (2) in Table 1) and a recent comprehensive study (Olszanski et al. 2019). Pulsar B0540 was also observed in 2019 (B) and has an identical shape to that during the 2016 (A) observation with a 10% reduction in the mean flux density, S mean , while the width W3 remained the same to within less than 2%. The width was determined by eye from the intersection of the indicated horizontal line with the profile data points and has a maximum error of 0.2% across all pulsars (B0540 yields the largest error owing to the smallest slope of the profile at the intersection). Emission duty cycles (W3) are between 2.8% and 9.7%, with B0823 having by far the shortest. The postcursor (PC) component of B0823 is trailing the main pulse (MP) by over 30°and hence is just outside our panel.

The Average Profiles
For B0823 we derive a mean flux density at 1730 MHz of 7.4 mJy, which compares well with the average of over 50 observations in a long-term study at 1700 MHz yielding 8.7 mJy (Daszuta et al. 2013). Figure 2 shows that pulsar B0540 has an asymmetric average profile, even though it belongs nominally to the core-single S t class. This pulsar has in all likelihood a second component in the profile at the trailing edge as was already pointed out in an earlier study at 430 MHz (Nowakowski 1991); another very recent multifrequency study also arrived at the conclusion that this pulsar has a more complicated emission pattern than originally assumed (Olszanski et al. 2019). On close inspection, the profile of the Dclass pulsar B0525 likely has an additional weak third component in the bridge between the two main conal components, which has also become apparent in a single-pulse multifrequency study at 327 MHz (Young & Rankin 2012).
The profile energy was determined by adding the fluxcalibrated data values between the 3% points and multiplying by the sampling time. The energy values are listed in Column (8) of Table 2, with the two D-class pulsars having the highest energy content. We also compare the full width at 10% maximum (W10) of the profiles from Figure 2 with the available W10 values from the literature in Table 3. The W10 width is seen to decrease with increasing frequency in the range 1.4-1.7 GHz, which indicates that there is an anticorrelation. Assuming a simple dipolar magnetic field geometry, this would indicate that higher frequencies are emitted at lower heights, a tendency that has been dubbed radius-to-frequency mapping (Cordes 1978 and references therein). However, it still has not been unambiguously demonstrated whether this is merely a height effect or other propagation effects in the magnetosphere play a role.

Single-pulse Peak Fluxes and Energies
For every single pulse in the time series we determined the phase location and flux density of its maximum and show this in Figure 3. The phases are labeled by data samples (phase bins), and the color bar gives the flux density on a logarithmic scale in units of mJy; the strongest single-pulse amplitude was recorded for B0823 with about 13.4 Jy. This plot, which we call the maxmap, gives a quick visual on the overall data quality and can also reveal several different emission features, such as the presence of RFI, nulling, or scintillation. If the data were contaminated by RFI, then there would appear scattered high-flux data points ("red") outside of the emission window. We did see several RFI-contaminated maxmaps (not shown) and decided not to include those pulsars in this discussion.
The maxmap shows that pulsar B0611 has the fewest nulls and mostly also weak emission. The interpulse (IP) pulsar B0823 has the strongest pulses and the narrowest emission window, and one can see around pulses 3000, 5000, and 7000 a significant increase in the flux density, which is due to interstellar scintillation. This pulsar has been shown to emit in two different modes, termed the bright (B) mode and the quiet (Q) mode (Sobey et al. 2015;Basu & Mitra 2019;Rankin et al. 2020). The Q-mode was found to persist for several hours at a time over the frequency interval 0.1-2.7 GHz and to be about 100× weaker than the B-mode. During our observing run, B0823 was emitting in its stronger B-mode. The two horizontal orange lines indicate the expected phase range for the IP, and it can be seen that some single pulses have their maximum value  in that range, but no specific pattern is obvious. Figure 4 shows the average profile and examples of a very bright IP and PC, which are stronger and narrower than the accompanying MP. In Figure 3 one can also clearly see the nulling sequences in B0301 as the vertical low-flux ("blue") lines, which is also corroborated by the corresponding missing pulses in the time series, hence not being RFI. We can visualize scintillations by plotting the single-pulse energies for the full time series in Figure 5. While the energy for the average profile was calculated by summing the data samples between the 3% points, for the single-pulse energies the phase range was extended symmetrically to either side of the 3% points to obtain a total phase range of two times the average profile width. With this criterion, care is taken of the widest pulses that are exceeding the width of the average profile as will be shown later. Radio emission from the pulsars has to pass through the inhomogeneous and turbulent interstellar medium, leading to scintillation, where the intensity reaching the observer is modulated in time and frequency. For our frequencies and dispersion measures (DMs), the strong scintillation regime is relevant, and we focus here on the diffractive branch with its scintillation bandwidth Δf DISS and timescale Δt DISS .
Usually, a Kolmogorov power spectrum for isotropic, incompressible turbulence is assumed, which leads to a frequency dependence according to Δf DISS ∼ f 4.4 and Δt DISS ∼ f 1.2 (Lorimer & Kramer 2005 and references therein). We obtained the scintillation bandwidth and timescale from the NE2001 Free Electron Density Model (Cordes & Lazio 2001) using the Python wrapper pyne2001 (pypi.org/project/ pyne2001/). Since both scintillation values are calculated at 1 GHz, they were scaled up to our observing frequency of 1.7 GHz with the Kolmogorov frequency dependencies and listed in Table 4.
The observing bandwidth is much wider than the scintillation bandwidth for B0525, B0540, and B0611, and hence no scintillation is expected. For B0301 the observing bandwidth is marginally smaller, but no scintillation on an 800-period timescale is visible. The refractive scintillation timescale for B0301 is Δt RISS = ( f/Δf DISS )Δt DISS ≈ 15,400 periods, which is much longer than the observing length of 3621 periods. However, it is suggestive from Figure 5 that some kind of energy modulation appears with a characteristic timescale on the order of 3000 periods (69 minutes), which then might be intrinsic to the pulsar. On the other hand, for B0823 one can clearly see three scintillation maxima in the single-pulse energies separated by about 2000 periods or 18 minutes, which is very close to the predicted 16.7 minutes based on the NE2001 model. A 2D autocorrelation analysis of the dynamic spectrum was published at coincidentally the same observing frequency of 1.7 GHz for B0823 (Daszuta et al. 2013). That study derives a diffractive scintillation timescale of 19.3 minutes, which is, within their reported errors, in good agreement with our visual estimate of 18 minutes.

Relations between Peak Fluxes, Energies, and Equivalent Widths
The single-pulse energies divided by the respective mean energy 〈E〉 are plotted as 50-bin histograms in Figure 6 and agree with the corresponding energies of the average profile as shown in Table 2. The histogram peak is listed in the main panels as "max," while the highest energy is listed as "high"; the mean energy is additionally included in units of mJy · s. We note that the two B0540 observations produce virtually an identical histogram, thereby proving the reproducibilty of our analysis and at the same time indicating that the pulsar B0540 did not noticeably change its emission pattern between 2016 and 2019. The insets show the off-pulse energy histograms that were obtained from the same number of data samples as the onpulse histograms and also divided by the on-pulse mean energy; the symmetry with respect to zero indicates, as expected, Gaussian noise.
Regarding the "max" and "high" energy values, the two Dclass pulsars (B0301 and B0525) do not show any particularly different behavior compared with the S t -class pulsars; however, they have a second sharp peak around zero mean energy resulting from a significant NF. This is a well-studied phenomenon and has first been reported as an additional peak in energy histograms for B1133+16 and B1237+25 (Backer 1970b) and since then in many other pulsars (e.g., Ritchings 1976;Vivekanand 1995;Rankin 2017;Basu & Mitra 2019). We are going to introduce a way to estimate the NF using cumulative probability distributions for single-pulse energies in a later section.
From Figure 6 we also see that the strongest pulses occur for B0823 with 9.5 〈E〉 and hence conclude that none of our observed pulsars were emitting GPs. As a comparison, a 9.5 〈E〉 pulse from B0823 has an energy of 38 mJy · s, while some of the most energetic Crab pulsar GPs detected at 1.4 GHz have about 4700 mJy · s (Bera & Chengalur 2019). It has been reported that B0823 emits single pulses at 0.15 GHz up to an energy of 25 〈E〉 (Sobey et al. 2015), which indicates that it can produce relatively more energetic pulses at low frequencies.
While the exact physics for the generation of GPs is currently unknown, several mechanisms were proposed that could lead to their generation. One of those considers induced Compton scattering of radio photons off ultrarelativistic secondary plasma particles, leading to a focusing of the radiation (Petrova 2004). According to this mechanism, there would be a constant maximum GP energy, since more focusing means stronger but narrower pulses. While this was originally applied to GPs, our results suggest that a constant maximum single-pulse energy might exist in many pulsars.
In order to extract single-pulse widths, the standard procedure is to assume a top-hat shape and divide the pulse area by the peak flux, thereby keeping the energy identical, but underestimating the true width. In Figure 7 we plot 50-bin histograms for these so-called equivalent widths. Pulses with a negative energy were removed, as well as weak pulses with a peak flux below 5 off-pulse standard deviations to reduce the number of outliers. The vertical magenta lines represent the equivalent width of the average profile, while the real W3 profile width is listed for comparison. All W3 profile widths turn out to be about 2-3 times broader than the equivalent profile widths, which was to be expected.
For B0823 there appear to be many single pulses that are wider than the average profile. We wanted to verify that pulses with large equivalent widths are actually wider than pulses with small equivalent widths and determine how these compare with the profile width. This is shown in Figure 8 for the two high-S/ N pulsars B0540 and B0823. In the left panel for B0540 single pulses in the equivalent width ranges   milliperiods. None of the shown profiles were scaled or shifted in phase. As can be seen, the pulses with large equivalent widths are indeed wider than the pulses with small equivalent widths, and hence the equivalent width turns out to be a useful proxy.
Several previous studies showed that averaging single pulses from intensity bins results in profiles of different shapes and widths. This was demonstrated for the Vela pulsar (Krishnamohan & Downs 1983); for B0540+23, B1133+16, and B1541+09 (Nowakowski 1991(Nowakowski , 1996; and for B0329+54 (McKinnon & Hankins 1993). In all cases the profiles turned out to be intensity dependent, with weaker profiles being generally wider. Invoking the dipolar magnetic field geometry, a natural explanation would be that weak profiles, being wider, are emitted at higher altitudes. Contrary to this, we find from Figure 8 that two profiles with about the same intensity can have a very different width (e.g., the [0-2] and [10-12] milliperiod width bins for B0823 both belong to roughly the [0-2] Jy intensity bin), thereby ruling out any relation between intensity and width or emission height. It seems that averaging in intensity bins will collect single pulses over a certain width range and produce an apparent intensity-height correlation.
In Figure 9 we show a single-pulse scatterplot with the equivalent width in milliperiods, peak intensity in Jy, and energy in mJy · s (color bar on top), which we refer to as a WIE plot. The horizontal magenta line gives the mean equivalent width for each pulsar, which can be compared with the histogram peaks of Figure 7. In each plot there are two pulses represented by gray circles (B0301 has only one), which belong to pulses that were not included in the color map. In this way, the visual representation is significantly enhanced without hiding any data. It turns out that each of the pulsars follows the same trend in the plots, where the maximum width keeps decreasing with increasing intensity, and we believe that it could indicate a general behavior. By choosing a certain "width-intensity region" from the plot and forming a profile, one can access a certain emission altitude with a given singlepulse intensity range and investigate, e.g., component evolution as shown in Figure 8.

Autocorrelation of Peak Fluxes and Energies
We employed autocorrelation analysis of single-pulse energies and peak fluxes to search for periodicities. The lag-1 autocorrelation coefficient, obtained from the normalized autocorrelation function (ACF), is a measure of nonrandomness, being close to 0 for a sequence of random numbers and approaching 1 for a highly periodic function. Early ACF studies yielded for the lag-1 coefficient of pulse amplitudes on B0833-45 about 0.4 (Jones & Wielebinski 1970) and on B0329+54 (CP 0328) about 0.3 (Shuter & McCutcheon 1970). Figure 10 shows our lag-1 coefficients for on-pulse and offpulse energies and peak fluxes. For the autocorrelation analysis no pulses were excluded and the lag-1 coefficient was computed for each 1000 consecutive single-pulse energies and peak fluxes. Strikingly, all lag-1 coefficients for each pulsar are roughly constant over the entire time series, with the lag-1 coefficients for the single-pulse energies always being larger than the ones for the peak fluxes. We believe that this could signify that single-pulse energies encode more fundamental information than peak fluxes, which is the reason why we decided to focus more on energies in this study. Table 5 lists all the mean values in Columns (2)-(4), with pulsars B0540 and B0301 having lag-1 coefficients for energies over 0.5, which indicates highly nonrandom behavior. For the off-pulse energies the means are all about zero, indicating that they constitute a sequence of random numbers as expected. The lag-1 mean for the peaks of  with the profiles computed for a narrow (orange dotted) and a wide (green solid) equivalent width bin, comprising 5 milliperiods for B0540 and 2 milliperiods for B0823. We note that the profile heights for the subsets are not scaled. B0611 is only marginally larger than the random number values, which could be due to the low S/N in this pulsar.
We now plot the full ACFs in order to obtain information on any periodicities in energy and intensity (i.e., pulse peak flux). Besides the ACF, another common technique is to compute longitude-resolved fluctuation spectra (LRFS), which are obtained by performing a 1D discrete Fourier transform (DFT) along a constant pulse longitude in a pulse stack. In this way one obtains fluctuation frequencies across the profile and can extract any P 3 modulation periods (e.g., Taylor & Huguenin 1971;Edwards & Stappers 2003;Weltevrede et al. 2006). In order to find any subpulse drifting features (denoted by P 2 ), the approach is to select a phase longitude range and to perform a 2D DFT, leading to the 2D fluctuation spectrum (2DFS). Since our main focus are single-pulse energies, we do not need the phase-resolved information and choose the ACF as the primary analysis tool. Both methods, however, yield consistent information.
We show in Figure 11 the normalized ACFs for energies and intensities obtained from a sequence of 100 consecutive pulses, with pulse numbers indicated in parentheses. The 50-pulse sequences displayed in Figure 1 are part of these 100 consecutive pulses. We chose these particular sequences of 100 pulses, because their ACF has well-defined oscillations, making the determination of the period very easy by eye. The main purpose here was to establish and compare ACF oscillations for energies and intensities and to confirm any published periodicities. ACF oscillations have been reported previously, but not all pulsars were found to produce them (e.g., Ritchings 1976). We notice that in all cases the intensity ACF traces well the energy ACF, except for B0611, where the intensities do not show clear oscillations. Table 5 lists the periods extracted from the energy ACFs by taking the center of the cusp of the first oscillation maximum; it is standard practice to identify the first ACF maximum with the modulation period (Kramer et al. 2002;Mitra et al. 2015;McKee et al. 2019). The periods for the two D-class pulsars are much longer than for the rest of the pulsars owing to extensive nulling. An early study had already shown oscillations in the ACF obtained from B1919+21 (CP 1919), with the first maximum at about 4 periods (Conklin et al. 1968). Later ACF and fluctuation spectra studies showed that many pulsars have periodic pulse amplitude modulations, typically being in the range of 2-15 periods (Sutton et al. 1970;Taylor et al. 1975;Rankin 1986

Maxmaps and Search for Drift Bands
To search for any drifting, we show in Figures 12 and 13 a close-up of the respective maxmaps from Figure 3 for the corresponding 100-pulse sequences. The top panels contain the average profile of all pulses (identical to Figure 2), with the 3% points indicated by the vertical lines. The bottom panels show the pulse stacks for the 100-pulse sequences, with the color bar indicating the flux density in units of Jy. Figure 9. Equivalent width as a function of peak flux with corresponding energy (WIE) plots with color bar in units of mJy · s on top. The horizontal line indicates the mean equivalent width. Single pulses with negative energies or peak flux below 5σ were not considered. Each subplot has one to two pulses represented by gray circles; these pulses were "shaved off" to better visualize the approximately constant energy envelope.
B0301+19. Figure 13 shows that there are two extended ranges of nulling. In the top panels the nulls appear as lowintensity, blue circles that extend to outside the vertical 3% lines. The 50-period modulation derived from the ACF is most likely resulting from the on-off nulling behavior and could be referred to as supermodulation. In the top panels of Figure 13 no organized drifting can be identified. In a previous 2DFS study this pulsar has been classified as a diffuse drifter with occasional drifting in the trailing component (Weltevrede et al. 2006).
B0525+21. In Figure 13 we can see three short nulling ranges, which likely also determine the 42-period supermodulation derived from the ACF, and no organized drifting is present. As can be seen from the average profile, the emission in the broad saddle region is suppressed to only 4% of the peak flux density at our frequency of 1.7 GHz. Interestingly, in the pulse range from 186 to 200 the two main components are alternatingly illuminated, creating the impression of a 2-period intensity modulation. Occasional diffuse drifting in both components has been previously reported (Weltevrede et al. 2006).
B0823+26. From the top panels of Figure 12 one can see that B0823 does have drift bands in both directions and only a few, mostly single nulls. While the pulse stack in the bottom panels shows very clearly the 5-to-6-period-long intensity modulations, it is the maxmap in the top panels that displays the drifting better. This pulsar has been classified as a diffuse drifter with a positive preferred drift direction, where Figure 10. Lag-1 autocorrelation coefficients for 1000 consecutive single pulses for on-pulse energies (blue circles), off-pulse energies (green diamonds), and peak fluxes (orange crosses). Off-pulse energies contain random noise, and the coefficients are around 0. The values on the abscissa refer to the single pulses used, e.g., the value 4 includes pulses 4000-4999.  Rankin (1986).

B0540+23(A) and (B).
Here we can see only few nulls and also short drift bands (about 5 periods) with no preferred direction. An intensity modulation is clearly visible in the bottom panels. Previous studies also found some short drift bands but do not report any modulation periods (Nowakowski 1991;Weltevrede et al. 2006). We confirm the short drift bands without preferred direction, but to the contrary, we do find a very clear modulation period P 3 of 15 periods (16 for the B observation).
B0611+22. For this pulsar we find almost no nulls, no drift bands, and no obvious intensity modulation in Figure 12, but we can still see oscillations in the energy ACF in Figure 11. Since no clear accompanying peak flux modulation is present, we have cautiously assigned a modulation of 10 periods as derived from the ACF first maximum. This pulsar has no previously reported modulation (Weltevrede et al. 2006;Seymour et al. 2014).

Cumulative Distribution Function for Energies
A pulsar survey has shown that the energy distribution for regular pulses in many pulsars obeys lognormal statistics (Burke-Spolaor et al. 2012); a lognormal distribution of energies was also found to exist in an MSP (Liu et al. 2016). Two previous reports specifically address some of our pulsars, and in both it was reported that the single-pulse energies follow lognormal statistics, for B0611 at 0.3 and 1.4 GHz (Seymour et al. 2014) and for B0823 at 0.15 GHz (Sobey et al. 2015). On the other hand, another large study finds that for many pulsars a lognormal energy distribution does not provide the best fit, but rather a newly suggested empirical "power-law divided by exponential" distribution (Mickaliger et al. 2018). As long as there is no unified theory for the emission mechanism, there will be no consensus on the distribution statistics. We are proposing below another distribution for single-pulse energies of regular pulses, which is based on the Weibull statistics and to the best of our knowledge has not been applied yet in this context. Figure 14 shows the cumulative probability for the singlepulse energies, including negative energies, but eliminating the largest energy value, which was generically assumed to be an outlier. The cumulative probability represents the probability that the single-pulse energy takes on a value less than or equal to x in units of mJy · s; the corresponding cumulative distribution function (CDF) can be expressed as F E (x) = P(E x). The insets expand the region around the origin and show that F E (0) ≠ 0 owing to the presence of negative energies representing pulse nulling. The curves for B0301, B0525, and (barely visible) B0823 have an inflection point, which we identify with the NF. For B0540 and B0611, where no inflection is visible, we take the ordinate intersection as the NF. The dashed lines in the insets show the intersection and inflection points with the corresponding cumulative probabilities, and we believe that this gives a useful NF estimate. Figure 11. Normalized ACFs for single-pulse energies (solid blue) and peak fluxes (dashed orange) of the 100 pulses indicated in parentheses. A 50-pulse sequence of those 100 pulses is shown in Figure 1. These particular 100-pulse sequences show very pronounced ACF oscillations, but many others do not. The period of the oscillations is given in Table 5. B0301+19 has an estimated NF of 15%, which agrees well with another study deriving 14% ± 4% (Herfindal & Rankin 2009). For B0525+21 we estimate 24%, which also agrees well with previous studies that are estimating 25% ± 5% (Ritchings 1976) and 28% ± 2% (Herfindal & Rankin 2009). In B0823+26 we estimate an NF of 4.9%. There is some spread in the published values, ranging from 1.8% (Sobey et al. 2015) to 3.9% (Basu & Mitra 2019) to 7% ± 2% (Herfindal & Figure 12. Top panels: average profile of all pulses (see Figure 2), with the vertical lines indicating the 3% points. Each circle in the scatterplot represents a single pulse with the corresponding peak flux and phase of the peak (see Figure 3); red colors indicate strong pulses and blue colors weak pulses. This same pulse sequence is also used for the ACFs in Figure 11. Bottom panels: corresponding pulse stack, with the color bar indicating flux density in Jy. Rankin 2009). B0540+23 does not have an inflection point, and so we estimate the NF from the ordinate intersection, which then presumably yields a lower limit. The A observation yields 1.5% and the B observation 2.6%, while a previous study reported no nulls (Nowakowski 1991). Apparently, there is a smallest NF below which no inflection point appears, and the only possibility then is to use the intersection. Finally, for B0611+22 we estimate a mere 0.6%. We realize that this approach is somewhat experimental; nevertheless, our values are generally in good agreement with the literature.
In Figure 14 we also show the outcome of a fit with the Weibull CDF to the cumulative probabilities, where the red overlaid curve is the fit and given by where λ > 0 is the scale parameter and k > 0 is the shape parameter. The corresponding shape and scale parameters, together with χ 2 and the two-sided Kolmogorov-Smirnov (K-S) test statistic rejecting the null hypothesis, appear in the main panels of the figure. Since the cumulative distribution crosses into negative energies, the fit was limited to the indicated range but extrapolated to the origin in the plot. Fits of the cumulative probabilities with the lognormal CDF according to where μ is the mean, σ > 0 the standard deviation of the variable's natural logarithm, and erf the error function were also performed as a comparison (not shown). The fit range was kept the same in both cases to obtain a fair comparison. It turns out that for all pulsars the χ 2 and the K-S statistic for the lognormal fit are larger than for the corresponding Weibull fit. The lognormal CDFs have a tendency to slightly overestimate the number of higher-energy pulses, and hence we conclude that at least for these pulsars the Weibull distribution yields the better fit. Figure 15 presents the average profile of B0823 normalized by the MP peak to show details of the IP and PC. As can be seen, the PC is narrower and brighter than the IP and connected by a bridge emission with the MP. The left inset of Figure 15 shows the lag-1 autocorrelation coefficients for the single-pulse Figure 14. Cumulative probability for the energies (blue circles) and fit with the Weibull CDF (red line). The fit range, shape and scale parameters, χ 2 , and K-S statistic are given inside the plot. The inset expands the region around the origin for estimating the NF. The shape parameters for both B0540 observations are identical. Figure 15. Expanded view of the average profile of B0823 with the IP and PC in detail; the ordinate is in units of normalized flux density. The left inset shows the lag-1 autocorrelation coefficients for the PC energies computed for sequences of successive 1000 pulses. The right inset shows the PC over the range used to calculate the energies (samples 3275-3475); note the doublepeaked structure of the PC. PC energies obtained for phase bins 3275-3475. Each data point corresponds to 1000 consecutive PC energies, analogous to the presentation in Figure 10. All lag-1 values are close to zero, suggesting that there are no correlations, just like the offpulse energies, which might, however, be due to the obviously very low S/N level. One recent study is deriving power spectra from the LRFS for the MP and the PC and finds for both a 5.5period modulation (Basu & Mitra 2019), albeit at a much lower frequency of 0.3 GHz. Table 6 summarizes our IP and PC heights relative to the MP peak flux and the IP-MP and MP-PC separations and lists for comparison the corresponding values of several previous studies. As expected, the emission becomes weaker at higher frequencies. The MP-PC separation increases by almost 20% at higher frequencies, while the IP-MP separation shows no frequency dependence.

Shape and Energy at Different Emission Heights
Assuming a dipolar magnetic field with a filled beam, a pulse with a narrower width is emitted closer to the neutron star surface, and a wider pulse is emitted at higher altitudes (Lorimer & Kramer 2005 and references therein). While the assumption of a dipolar field geometry is probably valid for normal pulsars, Neutron Star Interior Composition Explorer (NICER) observations of thermal X-ray pulses from the surface of an MSP have a more complicated, multipole magnetic geometry suggested (Bilous et al. 2019 and references therein). Furthermore, there are indications that the emission of some pulsars does not completely fill the flux tube, being limited by the pair plasma supply (Rankin et al. 2020), and so it is not straightforward to calculate absolute emission heights for a given profile width. Nevertheless, for normal pulsars over a narrow emission bandwidth, it is to be expected that a wider profile generally corresponds to a higher emission region.
Hence, under this assumption, selecting a certain equivalent width range permits accessing the emission cone at different altitudes and evaluating specific emission properties, such as number of components. In Figure 8 it becomes evident for B0540 that the wider profile (solid green line) has a pronounced shoulder at the trailing edge, indicating the presence of a second component, as was already reported for this pulsar at 0.43 GHz (Nowakowski 1991). On the other hand, the narrow profile (dotted orange line) has a singlecomponent shape. One might speculate whether the emission is evolving from a core-single S t profile at lower heights to a conal-double D profile at larger heights.
The emission peak for B0823 shifts by 1°.1 toward later phases at larger emission heights, while maintaining its nominal S t class. This is reminiscent of the shift toward later phases of the peak of the wider Q-mode with respect to the narrower B-mode. A Q-mode shift of about 3°at 0.3 GHz (Rankin et al. 2020) and of just over 1°at 1.5 GHz (Sobey et al. 2015) has been reported. While the amount of the shift in our case is equal to the reported Q-mode shift, the peak intensity remains the same. Since the Q-mode is conal and 100× weaker, suggesting any relation between the two cases is tempting but currently speculative at best.
The width-intensity-energy plots in Figure 9 suggest that the falling-off envelope for the maximum width represents a constant maximum energy for a single pulse. Several previous studies of MSPs show similar width-versus-peak flux plots, which demonstrate that stronger single pulses are generally narrower than weaker single pulses, but without discussing corresponding energies (Jenet et al. 1998;Vivekanand 2000;Liu et al. 2016). Two other works examine GPs from the Crab pulsar (Popov & Stappers 2007;Majid et al. 2011) and conclude that stronger GPs are generally narrower than weaker GPs. All those results are either for GPs or for pulses from MSPs, while we show here that such a relationship is also valid for regular pulses from normal pulsars. Hence, it appears to be a general phenomenon across all pulsar categories.
As mentioned above, an energy-limiting mechanism was proposed in the context of the induced scattering model to explain the generation of GPs, where an increase in intensity of the GPs was due to focusing. According to this model, the narrowing of the pulses leads to a maximum constant energy (Petrova 2004). Due to a varying scattering efficiency, not all GPs will be emitted with the maximum possible energy and hence populate the plot below the envelope. Another large study of single-pulse energetics suggests that there may exist a self-balancing effect that regulates the net energy contained across the full emission phase (Burke-Spolaor et al. 2012). This energy regulation would limit the emission at certain phases if there are outbursts at other adjacent phases, thereby maintaining a fairly steady maximum total energy output per pulse. We can speculate whether some focusing mechanism or other type of energy-regulating mechanism perhaps is prevalent in all pulsars and so might also explain our WIE plots.

Energy, Intensity, and Phase Modulations
ACFs can provide information on any underlying periodicity in single-pulse energies and intensities from the lag of their first maxima. In Figure 11 we find ACF oscillations in all our pulsars, including two without previous reports. In all but one Note. The heights are given relative to the MP; our separations consider the peak-to-peak distance; the values for Rankin+2020 and Young+2012 were estimated from their figures.
case do the energy and intensity ACF oscillations trace each other, however, with the energy oscillations showing clearer periodicities. Following is a brief discussion broken down by pulsar.
B0301+19. This pulsar was reported to have a short-period modulation of 6.4 periods from the leading component, while the saddle region has a modulation of about 20 periods (Rankin 1986 and references therein). LRFS/2DFS studies at 1.4 GHz classify this pulsar as a diffuse drifter with a broad feature in the power spectrum and assign a 5.2-period modulation (Weltevrede et al. 2006 and references therein). Our analysis identifies a modulation period of about 50 periods, which we interpret as a supermodulation resulting from the extensive nulling. Interestingly, there is one study that reports a 51-period modulation feature in LRFS (Herfindal & Rankin 2009), which is almost exactly the same periodicity we identify.
B0525+21. A 5.6-period modulation and a 42-period modulation produced by the saddle region between the two components have been reported for this pulsar (Rankin 1986 and references therein), as well as a 3.8-period modulation (Weltevrede et al. 2006 and references therein). Our analysis finds also a 42-period feature in the ACF, even though this resulted from single-pulse energies, which of course cannot distinguish between component and saddle modulation. Possibly, longer modulation periods in D-class pulsars describe supermodulation from nulling.
B0823+26. Many studies focused on this pulsar, and the P 3 modulation periods found at frequencies below 0.5 GHz are in the range 5-5.6 (Lang 1969;Taylor et al. 1969;Nowakowski et al. 1982;Rankin 1986;Sobey et al. 2015;Basu & Mitra 2019), while at 1-2 GHz frequencies the reported period is in the range 6.7-7 (Nowakowski et al. 1982;Weltevrede et al. 2006). All those quoted studies were applying LRFS analysis, while we use the ACF method of pulse energies and derive a 5.3-period modulation, being closer to the lower frequency values.
B0540+23 and B0611+22. In both pulsars no P 3 modulation has been previously observed (Rankin 1986;Nowakowski 1991;Weltevrede et al. 2006). Our results for B0540 show a very clear 15-period energy and intensity modulation in 2016 and a 16-period modulation in 2019. For B0611 we find a possible 10-period energy modulation, but no clear amplitude modulations are visible. Since P 3 has been usually associated with drift features in conal pulsars related to spark circulation, this could indicate that energy modulations are a common phenomenon in all pulsars. Modulated emission has been observed recently in the core component of the pulsar B1946 +35, which cannot be explained by drift resulting from spark circulation (Mitra & Rankin 2017). This modulation was attributed to an emission-pattern change that hitherto has not been reported and might be responsible for modulations in other core-single pulsars as well.
Phase modulations and drift bands, on the other hand, can be examined in the maxmaps from Figures 12 and 13. We find no organized drift features in the two D-class pulsars B0301 and B0525, which is consistent with the spark circulation and conal emission models, assuming a line of sight cutting through at large angles. In this model only a tangential cut through the cone could produce appreciable drifting. Hence, we can speculate whether the single component in B0823 and B0540 actually might represent such a tangential cut, since among our pulsars only these two have occasional short drift bands. Alternatively, it could also indicate that some pulsars show a mixture of core and conal emission, such as was found for B0823 regarding its B/Q-modes (Rankin et al. 2020).

Are Energies Weibull Distributed?
In Figure 14 we suggest to fit the distribution of energies with the Weibull CDF (Equation (1)), instead of using the usually assumed lognormal CDF. With a shape parameter of k = 1, the Weibull CDF becomes the exponential distribution, which describes the time between independent events of a Poisson point process, such as the time between radioactive decays. With k = 2 it becomes the Rayleigh distribution, which can describe, e.g., a wind speed distribution in two dimensions, assuming that each velocity component is independent and normally distributed. The likely total power output of a wind turbine per year is usually assessed using a Weibull distribution with a variable shape parameter. A shape value approaching one would indicate many calm days and occasional days with strong gusts, while a shape value approaching 3 (and above) means many days with a fairly steady breeze and wind speeds around the average (e.g., Serban et al. 2020;Shu & Jesson 2021). A Weibull distribution with shape values above 3 resembles a normal distribution.
The Weibull distribution has also been applied in pulsar and fast radio burst (FRB) research. One study investigates the statistics of null and burst durations in pulsar B0031-07 and finds that the processes can be fit with an exponential (k = 1) and a Rayleigh (k = 2) distribution (Vivekanand 1995). Another study describes the intervals between bursts of a repeating FRB using a Weibull distribution with shape parameter k = 0.34, allowing for non-Poissonian statistics favoring clustering (Oppermann et al. 2018). Also, one study applies the Weibull distribution to study the timescale of mode changes in pulsar J1326-6700 and finds a shape parameter of 0.8 (Wen et al. 2020).
Our original idea was to search for a Rayleigh energy distribution, which was based on the geometry of the emission cone. Some observational analyses suggest that the emission cone might contain beamlets that sweep across the line of sight, thereby giving rise to single-pulse microstructure (Kramer et al. 2002;De et al. 2016). Among the pulsars investigated in this work, microstructure has been previously observed in B0301, B0525, B0540, and B0823 (Lange et al. 1998;Kramer et al. 2002;Mitra et al. 2015), raising the question whether those presumed beamlets could possibly have a lateral random jitter motion, superimposed on the periodic rotational motion. If so, the beamlets could only be moving in two orthogonal directions, and we wanted to see whether this might be producing a Rayleigh distribution in the single-pulse energies. We found that the cumulative distribution of single-pulse energies can indeed be fit very well with a Rayleigh distribution, but only the upper higher-energy part. Therefore, we decided to fit the more general Weibull CDF with a variable shape parameter, which turns out to closely fit all energies.
The shape parameters were found to be between 0.94 and 2.45, with the two B0540 observations yielding exactly the same value of 1.30, albeit with a slightly different scale. The fact that both B0540 observations yield the same shape parameter indicates that this pulsar did not change its emission characteristics between 2016 and 2019 and also demonstrates that all our analysis steps are fully reproducible. The 9% smaller scale parameter compares well with the 10% smaller average profile mean flux in Figure 2. We attribute the low shape parameter for B0823 to strong scintillations, which amplify some pulses to high energy values, thereby producing the "occasional strong gusts," metaphorically borrowed from the wind discussion above. Conversely, for B0611 we forecast a "fairly steady breeze."

Some Speculations on the IP and PC in B0823+26
Even though at 1.7 GHz the PC height is 0.48% of the average MP and the IP only 0.19%, we can still extract some characteristics and make a few speculations. The right inset of Figure 15 shows a close-up of the PC, where we can see that it has a well-resolved double-peaked structure with nonzero bridge emission at earlier phases. To our knowledge this is the first observation of a double-peaked PC, which almost resembles a perfectly symmetric conal-double D profile. On the other hand, the IP resembles a slightly asymmetric coresingle S t profile. One study at 0.3 GHz shows close-ups of the IP and PC, and it is interesting to note that their IP has exactly the same shape as in our case, while their PC seems to have multiple components (Rankin et al. 2020). This implies that the PC has a rather strong frequency dependence, unlike the IP.
We show in Figure 4 a particularly bright single-pulse IP (pulse no. 447) and single-pulse PC (pulse no. 7931). The peak location of that particular single-pulse PC is at phase sample 3355, which is exactly at the center of the left peak of the average PC, thereby demonstrating that the double-peaked PC profile constitutes indeed two separate components. The intensity of the single-pulse IP displayed in Figure 4 is 340 mJy, about 230 times stronger than the average IP, while the energy is only 27 times higher owing to the pulse being extremely narrow. It might be justified to call this particular single-pulse IP a giant IP. The fact that the giant IP is extremely narrow supports the conclusion drawn from Figure 9, where the brightest single pulses were shown to be the narrowest, thereby limiting the maximum possible single-pulse energy output. This same energy-limiting mechanism obviously also applies to IPs.
In conclusion, we observed five normal bright pulsars at 1.7 GHz with 100 MHz bandwidth and analyzed single-pulse energies, peak fluxes, and equivalent widths. Our results suggest that there is a pulsar-specific upper limit for singlepulse energies, limiting the maximum width for the pulses and appearing as an envelope in a width-intensity plot. Assuming that a wider pulse generally corresponds to a higher emission altitude, selecting pulses with certain widths could allow, e.g., investigating the component evolution of the average profile. An autocorrelation analysis of single-pulse energies shows that all observed pulsars have a periodic energy modulation. We find a so-far-unreported energy modulation of 15 periods in B0540+23 and 10 periods in B0611+22. From the cumulative distribution of the energies it is possible to estimate the pulse NF, which for B0611+22 is a very low 0.6%. A fit of the cumulative energy distribution shows that a Weibull function yields a better fit than the usually assumed lognormal function. Finally, for B0823+26 we find that the PC has a symmetric double-peaked structure and we recorded a giant IP.