A Pilot Study of Nulling in 22 Pulsars Using Mixture Modeling

The phenomenon of pulsar nulling, observed as the temporary inactivity of a pulsar, remains poorly understood both observationally and theoretically. Most observational studies that quantify nulling employ a variant of Ritchings (1976)'s algorithm which can suffer significant biases for pulsars where the emission is weak. Using a more robust mixture model method, we study pulsar nulling in a sample of 22 recently discovered pulsars, for which we publish the nulling fractions for the first time. These data clearly demonstrate biases of the former approach and show how an otherwise non-nulling pulsar can be classified as having significant nulls. We show that the population-wide studies that find a positive correlation of nulling with pulsar period/characteristic age can similarly be biased because of the bias in estimating the nulling fraction. We use our probabilistic approach to find the evidence for periodicity in the nulls in a subset of three pulsars in our sample. In addition, we also provide improved timing parameters for 17 of the 22 pulsars that had no prior follow-up.


INTRODUCTION
Pulsar nulling, initially observed by Backer (1970a), is the absence of observed emission from a pulsar for one or more pulse periods. Observationally, the phenomenon of pulsar nulling remains poorly understood. It is clear that nulling is a broadband phenomenon, observed from 102 MHz (Davies et al. 1984) to 8.35 GHz (Honnappa et al. 2012). However, it is not firmly established whether nulling is simultaneous over this frequency range using a large sample of nulling pulsars. Prior studies found contradictory conclusions. For example, observing over two frequency ranges, 50-140 MHz and 275-430 MHz, Taylor et al. (1975) found that nulls are simultaneous in two different pulsars (PSR B0031−07, PSR B0809+74), while Davies et al. (1984) found the evidence for excessive nulls in single pulses at 102 MHz compared to 406 MHz in PSR B0809+74. A more recent study by Gajjar et al. (2014a) found that the nulls are highly coherent in three pulsars at four different frequencies -313, 607, 1380, and 4850 MHz. In addi-Corresponding author: Akash Anumarlapudi aakash@uwm.edu tion, it is also not clear whether pulsars null randomly. Redman & Rankin (2009) and Gajjar et al. (2012) found that nulls might not occur randomly but might be clustered, where nulls and bursts tend to occur in groups, but the latter found that the null durations can be random. However, for many of these results the dependency of the nulling inferences on signal-to-noise ratio makes it hard to robustly interpret their findings.
Although the formation of a pair cascade and the radiation from these accelerated pairs in the pulsar magnetosphere is often invoked to explain the observed emission from a pulsar (Ruderman & Sutherland 1975), a full theory of pulsar magnetospheres and its emission to explain the diverse morphology in pulse profiles and phenomenology is yet to be developed. As such, the theory of pulsar nulling remains largely speculative, though it is often attributed to one of two classes: i) inherent to the magnetosphere itself such as loss of coherence condition required for radio emission, e.g., Filippenko & Radhakrishnan (1982), or the depletion of pairs in the magnetosphere themselves, e.g., Kramer et al. (2006) or ii) geometrical factors external to the magnetosphere such as the line of sight traversing through the 'empty' region between rotating emission carousels, e.g., Herfindal & Rankin (2007, 2009. Further progress may require ad-arXiv:2301.13258v2 [astro-ph.HE] 13 Sep 2023 ditional observational data to understand how the properties of nulling relate to the properties of the pulsars themselves. Nulling as a phenomenon may be related to other more extreme forms of intensity modulation, where the pulses can disappear for hours to months in the cases of rotating radio transients (RRATs; McLaughlin et al. 2006) or intermittent pulsars (Kramer et al. 2006;Lyne 2009). However, the connection between these populations is not clear. Furthermore, pulsar nulling is often discussed in tandem with two other forms of single pulse variations: mode changing -a phenomenon in which an otherwise stable pulse profile switches between multiple shapes (or modes) (Backer 1970b) and sub-pulse drifting -a phenomenon in which the single pulse phase shows a uniform periodic drift (Drake & Craft 1968). Regardless, in all of these cases the appearance of these phenomena can be limited by instrumental sensitivity: without enough sensitivity to probe single pulses at high significance, one cannot be certain whether the pulsar emission is truly missing during the nulls or the pulsar switches to an alternate mode with lower intensity. Together all three are often thought of as different representatives of a larger underlying phenomenon of sub-pulse intensity variations (Lorimer & Kramer 2004).
Nulling is usually quantified by the fraction of pulses where there is no discernible emission, called the Nulling Fraction (NF). NF can vary from 0 -in the case of standard emission picture that shows no nulls -to 1, in the extreme case where the pulsar emission is visible only between long nulls (intermittent pulsars and RRATs). NF has been measured in roughly 8% of pulsars, but this has more to do with the lack of single pulse studies as opposed to nulling being restricted to a small subset of pulsars. This smaller data set of nulling pulsars is entirely restricted to normal (not recycled) pulsars, owing to the high sensitivity demands that would be needed to observe single pulses of millisecond pulsars (MSPs), although some recent studies  have been conducted in a sample of bright MSPs which did not find a signature of nulling with high confidence. In addition, there can be a bias against discovering normal pulsars which tend to have a high NF, or are intermittent. Hence the fraction (8%), can only be considered as a conservative lower limit.
Such a small data set restricts our ability to infer population-wide properties, which might give clues to the origin of the phenomenon, and hence studies done thus far have not reached a consensus. An initial study done by Ritchings (1976) claimed a correlation between NF and pulsar period (with longer period pulsars experiencing higher NF) and also a stronger correlation with the characteristic age. Wang et al. (2007) also suggested a correlation with spin-down age, albeit qualitatively, with older pulsars experiencing higher NF, before eventually crossing the death line. Konar & Deka (2019) found that there may be two different populations of pulsars separated by a NF of ∼ 40% but did not find correlations with any intrinsic pulsar properties, while Sheikh & MacDonald (2021) claimed that there is no strong evidence for the existence of two sub-populations. All of these studies may be significantly biased since the samples used are restricted to the pulsars that explicitly showed nulling.
In general, most studies (Wang et al. 2007;Gajjar et al. 2012Gajjar et al. , 2014bHerfindal & Rankin 2009) estimate NF using the methodology (or a variant) proposed by Ritchings (1976). But as Kaplan et al. (2018) demonstrated, this method can suffer strong biases in the case of weaker pulsars which can lead to overestimating the NF and classifying an otherwise standard weak pulsar as a nulling pulsar. This can also lead to systematic biases in population inferences. In addition, Kaplan et al. (2018) proposed an alternate method in which they use Gaussian Mixtures to model the single pulse intensities and estimate the NF , and demonstrate the reliability of this method in accurately measuring the NF in weaker pulsars. In this study, we expand on the Gaussian Mixture Model (GMM) of Kaplan et al. (2018) 1 to generalize their method and apply it to a larger sample of 22 pulsars 2 .
Pulsars selected for this study were discovered as a part of the Green Bank North Celestial Cap (GBNCC) pulsar survey (Stovall et al. 2014) in 2-min drift scans at 350 MHz with a 100 MHz bandwidth and with data sampled every 81.92 µs. At 350 MHz the beam size is 36 ′ (Full Width at Half Maximum; FWHM) and hence the astrometric precision prior to a coherent timing solution is limited by the beam size depending on the Signal-to-noise ratio (SNR) of the discovery candidate. These were later followed up at the Green Bank Telescope (GBT) and Arecibo Observatory (AO) to improve their timing solutions and establish their nulling characteristics.
The structure of this paper is as follows: In Section 2, we detail our data acquisition and reduction methods, and provide updated timing solutions for the pulsars in this study. We then describe the mixture model and provide our basic results in Section 3. Finally, we present the implications of the results in Section 4 and conclude in Section 5.

Observations and Data Reduction
A sample of 22 recently discovered pulsars was selected for this pilot study if they showed any signs of intermittency in their discovery plots 3 . Data for 15 out of 22 pulsars were collected using the 100-m Robert C. Byrd Green Bank Telescope (GBT) (hereafter referred to as the GBT sample), operating at 820 MHz with a bandwidth of 200 MHz, in 2 hr contiguous scans, with the primary aim of determining the pulsars' nulling characteristics (project code 18A−436; PI: J. Swiggum). Data for another nine pulsars were collected at the 300-m William E. Gordon Arecibo Observatory (AO) operating at 430 MHz over a bandwidth of 24 MHz, with the goals to both establish coherent timing solutions and determine nulling characteristics (project code P3436; PI: J. Swiggum) (hereafter referred to as the AO sample). Two pulsars in our sample, PSR J0414+31, and PSR J1829+25, were observed at both observatories.
Six of the 15 pulsars in the GBT sample already had coherent timing solutions (Lynch et al. 2018) and the data for these were collected in coherent search mode using the Green Bank Ultimate Pulsar Processing Instrument (GUPPI; Ransom et al. 2009) with 128 frequency channels sampled at 10.24 µs and retaining full polarization information. The remaining nine pulsars had no prior follow-up campaigns and so we first improved their positions using gridding observations and then observed them in incoherent search mode with 2048 frequency channels sampled at 40.96 µs. Data for the AO sample were collected in coherent search mode using the Puerto Rico Ultimate Pulsar Processing Instrument 4 (PUPPI), with 64 channels sampled at 40.96 µs, over a span of ∼ six months to establish coherent timing solutions in addition to studying the nulling properties. A summary of observations for each pulsar is provided in Tables 1 and  2. Starting with the raw search mode data, we used dspsr (van Straten & Bailes 2011) to fold the data. We then used pazi, the interactive zapping routine in psrchive ) to remove radio frequency interference (RFI)-affected frequency channels and single pulses. For GBT data, we also made use of RFI scans taken at the observatory 5 , when available, to 3 See the GBNCC discovery page: http://astro.phys.wvu.edu/ GBNCC. 4 http://www.naic.edu/puppi-observing/ 5 https://greenbankobservatory.org/rfi-gui-user-guide/ identify the frequency bands that are affected by RFI, which are otherwise not obvious visually. In some cases, we found that one of the polarization channels was persistently affected by RFI, and in such cases we excluded data from that polarization channel at the cost of SNR. Fortunately, this did not have a significant impact on the determination of the nulling fractions. Some of the AO data had periodic "drop-outs" in the data with submillisecond periodicity at zero dispersion measure (DM), caused by data rate overflow during the observations. We cleaned these "drop-outs" by replacing the data with NaN values and being careful to exclude those when folding/averaging. After cleaning the RFI, both for timing and estimating nulling, we averaged polarizations to measure the total intensity.

ON/OFF histograms
Once we had improved the timing solution, we used dspsr in single pulse mode to generate single pulses for all scans and used psradd, from psrchive, to phase align pulses from different scans after cleaning the data for RFI. We then averaged the data along the polarization and frequency axes to obtain the pulse intensity of the single pulses as a function of the rotational phase and generated single pulse stacks such as that shown in Figure 2.
The most important aspect in estimating the nulling fraction is determining the "ON"-pulse and "OFF"- Note-Quantities in parentheses are 1σ uncertainties on the last digit. a Coherent timing solutions are given in Lynch et al. (2018) b Timing solution is obtained by combining AO and GBT data.   Table 3. pulse phase windows. The single pulse intensities in the "OFF"-pulse window should be entirely due to radiometer noise, while the intensities in the "ON"-pulse window should be the sum of the radiometer noise component (same as the "OFF"-pulse window) and the pulsar emission component. We first generated the average pulse profile to visually select on and off windows of the same widths. We then fit a sixth-order polynomial as a function of pulse phase to each single pulse (similar to Rosen et al. 2013;Lynch et al. 2013;Kaplan et al. 2018) after masking the ON/OFF windows to remove any trends and construct a flat baseline. We recorded the ON/OFF intensities as the sum of the baseline-subtracted intensities across the windows. Finally, we constructed histograms of the ON/OFF intensities which we used to determine the nulling properties. Figure 2 shows the single pulse intensity distribution in the ON/OFF window. The OFF histogram can be accurately described by a single component (Gaussian noise), but the ON histogram can have multiple components -"null" and "emission" components. The presence of nulling manifests in the ON histogram as an excess of samples at levels consistent with the OFF component, which we refer to as the null component. The residual distribution, after removing the null component, is supposed to be a realization of pulsar's emission distribution (hereafter referred to as 'emission' component). The emission component can be a single distribution or a combination of multiple distributions. The ON distribution can be thought of as the sum of the null and the emission components.

Determining Nulling Frations
As demonstrated by Kaplan et al. (2018), Ritchings' method can give biased estimates for NF (hereafter referred as NF r ) in pulsars where the emission component is close to the noise level. Therefore, following Kaplan et al. (2018) we adopt a method which models the ON/OFF histograms using a mixture model (MM). This means that the intensities x can be considered as random draws from the probability density function (PDF) where the F n functions are the individual probability density functions parameterized by the set {θ n }, c n are the weights. In the case where all the F n functions are the same and are normal distributions where {µ n } and {σ n } are the means and standard deviations of component n, this reduces to a Gaussian mixture model (GMM), but more general models are considered.
There is an additional constraint that the weights c n add to one: m n=1 c n = 1, which comes from the normalization of the PDF, which leaves the total number of free parameters to be determined as m n=1 dim({θ n }) model parameters, and m − 1 latent parameters.
In general, the OFF histogram can be well-described by a Gaussian as expected of radiometer noise (assuming that RFI has been sufficiently removed), and this is what we observe in our data. The emission component usually can be described by a single Gaussian as well. However, there are cases when it deviates from a single Gaussian component. More than one component is a possibility considered in Kaplan et al. (2018), which can be tested against the single-component model through a model comparison test. However, we also consider non-Gaussian models here. Specifically, multi-path propagation of the pulses through the interstellar medium (ISM) (Smith 1973;Bhat et al. 2003;Lorimer & Kramer 2004), can result in the emission distribution having long tails towards higher intensities. This effect can be reasonably well described by the intensity distribution which is a convolution of a Gaussian N (x; µ, σ) and a one-sided exponential 1 τ exp(−x/τ )U(x), where U (x) is the Heaviside or step function, erfc(x) is the complementary error function, and τ is the decay time of the exponential (McKinnon 2014). Hence we try to model the emission component using multi-component Gaussians and Gaussians with exponential tails and rank them using their Bayesian Information Criterion (BIC) values to choose the best-fit model.
We employ the scikit-learn Gaussian mixture model (Pedregosa et al. 2011) to derive an initial fit for the ON and OFF histograms. This is based on the expectation-maximization (EM) algorithm, in which parameters are estimated by maximizing the likelihood function L(data |θ) (see Ivezić et al. 2020, for details). This produces a very good fit for the OFF histogram. However, in the case of weaker pulsars where the emission can be confused with the background, Kaplan et al. (2018) showed that this method can still fail in producing a reliable fit for the null and emission components of the ON histogram simultaneously, although this bias can be small compared to the Ritchings' algorithm. As such, a refined fit for the null and emission components can be obtained by performing a Markov-Chain Monte Carlo (MCMC) analysis.
For MCMC analysis, the likelihood function is given by following p(x i |θ) from Equation 1. The priors chosen are: • Initial Gaussian fit from the EM algorithm for the off-pulse mean and standard deviation • Uniform between the bounds dictated by the onpulse intensities for the parameters governing the pulsar emission component We use the emcee (Foreman-Mackey et al. 2013) ensemble sampler to sample the posterior. We initialize 32 walkers within a ±5σ range of the initial fit values of the parameters. To account for the finite correlation length of the chains and produce independent samples, we first let the walkers "burn-in" to erase their starting conditions, and we then let the walkers explore the parameter space until we have at least 100 independent samples for each walker. Figure (3, left column) shows the pulse intensity histograms for PSR J0325+6744: a pulsar in which the emission component is easily discernible from the noise; and PSR J1529−26: a pulsar where these two start to blend into each other. Looking at the null component in the ON histogram for the two pulsars, the evidence for nulling is clear in J0325+6744 while J1529−26 behaves like a non-nulling pulsar whose emission is weak. The blue, green and orange-filled regions show the fit for the OFF, null, and emission components respectively, and the black dotted line shows the overall fit for the ON component. The posteriors for the model parameters are presented in Figure (3, right column) with the point estimates (median 6 ) of the NF from MM given in Table 4.
For PSR J0325+6744, where the null and emission components are well separated (bright pulsars), our method yields a NF = 53.92 ± 0.81% while Ritchings' Note-Naming convention for the model represents the model used to describe the emission histogram (G=Gaussian, Eg=Exponentially modified Gaussian) followed by the number of components in the ON histogram. a We find that in extreme cases (non-nulling/highly-nulling), one of the distributions is confined to very few bins and so we quote this range rather than fitting for it. b We find that there are no single pulses with NP>0.5.
c We observe quasi-periodicity in these cases.
method (see Ritchings 1976;Wang et al. 2007;Kaplan et al. 2018, for implementation) gives a comparable estimate of 55.01%. However in the case of a weaker pulsar, PSR J1529−26, where the emission component is closer to the background noise, our method gives a best-fit value of NF = 5.55 ± 4.4% compared to 48.1% given by the Ritchings' method. The latter is significantly overestimated and can easily lead to (mis)classifying the source as a nulling pulsar, further illuminating the bias of Ritchings' method in weaker pulsars.
Full results for all the 23 pulsars, including the single pulse stacks, posteriors from the MCMC run and the resultant ON/OFF histogram model fits are shown in Appendix A.

Nulling Correlations
After determining the nulling properties we wish to know whether the locations and durations of nulls are completely random, or if there is any correlation between different nulling and emission episodes in a pulsar. Specifically, given a single pulse that shows emission (or that nulls), how likely are we to see emission for the next pulse, and are there any patterns of longer duration?
We test this using the probability of a null (the nulling "responsibility") evaluated for each individual pulse, given by We divided the data into stacks of 256 pulses (similar to Ritchings 1976; Herfindal & Rankin 2009) to calculate more robust estimates and to be less sensitive to longterm variations like scintillation and system temperature changes, and use equation 4 to calculate the probability of a given single pulse being a null. We then looked for periodic signature by taking the Fourier transform (FT) within each stack and co-adding the power from all stacks incoherently. Figure 4 shows the resultant spectrum for PSR J0414+31, in which a certain pattern of combination of emission and nulls seems to be periodic over ∼28 pulse periods. We estimate the significance of peaks in the stacked power spectra assuming that the null distribution from n stacks follows a χ 2 distribution with 2n degrees of freedom (this assumes white noise). We see significant periodic or quasi-periodic (a significant broad peak in the power spectrum) signatures in a few other pulsars, and tabulate their periods in Table 4. In the case of precise period measurements, we estimate the uncertainty as described in Ransom et al. (2002). However, this only points to the periodic nature of a certain pattern of emission and nulls. To find how emissions and nulls are 'bunched', we look for the distribution of continuous emissions and nulls, where we use NP I =0.5 to be the boundary between an emission and a null. Figure 5 shows the emission and null length distri-  butions for the single pulses of PSR J0414+31. We find that these distributions can be well described by an exponential distribution (p(x) = τ −1 exp(−x/λ)), where x is the null or emission length and the mean duration of the episode is λ. We find that for PSR J0414+31, the emission episodes have a characteristic period of four periods, whereas the nulls are two periods long, which is consistent with the observed nulling fraction of ∼ 33% (see Table 4). We repeat this for all the pulsars and the results are tabulated in Table 4.

Sub-pulse Drifting
Beyond nulling, we also look for any correlations between nulling and sub-pulse drifting. Drifting is usually characterized by two periods: the drifting period P 3 , defined as the period for which the pulse is seen at the same longitude (phase), and P 2 , the spacing between two sub-pluses within the same single pulse (see Figure  6). To estimate both, we prepared the data by selecting only the on-pulse window of data (n p phase bins) for all the single pulses (n s single pulses). We then calculated Longitude Resolved Fluctuation Spectra (LRFS, Backer 1970c), where we take a 1-D Fourier transform of the (n s × n p ) data along the n s axis. Figure 6 shows one of the two pulsars in our sample, J1822+02, that shows clear signs of drifting. A period P 3 of ∼ 28 pulse periods and P 2 of ∼ 35/1024 pulse periods can be clearly seen. We also find the evidence for drifting in PSR J1829+25 (see figure 7), with a P 3 of ∼ three pulse periods and a P 2 of 1/128 pulse periods, with similar inferences in the data from both AO and GBT.

Kaplan et al. (2018) demonstrated the bias of Ritch-
ings' method for weaker pulsars through simulated data, where the mixture model was able to recover the true injected nulling fraction. They also showed that for Gaussian mixtures, an analytical correction can correct the biased estimate of Ritchings' method to find the true value. We extend the same technique using our sample of 22 pulsars. Figure 8 shows the comparison of the NF estimates derived using both methods. The blue points show the NF estimate derived using Ritchings' algorithm (NF r ), the orange points show NF r estimate corrected for the bias (as in Kaplan et al. 2018), and the green points show the NF derived using mixture modeling. In the case of highly nulling pulsars, the contamination of the null component from the emission component can be small, and both methods perform comparably. However, in the case of pulsars with small NF a systematic bias can be seen as the pulsar emission component becomes blended with the background noise, and the fact that the green and orange points agree quite well demonstrates our confidence in estimating the bias in the Ritchings method and the utility of mixture models.

Is the Nulling Fraction Correlated with Pulsar
Properties?
Comparing the nulling estimates from the mixture modelling and Ritchings' method in Table 4, it can be seen that there can be significant differences between these estimates. Such a scenario can lead to significant biases in population-wide studies that look for correlation between nulling fraction and pulsar properties. Figure 9 shows the most complete list of nulling pulsars, extended from Konar & Deka (2019), on the P −Ṗ diagram. We do not find any clear visual trends of NF with respect to period (P ), spin-down rate (Ṗ ), characteristic age (τ c ), or surface magnetic field (B surf ), although we emphasize that most of the pulsars here (142/164) have their NF estimates derived using some variant of the Ritchings method.
Our sample size of 22 pulsars is too small to derive reliable correlations. However, we can test the similarity/disparity in the correlations obtained using nulling estimates derived with mixture models versus the Ritchings algorihtm. We use the Spearman correlation test, a non-parametric correlation test to quantify any correlations between the relevant parameters (P /Ṗ /B surf /τ c ) and NF. Table 5 shows the correlation coefficients of nulling fraction with parameters of interest (P ,Ṗ , B surf , τ c ). In no case do we see an evidence for strong correlations but we can see large differences between these coefficients obtained using the NF derived using the two methods. We emphasize that the values of these have to be taken with a high degree of caution, given the relative sample size under study and the presence of outliers. In particular we find that PSR J2310+6706 turns out to be a strong outlier, especially in the τ c and B surf space and this significantly affects the results (see Table 5), further illustrating the limitations of a small sample size. Previously, using a sample size (23) comparable to ours, Wang et al. (2007) qualitatively found that NF is related to age with older population experiencing larger nulling fractions. Ritchings (1976) found a positive correlation both with the pulsar period and age in a sample (32) slightly larger than the one in this study. However, as mentioned above those and most other nulling estimates in the literature are derived using some variant of Ritchings' algorithm. Computing the Spearman coefficient for all of the archival sources we cannot confirm either correlation and suggest caution in interpreting results using Ritchings' algorithm.
However, we also note that the source of this disparity does not seem to be straightforward: For a sample of pulsars with a given SNR, the energy per single pulse will be lower for pulsars with shorter periods, which means that the NF estimates for the shortperiod pulsars should experience larger biases and have higher nulling fractions measured with the Richtings' method. Under the (overly simplistic) assumption of a uniform distribution of luminosity with period (cf. Faucher-Giguère & Kaspi 2006;Bates et al. 2014), the correlation of inferred nulling fraction with period will then be negative which is contrary to the previous studies. This suggests that the source of this bias is not simple and needs careful understanding of the under- Note-Not all the pulsars in the sample havė P measurements. Hence the sample size used for period is larger. The two rows for each parameter correspond to the rank coefficients including and excluding PSR J2310+6706 (see Figure 10).
lying distribution of NF with pulsar properties and a larger sample of pulsars with more robust and unbiased NF estimates.

Is Nulling Periodic?
As shown in Section 3.2, we find that nulling appears periodic/quasi-periodic in a subset of pulsars, with their periods noted in Table 4. Herfindal & Rankin (2007, 2009) also find evidence for such signatures and attributd this to the line of sight passing through a structured rotating carousel. In addition we also find that in PSR J0414+31, which was observed at two different frequencies with different instruments, this period is the same. It should be noted that the frequency resolution here is ∼ 0.004 pulse period −1 (from the stacks of 256 pulses) and so we will be insensitive to any changes that are finer than this. Although significant correlations can not be drawn from these periodicities given our sample size and the number of pulsars that show periodic nulling, the occurrence of such a phenomenon in modest set of pulsars in our sample suggests that this might not be uncommon and should be searched for in future data.

CONCLUSIONS
In this study, we have extended the Gaussian mixture model of Kaplan et al. (2018) to study nulling behavior in 22 pulsars, spanning a wider range of properties than in the initial paper but still not selected independent of nulling behavior. We find that all pulsars can be well-represented by mixture model, but we find that a single Gaussian is not sufficient to describe the emis- Nulling Fraction Figure 9. Period-period derivative (P −Ṗ ) diagram highlighting nulling pulsars. Shown in grey circles are all the pulsar from the ATNF catalog (Manchester et al. 2005), in colored circles are the archival nulling pulsars from Konar & Deka (2019) and in diamonds are the pulsars from this study. The contours represent lines of constant characteristic age τc and dipolar surface magnetic field (B surf ). The color bar shows the nulling fraction which ranges from 0 to 1. No clear discernible trend of NF with any of P /Ṗ /B surf /τc is visible.
sion component in some pulsars 7 . Similar to Kaplan et al. (2018), we find that previous methods used to estimate NF can suffer significant biases when the pulsar emission is weak compared to the background noise. Such biases may lead to misinterpreting weak pulsars as nulling pulsars. We also show that these biases may lead to spurious correlations between the NF and pulsar properties in population-wide studies. Drawing on the more robust statistics that we calculate, we find that nulling can appear periodic, with three pulsars in our sample showing this behavior. Two pulsars in our sample, PSR J1822+02 and PSR J1829+25, shows clear signs of sub-pulse drifting, and they have an inferred nulling fraction consistent with 0. In contrast, studies like Gajjar et al. (2014a); Davies et al. (1984) find sub-pulse drifting in pulsars that exhibit moderate nulling, indicating that sub-pulse drifting and nulling might be two independent manifestations of sub-pulse intensity variations. In all cases we look forward to using larger, less-biased samples to more robustly explore the nulling population and seeing if it is related to other phenomenology.
Two pulsars in our sample, PSR J0414+31 and PSR J1829+25, were observed at two different frequencies (430 MHz and 820 MHz), albeit not simultaneously. PSR J1829+25 has nulling estimates that agree at both frequencies, consistent with 0, but we find that PSR J0414+31, has NF estimates in tension at the ∼ 2σ level, with the NF higher at lower frequencies. Although it is hard to draw definite conclusions from these two pulsars since the observations are not simultaneous, it emphasizes the need for simultaneous observations at multiple frequencies (or across a larger bandwidth). Observing at 4 different frequencies (325, 610, 1400, 4850 MHz), Gajjar et al. (2014a) find coherent nulling in three different pulsars whereas Bhat et al. (2007) find the evidence for null excess at lower frequencies in PSR B1133+16 further emphasizing the need for multi-frequency observations in a larger sample to find whether nulling is universally broadband.
One of the pulsars in our sample (PSR J2310+6706) has a two-component profile with a faint leading peak in addition to the primary peak. The very low SNR of the leading component limits our ability to find a stringent estimate of the NF independent of the primary component, but we find that the NF values obtained from each component is consistent. Analyzing nulling characteristics in pulsars with multi-component pulse profiles  Figure 10. Scatter plot showing the NF of the pulsars in this study vs their properties. It can be seen that the pulsars appear scattered in the P /Ṗ space. However, with the exclusion of PSR J2310+6706 which appears as an outlier in the τc/B surf space, a rough trend can be seen that of NF decreasing with the age τc and increasing with the surface magnetic field B surf . The correlation coefficients are given in Table 5.
with a robust method like mixture modeling can provide insights into the simultaneous nulling in the different regions of the pulsar's magnetosphere.
So far we have only analyzed normal, non-recycled pulsars. Current sensitivity limitations restrict the sample of nulling pulsars to normal pulsars (as is evident from Figure 9), while MSPs are largely unexplored. Initial single pulse studies done by Rajwade et al. (2014) do not find any compelling evidence for nulling in MSPs. Using the mixture model technique, which does not suffer from the same biases at low signal-to-noise, for MSPs, together with newer higher-sensitivity facilities may help explore whether the nulling phenomenon affects all pulsars, or is limited to a sub-population.
We thank an anonymous referee for helpful suggestions that clarified this work. AA, JS, and DK receive support from National Science Foundation (NSF) Physics Frontiers Center award numbers 1430284 and 2020265. AA thanks Alex McEwen for helpful discussions during the data reduction stage. The Arecibo Observatory is a facility of the NSF operated under cooperative agreement (#AST-1744119) by the University of Central Florida (UCF) in alliance with Universidad Ana G. Méndez (UAGM) and Yang Enterprises (YEI), Inc. The Green Bank Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc.