A Statistical Analysis of the Nulling Pulsar Population
MMNRAS , 1–12 (2020) Preprint 2 February 2021 Compiled using MNRAS L A TEX style file v3.0
A Statistical Analysis of the Nulling Pulsar Population
Sofia Z. Sheikh, , ★ Mariah G. MacDonald, , Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA Center for Exoplanets and Habitable Worlds, The Pennsylvania State University, University Park, PA 16802, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Approximately 8% of the ∼ nulling fraction , which describes the percentage of time a given pulsar will be found in a nullstate. In this paper, we perform the most thorough statistical analysis thus far of the properties of141 known nulling pulsars. We find weak, non-linear correlations between nulling fraction andpulse width, as well as nulling fraction and spin period which could be attributed to selectioneffects. We also further investigate a recently-hypothesized gap at 40% nulling fraction. Whilea local minimum does exist in the distribution, we cannot confirm a consistent and unique breakin the distribution when we investigate with univariate and multivariate clustering methods,nor can we prove the existence of two statistically distinct populations about this minimum.Using the same methods, we find that nulling pulsars are a statistically different populationfrom normal, radio, non-nulling pulsars, which has never been quantitatively verified. Inaddition, we summarize the findings of the prior nulling pulsar statistics literature, which arenotoriously contradictory. This study, in context, furthers the idea that nulling fraction alonedoes not contain enough information to describe the behaviour of a nulling pulsar and that otherparameters such as null lengths and null randomness, in addition to a better understanding ofselection effects, are required to fully understand this phenomenon. Key words: radio pulsars — astrostatistics — galactic radio sources
Pulsar “nulling” (first reported in Backer 1970) is a relatively-common phenomenon where otherwise strong and steady radiationfrom a pulsar decreases by at least a factor of 100, and possiblyentirely, over all frequencies (Ritchings 1976; Wang et al. 2007a).These nulls occur for durations of a single pulse period all the wayup to days, and affect all components of a pulse simultaneously,(e.g. both core and conal components (Rankin 1986), interpulses(Hankins 1984)).Nulling pulsars can be identified by looking for gaps in emis-sion in phase-time diagrams (Wang et al. 2007a). The nulling frac-tions (NFs) can then be measured with pulse energy distributions,which will have an overabundance of pulses near zero energy, some-times leading to a separate well-resolved peak in pulsars with largeNF. These overabundances can be measured and reported as a per-centage of the total number of pulses that are nulled (Ritchings1976).There are some caveats to this method. Firstly, nulling fractionis only one of the parameters that characterizes nulling behaviour.Other parameters, including the duration of a nulling episode and ★ E-mail: [email protected] (SZS) Factors up to 1420 have been found (e.g., B1706-16, Naidu et al. 2018) distribution of these durations, are needed to fully describe thenulling behaviour of a pulsar (Gajjar et al. 2012). Secondly, sensi-tivity limits make it difficult to identify whether pulsars are entirelydormant during a null, or whether there is some residual low-levelemission (Wang et al. 2007a).Some nulling pulsars are reported to show changes in theirenergetics before, after, or during a pulse (e.g. an increase in pulseenergy after a null or a slow-down rate that drops during a null,Backer 1970; Kramer et al. 2006), suggesting that nulls, like otherphenomena such as mode-changing , occur due to changes in a pul-sar’s entire magnetospheric environment (Lyne & Ashworth 1982;Lyne et al. 2010; Gajjar et al. 2014a). This, when considered withthe simultaneous broadband cessation during a null, implies a non-uniform flow of particles into the emission region or a global break-down of the emission mechanism (in a “tipping point” like insta-bility, Wang et al. 2007a; Gajjar et al. 2014b), and not changes inrelative geometry.Determining correlations between nulls and other pulsar prop-erties and behaviour allows for characterization and investigationsof the pulsar emission mechanism and its evolution over time. Pul- The likely instability-induced switching from one primary pulse profileto other secondary ones (Wang et al. 2007a) © a r X i v : . [ a s t r o - ph . H E ] J a n Sheikh and MacDonald sar characteristic age, pulse morphology, and pulse periods havebeen closely investigated in relation to NF, but there is currently noconsensus about their impact on NF. Some studies have reporteda positive correlation between characteristic age and nulling frac-tion (Ritchings 1976; Wang et al. 2007a), interpreting this as a signthat pulsars null more and more as they age until they “die” andstop emitting entirely (Ritchings 1976). Other studies find no cor-relation with age at all (Biggs 1992), or no correlation with agewithin a morphological group (Rankin 1986). Using classes fromRankin (1983a), Rankin (1986) suggested that different pulse mor-phologies lead to different NFs (with triple and multiple profileshaving the highest NFs). However, even using the same classes,Wang et al. (2007a) did not find such a correlation and suggestedthat the known relationship between complex pulse morphologyand older age (Huguenin et al. 1971) could have caused the previ-ously observed correlated behaviour. Suggestions of a correlationbetween NF and period (Ritchings 1976; Biggs 1992) are hard toentangle from the previous discussion of correlations with age andmorphology, and are also debated in the literature (Rankin 1986).The statistical methods used in the majority of this previousliterature have been preliminary at best. Most studies graphed rel-evant variables against each other and looked for trends via visualinspection, though they were working with much less data and oftenadmitted the limitations to claiming a correlation. Rankin (1986)calculated means and standard deviations of different pulsar popu-lations to claim that they were distinct. Some studies, specificallyBiggs (1992) and Li & Wang (1995), attempted to account for theupper limits and uncertainties on the nulling fraction measurementsby using survival analysis methods such as Cox’s proportion hazardmodel (Cox 1972) and Kendall’s rank correlations (Kendall 1938).Most recently, Konar & Deka (2019) quantified correlations withnulling fraction using Pearson’s R.Konar & Deka (2019) provided the results of a thorough liter-ature survey for NF values and the most recent statistical survey ofnulling pulsar properties. The authors tentatively propose a gap inNF values around 40% and claim that a Kolmogorov-Smirnov teston the spin-period of nulling pulsars with NF values greater thanand less than 40% imply that there are two distinct populations ofpulsars (Konar & Deka 2019). The authors also did not find corre-lations between NF and any intrinsic pulsar parameters, in contrastwith previous studies.In this paper, we will expand upon the table of pulsar param-eters provided by Konar & Deka (2019), investigate correlationsbetween intrinsic pulsar parameters and null fractions, compare thenulling and non-nulling pulsar populations, and statistically char-acterize and attempt to confirm the two populations around thesuggested 40% NF gap from Konar & Deka (2019). In Section 2,we explain how we compile and prepare our data and define allof the parameters that will be used in the statistical analysis de-scribed in Section 3. We then explore the existence of the twodistinct populations in null fraction suggested by Konar & Deka(2019) in Section 4, and look at trends between nulling fraction andmorphological class in Section 5. We provide the results from thecorrelation analyses between the nulling fraction and between otherparameters in the nulling pulsar population in Section 6. In Section7, we look at the differences between the nulling and non-nullingpulsar populations before discussing our results in Section 8 andconcluding in Section 9.
The data used in the analyses in Sections 4 – 7 are comprised of15 pulsar parameters and derived quantities. We show 14 of the 15parameters used in the study in Table 1. We discuss the remainingparameter, morphology, below. We gathered these parameters fora set of 141 pulsars contained in Tables 1–5 of Konar & Deka(2019), excluding pulsars in their sample without a measured nullingfraction. Pulsars with 2 or 3 measurements of nulling fraction in theliterature have 2 or 3 rows in the table, respectively, for a total of162 rows . In our analyses, we use all 162 measurements of nullingfraction, as if each measurement were from a separate pulsar. Table1 includes a description of each parameter and an explanation ofhow it was collected.In addition to the parameters shown in Table 1, we compileRankin (1983a) pulse morphology classes for 141 of the 162 pulsarentries in our sample. 120 of these have been previously classified,and we assigned classifications to 22 additional pulsar entries basedon published pulse profiles and analyses. These classifications are:S 𝑑 (conal-component single-peaked profile), S 𝑡 (core-componentsingle-peaked profile), D (double-peaked profile), T (triple-peakedprofile), and M (multiple profile containing more than three peaks instructure). We added an additional letter U (“undetermined”) for theremaining 20 pulsars and did not consider them in the morphologyanalysis in Section 5. Cluster analyses are a set of methods which can determine howmany distinct subsets lie within a single dataset, which points be-long to each of the subsets, and whether the subsets themselvesare statistically distinct. We employ multiple different clusteringtechniques, described below.
A kernel density estimation (KDE) is a nonparametric estimationof probability density functions, defined by:ˆ 𝑓 ( 𝑥 ; 𝐻 ) = 𝑛 − 𝑛 ∑︁ 𝑖 = 𝐾 𝐻 ( 𝑥 − 𝑋 𝑖 ) (1)where 𝑥 = ( 𝑥 , 𝑥 ) 𝑇 and 𝑋 𝑖 = ( 𝑋 𝑖 , 𝑋 𝑖 ) 𝑇 for 𝑖 =1,2,...,n a randomsample drawn from a density 𝑓 , 𝐾 ( 𝑥 ) is the kernel which is a sym-metric probability density function, and 𝐻 is the positive-definiteand symmetric bandwidth matrix. We use an axis-aligned normalkernel. When estimating the location of the gap in nulling fraction,we use a variety of bandwidths from U[1.0,100.0] (see Section 4for full analysis). Jenks’ Natural Breaks Optimization (JNBO) is a univariate classinterval method originally developed for the classifying geographicdata for cartography (Jenks 1977). Now it is used in many fields There are actually 163 individual estimates of nulling fractions, but theestimate for B1713-40 from Wang et al. (2007b) was accidentally omittedin our analyses. MNRAS000
A kernel density estimation (KDE) is a nonparametric estimationof probability density functions, defined by:ˆ 𝑓 ( 𝑥 ; 𝐻 ) = 𝑛 − 𝑛 ∑︁ 𝑖 = 𝐾 𝐻 ( 𝑥 − 𝑋 𝑖 ) (1)where 𝑥 = ( 𝑥 , 𝑥 ) 𝑇 and 𝑋 𝑖 = ( 𝑋 𝑖 , 𝑋 𝑖 ) 𝑇 for 𝑖 =1,2,...,n a randomsample drawn from a density 𝑓 , 𝐾 ( 𝑥 ) is the kernel which is a sym-metric probability density function, and 𝐻 is the positive-definiteand symmetric bandwidth matrix. We use an axis-aligned normalkernel. When estimating the location of the gap in nulling fraction,we use a variety of bandwidths from U[1.0,100.0] (see Section 4for full analysis). Jenks’ Natural Breaks Optimization (JNBO) is a univariate classinterval method originally developed for the classifying geographicdata for cartography (Jenks 1977). Now it is used in many fields There are actually 163 individual estimates of nulling fractions, but theestimate for B1713-40 from Wang et al. (2007b) was accidentally omittedin our analyses. MNRAS000 , 1–12 (2020) ulling Pulsar Statistics Name Symbol Units Equation Data Reference DescriptionPulse Period 𝑃 s N/A Konar & Deka (2019),from Manchester et al.(2005) Barycentric spin periodSpindown Rate (cid:164) 𝑃 s/s N/A Manchester et al. (2005) First time derivative ofpulse periodMagnetic Flux 𝐵 𝑠 G 𝐵 𝑠 ∝ ( 𝑃 (cid:164) 𝑃 ) / Konar & Deka (2019),from Manchester et al.(2005), Ostriker & Gunn(1969)* Surface magnetic flux den-sityDispersion Measure DM cm / pc N/A Manchester et al. (2005) Indirect measure of pulsardistance from frequency-dependent arrival time de-lay (e.g., Jokipii & Lerche1969)Characteristic Age 𝜏 𝑐 yr 𝜏 𝑐 ∼ 𝑃 (cid:164) 𝑃 Existing table values Measure of pulsar age fromperiod and spindown rateKinetic Age 𝜏 𝑘 yr 𝜏 𝑘 = 𝑧𝑣 𝑧 Derived with 𝑧 and 𝑣 𝑧 from Manchester et al.(2005), using AstroPy (Robitaille et al. 2013) Period-independent pulsarage derived from pulsar’sdistance from the galacticplane and z-velocity awayfrom the plane † Radio Luminosity 𝐿 𝑅 mJy/kpc 𝐿 𝑅 ∝ 𝑆 𝑑 Derived with 𝑑 and 𝑆 from Manchester et al.(2005) Radio luminosity at 400MHz, quantifies intrinsicpulsar powerSpindown Luminosity (cid:164) 𝐸 mJy/kpc (cid:164) 𝐸 ∝ (cid:164) 𝑃𝑃 − Manchester et al. (2005),Gold (1969)* Spindown energy loss rateMaximum MagneticFlux Density 𝐵 𝑙𝑐 G 𝐵 𝑙𝑐 ∝ 𝑃 (cid:164) 𝑃 − Existing table values, Gol-dreich (1969); Lyne et al.(1975)* Magnetic field in the “lightcylinder” region of the pul-sar magnetosphere, whereparticle velocities are nearcPlasma Flow Parameter 𝑄 None 𝑄 ∝ 𝑃 . (cid:164) 𝑃 − . Existing table values, Be-skin et al. (1984)* Derived parameter empir-ically related to pulsarmorphology, with
𝑄 << 𝑊 ms N/A Manchester et al. (2005) Width of the average pulseprofile measured at 50%and 10% of the peakbrightnessMagnetic InclinationAngle 𝛼 ◦ N/A Malov (1990); Malov &Nikitina (2011); Nikitina& Malov (2017) ‡ Angle between spin axisof the pulsar and magneticmomentNulling Fraction NF % N/A Tables 1–5 of Konar &Deka (2019) Percentage of time that apulsar is in a nulled state
Table 1.
14 of the 15 pulsar parameters used in the following statistical analyses. Equations are only given for derived parameters. References marked with *give the equations used for derived parameters. † These data were only available for 96 pulsars in the sample. It should be noted that pulse widths change significantly when measured at differentfrequencies, or with different time resolutions, and that any use of these data must take into account the large additional uncertainties introduced to a potentialcorrelation. ‡ We compiled 𝛼 values (reported as 𝛽 in some references) from the literature for 107 pulsars in our sample. As noted by Nikitina & Malov (2017), 𝛼 can be calculated from pulse widths (assuming that the line of sight passes through the center of the radiation cone as in Rankin 1990), polarizationinformation from the main cone, or position angle and mean profile shape, as in Lyne & Manchester (1988). The reported 𝛼 values used here are averagesfrom the different methods, where possible. It can be difficult to measure 𝛼 values with the pulse width method for pulsars with extremely narrow pulse widths(e.g., J1717-4054, Kerr et al. 2014). beyond geography to find a user-specified number of homogeneousclassification groups in one-dimensional data (e.g., Rahadiantoet al. 2015). Once the number of classes has been defined, the JNBOalgorithm tries to iteratively minimize the sum of the variances fromthe class means (Coulson 1987). Jenks-Fisher is a more runtime- efficient version of this same algorithm, incorporating the methodput forth in Fisher (1958) (Bivand 2019). MNRAS , 1–12 (2020)
Sheikh and MacDonald
Hierarchical clustering is a clustering method which does not imme-diately produce a user-specified number of clusters. Agglomerativealgorithms such as hclust in R instead look for the two “clos-est” points (by some linkage metric) and then combine them into acluster, iteratively repeating until all points are incorporated into asingle cluster (Leisch 1999; Müllner et al. 2013). The results for anynumber of clusters are thus pre-calculated. The main drawback ofhierarchical clustering is its large computational cost, but its abil-ity to use different linkage metrics makes it a particularly flexiblesolution for clustering (Leisch 1999).
Bagged (“ b ootstrap agg regat ed ”) clustering is a clustering methodthat is both a partitioning method, like k-means clustering (Mac-Queen et al. 1967), and a hierarchical method (as described inSection 3.1.3) (Leisch 1999). Partitioning methods such as k-meansare less effective at recovering structure if the clusters are not convexand, as mentioned in the previous section, hierarchical clusteringcan be computationally intensive; bagged clustering algorithms suchas bclust in R avoid both drawbacks (Leisch 1999). Geometrical interval classification classifies data into intervals thathave a geometric series. The algorithm forms these intervals byminimizing the sum of the squares of the number of points in eachclass. Interval classifications also aim for each class range to beapproximately the same size and for the change between intervals tobe consistent (Jianhua 2002). This algorithm is designed to accom-modate continuous data and looks to highlight changes in the dataat extreme values as well as at moderate values, making it ideal forheavily skewed and non-normal distributions.
The Kolmogorov-Smirnov (K-S) test is a nonparametric test of one-dimensional probability distributions that can be used to comparetwo samples. The K-S statistic quantifies the distance between theempirical distribution functions of the two samples and is 𝐷 𝑛,𝑚 = 𝑠𝑢𝑝 | 𝐹 ,𝑛 ( 𝑥 ) − 𝐹 ,𝑚 ( 𝑥 )| (2)where 𝐹 ,𝑛 and 𝐹 ,𝑚 are the empirical distribution functions of thetwo samples, and 𝑠𝑢𝑝 is the supremum function. The null hypothesisthat the two samples are from the same underlying population isrejected at a level 𝛼 if 𝐷 𝑛,𝑚 > 𝑐 ( 𝛼 ) √︃ 𝑛 + 𝑚𝑛𝑚 ,where 𝑐 ( 𝛼 ) = √︂ − ln (cid:16) 𝛼 (cid:17) Astronomers have been using the K-S test, for two-samples aswell as for goodness-of-fit, for decades now (e.g., Peacock 1983).This is mainly because it easy to compute and quite simple to un-derstand, but the test is also distribution-free and can be universallyapplied without restrictions to sample size. Unfortunately, the K-Stest is insensitive when the differences between the empirical dis-tribution functions are not in the middle, but rather at the beginningand the end. The Anderson-Darling (A-D) test oftentimes tests for normal-ity, or compares a sample’s empirical distribution function to that ofanother known form (e.g. log-normal, exponential). Using the A-Dmeasure of agreement between distributions, though, it is possibleto nonparametrically test a sample, or multiple samples, to see ifthey come from the same distribution.The test is based on the distance between the two empiricaldistribution functions: 𝐴 = 𝑛 ∫ ∞−∞ ( 𝐹 𝑛 ( 𝑥 ) − 𝐹 ( 𝑥 )) 𝐹 ( 𝑥 )( − 𝐹 ( 𝑥 )) 𝑑𝐹 ( 𝑥 ) (3)Like the K-S test, the A-D test is both nonparametric anddistribution free, but it does not suffer from the same insensitivityissues as the K-S test since it compares the two functions over theentire interval. Also like the K-S test, the null hypothesis states thatthe two tested samples are from the same underlying population. Incomparing two populations, we employ both K-S ad A-D tests. The Pearson correlation coefficient, also known as “Pearson’s R” orbivariate correlation, measures the linear correlation between twodifferent variables (Pearson 1895). The coefficient ranges between-1 and 1, where 0 indicates no correlation, 0.3–0.5 indicates a weakcorrelation, 0.5–0.7 indicates a moderate correlation, and greaterthan 0.7 indicates a strong correlation; negative values of the sameranges indicate anti-correlations. The Pearson correlation coeffi-cient is used widely across scientific fields when analyzing data.The coefficient is the covariance of the two variables being tested,divided by the product of their standard deviations. Since we oftendo not know the entire population and merely study samples of thepopulation, the coefficient is of the sample and calculated by: 𝑟 = 𝑟 𝑥𝑦 = Σ 𝑥 𝑖 𝑦 𝑖 − 𝑛 ¯ 𝑥 ¯ 𝑦 ( 𝑛 − ) 𝑠 𝑥 𝑠 𝑦 (4)where 𝑛 is the number of data points in the sample, 𝑥 𝑖 and 𝑦 𝑖 areeach data point, ¯ 𝑥 and ¯ 𝑦 are the sample means for x and y, and 𝑠 𝑥 and 𝑠 𝑦 are the sample standard deviations for x and y.Although Pearson’s correlation coefficient is widely used, thecoefficient only measures how linear a correlation is; it is thereforesometimes insufficient since two variables can be highly corre-lated without their correlation being linear. Many correlations inastronomy are power laws or exponential, and Pearson’s correlationcoefficient is insensitive in these situations. As such, we also ana-lyze the data by looking at the Maximal Information Criteria (MIC,Reshef et al. 2011). Similar to Pearson’s correlation coefficient, theMIC ranges between 0 and 1, with 0 indicating no correlation, andhigher values corresponding with stronger correlations. The MIC isdefined as:MIC(D) = max 𝑋𝑌 < 𝐵 ( 𝑛 ) 𝐼 ∗ ( 𝐷, 𝑋, 𝑌 ) log ( min ( 𝑋, 𝑌 )) (5)where 𝐷 = 𝐷 ( 𝑥, 𝑦 ) is the set of 𝑛 ordered pairs of elements of 𝑥 and 𝑦 in a space partitioned by an 𝑋 × 𝑌 grid, grouping the 𝑥 and 𝑦 values in 𝑋 and 𝑌 bins, respectively, 𝐵 ( 𝑛 ) = 𝑛 𝛼 is the search-grid size, with the 𝛼 parameter defined by the sample size (here, 𝛼 = .
75 Albanese et al. 2018), and 𝐼 ∗ ( 𝐷, 𝑋, 𝑌 ) is the maximummutual information over all grids 𝑋 × 𝑌 of the distribution inducedby 𝐷 on a grid having 𝑋 and 𝑌 bins.The MIC and other maximal information-based nonparametric MNRAS000
75 Albanese et al. 2018), and 𝐼 ∗ ( 𝐷, 𝑋, 𝑌 ) is the maximummutual information over all grids 𝑋 × 𝑌 of the distribution inducedby 𝐷 on a grid having 𝑋 and 𝑌 bins.The MIC and other maximal information-based nonparametric MNRAS000 , 1–12 (2020) ulling Pulsar Statistics exploration methods were introduced by Reshef et al. (2011), whoshowed that the MIC for a random relationship was approximately0.18, with anything near or below that number indicating no corre-lation. For this work, we do not use an MIC cut-off, but instead relyon p-values generated as described in Section 6.1. NF is often reported in the literature with associated upper andlower limits. In addition, some pulsars have only been seen to nullonce, leading to a NF that is only a moving upper limit until thepulsar is observed to null again. These uncertainties and upperlimits must be addressed to ensure that any discovered correlationsare statistically significant, so we account for them in the followingmanner. For pulsars with an uncertainty and a nominal value, wedraw NF values for each pulsar from a normal distribution centeredon the nominal value, with a standard deviation of the publisheduncertainty. If the published NF is an upper limit, we instead drawthe NF from a log uniform distribution of U[0, 𝑁𝐹 𝑢 ], where 𝑁𝐹 𝑢 isthe upper limit of the NF. We chose the log uniform distribution toaccount for the skewness of the NF distribution. We test the uniformdistribution as well to ensure that the choice of distribution is notsignificantly affecting the results, as discussed in Section 6.1.We produce an uncertainty-adjusted value for each NF as de-scribed above. When exploring correlations between NF and otherparameters, we produce 1000 of these sets of corrected NFs andcalculate the median of the correlation measurements from eachof these 1000 sets. To evaluate the significance of that median, wepermute our data as explained in Section 6.1. ∼ Recently, Konar & Deka (2019) suggested a gap in nulling fraction(NF) at ∼ classInt (Bivand 2019) to run a series of univariate class interval algorithmswith two classes on 1000 sets of corrected NF (see Section 3.1 for adescription of the methods). Hierarchical clustering, Jenks-Fisher,bagged clustering, and Jenks natural breaks opimization methodsall identified the natural breakpoint in the data at 47 . ± . . ± . . ± . . ± . cartography , did not identify a break orgap near 40%. In fact, the algorithm returned that all 1000 sets ofcorrected NF were continuous distributions, suggesting that the gapthat we see and the breaks that are identified by the other methodsare artifacts of an under-sampled population.In order to directly compare to Konar & Deka (2019), we alsotest this method on the uncorrected NF, recovering gaps at 37% for All statistics in this paper were performed using R 3.4.1 (R Core Team2019). Nulling Fraction I nde x ( S o r t ed ) Figure 1.
Results of Jenks Natural Breaks Optimization (Sections 3.1.2)with two classes run on the uncorrected nulling fraction distribution of 162pulsar entries. All methods identified the class breakpoint to be in the gapbetween 30% and 45% for the uncorrected NF distribution, as well as for1000 sets of corrected NF, in accordance with Konar & Deka (2019). Baggedclustering and hierarchical clustering also identified a secondary breakpointat 58%. Geometrical interval classification, which does not assume a gapas the other methods do, concludes that both the corrected and uncorrectedNF distributions are continuous and absent of true breaks. This suggests thatthe gap seen at 40% is an artifact of small sample size. each method, except geometric interval classification which againdoes not identify any breaks in the data.We then estimate the distribution of NF using a 1D KernelDensity Estimator (KDE), using the uncorrected NF distribution.Using 100 bandwidths drawn from a uniform distribution U[1.0,100.0] , we confirm a minimum in the KDE of 37 . ± .
19. Wedraw multiple bandwidths such that the location of the minimum isnot dependent on our bandwidth.We then perform a KDE on 16000 samples of the correctedNF, again using 100 bandwidths drawn from a uniform distributionU[1.0, 100.0]. Regardless of bandwidth, there does exist a minimumin the KDE around 37%, but there also exist minima at ∼
20% andat ∼
60% (see Figure 2). We thus find that the gap seen by Konar& Deka (2019) is not a boundary between two distinct populations,but is simply an artifact of the uncertainties on NF estimates andan under-sampled population. Once uncertainties are accounted for,there does appear to be a local minimum at around 40%, but otherminima exist, suggesting that this minimum is solely due to thesmall sample size.
We run both Kolmogorov-Smirnov and Anderson-Darling two-sample tests (see Section 3.2) on 1000 sets of corrected nullingfraction measurements. We define our two populations as pulsarswith
𝑁𝐹 < .
5% and pulsars with
𝑁𝐹 > . (cid:164) 𝑃 , 𝐵 𝑠 , DM, log 𝜏 𝑐 , log 𝜏 𝐾 , 𝐿 𝑅 , (cid:164) 𝐸 , 𝐵 𝑙𝑐 , 𝑄 , W50, W10, and 𝛼 ) and report the median p-valuesfrom the 1000 sets in Table 2.Given the large p-values resulting from both of these tests, we Bandwidths drawn from U[1.0, 100.0] span the entire range of NFMNRAS , 1–12 (2020)
Sheikh and MacDonald . . . . . . Nulling Fraction d N Figure 2.
KDE of 16000 samples of corrected nulling fraction, with abandwidth of 15.0. There does exist a local minimum in the distributionat around 40%, as noted by Konar & Deka (2019), but this minimum isaccompanied by two other local minima, at 20% and at 60%. This suggeststhat the gap seen in the data is due to small sample size.Parameter K-S A-D 𝑃 𝐵 𝑠 (cid:164) 𝑃 𝜏 𝑐 𝜏 𝐾 𝐿 𝑅 (cid:164) 𝐸 𝐵 𝑙𝑐 𝑄 𝛼 Table 2.
P-values resulting from Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) two-sample tests comparing the proposed populations ofpulsars with
𝑁 𝐹 < .
5% and pulsars with
𝑁 𝐹 > .
5% for spin periodin seconds ( 𝑃 ), magnetic field in Gauss ( 𝐵 𝑠 ), spindown rate in seconds persecond ( (cid:164) 𝑃 ), dispersion measure (DM), characteristic age ( 𝜏 𝑐 ), kinematicage ( 𝜏 𝑘 ), radio luminosity at 400 MHz ( 𝐿 𝑅 ), spindown luminosity ( (cid:164) 𝐸 ),pulse widths at 50% (W50) and 10% (W10), and magnetic inclination angle 𝛼 in 162 nulling pulsar entries with measured values of 𝑁 𝐹 ( % ) taken fromKonar & Deka (2019). We run these tests with 1000 sets of corrected NFmeasurements and report the median p-value. Given the above values, wefail to reject the null hypothesis that the two samples are drawn from thesame population. fail to reject the null hypothesis that nulling pulsars above and belowa nulling fraction of ∼
40% are drawn from the same population. Wedo find that the two samples are distinct ( 𝑝 < 0.05) in DM and 𝛼 and potentially distinct in (cid:164) 𝐸 , 𝐵 𝑙𝑐 , and 𝑄 . Although these factorsindicate separate populations, the other parameters, such as spinperiod and pulse width, have large p-values even though they are correlated with nulling fraction (see Section 6). This indicates thatthe apparent difference in populations is an artifact of selectionbias or a noise limit. We additionally run these two-sample tests onthe uncorrected NF measurements. We find that the two uncorrectedsamples are distinct in DM, log 𝜏 𝑐 , 𝐿 𝑅 , (cid:164) 𝐸 , 𝑄 , and 𝛼 , but again notdistinct in the more characteristic parameters. We therefore concludethat the two sets of pulsars are drawn from the same population eventhough they are statistically distinct in some parameters. We divide our sample into the morphological classifications of pulseprofile (Rankin 1983a). Using Kolmogorov-Smirnov and Anderson-Darling tests, we investigate whether the populations are statisticallydistinct in nulling fraction and characteristic age, both of which havebeen claimed in the literature (e.g., Hesse & Wielebinski 1974;Rankin 1986).We find that the relation between NF and morphological classis equivocal at best. Pulsars with double profiles (class “D”) are astatistically distinct population in NF when compared to all otherclasses, with p-values < .
05, except multiple profiles (class “M”),but no other classes are statistically distinct from each other.The relation between characteristic age and morphologicalclass is extant in our sample, but subtle. We find that S 𝑡 nullingpulsars are the youngest morphological class by median and mean,followed by 𝑇 , then S 𝑑 , then D, then M, which is consistent withthe age-morphology findings in Rankin (1986). However, in char-acteristic age, we cannot prove that neighboring classifications aredistinct (e.g., S 𝑡 and T), but any given class is statistically distinctfrom all non-neighboring classes (e.g., S 𝑡 is distinct from S 𝑑 , D,and M). The gradient of age across profile classes can be seen inFigure 3, but neighboring CDFs are too similar to be distinguishedby the K-S and A-D tests. We calculate the Pearson-R correlation value and the Maximal In-formation Coefficient (MIC) between nulling fraction and each ofthe fifteen parameters described in Section 2.We determine the significance of the MIC for a parameter 𝑥 by permuting the NF and 𝑥 values within 1,000 bootstrappedsub-samples of the dataset. We then calculate the p-value as thepercentage of permuted outcomes that meet or exceed the medianMIC value of the unscrambled and uncertainty-adjusted (correctedNF) data (see Section 3.4). We consider correlations with p-valuesless than 0.05 as statistically significant.In the past, only linear correlations have been investigated fornulling pulsars, but in this paper we are also able to investigate non-linear correlations by using the MIC (Section 3.3). The differencebetween the two correlation measures (MIC - 𝑅 ) describes thenon-linearity of a correlation.We find no correlations shown by Pearson-R between NF andany other variable. However, MIC reveals weak non-linear correla-tions between NF and: 𝑃 (MIC = 0.328, p-value = 0.009), (cid:164) 𝐸 (MIC= 0.334, p-value = 0.005), 𝐵 𝑙𝑐 (MIC = 0.326, p-value = 0.007), 𝑄 (MIC = 0.306, p-value = 0.012), and W10 (MIC = 0.337, p-value =0.003). MNRAS000
05, except multiple profiles (class “M”),but no other classes are statistically distinct from each other.The relation between characteristic age and morphologicalclass is extant in our sample, but subtle. We find that S 𝑡 nullingpulsars are the youngest morphological class by median and mean,followed by 𝑇 , then S 𝑑 , then D, then M, which is consistent withthe age-morphology findings in Rankin (1986). However, in char-acteristic age, we cannot prove that neighboring classifications aredistinct (e.g., S 𝑡 and T), but any given class is statistically distinctfrom all non-neighboring classes (e.g., S 𝑡 is distinct from S 𝑑 , D,and M). The gradient of age across profile classes can be seen inFigure 3, but neighboring CDFs are too similar to be distinguishedby the K-S and A-D tests. We calculate the Pearson-R correlation value and the Maximal In-formation Coefficient (MIC) between nulling fraction and each ofthe fifteen parameters described in Section 2.We determine the significance of the MIC for a parameter 𝑥 by permuting the NF and 𝑥 values within 1,000 bootstrappedsub-samples of the dataset. We then calculate the p-value as thepercentage of permuted outcomes that meet or exceed the medianMIC value of the unscrambled and uncertainty-adjusted (correctedNF) data (see Section 3.4). We consider correlations with p-valuesless than 0.05 as statistically significant.In the past, only linear correlations have been investigated fornulling pulsars, but in this paper we are also able to investigate non-linear correlations by using the MIC (Section 3.3). The differencebetween the two correlation measures (MIC - 𝑅 ) describes thenon-linearity of a correlation.We find no correlations shown by Pearson-R between NF andany other variable. However, MIC reveals weak non-linear correla-tions between NF and: 𝑃 (MIC = 0.328, p-value = 0.009), (cid:164) 𝐸 (MIC= 0.334, p-value = 0.005), 𝐵 𝑙𝑐 (MIC = 0.326, p-value = 0.007), 𝑄 (MIC = 0.306, p-value = 0.012), and W10 (MIC = 0.337, p-value =0.003). MNRAS000 , 1–12 (2020) ulling Pulsar Statistics . . . . . log τ c F n ( x ) StTSdDM
Figure 3.
Cumulative Distribution Functions of characteristic age for thefive morphological classifications S 𝑡 , T, S 𝑑 , D, and M. We find classesto be distinct from other classifications, but cannot prove that neighboringclassifications (e.g. S 𝑡 and T) are distinct from one another via K-S or A-Dtests. The correlations with (cid:164) 𝐸 , 𝐵 𝑙𝑐 , and 𝑄 are not particularly in-teresting because all three parameters probably take their correla-tion from their heavy 𝑃 -dependence. However, the weak non-lineartrends with 𝑃 and pulse width align with reports in the previousliterature (Biggs 1992; Li & Wang 1995) and could point to funda-mental physical nulling behaviour. We show the median NF values(for those with upper limits) in Figures 4 and 5. The measured trendsmight be driven by selection effects, as illustrated by the formationof very similar distributions by the nulling pulsar population andrandomly scrambled samples from the non-nulling pulsar popula-tion in Figures 4 and 5. We therefore conclude that, while thesecorrelations are statistically significant, they are not necessarily sci-entifically significant and caution against drawing conclusions fromthem until we have more data. In the meantime, future theories ofpulsar nulling should account for these trends with period and pulsewidth.To determine whether our method for accounting for the uncer-tainties was affected by the choice of distribution, we also performthe same test using a uniform distribution for the uncertainty ad-justment instead of a log uniform distribution for NFs reported asan upper limit. We find that two additional correlations, with 𝐿 𝑅 and 𝑊
50, are only significant in the uniform and log uniform dis-tributions, respectively. Because these correlations are dependenton the uncertainties, and there is no clear scientific basis to chooseone distribution over the other, we do not report these parametersas correlated as to NF.We find numerous correlations between NF and the other pa-rameters within morphological classes. However, as the sample sizesare extremely small (13–36 pulsars), we note only that all classesuniformly show statistically significant correlations with the param-eters found for the entire population: 𝑃 , (cid:164) 𝐸 , 𝑄 , W10 and W50, withthe exception of 𝐵 𝑙𝑐 which is correlated with NF in all classes ex-cept D. We also recover two potentially interesting correlations notseen in the overall population: 𝐵 𝑠 and (cid:164) 𝑃 are non-linearly correlated N u lli ng F r a c t i on 0 . . . Figure 4.
Nulling fraction versus the pulse width (in milliseconds) for 162nulling pulsar measurements in red triangles. Nulling fractions that haveonly been measured as upper limits are displayed as this upper limit. Forcontext, we plot the periods of the non-nulling populations against randomlyassigned nulling fractions from the nulling population in semi-transparentgrey squares to show the expected density with no correlation. Given thatthe randomly sampled non-nulling pulsar population shows a similar distri-bution to the nulling pulsar distribution, the correlation is most likely dueto selection effects. However, the paucity of points in the lower right quad-rant is different between the two distributions and could be suggestive of ascientifically interesting trend. with NF in all classes except S 𝑡 . However, because 𝐵 𝑠 is derivedfrom 𝑃 and (cid:164) 𝑃 , this correlation is probably just reflective of a trendwith (cid:164) 𝑃 . We calculate the Pearson-R correlation value and the Maximal In-formation Coefficient (MIC) for every combination of the fifteenparameters described in Section 2, with the exception of morphol-ogy. Many correlations between pulsar parameters within thenulling pulsar dataset are expected because the quantities are de-rived from one another. For example, we expect to see (and do see)a correlation between period 𝑃 and characteristic age 𝜏 𝑐 . This classof derived correlations accounts for the majority of correlationsthat we find, and while scientifically unoriginal, is still valuable toconfirm our methods.We find a weak but significant non-linear correlation betweendispersion measure DM and surface magnetic field 𝐵 𝑠 (MIC =0.372, p-value = 4 × − ), however we do not see this correlationin the overall non-nulling population. There is a paucity of pulsarswith low surface magnetic fields at large DMs, suggesting that thecorrelation we find is due to a detection bias against distant, lowenergy pulsars. The nulling population is most likely more sensitiveto this bias because nulling pulsars must be observed for a longerspan of time to reach the same signal-to-noise ratio as non-nullingpulsars during pulsar discovery surveys.There is also a weak but significant non-linear correlation be-tween spindown rate (cid:164) 𝑃 and width of the pulse W10 (MIC = 0.318, MNRAS , 1–12 (2020)
Sheikh and MacDonald N u lli ng F r a c t i on 0 . . . Figure 5.
Nulling fraction versus the pulse period for 162 nulling pulsarmeasurements in red triangles. Nulling fractions that have only been mea-sured as upper limits are displayed as this upper limit. For context, we plot theperiods of the non-nulling populations against randomly assigned nullingfractions from the nulling population in semi-transparent grey squares toshow the expected density with no correlation. Since the randomly sampledand scrambled non-nulling pulsar population shows a similar distribution tothe nulling pulsar distribution, the correlation is most likely due to selectioneffects. p-value = 0.006), a correlation that we also see in the non-nullingpopulation. This relationship has been previously reported withinalternating spin-down states in individual pulsars (e.g., Lyne et al.2010; Perera et al. 2016).
To verify that nulling and normal, radio pulsars are two distinct pop-ulations, we perform Kolmogorov–Smirnov (K-S) and Anderson-Darling (A-D) two-sample tests, comparing 162 nulling pulsar en-tries (see Section 2) with 2307 radio pulsars in nine parameters thatare available via ATNF (Manchester et al. 2005). We present theresulting p-values from both tests in Table 3. We find that p-valuesfrom all parameters except 𝑊
10 and 𝑊
50 are less than 𝛼 = . To put this work in context, we compile the correlations that havebeen claimed or suggested by all nulling pulsar statistical studies,including this one, in Table 4. Table 4 highlights the lack of aconsensus in the literature about correlations between nulling frac-tion (NF) and intrinsic pulsar parameters. We are the only studyto have searched for correlations in every previously suggestedparameter (see Section 2 for a more thorough description of the included parameters). Only recently have enough nulling pulsarsbeen discovered to draw statistically significant conclusions aboutthe population — this can be seen in the chronological ordering ofTable 4.In Section 4.2, we conclude that pulsars with NF below 37.5%and pulsars with NF above 37.5% are drawn from the same un-derlying population. In some variables, however, we find the twopopulations to be distinct such as in Dispersion Measure (DM). DMis dependent on distance which reveals a potential observationalbias: pulsars that null for longer are less likely to appear in e.g.a drift scan survey, so we would expect discovered pulsars withhigher nulling fractions to be those that are brighter and/or closer toEarth. To confirm that our results are not being affected by nullingpulsars with multiple NF measurements, , we re-run our analysis afew times by randomly choosing a singular measurement of NF foreach pulsar and by omitting these pulsars entirely. We find this doesnot change our results.In Figure 7, we show a 𝑃 – (cid:164) 𝑃 plot for all of the pulsars inthe ATNF database (Manchester et al. 2005), with nulling pulsarshighlighted in red triangles. The “tail” in the plot towards the bottomleft corner is comprised of the millisecond pulsar population, whichhas extremely short spin periods and very small (cid:164) 𝑃 values. Rajwadeet al. (2016) suggested that millisecond pulsars might not null atall. The current nulling pulsar population, as shown in Figure 7supports this conclusion.Given the non-linear correlation between NF and 𝑃 that wefound in this work, the absence of nulling millisecond pulsars pointsto a fundamental relationship between spin period and nulling be-haviour. This motivated the decision to remove both millisecondpulsars and all other non-normal pulsars from the comparison be-tween nulling and non-nulling pulsars in Section 7.To ensure that the uncertainties resulting from pulse widthsmeasured at different frequencies are not obscuring our results, weperform the following. First, we find that our pulse widths weremeasured at the following frequencies: roughly one third of themeasurements (both W10 and W50) were made at 1400 MHz, onethird at 408 MHz, and the rest between those frequencies, withonly a single outlier at 5000 MHz (the W10 value for B1641-45).This factor of 3.5 in observing frequency translates to differencesin measured pulse widths. Previous literature shows that, becauseobserving frequency is related to pulse width by a power law, dif-ferences in observing frequency of factors <
10 do not greatly affectpulse widths in most cases (e.g., D’Amico et al. 1998; Chen &Wang 2014). Rankin (1983b) found that a change of pulse width bya factor of 2 requires multiple orders-of-magnitude changes in ob-serving frequency. To prove that this potential source of error doesnot affect our results, we add in errors of 60% and perform a boot-strapped analysis. We find that, while the strength of the correlationdecreases in some of our trials, our conclusions do not change. Westill find a nonlinear, weak, but statistically significant correlationbetween pulse width and NF.While we searched for correlations between NF and other pa-rameters within pulse profile morphological classes in Section 6.2,we were limited by small-number statistics. In order to determinewhether nulling behaviour depends on morphological class, a much Our sample contains 141 individual pulsars, but numerous pulsars havedifferent estimates of nulling fraction. We therefore typically treat eachindividual nulling fraction measurement as a different pulsar
60% corresponds to the largest errors from Chen & Wang (2014), assum-ing pulse widths varied by a factor of 10, and ∼000
60% corresponds to the largest errors from Chen & Wang (2014), assum-ing pulse widths varied by a factor of 10, and ∼000 , 1–12 (2020) ulling Pulsar Statistics Parameter No. Nulling No. Radio K-S A-D 𝑃
162 2307 1.313e-05 4.196e-09 (cid:164) 𝑃
162 1958 0.03879 0.02330 𝐵 𝑠
162 1958 0.007411 0.005287DM 162 2306 4.408e-14 7.058e-23 𝜏 𝑐
162 1958 0.002137 0.001875 𝐿 𝑅
131 700 0.0293 0.04452 (cid:164) 𝐸
161 1958 0.0002043 1.714e-06 𝑊
50 158 1971 0.0363 0.09156 𝑊
10 145 1154 0.1802 0.08264
Table 3.
Resulting p-values from Kolmogorov–Smirnov (K-S) and Anderson-Darling (A-D) two-sample tests. Here, 𝑃 is spin period, (cid:164) 𝑃 is spin-down rate, 𝐵 𝑠 is the surface magnetic field strength, DM is the dispersion measure, 𝜏 𝑐 is characteristic age, 𝐿 𝑅 is the radio luminosity at 400 MHz, (cid:164) 𝐸 is the spin-downluminosity, 𝑊
50 is the pulse width at 50% peak brightness, and 𝑊
10 is the pulse width at 10% brightness. Given our small p-values, we reject the nullhypothesis that these two samples are from the same population, even though we cannot claim they are distinct in 𝑊
10 or 𝑊
50. We instead conclude thatnulling pulsars are not a subset of radio pulsars and are a distinct and separate population. See Figure 6 for the cumulative distribution functions of the twopopulations in all nine parameters.Study Sample Size 𝑃 (cid:164) 𝑃 𝐵 𝑠 𝜏 𝐶 𝛼 𝐿 𝑅 (cid:164) 𝐸 𝐵 𝑙𝑐 𝑄 𝑀 𝜏 𝑘 𝑊 𝑃 DMHesse & Wielebinski (1974) 15 (cid:88) (cid:88)
Ritchings (1976) 32 (cid:88) X (cid:88) (cid:88) Zhen-ru & Yi (1981) 32 (cid:88) (cid:88) (cid:88)
Rankin (1986) 59 X X X X X (cid:88)
Biggs (1992) 72 (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88)
Li & Wang (1995) 72 (cid:88)
Lazaridis & Seiradakis (2006) 19 (cid:88)
Wang et al. (2007a) 64 † X (cid:88) XCordes & Shannon (2008) 69 † (cid:88) Konar & Deka (2019) 204* X XThis Work 162 (cid:88)
X X X X X (cid:88) (cid:88) (cid:88)
X X (cid:88) X Table 4.
Pulsar parameters studied for correlation with nulling fraction in the literature, along with the results of this work. 𝑃 is spin period, (cid:164) 𝑃 is spin-downrate, 𝐵 𝑠 is the surface magnetic field strength, 𝜏 𝐶 is characteristic age, 𝛼 is the inclination angle between the magnetic and spin axes, 𝐿 𝑅 is the radio luminosityat 400 MHz, (cid:164) 𝐸 is the spin-down luminosity, 𝐵 𝑙𝑐 is the magnetic field strength at the light cylinder, 𝑄 is a dimensionless plasma flow parameter from Beskinet al. (1988), 𝑀 is a categorical variable representing the morphological class of the pulse profile (as per Rankin (1983a)), 𝜏 𝑘 is the kinetic age of the pulsar,and 𝑊 𝑝 is the total pulse width. Cells with an "X" indicate that the authors investigated that particular variable in their study and found no correlation withnulling fraction. Shaded cells with a checkmark indicate that a study suggested a correlation between the given variable and nulling fraction. Unshaded cellsindicate that a study did not use that particular parameter. We investigated correlations with MIC and 𝑅 and found that the five correlations that we foundwere non-linear. Morphology is a categorical variable so we used K-S and A-D tests to determine the association between morphology and nulling fraction.Li & Wang (1995) used equivalent width for the pulse width parameter, while we saw a correlation with W10 but not W50. † The authors did not report thenumber of nulling pulsars, so we estimated the number from provided figures. *This study also included pulsars that were known to null but did not have anestimated nulling fraction. Only 162 values were used in the correlation analyses. larger nulling population will be needed. One indirect way to in-crease the sample of nulling pulsars is to acknowledge the link be-tween mode-changing behaviour and nulling behaviour (e.g., Timo-khin 2010). The fraction of time that a mode-changing pulsar spendsin a less-energetic mode can be used as a de facto nulling fraction.Selection effects influence the population of nulling pulsarsmore than the overall pulsar population, as discussed in Biggs(1992). For example, an absence of emission from a pulsar is eas-ier to detect if the pulsar has a higher flux density (i.e., is nearerto Earth or more luminous). The issue is complicated by an addi-tional selection bias due to lack of sensitivity towards single pulses.For the majority of pulsars, multiple pulses must be integrated inorder to recover a pulse profile. Single pulse nulls can thus eas-ily be missed. Nearer pulsars will have lower DMs, explaining thecenter-left plot in Figure 6. Thus it is hard to characterize intrinsicdifferences between nulling and non-nulling pulsar parameters untilthe selection effects have been characterized. This could perhaps beaccomplished via simulated injection-recovery testing.Nulling fraction is only a single, and admittedly the crudest,measure of pulsar nulling behaviour. Gajjar (2017) suggested thatthe lack of a consistent trend of correlation between pulsar pa- rameters and NF was indicative of NF being a poor metric ratherthan any inherent behavioural cause. Nulling behaviour can also becharacterized with nulling length, which has been reported to be cor-related with age (Lazaridis & Seiradakis 2006). Yang et al. (2014)defines more complex parameters to quantify nulling that relate thelengths of the on and off states, and finds that spin-down rate andspin-down luminosity are correlated with these parameters. Nullingrandomness is yet another way to parameterize nulling behaviour(Gajjar 2017). Polarization parameters and timing noise, as of yetuncompiled for these pulsars, could also be important to character-izing nulling behaviour. Future statistical studies should incorporatemore nulling-related parameters once the data are available to doso.
MNRAS , 1–12 (2020) Sheikh and MacDonald . . . P (s) F n ( x ) − − . . . log P ⋅ ( s s ) . . . log B s . . . DM F n ( x ) . . . log τ c . . . log L R400
31 33 35 37 . . . log E ⋅ F n ( x ) nullingradio . . . log W50 1.0 2.0 3.0 . . . log W10 Figure 6.
Cumulative distribution functions comparing nulling pulsars in black to normal, radio pulsars in blue. Here, 𝑃 is spin period, (cid:164) 𝑃 is spin-down rate, 𝐵 𝑠 is the surface magnetic field strength, DM is the dispersion measure, 𝜏 𝑐 is characteristic age, 𝐿 𝑅
400 is the radio luminosity at 400 MHz, (cid:164) 𝐸 is the spin-downluminosity, 𝑊
50 is the pulse width at 50% brightness, and 𝑊
10 is the pulse width at 10% brightness. We conclude that these two populations are indeed twopopulations and not samples from a singular population.
About a tenth of all known pulsars temporarily cease their normalemission in a phenomenon known as “nulling.” Here, we investi-gate this population of 141 pulsars, using a variety of techniquesdescribed in Section 3. We first explore whether there is a gap inthe nulling fraction (NF) distribution around 40%, as suggested byKonar & Deka (2019), with cluster analyses. Then we determine ifthis gap does indeed separate two distinct pulsar populations. Us- ing both Kolmogorov-Smirnov and Anderson-Darling two-sampletests, we then investigate the statistical distinctness of different mor-phological classes within the overall nulling population. Followingthe suggestion in past studies that NF is correlated with a varietyof pulsar parameters, we test for such correlations between funda-mental pulsar parameters (see Section 2) and NF using both theMaximal Information Criterion and Pearson’s R. We also investi-gate correlations between all of the parameters in the nulling pulsarpopulation, independent of their behaviour with NF. Lastly, we ex-
MNRAS000
MNRAS000 , 1–12 (2020) ulling Pulsar Statistics P (s) P ⋅ ( ss ) − − Figure 7.
Period vs. Period Derivative plot for non-nulling pulsars (blackdots) and nulling pulsars (red triangles). The absence of red triangles in thebottom left of the plot indicates that there are no known millisecond pulsarsthat null. plore whether nulling pulsars and non-nulling pulsars are distinctpopulations. Following these analyses, we conclude the following: • There is a local minimum in the nulling fraction of the nullingpulsar population at about 40%, but we can not identify a natu-ral breakpoint, and this minimum does not divide two statisticallydistinct populations (Section 4). • In terms of nulling fraction, pulsars in morphological class Dare statistically distinct from all other classes except M. No otherclasses are statistically distinct from one another (Section 5). • Morphological class is related to age in the nulling pulsarpopulation, with the classes ordered, youngest to oldest, as S 𝑡 , T,S 𝑑 , D, M, which is consistent with Rankin (1986). This age gradientis statistically distinguishable, but also subtle with large scatter(Section 5). • We find nulling fraction in the nulling pulsar population to beweakly correlated with pulse period 𝑃 and full pulse width 𝑊
10, aswell as with parameters that are derived from 𝑃 ( (cid:164) 𝐸 , 𝐵 𝑙𝑐 , 𝑄 ), howeverthese trends should be taken with caution because of strong selectioneffects (Section 6.1). • We find dispersion measure DM and surface magnetic field 𝐵 𝑠 to be weakly, non-linearly correlated in the nulling pulsar populationbut not in the non-nulling population, likely due to selection bias(Section 6.2). • We find spindown rate (cid:164) 𝑃 and pulse width W10 to be weakly,non-linearly correlated in the nulling pulsar population as well as inthe non-nulling population, consistent with previous studies (Sec-tion 6.2). • Nulling and normal, radio pulsar populations are statisticallydistinct in nine fundamental parameters: period, spindown rate,surface magnetic field, dispersion measure, characteristic age, radioluminosity, spindown luminosity, and 10% pulse width and 50%pulse width (Section 7).This list of statistical outcomes might feel disjointed and de- tached from the physical sample, but each result gives us a clueabout either 1) the fundamental physics behind the pulsar nullingphenomenon or 2) the biases that must be accounted for when suchanalyses are performed.Future work must thoroughly investigate the selection effectspresent in the nulling pulsar sample. This can be achieved via creat-ing synthetic populations and simulating detection rates for pulsarswith different nulling fractions, dispersion measures, radio lumi-nosities, etc. for a range of instruments and survey parameters. Astudy of this magnitude would allow us to disentangle the connec-tion between nulling fraction and inherent pulsar properties fromthe effect of detection biases. As always, over time, the nulling pul-sar sample will become better characterized by the acquisition ofpolarizations, nulling lengths, null distributions, timing noise, andother data, which will lead to a better understanding of the emissionmechanism and its relationship to nulling behaviour.
ACKNOWLEDGEMENTS
We would like to thank the referee for the helpful commentsand suggestions that have improved this manuscript. We wouldlike to thank Vishal Gajjar, Alex Wolszczan, and Jason Wrightfor the helpful discussions. MGM acknowledges that this mate-rial is based upon work supported by the National Science Foun-dation Graduate Research Fellowship Program under Grant No.DGE1255832. Any opinions, findings, and conclusions or rec-ommendations expressed in this material are those of the au-thor and do not necessarily reflect the views of the NationalScience Foundation. We acknowledge use of the ATNF catalogat .The Center for Exoplanets and Habitable Worlds is supported bythe Pennsylvania State University, the Eberly College of Science,and the Pennsylvania Space Grant Consortium.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Albanese D., Riccadonna S., Donati C., Franceschi P., 2018, GigaScience,7, giy032Backer D., 1970, Nature, 228, 42Beskin V. S., 2009, MHD flows in compact astrophysical objects: accretion,winds and jets. Springer Science & Business MediaBeskin V., Gurevich A., Istomin Y. N., 1984, Astrophysics and Space Sci-ence, 102, 301Beskin V. S., Gurevich A. V., Istomin Y. N., 1988, Astrophysics and SpaceScience, 146, 205Biggs J. D., 1992, The Astrophysical Journal, 394, 574Bivand R., 2019, classInt: Choose Univariate Class Intervals. https://CRAN.R-project.org/package=classInt
Chen J. L., Wang H. G., 2014, ApJS, 215, 11Cordes J. M., Shannon R., 2008, The Astrophysical Journal, 682, 1152Coulson M. R., 1987, Cartographica: The International Journal for Geo-graphic Information and Geovisualization, 24, 16Cox D. R., 1972, Journal of the Royal Statistical Society: Series B (Method-ological), 34, 187D’Amico N., Stappers B. W., Bailes M., Martin C. E., Bell J. F., Lyne A. G.,Manchester R. N., 1998, MNRAS, 297, 28MNRAS , 1–12 (2020) Sheikh and MacDonald
Fisher W. D., 1958, Journal of the American statistical Association, 53, 789Gajjar V., 2017, arXiv preprint arXiv:1706.05407Gajjar V., Joshi B. C., Kramer M., 2012, Monthly Notices of the RoyalAstronomical Society, 424, 1197Gajjar V., Joshi B. C., Wright G., 2014a, Monthly Notices of the RoyalAstronomical Society, 439, 221Gajjar V., Joshi B. C., Kramer M., Karuppusamy R., Smits R., 2014b, TheAstrophysical Journal, 797, 18Gold T., 1969, Nature, 221, 25Goldreich P., 1969, Publications of the Astronomical Society of Australia,1, 227Hankins T., 1984, in Bulletin of the American Astronomical Society. pp468–469Hesse K., Wielebinski R., 1974, Astronomy and Astrophysics, 31, 409Huguenin G., Manchester R., Taylor J., 1971, The Astrophysical Journal,169, 97Jenks G. F., 1977, Department of Geographiy, University of Kansas Occa-sional PaperJianhua X., 2002, Beijing: Higher Education Press, 2002, 37Jokipii J., Lerche I., 1969, The Astrophysical Journal, 157, 1137Kendall M. G., 1938, Biometrika, 30, 81Kerr M., Hobbs G., Shannon R. M., Kiczynski M., Hollow R., Johnston S.,2014, Monthly Notices of the Royal Astronomical Society, 445, 320Konar S., Deka U., 2019, Journal of Astrophysics and Astronomy, 40, 42Kramer M., Lyne A. G., O’Brien J. T., Jordan C. A., Lorimer D. R., 2006,Science, 312, 549Lazaridis K., Seiradakis J. H., 2006, in AIP Conference Proceedings. pp309–315, doi:10.1063/1.2347995Leisch F., 1999, Technical Report 51, Bagged clustering. Vienna Universityof Economics and Business AdministrationLi X., Wang Z., 1995, Chinese Astronomy and Astrophysics, 19, 302Lyne A. G., Ashworth M., 1982, Monthly Notices of the Royal AstronomicalSociety, 204, 519Lyne A., Manchester R., 1988, Monthly Notices of the Royal AstronomicalSociety, 234, 477Lyne A., Ritchings R., Smith F., 1975, Monthly Notices of the Royal Astro-nomical Society, 171, 579Lyne A., Hobbs G., Kramer M., Stairs I., Stappers B., 2010, Science, 329,408MacQueen J., et al., 1967, Proceedings of the fifth Berkeley symposium onmathematical statistics and probability, 1, 281Malov I. F., 1990, Astronomicheskii Zhurnal, 67, 377Malov I., Nikitina E., 2011, Astronomy reports, 55, 19Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, The AstronomicalJournal, 129, 1993–2006Müllner D., et al., 2013, Journal of Statistical Software, 53, 1Naidu A., Joshi B. C., Manoharan P., Krishnakumar M., 2018, MonthlyNotices of the Royal Astronomical Society, 475, 2375Nikitina E. B., Malov I. F., 2017, Astronomy Reports, 61, 591Ostriker J., Gunn J., 1969, The Astrophysical Journal, 157, 1395Peacock J., 1983, Monthly Notices of the Royal Astronomical Society, 202,615Pearson K., 1895, proceedings of the royal society of London, 58, 240Perera B. B., Stappers B. W., Weltevrede P., Lyne A. G., Rankin J. M., 2016,Monthly Notices of the Royal Astronomical Society, 455, 1071R Core Team 2019, R: A Language and Environment for StatisticalComputing. R Foundation for Statistical Computing, Vienna, Austria,
Rahadianto H., Fariza A., Hasim J. A. N., 2015, in 2015 InternationalConference on Data and Software Engineering (ICoDSE). pp 195–200Rajwade K., Arjunwadkar M., Gupta Y., Kumar U., 2016, Under PreparationRankin J. M., 1983a, Toward an Empirical Theory of Pulsar Emission.I. Morphological Taxonomy, https://patentimages.storage.googleapis.com/41/b5/ca/69ffeea861af61/US8949899.pdf
Rankin J. M., 1983b, ApJ, 274, 359Rankin J. M., 1986, The Astrophysical Journal, 301, 901Rankin J. M., 1990, The Astrophysical Journal, 352, 247Reshef D. N., et al., 2011, science, 334, 1518 Ritchings R., 1976, Monthly Notices of the Royal Astronomical Society,176, 249Robitaille T. P., et al., 2013, Astronomy & Astrophysics, 558, A33Timokhin A. N., 2010, A model for nulling and mode chang-ing in pulsars, doi:10.1111/j.1745-3933.2010.00924.x, https://academic.oup.com/mnrasl/article-lookup/doi/10.1111/j.1745-3933.2010.00924.x
Wang N., Manchester R., Johnston S., 2007a, Monthly Notices of the RoyalAstronomical Society, 377, 1383Wang N., Manchester R., Johnston S., 2007b, Monthly Notices of the RoyalAstronomical Society, 377, 1383Yang A., Han J., Wang N., 2014, Science China: Physics, Mechanics andAstronomy, 57, 1600Zhen-ru W., Yi C., 1981, Symposium - International Astronomical Union,95, 215This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000