# Constraining the selection corrected luminosity function and total pulse
count for radio transients

Fengqiu Adam Dong Department of Physics and Astronomy, University of British
Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1 Canada  Antonio
Herrera-Martin David A. Dunlap Department of Astronomy & Astrophysics,
University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada
Department of Statistical Science, University of Toronto, Ontario Power
Building, 700 University Avenue, 9th Floor, Toronto, ON, Canada M5G 1Z5
Ingrid Stairs  Department of Physics and Astronomy, University of British
Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1 Canada  Radu V. Craiu
Department of Statistical Science, University of Toronto, Ontario Power
Building, 700 University Avenue, 9th Floor, Toronto, ON, Canada M5G 1Z5
Kathryn Crowter Department of Physics and Astronomy, University of British
Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1 Canada  Gwendolyn M.
Eadie David A. Dunlap Department of Astronomy & Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada  Department of
Statistical Science, University of Toronto, Ontario Power Building, 700
University Avenue, 9th Floor, Toronto, ON, Canada M5G 1Z5  Emmanuel Fonseca
Department of Physics and Astronomy, West Virginia University, Morgantown, WV
26506-6315, USA  Center for Gravitational Waves and Cosmology, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Deborah Good Physics Department,
University of Montana, Charles H Clapp, Building Rm 126, 32 Campus Drive,
Missoula, MT 59812  James W. Mckee E.A. Milne Centre for Astrophysics,
University of Hull, Cottingham Road, Hull, HU6 7RX, UK  Centre of Excellence
for Data Science, Artificial Intelligence and Modelling, University of Hull,
Cottingham Road, Hull, HU6 7RX, UK  Bradley W. Meyers International Centre for
Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia
Aaron B. Pearlman Department of Physics, McGill University, 3600 rue
University, Montréal, QC H3A 2T8, Canada  Trottier Space Institute, McGill
University, 3550 rue University, Montréal, QC H3A 2A7, Canada  David C.
Stenning Department of Statistics and Actuarial Science, Simon Fraser
University, 8888 University Drive, Burnaby, B.C., V5A 1S6 Canada

###### Abstract

Studying transient phenomena, such as individual pulses from pulsars, has
garnered considerable attention in the era of astronomical big data. Of
specific interest to this study are Rotating Radio Transients (RRATs),
nulling, and intermittent pulsars. This study introduces a new algorithm named
LuNfit, tailored to correct the selection biases originating from the
telescope and detection pipelines. Ultimately LuNfit estimates the intrinsic
luminosity distribution and nulling fraction of the single pulses emitted by
pulsars. LuNfit relies on Bayesian nested sampling so that the parameter space
can be fully explored. Bayesian nested sampling also provides the additional
benefit of simplifying model comparisons through the Bayes ratio. The
robustness of LuNfit is shown through simulations and applying LuNfit onto
pulsars with known nulling fractions. LuNfit is then applied to three RRATs,
J0012+5431, J1538+1523, and J2355+1523, extracting their intrinsic luminosity
distribution and burst rates. We find that their nulling fraction is 0.4(2),
0.749(5) and 0.995(2) respectively. We further find that a log-normal
distribution likely describes the single pulse luminosity distribution of
J0012+5431 and J1538+1523, while the Bayes ratio for J2355+1523 slightly
favors an exponential distribution. We show the conventional method of
correcting selection effects by “scaling up” the missed fraction of radio
transients can be unreliable when the mean luminosity of the source is faint
relative to the telescope sensitivity. Finally, we discuss the limitations of
the current implementation of LuNfit while also delving into potential
enhancements that would enable LuNfit to be applied to sources with complex
pulse morphologies.

pulsars

††facilities: CHIME††software: PRESTO (Ransom, 2001), CHIPSPIPE (Dong et al.,
2023), Numpy (Harris et al., 2020), Scipy (Virtanen et al., 2020), sigpyproc3

##  1 Introduction

Neutron stars result from the remnants of massive stars after they have
undergone core-collapse supernova. Rapid rotation of neutron stars likely
emerges due to angular momentum conservation following stellar collapse. The
rapid rotation combined with beamed radiation creates a lighthouse effect and
results in pulsed signals at observatories on Earth. The observed emission
exhibits a frequency-dependent delay due to cold-plasma dispersion, imparted
onto the radio signal as it traverses the interstellar medium. The dispersion
is dependent on the integrated line-of-sight electron density and is inversely
proportional to the square of the observing frequency (i.e. $\Delta
t=\frac{1}{2.41\times 10^{-4}}~{}\text{s}~{}\text{DM}$
pc/cm${}^{3}~{}((\nu_{1}/\text{GHz})^{2}-(\nu_{2}/\text{GHz})^{2})$ (Lorimer &
Kramer, 2004) where $\Delta t$ is the time delay and $\nu_{1}$ and $\nu_{2}$
are the upper and lower observing frequencies. DM is the dispersion measure).
The conventional technique to detect pulsars is to first correct for the
dispersion and then fold the observation on the pulsar period to build signal-
to-noise (S/N). Therefore, pulsar detection sensitivity can be significantly
boosted by longer integration times. Pulsar survey strategies must strike a
balance between integration time and sky coverage when using conventional
single-dish radio telescopes (Manchester et al., 2001; Cordes et al., 2006;
Stovall et al., 2014). Often, but not always, this leads to surveys focussing
on the Galactic plane. A counter-example is the Green Bank North Celestial Cap
survey which surveyed the whole northern hemisphere (Stovall et al., 2014).
Because of the large sky coverage, the survey only dwelled for 120 s per
pointing. The contemporary pulsar astronomy landscape operates in a distinctly
different parameter space such as the Canadian Hydrogen Intensity Mapping
Experiment (CHIME) and the Low-Frequency Array (LOFAR) (CHIME/FRB
Collaboration et al., 2018; Sanidas et al., 2019). Their wide field of view
allows for simultaneous significant sky coverage and long dwell times. With
these new capabilities, intermittent pulsars, nulling pulsars, and RRATs
become easier to find and characterize.

It has long been known that some pulsars do not emit on every rotation
(Backer, 1970). This phenomenon is called nulling, and up until recently, the
population of known nulling pulsars has been small compared to the total
pulsar population, likely due to the aforementioned survey methods (Konar &
Deka, 2019). With access to instruments with wide fields of view, this may be
changing. For example, Ng et al. (2020) showed that even long-known pulsars
thought to be persistent emitters could exhibit nulling behavior. It is
important to understand the nulling phenomenon as this can lead to better
models for pulsar emission mechanisms (van Leeuwen et al., 2003). On the
extreme end of the nulling pulsar spectrum are RRATs, highly intermittent
pulsars that can only be detected via their single pulses. As a large portion
of the pulsar population likely consists of RRATs (McLaughlin et al., 2006;
Keane & McLaughlin, 2011), understanding their intrinsic characteristics is
vital to understanding the pulsar population. However, RRATs can not be folded
to increase S/N and thus may require long observations to detect appreciable
amounts of single pulses.

A critical aspect of population analysis lacking in the current literature is
probing the intrinsic luminosity function and burst rate of highly
intermittent sources using only the detected pulsar single pulses. Many
studies of RRATs only describe an observed burst rate (e.g. Karako-Argaman et
al., 2015; Deneva et al., 2016; Good et al., 2021; Dong et al., 2023) rather
than correcting for selection effects to probe the intrinsic values. There
have been attempts to probe the intrinsic parameters for nulling pulsars
(Ritchings, 1976; Kaplan et al., 2018; Anumarlapudi et al., 2023) and also
another similar class of astrophysical phenomena, repeating Fast Radio Bursts,
luminous extragalactic radio bursts with similar burst characteristics to
RRATs (Li et al., 2021). The most established method described in Ritchings
(1976) relies on forming histogram bins and calculating
$\Delta=\rm{ON}-NF\times\rm{OFF}$ for each histogram bin, where $\rm{ON}$ and
$\rm{OFF}$ are the histogram bins of the integrated intensities during the on-
pulse window and off-pulse window respectively and $NF$ is the nulling
fraction. One would then minimize $|\Delta|$ by changing $NF$ to find the
preferred nulling fraction. However, this method and the Gaussian mixture
model in Kaplan et al. (2018) operate on the folded time series of the pulsar
observations and necessarily requires two conditions: the first of these
conditions is a period and a correctly identified pulse arrival phase. This is
not always possible with RRATs due to their extreme sporadicity (i.e. Dong et
al., 2023). The second requirement is recording the full-time series for the
observation, even during times of no emission. This is problematic for some
contemporary commensal instruments such as CHIME/FRB (CHIME/FRB Collaboration
et al., 2018), which only save data segments around detected single pulses.
While the method described in Li et al. (2021) looks at only single pulses,
they account for the missed pulses by “scaling up” by the measure detection
fraction. However, we show in section 5.1 that this method can be unreliable,
especially when the source has emitted many bursts below the observing
telescope detection threshold. In an attempt to avoid the limitations of
earlier works, we implement a new method, Luminosity and N fit (LuNfit), that
utilizes Bayesian inference using the nested sampling algorithm (Skilling,
2004) to constrain the intrinsic single pulse luminosity distribution and the
burst rates simultaneously. LuNfit is designed to remove any systematic biases
in the single pulse luminosity distribution and burst rate by empirically
measuring the selection effects of the telescope and detection pipeline.

This study is presented in the following way. Section 2 describes the CHIME
instrument and the implementation of LuNfit. We discuss the likelihood and
derivation of LuNfit. Section 2.1 described the underlying luminosity
distributions used. Section 2.2 details the priors used in LuNfit. Section 2.3
addresses the issue of discrete parameters. Section 2.4 details the
CHIME/Pulsar backend, the detection of pulses, and how single pulse
characteristics are retrieved following the detection of pulses. Section 2.5
details the process of injecting simulated pulses into the CHIME/Pulsar
detection pipeline to find the selection biases of the telescope. Section 3
details the simulations we performed to show the efficacy and robustness of
LuNfit. Section 3.1 compares LuNfit results with the known nulling fractions
of two pulsars, B1905+39 and J2044+4614. Section 3.2 describes the conversion
to flux units. Section 4 applies LuNfit to three RRATs, J0012+54, J1528+2345,
and J2355+1523. We present their intrinsic luminosity distributions and burst
rates. Finally, Section 5 compares LuNfit and conventional methods of
correcting selection biases and discusses LuNfit’s potential future
applications to pulsars, RRATs, and FRBs.

##  2 Methods

Table 1:  | CHIME/Pulsar
---|---
Receiver noise temperature |  $\sim$50K
Frequency Range | 400-800MHz
Number of beams | 10 (Tracking)
Beam width (FWHM) | 30’(400MHz)-15’(800MHz)
Time resolution | 327.68$\mu$s
Search Frequency Resolution | 390.625kHz
Coherent Dedispersion | Yes

CHIME is a commensal transit radio telescope situated in British Columbia,
Canada. It comprises three backend instruments: CHIME/cosmology, CHIME/FRB,
and CHIME/Pulsar (CHIME/FRB Collaboration et al., 2018; CHIME/Pulsar
Collaboration et al., 2021; CHIME Collaboration et al., 2022). Our study only
utilizes the CHIME/Pulsar backend, which processes up to 10 simultaneous
steerable tracking beams as sources transit the CHIME sky. The properties of
CHIME/Pulsar are provided in Table 1. The transit nature of the CHIME
telescope allows it to cover the whole northern sky every day, and while each
observation is only around $\sim 10$ minutes CHIME/Pulsar can achieve daily
cadence on astrophysical sources. Every day CHIME/Pulsar can observe $\sim
500$ individual sources. For this study, CHIME/Pulsar collects high-time
resolution spectra in the form of Sigproc111https://sigproc.sourceforge.net/
style filterbank data. CHIME’s sensitivity changes relative to the zenith
angle of the source; therefore, to maintain consistency within each
observation, we use only the data in which pulses from the nominal pulsar
position would be visible through the entire frequency range.

For this study, we are focussed on RRATs, which are only detectable via single
pulses. A Bernoulli process is sequence of binary random variables i.e. they
take only values of 0 or 1. Each pulse from an RRAT can be described as one of
these binary variables where 0 represents non-detection, and 1 represents
detection. While RRAT pulses have not been modelled as Bernoulli process
previously, RRATs are often assumed to follow a Poisson distribution, the
extreme case for the Bernoulli process (e.g. Chen et al., 2022; McKenna et
al., 2024). During any observation of a RRAT, any radio telescope will detect
only a portion of the emitted pulses due to limited sensitivity. This is
called the “selection’ of the telescope. Therefore, Assuming the selection
effects of the telescope are understood, we forward model the intrinsic
luminosity function of the RRAT. For N pulses emitted by a RRAT, the
likelihood is a series of Bernoulli trials for each pulse emitted. This can be
written as

$L(N,\phi_{1}...\phi_{N})\propto\prod_{i=1}^{N}\phi_{i}^{d}(1-\phi_{i})^{1-d}$
(1)

where $L$ is the likelihood, $i$ is the pulse number, $N$ is the total (but
unknown) number of pulses emitted by the RRAT, $\phi_{i}$ is the probability
of detecting pulse $i$ and $d=1$ if a detection was made, or $d=0$ otherwise.
Removing the terms in the product that go to one (e.g., when $d=1$, then the
terms $(1-\phi_{i})^{1-1}=1)$, we can rewrite this as

$L(N,\phi_{1}...\phi_{N})\propto\bigg{(}\prod_{i=1}^{n}\phi_{i}\bigg{)}\bigg{(}\prod_{j=n+1}^{N}(1-\phi_{j})\bigg{)}$
(2)

where $n$ is the number of detected pulses. The first product results from all
the pulses detected by the telescope, and the second product results from all
those not detected by the telescope. In practice, $\phi_{i}$ depends on the
pulse’s brightness and the intrinsic luminosity function of the RRAT.
Therefore,

$\phi_{i}(S_{i,det},\vec{\alpha})=P(det|S_{i,det})\int
P(S_{i,det}|S_{T},\sigma_{d})P(S_{T}|\vec{\alpha})\text{d}S_{T}$ (3)

where $S_{i,det}$ is the detected brightness of pulse $i$. $S_{T}$ is the true
intrinsic brightness as emitted by the RRAT. We emphasize the difference
between $S_{T}$ and $S_{det}$ where $S_{det}$ is the detected brightness and
comes with the errors due to the stochastic background noise of the telescope.
$S_{T}$ is the intrinsic brightness of the pulse emitted from the source
without any noise applied to it. The latter, even if inherently low, can be
amplified by random noise fluctuations, rendering it detectable and vice-
versa. Expressed differently, $S_{det}=S_{T}+S_{\sigma}$ where $S_{\sigma}$ is
the stochastic noise of the data. Depending on the chosen $S_{i,det}$ cutoff
this effect can be non-negligible. Figure 2 shows that this effect is most
apparent at the turn-off locations, near detection fractions of 0 and 1.
$P(det|S_{i,det})$ encodes all the selection biases of the telescope and
detection pipeline. $P(S_{T}|\vec{\alpha})$ is the probability density of the
intrinsic luminosity distribution of the RRAT where $\vec{\alpha}$ is a vector
containing the parameters of the luminosity model, that is,
$\vec{\alpha}=[\mu,\sigma]$, or $\vec{\alpha}=[k]$ for a log-normal or
exponential distribution, respectively. This is discussed in further in
section 2.1. $P(S_{i,det}|S_{T},\sigma_{d})$ designates the Gaussian detection
error with $\sigma_{d}$ being the known and fixed error. $\sigma_{d}$ is
measured empirically by performing injections, as discussed in Section 2.5.
Due to $S_{T}$ being an unmeasurable parameter, we marginalize over it. A
computational efficiency boost is achieved by realizing that
$P(S_{i,det}|S_{T},\sigma_{d})$ is a Gaussian and
$P(S_{det}|S_{T},\sigma_{d})=P(S_{det}-S_{T}|0,\sigma_{d})$ allowing the
right-hand integral of equation 3 to be written as a convolution.

For pulses that aren’t detected, we assume that the probability of non-
detection is the same for every pulse. That is,

$1-\phi_{j}(\vec{\alpha})=1-\int
P(det|S_{T})P(S_{T}|\vec{\alpha})\text{d}S_{T}.$ (4)

This assumption must be made as there is no way to know the brightness of a
pulse that is not detected. Finally, as the pulses are assumed to be
independent, the order does not matter. The independence assumption does not
consider refractive scintillation or the clustering of RRAT pulses. Refractive
scintillation can change the overall shape of $P(S_{T}|\vec{\alpha})$ and
occurs on timescales of weeks (Narayan, 1992). Therefore, when applying LuNfit
to data that spans longer than a few weeks, we point out that the fitted
luminosity function will be the intrinsic luminosity function modulated by
refractive scintillation. However, carefully accounting for the effects of
refractive scintillation with different classes of luminosity functions is out
of the scope of this study. We note that in this study, the timescales for
observations in Section 4 span months and therefore, the luminosity function
that has been fit is likely modulated by refractive scintillation. We do not
make any attempts to remove these effects. Despite this, as seen in Section 4,
it is unlikely that the shape of the luminosity function is much different
from a log-normal or exponential distribution. Any clustering of RRAT pulses
should not cause a problem for LuNfit as there is no dependence on each
pulse’s arrival time. We assume that all pulses from each cluster conform to a
global luminosity distribution and do not consider the case where each cluster
possesses an independent luminosity distribution. Thus, the likelihood is
multiplied by ${N\choose n}$ and becomes

$L(\vec{\alpha},N|n,S_{det})={N\choose n}\bigg{(}1-\int
P(det|S_{T})P(S_{T}|\vec{\alpha})\text{d}S_{T}\bigg{)}^{N-n}\prod_{i=1}^{n}P(det|S_{i,det})\int
P(S_{i,det}|S_{T},\sigma_{d})P(S_{T}|\vec{\alpha})\text{d}S_{T}.$ (5)

For this analysis, we take a Bayesian approach utilizing nested sampling
(Skilling, 2004) through the dynesty package (Speagle,
2020)222https://dynesty.readthedocs.io/en/stable/. Bayesian nested sampling is
a method which breaks up the posterior into many nested “slices”, assigns
weights and volumes to each slice and recombines it again to form the
posterior and evidence. It is a distinct way to approach Bayesian analysis
compared to Markov Chain Monte Carlo and possesses multiple benefits. For
example, nested sampling can sample the whole parameter space and avoid
problems with multimodal distributions. By extension, obtaining the model
evidence (marginal likelihood) for model comparisons through the Bayes ratio
is trivial. Both these features are crucial for LuNfit.

We use the detected signal to noise (S/N) for $S_{i,det}$. The strict
definitions are given in section 2.4. Nevertheless, other brightness metrics,
like detected flux density or fluence, can be employed instead. Details of
changing to flux density units are given in section 3.2.

###  2.1 Underlying Luminosity Distributions

One caveat in this formulation is the luminosity distribution chosen to
describe the detected single pulses, where choosing the incorrect model can
significantly alter the luminosity and nulling characteristics of the source
in question. However, as we employ a Bayesian nested sampling approach, we use
the Bayes ratio to differentiate between different models. In section 3, we
demonstrate with simulations that the Bayes ratio is robust in general. For
this reason, we have selected two models to describe pulsar emission
characteristics. The first model is the log-normal distribution parameterized
by the location, $\mu$, and scale, $\sigma$, described by

$P(S_{T}|\mu,\sigma)=\frac{1}{S_{T}\sigma\sqrt{2\pi}}\exp\left(-\frac{(\ln
S_{T}-\mu)^{2}}{2\sigma^{2}}\right).$ (6)

The second model is the exponential distribution parameterized by $k$, given
by:

$P(S_{T}|k)=k\exp(-kS_{T}).$ (7)

We chose these two models because they have demonstrated their ability to
explain various neutron star single pulse emission phenomena, such as RRAT
pulses (Burke-Spolaor et al., 2011; Mickaliger et al., 2018; Meyers et al.,
2018, 2019; Good et al., 2021; Dong et al., 2023), giant pulses (Bera &
Chengalur, 2019; Karuppusamy et al., 2010; McKee et al., 2019), and regular
single pulses (Burke-Spolaor et al., 2012; Liu et al., 2015; Bilous et al.,
2022). Repeating FRBs have shown double-peaked Gaussian distributions which
could also be considered in future studies (Li et al., 2021) .

###  2.2 Priors

Table 2: Prior types and parameters Parameter | Prior Type | $p_{1}$ | $p_{2}$
---|---|---|---
$\mu$ | Gaussian | 0 | 4
$\sigma$ | Inverse Gamma | 1.938 | 1
k | Inverse Gamma | 1 | 1
N | Uniform | n | Nrot = ($T_{obs}$/Period)

We assign general priors for both the luminosity function parameters and the
number of pulses. These priors are given in Table 2. Unlike Markov Chain Monte
Carlo methods that do not necessarily require normalized priors, we need the
Bayes Ratio for model comparisons and thus need robust evidence estimates. The
priors are chosen according to expected behavior. However, they do not
restrict unexpected behavior. To achieve this, all priors do not have strict
limits unless those limits are unphysical, such as in the case of $N$. All
chosen values on the priors have been chosen based on simulations of pulsars
for CHIME, where the underlying unit is S/N. Users of LuNfit should adjust
these values if they are using flux, fluence, or other brightness metrics.
For $\mu$, we choose a Gaussian prior described by the following equation

$P(\mu)=\frac{1}{p_{2}\sqrt{2\pi}}\exp{\bigg{(}-\frac{(\mu-
p_{1})^{2}}{2p_{2}^{2}}\bigg{)}}$ (8)

such that the prior contains $\sim 68$% of the probability mass between $-4$
and $4$. This ensures little bias towards particularly faint or bright
pulsars. Note that $\mu$ has units of ln(S/N), such that -4 corresponds with
S/N$\approx$0.02 and 4 corresponds to S/N$\approx$54.60. Therefore, we allow a
wide range of median values. If the user knows that the pulsar of interest is
particularly bright/faint relative to the telescope sensitivity, the prior on
$\mu$ should be adjusted to reflect that prior knowledge.

For $\sigma$ and $k$ we choose an inverse gamma distribution as the prior.
They are described by

$P(\alpha)=\frac{p_{2}^{p_{1}}}{\Gamma(p_{1})}(1/\alpha)^{p_{1}+1}\exp({p_{2}/\alpha})$
(9)

These priors ensure support for $\sigma$ and k across the viable parameter
range between 0 and $\infty$. However, we choose $p_{1}$ and $p_{2}$ such that
for $\sigma$, $\sim$90% of the probability mass is between 0 and 2, and for k,
$\sim$90% of the probability mass is between 0 and 10. The ranges represent
the most likely parameter space for detectable pulsars. Since $\sigma$ is
unitless, we believe that the priors chosen for $\sigma$ should be near
universal for all pulsars. We note that the prior on $\sigma$ did not need to
be altered in any set of simulated parameters. Although k is also unitless, k
determines the mean of an exponential distribution. Like $\mu$, the chosen
prior for k can differ if the pulsar is particularly bright/faint relative to
the telescope sensitivity.

The prior for N is determined by the number of detections made and the amount
of time we have observed the source. We assign a uniform distribution between
$n$ and $N_{rot}=\lfloor T_{obs}/P\rfloor$ where $\lfloor\rfloor$ is the floor
function, $T_{obs}$ is the amount of observation time, and $P$ is the pulsar
period. If $P$ is unknown, we set $N_{rot}=\lfloor T_{obs}/1\text{ms}\rfloor$.

###  2.3 Discrete N

Discrete parameters are non-trivial to implement in Bayesian nested sampling.
As N is necessarily discrete, we implement a pragmatic approach approximating
its discreteness. Dynesty draws samples from the inverse cumulative
distribution function (CDF). Therefore, to ensure discrete values of N are
drawn, we use the
scipy.stats.randint.ppf333https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.randint.html
function. This applies the following CDF:

${\rm CDF}=\frac{\lfloor N\rfloor-n+1}{N_{rot}-n}$ (10)

While this implementation may cause errors if n and N are $\mathcal{O}(1)$, in
the case of this study, both n and N and their associated errors are much
greater than 1. Therefore, we do not expect significant effects on the
inference of N.

###  2.4 CHIME and detection of pulses

Here, we describe how we obtain $S_{i,det}$ for each of the detected pulses.
To identify single pulses from known RRATs and pulsars, we first search the
CHIME/Pulsar filterbank data using the CHIME/Pulsar Single-pulse PIPEline
(CHIPSPIPE) (Dong et al., 2023). Once we detect a pulse, we extract 2 seconds
of data around the Time of Arrival (TOA) and downsample it by a factor of 3,
achieving a time resolution of approximately 1 ms. Next, we mask the channels
contaminated by Radio Frequency Interference (RFI), dedisperse the data to the
nominal DM of the pulsar, and average the data in frequency, producing a
single time series per pulse. To remove the underlying baseline, we remove the
pulse from the data and fit polynomials up to 10th order. For each fit, we
calculate a reduced $\chi^{2}$ value and accept the fit where the reduced
$\chi^{2}$ is closest to 1. This process is important as the CHIME/Pulsar
baseline can be highly variable. Once the baseline is fit, we subtract it from
the frequency collapsed time series, and the remaining background noise is
described by the standard deviation of $X_{n}$, the baseline-subtracted, pulse
removed time series ($\sigma_{noise}$).
To determine the amplitude of the pulse, we subtract the fitted baseline from
the original 2-second chunk of data. Then, we perform a maximum likelihood fit
of a Gaussian pulse given by the equation:

$F(X)=A\exp\left(-\frac{1}{2}\frac{(X-X_{loc})^{2}}{W^{2}}\right)+B$ (11)

where $X_{loc}$ represents the location of the Gaussian peak, $X$ is the
frequency collapsed time series, $A$ is the peak amplitude, $W$ is the
Gaussian width, and $B$ accounts for any remaining DC offset in the time
series. Throughout our study, unless otherwise specified, we use
$S/N=S_{det}=\frac{A}{\sigma_{noise}}$. We perform the Gaussian pulse fit on
each pulse detected by CHIPSPIPE. After fitting, we discard all $S_{det}<2$ to
ensure that every pulse is a good fit to the Gaussian. This fiducial cut-off
is also reflected in the selection effects that are measured.

###  2.5 Injections

Figure 1: The width effects of CHIPSPIPE. The plot shows the probability of
retrieval for different $S_{det}$ at varying widths. For regular slow pulsars,
the Gaussian widths are generally greater than 10ms. Thus, they are not
affected much by CHIPSPIPE selection effects. For sources that exhibit high
levels of pulse-to-pulse width variability, one needs to model the width along
with $S_{det}$ in LuNfit. For the current iteration of LuNfit, we only chose
pulsars with Gaussian widths larger than 5 ms to avoid the sensitivity drop-
off.

Figure 2: The selection effects of CHIPSPIPE. This plot shows the probability
of detection given $S_{i,det}$ or $S_{T}$. The left shows the fitted
$P(det|S_{T})$ compared to the forward modeled $P(det|S_{i,det})$. The right
shows how the forward model $P(det|S_{i,det})$ compared to the measured
values. This study uses the LuNfit implementation. We show the full selection
effects here, but when performing LuNfit, we cut off the selection at
$S_{i,det}=2$ to ensure robust pulses and retrieval probability measurements.
The dotted line indicates this cut. The blue line is a best-fit polynomial for
$P(det|S_{T})$. Figure 3: The detected CHIPSPIPE single pulse search S/N vs
the injected $S_{T}$ value. Not all injected pulses are detected for
detections of $S_{T}<\sim 3$. Thus, the CHIPSPIPE S/N will be biased high as
higher S/N pulses are being preferentially detected. This causes the
flattening of the dependence at low $S_{T}$.

In this section, we describe the process of obtaining $P(det|S_{i,det})$,
$P(det|S_{T})$, and $P(S_{i,det}|S_{T},\sigma_{d})$, with the primary goal of
empirically determining the retrieval probability (or detection fraction) at
different $S_{T}$ and $S_{i,det}$. These values are obtained empirically since
the complexity of CHIME/Pulsar and CHIPSPIPE makes it impossible to derive
them analytically. We inject simulated Gaussian pulses into real CHIME/Pulsar
data. The simulated pulses are broadband with flat spectral indices and widths
equal to the mean of the observed widths of the source in question. CHIPSPIPE
is not sensitive to width changes for regular slow pulsars or RRATs, so
injecting with one width per pulsar is sufficient. Figure 1 shows the effects
of width on pulse retrieval for CHIPSPIPE. We find that as long as most pulse
widths remain above 10 ms then the width effect of CHIPSPIPE is minimal.
Discussions of situations where one must include width effects in the
likelihood are provided in more detail in section 5.3. Since there are
hundreds of hours of CHIME/Pulsar observations for each source spread across
hundreds of observations, we spread out the injections over multiple
observations to perform about 20,000 injections per pulsar. Spreading the
injections across multiple observation epochs allows for an accurate
description of day-to-day RFI variations.

For a given injection we find the amplitude of the injected pulse in the
following way. We subtract the baseline by fitting a polynomial around each
injection TOA, and we measure $\sigma_{noise}$ in the same way that is
described in section 2.4. One can only control $S_{T}$ and not $S_{det}$ when
injecting. Thus, we inject the desired $S_{T}$ based on the $\sigma_{noise}$.
For each $S_{T}$, we inject multiple copies at different TOAs in different
observation files. The injected “observations” are saved and processed using
the same pipelines for real pulsar data. An example of the recovery rate is
shown in Figure 2, and binomial errors are assumed for each data point
$\bigg{(}\sqrt{\frac{p(1-p)}{m}}\bigg{)}$, where $p$ is the measured recovery
fraction and $m$ is the total number of injections at a specific $S_{T}$. The
low $S_{T}$ values observed here are due to a matter of S/N definition, and
Figure 3 illustrates the conversion between CHIPSPIPE/PRESTO detection
significance and $S_{T}$.

We assume that $P(S_{i,det}|S_{T},\sigma_{d})$ is Gaussian, and as $S_{T}$ is
a controlled parameter, $\sigma_{d}$ is the only unknown. To measure
$\sigma_{d}$, we fit each injected pulse as described in section 2.4.
Therefore for each group of injected $S_{T}$ we obtain a group of detected
$S_{i,det}$. We then take the standard deviation of the retrieved $S_{i,det}$
for each of the three brightest values of $S_{T}$ to use for $\sigma_{d}$.
Since $P(S_{i,det}|S_{T},\sigma_{d})$ is dominated by the stochastic noise in
the data, we assume $\sigma_{d}$ is constant across all injected $S_{T}$.

Two disparate methods are provided to measure $P(det|S_{i,det})$. The first
approach involves measuring $S_{i,det}$ directly from the injected
“observations”. For each injected pulse, we measure $S_{i,det}$ by the method
outlined in section 2.4. We bin them and calculate the recovery fraction at
each bin by comparing them with the detected injections. This approach is the
most accurate but may not be feasible for systems like CHIME/FRB, where one
cannot go back through the data stream to calculate $S_{det}$ for non-
detections. However, we injected into filterbank data, and even non-detections
can be fitted as we know the data segment where the injections were placed.
This is shown by red circles in the right panel of Figure 2.

In cases where the first approach is not feasible, a second approach is
provided, where one can forward model based on $P(S_{i,det}|S_{T},\sigma_{d})$
and $P(det|S_{T})$. It is found that for CHIPSPIPE, $P(det|S_{i,det})$ can be
well fit by a piecewise nth-order polynomial.

$\displaystyle\begin{aligned}
P(det|S_{det})=\begin{cases}0\qquad\qquad&S_{d}<S_{d0}\\\
c_{0}+c_{1}x_{1}+c2x_{2}^{2}...+c_{n}x_{n}^{n}\qquad\qquad&S_{d0}<S_{d}<S_{d1}\\\
1&S_{d}>S_{d1}\end{cases}\end{aligned}$ (12)

where $S_{d0}$ is determined by when $P(det|S_{det})$ first reaches 0 and
$S_{d1}$ is when $P(det|S_{det})$ first reaches 1, enforcing the function to
be continuous.

$P(det|S_{T})=\int P(det|S_{det})P(S_{det}|S_{T})dS_{det}$ (13)

Equation 13 shows the relationship between the measured probabilities,
$P(det|S_{T})$, and the derived property $P(det|S_{det})$. Forward modeling
can be done by using the user’s preferred fitting method. In this case, we
show a maximum likelihood fit for $P(det|S_{det})$ in Figure 2. We perform
many maximum likelihood fits with the polynomials of equation 12 ascending
from first order to tenth order. The polynomial with the closest reduced
$\chi^{2}_{r}$ to 1 is accepted. Therefore, the polynomial order can change
for each set of injections depending on the best fit.

##  3 Simulations

To assess the limitations of LuNfit, we simulate various sets of pulsars with
different luminosity parameters and numbers of detected pulses. For each set
of luminosity parameters, we simulate a pulse, add Gaussian noise based on the
measured $P(S_{i,det}|S_{T},\sigma_{d})$, and then determine whether it will
be detected according to the measured $P(det|S_{i,det})$. This process is
repeated until the desired number of pulses is obtained. We then employ LuNfit
to characterize each set of simulations. In total, 20 different sets of
parameters are simulated for each luminosity distribution. We simulate one
realization for each set of $\vec{\alpha}$ parameters for the following
scenarios:

  1. 1.

150 detections, which represents the minimum required for a reliable fit for
LuNfit.

  2. 2.

500 detections, a reasonable number typically observed for the most prolific
RRATs discovered by CHIME (Dong et al., 2023).

  3. 3.

1000 detections, approximately the number of detections made by the Five-
hundred-meter Aperture Spherical Telescope (FAST) for the most prolific FRBs
(Zhou et al., 2022; Li et al., 2022; Zhang et al., 2023).

  4. 4.

5000 detections, The best case scenario that can be achieved for the most
prolific RRATs.

These simulation results are presented in Appendix Figure 10 and 11. Note that
all simulations of the log-normal luminosity model were made with
$\sigma=0.5$. In summary, we recover the correct parameters in almost all
situations, with the exponential model requiring fewer detections for the same
error region compared to the log-normal model. This is due to the fact that
the log-normal distribution is characterized by two parameters (as opposed to
one for an exponential distribution) which allows more flexibility when
modelling the data. The recovered parameters are generally more accurate when
there are more detected pulses and (or) the median luminosity is higher. The
true N will be lower for a higher median luminosity, and LuNfit will not need
to account for as many missed pulses, leading to smaller error regions.
Conversely, for lower median luminosities, LuNfit accounts for more missed
bursts, and the error regions are larger. This characteristic is universal
across luminosity models. These simulations allow us to understand the
limitations of LuNfit and provide valuable insights into its performance under
various scenarios.

We compare the exponential and log-normal luminosity distributions via the
Bayes ratio to see if LuNFit can correctly identify the correct true
underlying distribution. We simulate four scenarios for both the log-normal
and exponential distributions: 150, 500, 1,000, and 5,000 detections (that is,
the true underlying distribution is log-normal or exponential). We simulate 20
parameter sets for each scenario. We vary $\mu$ from -0.5 to 1 for the log-
normal model while keeping $\sigma=0.5$. For the exponential distribution, we
vary k from 0.5 to 2. Appendix Figures 12 and 13 show the Bayes ratios between
the log-normal and exponential models for all of our simulations. LuNfit can
retrieve the correct luminosity distribution for most cases. In general, when
a pulsar follows an exponential distribution, LuNfit can always choose the
correct model. On the other hand, with only 150 detections, the Bayes ratio is
always low for the log-normal distribution. The Bayes ratio for log-normal
distributions increases as more detections are made or $\mu$ increases. The
reason for the $\mu$ dependence is that log-normal distributions approach
exponential distributions for low values of $\mu$. Therefore, it becomes more
difficult to parse the two distributions apart.

###  3.1 Validation on known pulsars

Figure 4: (a) LuNfit for B1905+39, as expected, the nulling fraction peaks at
$\sim 0$. The shaded region is forbidden due to negative nulling fractions. We
show this region to demonstrate that even when physical constraints are not
set, LuNfit converges on reasonable values. (b) LuNfit for B1905+39 again,
with the physical limits for nulling fraction. (c) LuNfit for J2044+4614. The
predicted nulling fraction according to the method outlined in Ritchings
(1976) is $\sim$0.41. This is shown by the red dashed line. The red-shaded
region shows the allowed nulling fractions provided by Ng et al. (2020). The
black dashed lines show the 68th percentile, 95th percentile confidence
intervals and the median.

We validate our method by performing LuNfit on two pulsars with known nulling
fractions, where the nulling fraction is given by $NF=1-\frac{N}{N_{rot}}$ and
$N_{rot}=\frac{T_{obs}}{P}$ is the total number of rotations that the
observation time allows. This is shown in Figure 4. B1905$+$39 (P=1.236 s) is
a bright, slow isolated pulsar and not known to exhibit any nulling. Indeed,
we find no evidence of nulling in the LuNfit results given in Figure 4. We
used one CHIME/Pulsar observation on MJD 59161, totaling 871 s, and detected
473 pulses. When LuNfit is applied to this pulsar, it accurately indicates a
nulling fraction of $-0.02(5)$, consistent with 0 and confirming its non-
nulling nature. In this calculation, we have lifted the prior on N and thus
allow $NF<0$ to show the efficacy of LuNfit. This is shown in Figure 4 (a) and
(b). In Figure 4 (a), we remove the constraints on the prior for N to show
that LuNfit converges on the correct value even when not constrained. The
constraints on N are then reapplied in Figure 4 (b).

J2044$+$4614 (P=1.393 s), on the other hand, is a known nulling pulsar, and
its nulling properties were previously discovered in Ng et al. (2020). We used
39 CHIME/Pulsar observations from MJD 58808 to MJD 58579 and detected 585
bursts with 34130 s of observations. We show the LuNfit results in Figure 4.
To compare our results with previously documented methods, we use the nulling
measurement method discussed in Ritchings (1976). The LuNfit derived nulling
fraction is 0.4(2) and agrees with the Ritchings (1976) method value of $\sim
0.41$. The two methods are shown in Figure 4 (c). In addition, the nulling
fraction as measured by Ng et al. (2020) is $>$9%, consistent with our
findings and shown as the shaded region in Figure 4 (c).

###  3.2 Conversion to flux density units

Conversion from the S/N based $S_{det}$ to Jy-based units can be done by
applying the single pulse radiometer equation (e.g Good et al., 2021),

$S_{peak}=S/N\frac{T_{sys}}{G\sqrt{\Delta\nu\Delta t~{}n_{p}}}$ (14)

to each detection. Then one would repeat the analysis detailed in section 2
with $S_{peak}$ in place of $S_{det}$, being careful to apply the same
radiometer conversion factor to both the $P(det|S_{i,det})$ and
$P(S_{i,det}|S_{T},\sigma_{d})$. We show the results of this in Table 3. Only
the luminosity function is altered in this change of variables, and thus, the
result for N does not change in Table 3. We caution, however, that due to
CHIME/Pulsar’s poor flux density calibration issues (Good et al., 2021; Dong
et al., 2023), these values should only be taken as a lower limit.

An alternative method for converting the flux is to transform the resultant
distribution. We can directly transform the units of the luminosity
distribution to flux density by

$P(S/N)\text{d}S/N=P(S_{peak})\text{d}S_{peak}.$ (15)

Where $P(S/N)$ is the probability distribution of the S/N, and $P(S_{peak})$
is the probability distribution of the peak flux density. Rearranging and
applying equation 14 gives

$P(S_{peak})=P(S/N)\bigg{(}\frac{T_{sys}}{G\sqrt{\Delta\nu\Delta
t~{}n_{p}}}\bigg{)}^{-1}.$ (16)

##  4 Results

Table 3: Table of results for J0012$+$5431, J1538$+$2345 and J2355$+$1523\. While both exponential and log-normal fit parameters are provided, Two of the three RRATs favor the log-normal model. Whereas J0012$+$5431 only weakly favors the log-normal distribution, J1538$+$2345 is has a strong preference for the log-normal distribution. The only RRAT to favor the exponential model is J2355+1523, but it is a weak Bayes ratio with large errors. In either case, both distributions for J2355+1523 strongly suggest a high nulling fraction. All the uncertainties provided are the 68th percentile confidence intervals. Pulsar | J0012+5431 | J1538+2345 | J2355+1523
---|---|---|---
# of detections | 194 | 6272 | 453
$\mu$ | -0.8(2) | 1.48(3) | 0.1(3)
$\sigma$ | 0.62(5) | 1.01(2) | 0.8(1)
$N_{logn}$ | 88000(36000) | 9400(200) | 2300(1900)
k | 2.1(1) | 0.141(2) | 0.62(3)
$N_{exp}$ | 39032(8700) | 8910(70) | 2010(160)
$\mu_{flux}$ | -4.0(2) | -1.94(3) | -3.1(6)
$\sigma_{flux}$ | 0.62(5) | 1.01(2) | 0.8(1)
$k_{flux}$ | 51(5) | 4.31(6) | 16(1)
Observation Time (s) | 442410 | 129196 | 352933
Period (s) | 3.025 | 3.449 | 0.913
Nulling Fraction | 0.4(2) | 0.749(5) | 0.995(3)
$\ln{Z_{ln}}-\ln{Z_{exp}}$ | 0.7(3) | 299.3(4) | -0.6(5)

We present our findings concerning the constrained luminosity function, the
intrinsic nulling fractions, and the total pulse number for three RRATs
discovered in Karako-Argaman et al. (2015) and Dong et al. (2023). They are
J0012+5431, J1538+2345 and J2355+1523. Detailed results are summarized in
Table 3, while visual representations of the dynesty nested sampling fits can
be found in Figures 5, 6, and 7.

Upon examination, a notable trend emerges; the fitted exponential distribution
consistently exhibits a higher median $S_{T}$ than the log-normal
distribution, encompassing all pulsars outlined in Table 3. This distinction
consequently results in a diminished total pulse count, N, for the exponential
distribution.

Figure 5: Corner plot and fit for J0012+5431. The log of the Bayes ratio is
0.7(3). (a) and (b) show the log-normal model, and (c) and (d) show the
exponential model. Figures (a) and (b) show the LuNfit corner plots for
J0012+5431. Figures (b) and (d) show the LuNfit fits as applied to the
detections of J0012+5431. The black dashed lines show the 68 percentile, 95
percentile, and the median. We also provide the selection effect for this
source as the blue line.

J0012+5431 is a RRAT first discovered in Dong et al. (2023). As a result, we
had an abundant amount of CHIME/Pulsar data available for this analysis. The
LuNfit results for J0012+5431 are presented in Figure 5. While it initially
appears to be extremely sporadic ($\sim 1.6$ pulses per hour), LuNfit shows
that this source is actually very faint and potentially less sporadic than
initially thought. Consequently, CHIME/Pulsar appears to be detecting only the
most intensely luminous tail of the distribution. We folded multiple
observations to find signs of continuous emission without success. However,
this is in line with expectations, considering the source’s faintness,
prolonged rotational period, and the brevity of CHIME/Pulsar observation
windows. The long rotational period dictates that within the brief $\sim
10$-minutes of a CHIME/Pulsar observation, accumulating a substantial number
of rotations to bolster the signal-to-noise ratio significantly is
unattainable. This suggests that follow-up with a more sensitive instrument
with longer integration times such as FAST or the GBT could confirm our
findings. With the Bayes Ratio only slightly leaning in favor of the log-
normal model, we present LuNfit’s results for both the log-normal and
exponential models in Figure 5.

Figure 6: Corner plot and fit for J1538+2345. As J1538+2345 highly prefers the
log-normal model, only that is shown. The left panel shows the log-normal
constraints on J1538+2345, and the right panel shows the fit given the mean
values. The black dashed lines show the 68 percentile, 95 percentile, and the
median.

J1538+2345 is a RRAT first discovered in Karako-Argaman et al. (2015). For
this RRAT, LuNfit heavily favors the log-normal distribution. Thus, only that
is provided in the LuNfit results in Figure 6. This source is extremely
prolific for a RRAT with a nulling fraction of 0.749(5). As many telescopes
have observed this source, we compare our results to that of the observational
burst rates reported by the Low Frequency Array (LOFAR) (Karako-Argaman et
al., 2015), the Green Bank Telescope (GBT) (Karako-Argaman et al., 2015) and
the Five Hundred Meter Aperture Synthesis Telescope (FAST) (Lu et al., 2019)
who reported observational nulling fractions of 0.94, 0.93, and 0.71
respectively. The FAST observations align the best with our selection
corrected nulling fraction of 0.749(5). This is likely due to the extreme
sensitivity of FAST, enabling it to probe the intrinsic properties of
J1538+2345. We obtained a much higher number of observations than Lu et al.
(2019). Therefore, we likely have a more accurate long-term average. The
slight discrepancy between the LuNfit results and the FAST results is likely
due to the length of observations. The study by Lu et al. (2019) was conducted
with only four observations, each lasting 30 minutes. The initial two
observations yielded significantly elevated levels of emitted pulses from
J1538+2345, which likely contributed to the lower nulling fractions reported
by FAST.

Figure 7: Fit and corner plots for J2355+1523. The log of the Bayes ratio is
2.7(4). (a) and (b) show the log-normal model, while (c) and (d) show the
exponential model. Figures (a) and (b) show the LuNfit corner plots for
J2355+1523. Figures (b) and (d) show the LuNfit fits as applied to the
detections of J2355+1523. The black dashed lines show the 68 percentile, 95
percentile, and the median. We also provide the selection effect for this
source as the blue line.

J2355+1523 is a RRAT first discovered in Dong et al. (2023). Among the RRATs
in that study, it was one of the most prolific sources identified, making it
an ideal candidate for LuNfit. The LuNfit results are shown in Figure 7. We
find that the Bayes ratio marginally favors the exponential distribution.
Thus, we provide LuNfit results for both models in Figure 7. In either case,
the source is bright, and the predicted nulling fraction is close to 1.

##  5 discussion

We have shown through this study that LuNfit is an effective tool for probing
the intrinsic luminosity function and burst rate of RRATs. This has been
validated through simulation and comparison with pulsars of known nulling
fractions. We then applied LuNfit to three known RRATs, finding their
intrinsic luminosity functions and nulling fractions. Below, we provide a
detailed discussion of LuNfit’s limitations, how it compares with other
techniques, and its use cases. Furthermore, we provide a discussion on how
LuNfit can be adapted to be used for exotic radio transients such as FRBs and
radio magnetars.

###  5.1 Comparisons with other methods

Figure 8: Comparison between a simple correction using the selection function
and LuNfit for an exponentially generated distribution with 1000 detections.
The ”simple correction” is the method described in Li et al. (2021), and a
full description is provided in section 5.1.

Figure 9: Comparison between a simple correction using the selection function
and LuNfit for a log-normal generated distribution with 1000 detections. The
red dashed line is the identity line. The ”simple correction” is the method
described in Li et al. (2021), and a full description is provided in section
5.1. All simulations are performed with a true $\sigma$ of 0.5.

The methods in the literature for correcting for the selection effects of
single pulses involve estimating the true burst count by scaling up by the
lost fraction, e.g., for FRB121102 (Li et al., 2021). This method involves
sampling the selection effects like LuNfit. Then, they place each detection
into a histogram bin. The histogram bin is divided by the detection fraction
to recover the “true count”. While this method can be effective when the
telescope’s sensitivity accurately captures the “turnover” of the luminosity
function, it becomes unreliable in cases where the source is dim, the
telescope lacks sensitivity or a combination of both.

Three factors contribute to the limitations of the Li et al. (2021) approach.
Firstly, accurately measuring the selection effects becomes extremely
challenging when the detection fraction is low. As the detection fraction
decreases, probing the selection function requires exponentially more
computing resources. To address this issue in our analysis, we cut off our
selection function at a minimum value of $S_{det}=2$. Secondly, even if the
selection effects are well measured, in parameter spaces where the detection
fraction is effectively zero, the simple correction used in conventional
methods leads to nonsensical results, as the corrected values go to infinity.
Lastly, the Li et al. (2021) correction requires binning of the detected
brightnesses, losing information in the process.

To illustrate the performance of LuNfit compared to the Li et al. (2021)
method outlined, we provide a comparison in Figure 8 and 9. For the Li et al.
(2021) method, we correct each histogram bin by dividing by the detection
fraction of that bin. We perform a maximum likelihood fit to the scaled bin
heights to get the luminosity function parameters. We find that LuNfit
exhibits much greater accuracy and precision compared to conventional methods
in all cases, especially when the true N is much higher than the detected n.
In summary, LuNfit offers an improved approach to handling selection effects
accurately and can provide more reliable results, especially in cases where
conventional methods may struggle due to low detection fractions.

###  5.2 Use for intermittent pulsars and comparisons with other tools

LuNfit has the potential to measure the nulling fraction of many intermittent
pulsars and RRATs using the vast observation capabilities of CHIME. While the
sample analyzed here is too small to say anything about the population, it
does show that there can exist a large range of nulling behavior between 0 and
1 nulling fractions. If raw filterbank data were saved, then the Gaussian
mixture model (Kaplan et al., 2018) or the Ritchings (1976) method also can
provide comparable results to LuNfit for measurements of the nulling fraction.
However, LuNfit has four benefits over the aforementioned tools. Firstly,
LuNfit can provide the single pulse luminosity function of intermittent
pulsars. Although, arguably, the Gaussian mixture model can also provide a
similar measure, this is not possible with the method from Ritchings (1976).
Secondly, as transient astronomy enters an era of big data, many observatories
and programs such as CHIME (CHIME/FRB Collaboration et al., 2018), UTMOST
(Farah et al., 2019), CRAFT (Macquart et al., 2010), MeerTRAP (Rajwade et al.,
2022) are not saving all observations, rather, only the segments that contain
identified pulses by their real-time pipelines. This poses a problem for the
Gaussian mixture model and the Ritchings (1976) method as they require data
measurements during non-detections. LuNfit avoids this problem by only
requiring a measurement of the selection function, which has already been
implemented in some experiments like CHIME/FRB (Merryfield et al., 2023).
Thirdly, the Gaussian mixture method and the Ritchings (1976) method require
identifying the pulse phase where the pulsar emits. If the source is extremely
intermittent, where there are only a few single pulses per observation, a few
pulses spread across many observations (i.e., some observations contain no
detectable emission), or the source is very faint such that even when folded,
no appreciable emission is seen, the analysis becomes near impossible with
both the Gaussian mixture method and the Ritchings (1976) method. This is not
a problem for LuNfit as long as enough single pulse detections were made on
the higher luminosity tail of the distribution. This is because LuNfit does
not require pulse phase information, only whether or not detections were made.
Finally, LuNfit gives the ability to select between different physically
driven models. In this study, we focused on the exponential and log-normal
distributions. However, other models, such as the double Gaussian, can be
used. LuNfit is based on nested sampling, so the Bayes ratio can always be
used to select the preferred model.

###  5.3 Potential for repeating FRB studies and limitations

We have outlined here the application of LuNfit for RRATs and intermittent
slow pulsars. LuNfit can also constrain other radio transients, such as
repeating FRBs. Tools like the Gaussian mixture model (Kaplan et al., 2018)
and the Ritchings (1976) method can not be used on repeating FRBs due to the
lack of strict periodicity. In our current implementation of LuNfit, we do not
include the width of the pulses in the analysis. This decision is based on the
observation that for slow pulsars, the width variations are generally small
enough to ignore their impact on the selection function. However, this
assumption may not hold for repeating FRBs that exhibit drastic width
variability (Pleunis et al., 2021). Therefore, LuNfit can only be applied to
repeating FRBs that show minimal width variation. Future iterations of LuNfit
will incorporate width dependence to address this issue. As more parameters
are introduced, the computational complexity of the analysis increases.
Dealing with selection effects in multiple parameters becomes more resource-
intensive. Therefore, optimizations will be required to retrieve selection
effects while including width variations effectively.

##  6 Conclusion and future work

In this study, we demonstrate the importance of accounting for selection
biases inherent to the telescope and detection pipeline. This allows the
accurate determination of the intrinsic pulse rate and luminosity function of
a slow pulsar or RRAT. To achieve this, we created an analysis framework named
LuNfit utilizing Bayesian nested sampling, leveraging the dynesty package.

We present detailed simulations and validation procedures for LuNfit in
section 2. This involved many sets of simulations of fake pulsars and
characterizing the nulling fraction of two pulsars, B1905+39 and J2044+4614,
where the nulling fraction is known. As a conclusive application of LuNfit, we
apply it to three known RRATs, successfully ascertaining their intrinsic
luminosity functions and burst rates. We show that the LuNfit nulling
fractions for J1538+2345 align most closely with the FAST observations. We
argue that this is likely due to the impressive sensitivity of FAST, enabling
them to probe the intrinsic properties of J1538+2345 without much selection
bias. The nulling fraction measured by FAST is likely biased slightly high due
to short observations. Notably, our findings in section 5.1 highlight the
limitations of conventional techniques in capturing the intrinsic luminosity
distribution and burst rate when compared with LuNfit.

Looking ahead, future work includes improving LuNfit by incorporating width
dependence and more model types, such as the single and double Gaussian. This
extension will render LuNfit more versatile, enabling its effective
utilization for sources emitting bursts with complex and varying morphologies.
This enhancement will allow LuNfit to be applied to all repeating Fast Radio
Bursts (FRBs), pulsars, and radio magnetars. Additionally, harnessing LuNfit’s
current capabilities, we intend to quantify the nulling fraction for numerous
intermittent pulsars and RRATs using CHIME/Pulsar observations. This will
provide a systematic understanding of the nulling slow pulsar population.
Finally, we plan to make the codebase publicly available with the next
iteration of LuNfit in the near future.

##  7 Acknowledgements

We acknowledge that CHIME is located on the traditional, ancestral, and
unceded territory of the Syilx/Okanagan people. We are grateful to the staff
of the Dominion Radio Astrophysical Observatory, which is operated by the
National Research Council of Canada. CHIME is funded by a grant from the
Canada Foundation for Innovation (CFI) 2012 Leading Edge Fund (Project 31170)
and by contributions from the provinces of British Columbia, Québec and
Ontario. The CHIME/FRB Project, which enabled development in common with the
CHIME/Pulsar instrument, is funded by a grant from the CFI 2015 Innovation
Fund (Project 33213) and by contributions from the provinces of British
Columbia and Québec, and by the Dunlap Institute for Astronomy and
Astrophysics at the University of Toronto. Additional support was provided by
the Canadian Institute for Advanced Research (CIFAR), McGill University and
the McGill Space Institute thanks to the Trottier Family Foundation, and the
University of British Columbia. The CHIME/Pulsar instrument hardware was
funded by NSERC RTI-1 grant EQPEQ 458893-2014. This research was enabled in
part by support provided by the BC Digital Research Infrastructure Group and
the Digital Research Alliance of Canada (alliancecan.ca).

F.A.D is supported by the UBC Four Year Fellowship

D.C.S. is supported by an NSERC Discovery Grant (RGPIN-2021-03985) and by a
Canadian Statistical Sciences Institute (CANSSI) Collaborative Research Team
Grant.

RVC is supported by an NSERC Discovery Grant RGPIN-2018-05663

A.B.P. is a Banting Fellow, a McGill Space Institute (MSI) Fellow, and a Fonds
de Recherche du Quebec – Nature et Technologies (FRQNT) postdoctoral fellow.

GME is supported by a CANSSI Collaborative Research Team Grant with support
from NSERC, and an NSERC Discovery Grant RGPIN-2020-04554.

Pulsar and FRB research at UBC are supported by an NSERC Discovery Grant and
by the Canadian Institute for Advanced Research

##  Appendix A Simulations

Figure 10: Simulations and recovery capabilities of LuNfit for the exponential
distribution. The axes show the LuNfit recovered quantities against the true
simulated quantities. N and k are the parameters constrained for an
exponential luminosity distribution. The red dashed line is the identify line.
All points should lie on this line for a perfect fitting algorithm. This
figure shows the recovered k (N) value given a simulated (true) k (N) value.
As expected for high true N, the error bars increase as this corresponds to an
intrinsically fainter pulsar. In other words, LuNfit needs to predict more
missed pulses for higher true N values. For an exponential distribution, a
fainter pulsar corresponds to a larger k. Therefore, we similarly find that
the errors for k increase for larger k. LuNfit performs significantly better
for an exponential distribution compared to a log-normal distribution due to
the smaller parameter space. Figure 11: Simulations and recovery capabilities
of LuNfit for the log-normal distribution. The axes show the LuNfit recovered
quantities against the true simulated quantities. $\mu$, $\sigma$, and N are
the parameters constrained for a log-normal luminosity distribution. All
simulations shown here have a $\sigma$ of 0.5. Again, the red dashed line is
the identity line. All points should fall on this line for a perfect fitting
algorithm. Like the exponential distribution shown in Figure 10, the fit is
worse for high true N and low true $\mu$ as this indicates a fainter pulsar.
The uncertainties of the log-normal luminosity model are markedly larger than
that of the exponential luminosity model. This is because of the larger
parameter space to explore. All simulations were performed with $\sigma=0.5$;
thus, all data points in the middle column should equal 1 (within
uncertainty). This figure shows that in most cases LuNfit is able to retrieve
the correct $\sigma$ value. Figure 12: Bayes ratio for a simulated exponential
distribution. The Bayes ratio is given in favor of the exponential
distribution for four sets of exponential luminosities distributions. That is,
positive values indicate a preference for the exponential distribution. There
are 150 (a), 500(b), 1000(c), and 5000(d) detected pulses in each simulation.
If the true underlying distribution is exponential, LuNfit always correctly
identifies this, as shown by the absence of negative values. In general, the
Bayes ratio is larger for low k and a larger number of detections. Figure 13:
Bayes ratio for simulated log-normal distributions. The Bayes ratio is given
in favor of the log-normal distribution for four sets of log-normal
luminosities. That is, positive values indicate a preference for a log-normal
distribution. There are 150 (a), 500(b), 1000(c) and 5000(d) detected pulses
in each simulation. The $\sigma$ for these simulations was set to 0.5. The
effects of a brighter median luminosity (higher $\mu$) are much more
pronounced for the log-normal distribution than the exponential distribution.
Higher $\mu$ values are clearly positively correlated with higher Bayes
ratios. This is due to LuNfit being able to completely fit the “turn-over” in
log-normal distributions. Additionally, the log-normal tends towards an
exponential for low $\mu$ values. Once again, a larger number of detections
aid the differentiation between log-normal and exponential distributions.

## References

  * Anumarlapudi et al. (2023) Anumarlapudi, A., Swiggum, J. K., Kaplan, D. L., & Fichtenbauer, T. D. J. 2023, ApJ, 948, 32, doi: 10.3847/1538-4357/acbb68
  * Backer (1970) Backer, D. C. 1970, Nature, 228, 42, doi: 10.1038/228042a0
  * Bera & Chengalur (2019) Bera, A., & Chengalur, J. N. 2019, MNRAS, 490, L12, doi: 10.1093/mnrasl/slz140
  * Bilous et al. (2022) Bilous, A. V., Grießmeier, J. M., Pennucci, T., et al. 2022, A&A, 658, A143, doi: 10.1051/0004-6361/202142242
  * Burke-Spolaor et al. (2011) Burke-Spolaor, S., Bailes, M., Johnston, S., et al. 2011, MNRAS, 416, 2465, doi: 10.1111/j.1365-2966.2011.18521.x
  * Burke-Spolaor et al. (2012) Burke-Spolaor, S., Johnston, S., Bailes, M., et al. 2012, MNRAS, 423, 1351, doi: 10.1111/j.1365-2966.2012.20998.x
  * Chen et al. (2022) Chen, J. L., Wen, Z. G., Yuan, J. P., et al. 2022, ApJ, 934, 24, doi: 10.3847/1538-4357/ac75d1
  * CHIME Collaboration et al. (2022) CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022, ApJS, 261, 29, doi: 10.3847/1538-4365/ac6fd9
  * CHIME/FRB Collaboration et al. (2018) CHIME/FRB Collaboration, Amiri, M., Bandura, K., et al. 2018, ApJ, 863, 48, doi: 10.3847/1538-4357/aad188
  * CHIME/Pulsar Collaboration et al. (2021) CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021, ApJS, 255, 5, doi: 10.3847/1538-4365/abfdcb
  * Cordes et al. (2006) Cordes, J. M., Freire, P. C. C., Lorimer, D. R., et al. 2006, ApJ, 637, 446, doi: 10.1086/498335
  * Deneva et al. (2016) Deneva, J. S., Stovall, K., McLaughlin, M. A., et al. 2016, ApJ, 821, 10, doi: 10.3847/0004-637X/821/1/10
  * Dong et al. (2023) Dong, F. A., Crowter, K., Meyers, B. W., et al. 2023, MNRAS, 524, 5132, doi: 10.1093/mnras/stad2012
  * Farah et al. (2019) Farah, W., Flynn, C., Bailes, M., et al. 2019, MNRAS, 488, 2989, doi: 10.1093/mnras/stz1748
  * Good et al. (2021) Good, D. C., Andersen, B. C., Chawla, P., et al. 2021, ApJ, 922, 43, doi: 10.3847/1538-4357/ac1da6
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Kaplan et al. (2018) Kaplan, D. L., Swiggum, J. K., Fichtenbauer, T. D. J., & Vallisneri, M. 2018, ApJ, 855, 14, doi: 10.3847/1538-4357/aaab62
  * Karako-Argaman et al. (2015) Karako-Argaman, C., Kaspi, V. M., Lynch, R. S., et al. 2015, ApJ, 809, 67, doi: 10.1088/0004-637X/809/1/67
  * Karuppusamy et al. (2010) Karuppusamy, R., Stappers, B. W., & van Straten, W. 2010, A&A, 515, A36, doi: 10.1051/0004-6361/200913729
  * Keane & McLaughlin (2011) Keane, E. F., & McLaughlin, M. A. 2011, Bulletin of the Astronomical Society of India, 39, 333, doi: 10.48550/arXiv.1109.6896
  * Konar & Deka (2019) Konar, S., & Deka, U. 2019, Journal of Astrophysics and Astronomy, 40, 42, doi: 10.1007/s12036-019-9608-z
  * Li et al. (2021) Li, D., Wang, P., Zhu, W. W., et al. 2021, Nature, 598, 267, doi: 10.1038/s41586-021-03878-5
  * Li et al. (2022) —. 2022, Nature, 601, E1, doi: 10.1038/s41586-021-04178-8
  * Liu et al. (2015) Liu, K., Karuppusamy, R., Lee, K. J., et al. 2015, MNRAS, 449, 1158, doi: 10.1093/mnras/stv397
  * Lorimer & Kramer (2004) Lorimer, D. R., & Kramer, M. 2004, Handbook of Pulsar Astronomy, Vol. 4
  * Lu et al. (2019) Lu, J., Peng, B., Liu, K., et al. 2019, Science China Physics, Mechanics, and Astronomy, 62, 959503, doi: 10.1007/s11433-018-9372-7
  * Macquart et al. (2010) Macquart, J.-P., Bailes, M., Bhat, N. D. R., et al. 2010, PASA, 27, 272, doi: 10.1071/AS09082
  * Manchester et al. (2001) Manchester, R. N., Lyne, A. G., Camilo, F., et al. 2001, MNRAS, 328, 17, doi: 10.1046/j.1365-8711.2001.04751.x
  * McKee et al. (2019) McKee, J. W., Stappers, B. W., Bassa, C. G., et al. 2019, MNRAS, 483, 4784, doi: 10.1093/mnras/sty3058
  * McKenna et al. (2024) McKenna, D. J., Keane, E. F., Gallagher, P. T., & McCauley, J. 2024, MNRAS, 527, 4397, doi: 10.1093/mnras/stad2900
  * McLaughlin et al. (2006) McLaughlin, M. A., Lyne, A. G., Lorimer, D. R., et al. 2006, Nature, 439, 817, doi: 10.1038/nature04440
  * Merryfield et al. (2023) Merryfield, M., Tendulkar, S. P., Shin, K., et al. 2023, AJ, 165, 152, doi: 10.3847/1538-3881/ac9ab5
  * Meyers et al. (2018) Meyers, B. W., Tremblay, S. E., Bhat, N. D. R., et al. 2018, ApJ, 869, 134, doi: 10.3847/1538-4357/aaee7b
  * Meyers et al. (2019) —. 2019, PASA, 36, e034, doi: 10.1017/pasa.2019.30
  * Mickaliger et al. (2018) Mickaliger, M. B., McEwen, A. E., McLaughlin, M. A., & Lorimer, D. R. 2018, MNRAS, 479, 5413, doi: 10.1093/mnras/sty1785
  * Narayan (1992) Narayan, R. 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 151, doi: 10.1098/rsta.1992.0090
  * Ng et al. (2020) Ng, C., Wu, B., Ma, M., et al. 2020, ApJ, 903, 81, doi: 10.3847/1538-4357/abb94f
  * Pleunis et al. (2021) Pleunis, Z., Good, D. C., Kaspi, V. M., et al. 2021, ApJ, 923, 1, doi: 10.3847/1538-4357/ac33ac
  * Rajwade et al. (2022) Rajwade, K. M., Bezuidenhout, M. C., Caleb, M., et al. 2022, MNRAS, 514, 1961, doi: 10.1093/mnras/stac1450
  * Ransom (2001) Ransom, S. M. 2001, PhD thesis, Harvard University, Massachusetts
  * Ritchings (1976) Ritchings, R. T. 1976, MNRAS, 176, 249, doi: 10.1093/mnras/176.2.249
  * Sanidas et al. (2019) Sanidas, S., Cooper, S., Bassa, C. G., et al. 2019, A&A, 626, A104, doi: 10.1051/0004-6361/201935609
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405, doi: 10.1063/1.1835238
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Stovall et al. (2014) Stovall, K., Lynch, R. S., Ransom, S. M., et al. 2014, ApJ, 791, 67, doi: 10.1088/0004-637X/791/1/67
  * van Leeuwen et al. (2003) van Leeuwen, A. G. J., Stappers, B. W., Ramachandran, R., & Rankin, J. M. 2003, A&A, 399, 223, doi: 10.1051/0004-6361:20021630
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Zhang et al. (2023) Zhang, Y.-K., Li, D., Zhang, B., et al. 2023, ApJ, 955, 142, doi: 10.3847/1538-4357/aced0b
  * Zhou et al. (2022) Zhou, D. J., Han, J. L., Zhang, B., et al. 2022, Research in Astronomy and Astrophysics, 22, 124001, doi: 10.1088/1674-4527/ac98f8