Bayesian Inference on the Radio-quietness of Gamma-ray Pulsars

The number of radio-quiet (RQ) and radio-loud (RL) gamma-ray pulsars that have been detected has increased rapidly during recent years thanks to the launch of the Fermi Gamma-Ray Space Telescope. The observed differences between RQ and RL pulsars are not expected to be intrinsic, but due to geometrical effects. Theoretically, the radio emissions are originated from the polar cap region, which is aligned with the magnetic axis of the pulsar (Ruderman & Sutherland 1975). On the other hand, the gamma rays are generally accepted to be originated from the outer gap (e.g., Cheng & Zhang 1998; Takata et al. 2006, 2008), which subtends a larger solid angle than the radio cones. Therefore, the probability of the line of sight (LOS) to intersect with the gamma-ray-emission region is higher. If the projection of the radio cone does not swap through the LOS, then the pulsar is seen as being RQ. Thus, the magnetic inclination angle, α, and the radio-cone opening half-angle, δ, play important roles in the modeling of gamma-ray pulsars.

Recently, Hui et al. (2017, hereafter H17) and Sokolova & Rubtsov (2016, hereafter S16) studied the observed and derived properties of RQ and RL gamma-ray pulsars. In our current study, we adopt the pulsar populations used in H17, which is extracted from the 117 gamma-ray-emitting pulsars detected in The Second Fermi Large Area Telescope (LAT, Atwood et al. 2009) Catalog of Gamma-ray Pulsars (2PC, Abdo et al. 2013) and The Fermi LAT Third Source Catalog (3FGL, Acero et al. 2015). Nevertheless, we also examine the reported discrepancy of the inferred fraction of RQ pulsars between these two studies.

Hui et al. (2017) performed a non-parametric analysis on the RQ and RL pulsar populations. Specifically, they applied the Anderson–Darling Test (the so-called A–D Test) to the unbinned data. They found that the magnetic field strength at the light cylinder ${B}_{{\rm{LC}}}$, the gamma-ray-to-X-ray flux ratio ${F}_{\gamma }/{F}_{{\rm{X}}}$, and the spectral curvature significance (CS) exhibit $\gt 3\sigma$ difference. In this paper, we try to verify or disprove their non-parametric frequentistic results using the advanced Bayesian formalism.

We apply a robust⁶ Bayesian statistical technique (Kruschke 2013) to the data to identify the key parameters that separate the groups of RQ and RL pulsars. We further obtain the posterior probability distribution of the fraction of RQ pulsars and the α-δ posterior distribution using Bayesian inference. Compared with conventional frequentist statistics, one of the advantages of Bayesian statistics is their ability to display the complete posterior probability distribution. The posterior provides information about the probability density, the uncertainty, modality, and correlation between parameters. Therefore, the uncertainties of the physical and statistical parameters can be readily obtained without the need to rely on complex error propagation; theoretical or numerical approximation; and sometimes point estimations that may not be representative if the distribution is highly skewed or multi-modal. It is also conceptually straightforward to integrate over the contributions of all but one parameter to get the marginal distribution, which displays complete information for the parameter in interest when the uncertainties of all other parameters are taken into account.

We detail the Bayesian analysis methods and results in respective subsections of Section 2. The conclusions and discussions are presented in Section 3, where the physical implications of the results are also discussed.

2.1. Difference between Radio-quiet and Radio-loud Gamma-ray Pulsars

We apply the Bayesian Estimation Supersedes the t-Test (BEST) technique, developed by Kruschke (2013), to the 12 observed distributions of the physical parameters listed in Tables 1 and 2 in H17. The analysis is performed using PyMC3 (Salvatier et al. 2016). The BEST technique employs a non-standardized t-distribution as the likelihood function

$$\begin{eqnarray}&&{\mathbb{P}}(x| \mu ,\sigma ,\nu )=\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{\nu +1}{2}\right)}{{\rm{\Gamma }}\left(\tfrac{\nu }{2}\right)}\left(\displaystyle \frac{1}{\sqrt{\pi \mu }\sigma }\right){\left[1+\displaystyle \frac{1}{\nu }{\left(\displaystyle \frac{x-\mu }{\sigma }\right)}^{2}\right]}^{-\tfrac{\nu +1}{2}},\end{eqnarray} \tag{ 1 }$$

where μ, σ, and ν are the statistical parameters of location, scale, and normality, respectively, and Γ is the gamma function. Because the t-distribution ${ \mathcal T }(\mu ,\sigma ,\nu )$ is broader than a Gaussian distribution such that its variance is undefined, ${\sigma }^{2}$ is not the variance. Note that the statistical parameter ν is a measure of how heavy the probability is in the tails. If ν is large ($\nu \gtrsim 30$), the t-distribution approximates a Gaussian distribution.

The t-distribution is used because it can better account for the contributions from the outliers in the population than a Gaussian distribution. Each physical parameter of RQ and RL pulsars is modeled with two t-distributions that are connected by the normality. As the observed distributions are highly skewed in linear space, we take the base-10 logarithm to all physical parameters, which means $\mathrm{log}{X}^{{\rm{RQ}}}\sim { \mathcal T }({\mu }_{X}^{{\rm{RQ}}},{\sigma }_{X}^{{\rm{RQ}}},{\nu }_{X})$ and $\mathrm{log}{X}^{{\rm{RL}}}\sim { \mathcal T }({\mu }_{X}^{{\rm{RL}}},{\sigma }_{X}^{{\rm{RL}}},{\nu }_{X})$. Therefore, we are essentially working with log-t-distributions.

The prior distributions of the statistical parameters are defined as

$$\begin{eqnarray}\left\{\begin{array}{rcl}{\mu }_{X}^{{\rm{RQ}}} & \sim & { \mathcal N }(\bar{X},2{S}_{X}),\\ {\mu }_{X}^{{\rm{RL}}} & \sim & { \mathcal N }(\bar{X},2{S}_{X}),\\ {\sigma }_{X}^{{\rm{RQ}}} & \sim & { \mathcal U }(0.1{S}_{X},10{S}_{X}),\\ {\sigma }_{X}^{{\rm{RL}}} & \sim & { \mathcal U }(0.1{S}_{X},10{S}_{X}),\\ {\nu }_{X}-1 & \sim & { \mathcal E }(1/29),\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where ${ \mathcal N }$, ${ \mathcal U }$, and ${ \mathcal E }$ are the normal, uniform, and exponential probability distributions, respectively. We use $\bar{X}$ and S_X to represent the mean and standard deviation of the physical parameter X, respectively, to distinguish from ${\mu }_{X}$ and ${\sigma }_{X}$, which have been defined earlier with different statistical meanings.

We perform a Markov chain Monte Carlo (MCMC) sampling using the above hierarchical Bayesian model. For each physical parameter, we use four chains and obtain 25,000 traces per chain, resulting in posterior distributions consisting of 100,000 traces. The inferred statistical parameters for individual physical parameter are listed in Table 1.

Table 1.  Bayesian Posterior Statistical Parameters

Physical Parameter	${\mu }^{{\rm{RQ}}}$	${\mu }^{{\rm{RL}}}$	${\sigma }^{{\rm{RQ}}}$	${\sigma }^{{\rm{RL}}}$	ν	${\mu }^{{\rm{RQ}}}-{\mu }^{{\rm{RL}}}$
$\mathrm{log}[p$ (ms)]	${2.227}_{-0.101}^{+0.103}$	${2.074}_{-0.107}^{+0.108}$	${0.219}_{-0.066}^{+0.087}$	${0.257}_{-0.074}^{+0.093}$	${40.031}_{-37.387}^{+107.513}$	${0.153}_{-0.148}^{+0.146}$
$\mathrm{log}[\dot{p}$ (s s−1)]	$-{13.58}_{-0.33}^{+0.33}$	$-{13.33}_{-0.28}^{+0.27}$	${0.70}_{-0.23}^{+0.29}$	${0.66}_{-0.22}^{+0.25}$	${32.42}_{-30.59}^{+102.61}$	$-{0.25}_{-0.42}^{+0.44}$
$\mathrm{log}[{B}_{{\rm{LC}}}$ (G)]	${3.606}_{-0.314}^{+0.320}$	${4.118}_{-0.288}^{+0.289}$	${0.686}_{-0.204}^{+0.269}$	${0.686}_{-0.200}^{+0.243}$	${41.088}_{-38.344}^{+106.775}$	$-{0.511}_{-0.428}^{+0.423}$
$\mathrm{log}[{B}_{{\rm{S}}}$ (G)]	${12.325}_{-0.163}^{+0.162}$	${12.370}_{-0.151}^{+0.153}$	${0.341}_{-0.123}^{+0.152}$	${0.346}_{-0.129}^{+0.144}$	${22.413}_{-20.764}^{+91.484}$	$-{0.045}_{-0.216}^{+0.220}$
$\mathrm{log}[\dot{E}$ (erg s−1)]	${35.335}_{-0.473}^{+0.488}$	${36.056}_{-0.403}^{+0.399}$	${1.027}_{-0.308}^{+0.395}$	${0.960}_{-0.269}^{+0.344}$	${42.299}_{-39.574}^{+108.965}$	$-{0.721}_{-0.622}^{+0.627}$
$\mathrm{log}[{E}_{{\rm{C}}}$ (GeV)]	${0.371}_{-0.097}^{+0.100}$	${0.328}_{-0.145}^{+0.147}$	${0.209}_{-0.076}^{+0.091}$	${0.314}_{-0.100}^{+0.123}$	${35.299}_{-33.382}^{+106.386}$	${0.043}_{-0.178}^{+0.170}$
$\mathrm{log}[{F}_{\gamma }/{F}_{{\rm{X}}}]$	${3.450}_{-0.301}^{+0.307}$	${2.489}_{-0.619}^{+0.598}$	${0.471}_{-0.222}^{+0.300}$	${1.100}_{-0.454}^{+0.577}$	${27.647}_{-26.357}^{+98.743}$	${0.962}_{-0.683}^{+0.684}$
log[CS]	${1.182}_{-0.110}^{+0.117}$	${0.903}_{-0.199}^{+0.201}$	${0.230}_{-0.101}^{+0.108}$	${0.417}_{-0.144}^{+0.173}$	${27.147}_{-25.653}^{+100.422}$	${0.279}_{-0.221}^{+0.232}$
log[VI]	${1.684}_{-0.047}^{+0.048}$	${1.700}_{-0.058}^{+0.055}$	${0.087}_{-0.038}^{+0.049}$	${0.105}_{-0.046}^{+0.059}$	${3.377}_{-2.247}^{+4.634}$	$-{0.016}_{-0.072}^{+0.074}$
$\mathrm{log}[{\rm{\Gamma }}]$	${0.168}_{-0.052}^{+0.045}$	${0.196}_{-0.052}^{+0.039}$	${0.086}_{-0.042}^{+0.055}$	${0.072}_{-0.040}^{+0.053}$	${6.822}_{-5.822}^{+53.288}$	$-{0.029}_{-0.060}^{+0.060}$
log[FWHM$\cup {{\rm{\Delta }}}_{\gamma }]$	$-{0.432}_{-0.101}^{+0.089}$	$-{0.491}_{-0.121}^{+0.105}$	${0.169}_{-0.082}^{+0.095}$	${0.214}_{-0.097}^{+0.119}$	${9.788}_{-8.777}^{+65.372}$	${0.059}_{-0.131}^{+0.135}$

Note. The statistical parameters obtained from the posterior distributions of the selected physical parameters are as follows: p is the spin period, $\dot{p}$ is the spin-down, ${B}_{{\rm{LC}}}$ is the magnetic field strength at the light cylinder, ${B}_{{\rm{S}}}$ is the surface magnetic field strength, $\dot{E}$ is the spin-down power, ${E}_{{\rm{C}}}$ is the spectral cutoff energy, ${F}_{\gamma }/{F}_{{\rm{X}}}$ is the gamma-ray-to-X-ray flux ratio, CS is the spectral curvature significance, VI is the variability index, Γ is the photon index, and FWHM$\cup {{\rm{\Delta }}}_{\gamma }$ is the union of FWHM and peak separation (${{\rm{\Delta }}}_{\gamma }$). The error bars are the 99% highest posterior density intervals (HPDIs). Because these values are obtained directly from the MCMC samplings, they represent the marginal distributions over the other statistical parameters. The values for ${\mu }_{X}^{{\rm{RQ}}}-{\mu }_{X}^{{\rm{RL}}}$ are directly measured from their marginal distributions.

Download table as:  ASCIITypeset image

We quantify the differences between the groups of RQ and RL pulsars by the distributions of ${\mu }_{X}^{{\rm{RQ}}}-{\mu }_{X}^{{\rm{RL}}}$. If ${\mu }_{X}^{{\rm{RQ}}}-{\mu }_{X}^{{\rm{RL}}}=0$ is excluded from the 99% highest posterior density interval (HPDI), then the two groups of pulsars are significantly different for this physical parameter. Among the 12 physical parameters listed in H17, we found that ${\mu }_{X}^{{\rm{RQ}}}-{\mu }_{X}^{{\rm{RL}}}$ for X being the spin period p, the magnetic field strength at the light cylinder ${B}_{{\rm{LC}}}$, the spin-down power $\dot{E}$, the gamma-ray-to-X-ray flux ratio ${F}_{\gamma }/{F}_{{\rm{X}}}$, and the spectral CS exhibit 99% difference. The marginal distributions of ${\mu }_{X}^{{\rm{RQ}}}-{\mu }_{X}^{{\rm{RL}}}$ for these 5 parameters are shown in Figure 1. Their observed and posterior distributions are plotted in Figures 2 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can check whether the Bayesian inference results make sense by the method of posterior predictive check (PPC). The PPC is a method to compare the predicted data generated by the posterior distributions with the observed data. We randomly choose 500 sets of the statistical parameters from the posteriors of the physical parameters, then we draw 10,000 values from each of the 500 predicted data sets. We then calculate the mean of these 500 data sets and compare the posterior predictive distribution with the mean of the observed data. We found that the means are all centered at the predictive distributions, which implies that they are consistent with each other.

We also verified our results with different prior combinations for Equation (2) by uniform or normal distribution. In the most conservative case of using all uniform distributions, the posteriors of ${\mu }_{X}^{{\rm{RQ}}}$, ${\mu }_{X}^{{\rm{RL}}}$, ${\sigma }_{X}^{{\rm{RQ}}}$, and ${\sigma }_{X}^{{\rm{RL}}}$ are unchanged. In the case of ${\nu }_{X}$, it is unchanged if its value was small (i.e., outliers are important). Its posterior favors large values for the cases that it was large (${\nu }_{X}\gtrsim 30$), which means that the distribution is preferably described by a normal distribution and the shape is unchanged because a t-distribution approximates a normal distribution for ${\nu }_{X}\gtrsim 30$.

2.2. Fraction of Radio-quiet Gamma-ray Pulsars

The detections of R RQ pulsars and $(N-R)$ RL pulsars, where N is the total number of pulsars, form a binomial distribution with success probability as the fraction of RQ pulsars ${f}_{{\rm{RQ}}}$. Thus, the likelihood of detecting R RQ pulsars among N pulsars is

$$\begin{eqnarray}&&{\mathbb{P}}(N,R| {f}_{{\rm{RQ}}})\propto {({f}_{{\rm{RQ}}})}^{R}\times {(1-{f}_{{\rm{RQ}}})}^{N-R}.\end{eqnarray} \tag{ 3 }$$

The equals sign is replaced by the proportional sign because the binomial coefficient can be absorbed into the normalization of the posterior distribution.

We use the non-informative Jeffreys prior for the binomial distribution,

$$\begin{eqnarray}&&{\mathbb{P}}({f}_{{\rm{RQ}}})=\displaystyle \frac{1}{\sqrt{{f}_{{\rm{RQ}}}(1-{f}_{{\rm{RQ}}})}},\end{eqnarray} \tag{ 4 }$$

which is indeed a beta distribution Beta $(\tfrac{1}{2},\tfrac{1}{2})$. The Jeffreys prior represents the fact that we are ignorant about the probability that the model is true (in our case, that ${f}_{{\rm{RQ}}}$ equals certain value), and it has the property that it is independent of the model parameterization.

The posterior distribution ${\mathbb{P}}({f}_{{\rm{RQ}}}| N,R)$ using the data from H17 is displayed in Figure 4. For comparison, the posterior using the data from S16 is also displayed. The frequentist point estimations, ${f}_{{\rm{RQ}}}=35/77=0.455$ for H17 and ${f}_{{\rm{RQ}}}=25/40=0.625$ for S16, are consistent with the 99% HPDIs of the posteriors, ${f}_{{\rm{RQ}}}={0.454}_{-0.159}^{+0.166}$ for H17 and ${f}_{{\rm{RQ}}}={0.627}_{-0.230}^{+0.195}$ for S16. Note that our result indicates that the inferred fraction of RQ pulsars from H17 and S16 are indeed consistent to each other because their HPDIs overlap.

Zoom In Zoom Out Reset image size

We can further calculate the expected ratio of RQ to RL pulsars to be detected given a value of p. The posterior distributions of ${p}^{{\rm{RQ}}}$ and ${p}^{{\rm{RL}}}$ obtained from Section 2.1 by taking the mean values of the distributions of ${\mu }_{p}^{{\rm{RQ}}}$ and ${\mu }_{p}^{{\rm{RL}}}$ are displayed in Figure 5(a). It is observed that RQ pulsars are usually found with larger p. This can be seen more clearly if we plot the relative numbers of RQ to RL pulsars expected to be detected given a value of p in Figure 5(b). The plot indicates that the probabilities of detecting a RQ and RL pulsar is 50/50 at $p\approx 120{\rm{-}}200$ ms.

Zoom In Zoom Out Reset image size

2.3. Constraining the Radio-cone Opening Half-angle and Magnetic Inclination Angle

If we assume that (1) the gamma rays are emitted from the outer gap (e.g., Cheng & Zhang 1998; Takata et al. 2006, 2008), (2) the radio emissions originate from the polar cap region, and (3) there is no intrinsic physical difference between RQ and RL pulsars, then the detection of radio signals depends completely on geometrical effects. Because all gamma-ray pulsars are detected in gamma rays by definition, we know that the LOS from Earth must cut through the outer gap. Therefore we can replace the random variable ${f}_{{\rm{RQ}}}$ in Section 2.2 with the radio-cone opening half-angle, δ, and the magnetic inclination angle, α. The fraction $1-{f}_{{\rm{RQ}}}$ is just the solid angle swapped by the radio cone over 2π. Using simple geometry, we have

$$\begin{eqnarray}{f}_{{\rm{RQ}}}(\alpha ,\delta )=\left\{\begin{array}{ll}\cos (\alpha +\delta ) & \mathrm{for}\ \alpha \lt \delta ,\\ 1-2\sin \alpha \sin \delta & \mathrm{for}\ \delta \leqslant \alpha \leqslant \tfrac{\pi }{2}-\delta ,\\ 1-\cos (\alpha -\delta ) & \mathrm{for}\ \alpha \geqslant \tfrac{\pi }{2}-\delta ,\end{array}\right.\end{eqnarray} \tag{ 5 }$$

and the likelihood of detecting R RQ pulsars among N pulsars is

$$\begin{eqnarray}{\mathbb{P}}(N,R| \alpha ,\delta ) & \propto & {[{f}_{{\rm{RQ}}}(\alpha ,\delta )]}^{R}\\ & & \times {[1-{f}_{{\rm{RQ}}}(\alpha ,\delta )]}^{N-R}.\end{eqnarray} \tag{ 6 }$$

It is expected that the spinning axis and the magnetic axis are randomly placed at the moment when the pulsar was born (e.g., Rookyard et al. 2015b). From symmetry, it is sufficient to only consider a single hemisphere, and the opening angle of the radio cone cannot be larger than 180°, so that we have $(\alpha ,\delta )\in [0,\tfrac{\pi }{2}]$. Therefore, we use a uniform prior for α:

$$\begin{eqnarray}&&\alpha \sim { \mathcal U }\left(0,\displaystyle \frac{\pi }{2}\right).\end{eqnarray} \tag{ 7 }$$

An empirical relation. $\delta \propto {p}^{-n}$, is found by past observations (e.g., Narayan & Vivekanand 1983; Lyne & Manchester 1988; Biggs 1990; Gil et al. 1993; Gil & Han 1996), where $n\approx 1/2$. We adopt the relation

$$\begin{eqnarray}&&\delta =6\buildrel{\circ}\over{.} 3\times {p}^{-1/2},\end{eqnarray} \tag{ 8 }$$

from (Gil & Han 1996) at 1.4 GHz for the outer beam. To obtain a reasonable δ prior, we repeat the BEST technique on the whole population, $p={p}^{{\rm{RQ}}}\cup {p}^{{\rm{RL}}}$, using the same procedure described in Section 2.1, and obtained the posterior distribution for $\mathrm{log}p$. Using Equation (8), we can construct the prior of δ using this t-distribution:

$$\begin{eqnarray}&&\mathrm{log}p\sim { \mathcal T }({\mu }_{p},{\sigma }_{p},{\nu }_{p}).\end{eqnarray} \tag{ 9 }$$

The posterior of this hierarchical Bayesian model is a function of α and p (or δ), conditioning on R and N. From the values of R and N from H17 and S16, we can compute the posteriors ${\mathbb{P}}(\alpha ,\delta | N,R)$ for both studies, which are plotted in Figure 6. The results suggest, based on the current detection statistics of RQ and RL pulsars, assuming Equation (8), that α is skewed toward values larger than 40°, and δ is between 10° and 35°.

Zoom In Zoom Out Reset image size

It may be interesting to remove the dependence of p in δ by inferring $\delta \times {p}^{1/2}$. However, p is not uniquely determined but is instead a random variable, so it is not possible to obtain the distribution for $\delta \times {p}^{1/2}$ because δ depends only on p. The distribution of δ is given once that of p is given. To verify this, we repeated the analysis using $\delta \times {p}^{1/2}$ instead of δ, which results in a point estimation corresponding to $6\buildrel{\circ}\over{.} 3$.

We demonstrated for the first time using Bayesian statistics to study the populations of RQ and RL gamma-ray pulsars. The analysis results are robust to outliers in the pulsar populations. The Bayesian approach has advantages over the conventional frequentist approach when dealing with scarce number of data points, as in the current case of RQ and RL gamma-ray pulsars. We obtained complete information of the posterior distributions for every physical and statistical parameters under study, in the sense that the uncertainties, modality, and correlations between parameters are readily visualized.

We showed that the spin period p, the magnetic field strength at the light cylinder ${B}_{{\rm{LC}}}$, the spin-down power $\dot{E}$, the gamma-ray-to-X-ray flux ratio ${F}_{\gamma }/{F}_{{\rm{X}}}$, and the spectral CS of the two groups of pulsars exhibit significant differences at the 99% level. In Figure 1, it can be seen that the RQ pulsars have p, ${F}_{\gamma }/{F}_{{\rm{X}}}$, and CS larger than those of the RL pulsars, while ${B}_{{\rm{LC}}}$ and $\dot{E}$ of RQ pulsars are smaller than RL pulsars. These behaviors can be explained, because both ${B}_{{\rm{LC}}}$ and $\dot{E}$ scale with ${p}^{-3}$ (Hui et al. 2017).

To further investigate the difference in ${F}_{\gamma }/{F}_{{\rm{X}}}$, we repeat the BEST analysis on the gamma-ray and X-ray luminosities, L_γ and ${L}_{{\rm{X}}}$, and the radiation efficiencies of gamma rays and X rays, ${L}_{\gamma }/\dot{E}$ and ${L}_{{\rm{X}}}/\dot{E}$, respectively, in the base-10 logarithmic scale. For L_γ and ${L}_{{\rm{X}}}$, we found that the ${L}_{{\rm{X}}}$ populations of RQ and RL pulsars are different at the 99% level, while that of L_γ has no significant difference. However, for ${L}_{\gamma }/\dot{E}$ and ${L}_{{\rm{X}}}/\dot{E}$, we found that ${L}_{\gamma }/\dot{E}$ of RQ and RL pulsars are different at the 95% level, while that of ${L}_{{\rm{X}}}/\dot{E}$ has no significant difference. Therefore, the observed difference in ${F}_{\gamma }/{F}_{{\rm{X}}}$ may be a mixed effect on the luminosities and spin-down powers in the gamma-ray and X-ray bands. Note that this effect is also subjected to the uncertainties of the adopted distances.

We showed that the averaged spin period p of the RL pulsars is smaller than that of the RQ pulsars. This result can be understood if the average width of the radio cone of RL pulsars is larger than that of the RQ pulsars. As we adopted in Equation (8), it has been considered that the width of the radio beam is related to the size of the polar cap, whose radius is typically ${R}_{{\rm{P}}}\sim {R}_{{\rm{NS}}}{({R}_{{\rm{NS}}}/{R}_{{\rm{LC}}})}^{1/2}\propto {p}^{-1/2}$, where ${R}_{{\rm{NS}}}$ is the radius of the neutron star, and ${R}_{{\rm{LC}}}={pc}/2\pi$ is the radius of the light cylinder. It has also been suggested that the radio emission altitude relative to the radius of the light cylinder increases with decreasing p, i.e., the emission altitude approaches to the light cylinder for smaller p (Kijak & Gil 2003). These empirical relations suggest that the width of the radio cone increases with decreasing p, hence the pulsar with a shorter period has a large chance to be detected as a RL pulsar. This explains that the average p of the RL pulsars is larger than that of the RQ pulsars. Since ${B}_{{\rm{LC}}}\propto {B}_{{\rm{S}}}{p}^{-3}$ and $\dot{E}\propto {B}_{{\rm{S}}}^{2}{p}^{-4}$, the average ${B}_{{\rm{LC}}}$ and $\dot{E}$ of the RL pulsars are larger than that of the RQ pulsars.

The gamma-ray spectral curvature measured by the Fermi-LAT rules out the classical polar cap scenario (e.g., Daugherty & Harding 1996) and supports the hypothesis that the gamma-ray emission is originated from the outer magnetosphere, e.g., Arons (1983) for the slot gap model, Cheng et al. (1986) for the outer gap model, and Spitkovsky (2006) for the current sheet. However, the spectral behavior above the cutoff energy at around 3 GeV, which decays slower than pure exponential function (Ackermann et al. 2012), has not been fully understood. It has been discussed that the spectral curvature is related to the non-stationary activity of the outer gap emission (Takata et al. 2016) or the magnetic field structure around the light cylinder (Viganò et al. 2015).

To account for the difference of the gamma-ray spectral curvature between the RL and RQ pulsars, we may expect another component that contributes to the high-energy emissions for the RL pulsars. For example, we speculate that inverse Compton (IC) process plays a role in the high-energy photon production of the RL pulsars. For the RL pulsars, which generally have wider radio cones and have emission regions closer to the light cylinder than their RL counterparts, it is highly plausible that part of the radio waves get into the acceleration region around the light cylinder. Because the typical Lorentz factor of the primary electrons/positrons is ${\rm{\Gamma }}\sim 3\times {10}^{7}$, IC scattering of the primary electrons/positrons with radio waves (∼1 GHz) will produce photons in the GeV regime. On the other hand, the probability of radio photons from the RQ pulsars getting into the gap is low. Hence, the spectral curvature of the RQ pulsars could be larger than that of the RL pulsars.

While the observed difference of the spectral curvature can be explained by the IC scenario, this can also be stemmed from observational bias. Different from the RL cases, RQ pulsars can only be detected through a blind search. This leads to the fact that the S/N of RQ pulsars are higher on average. For RL pulsars, with the help of the radio ephemeris, fainter pulsars can also be detected. For these fainter sources, the cutoff may not be constrained. This might lead to an apparently less curved spectra.

The differences between the results using the data from H17 and S16 arise from the fact that they used different pulsar samples. However, we found that their Bayesian inferred fractions of RQ pulsars are indeed consistent with each other. H17 adopted all non-recycled gamma-ray pulsars from 2PC and 3FGL, while S16 employed a blind search for gamma-ray pulsars using the Fermi data alone, resulting in fewer number of data points but plausibly bias-free.

From Figure 6(b), it is observed that the α-δ posterior distribution of S16 is multi-modal, and extends to lower vales of α. This is because S16 has fewer data points but a higher fraction of RQ pulsars. Because the ranges of plausible values of δ are similar for H17 and S16, smaller α results in smaller solid angle on the sky, i.e., more RQ pulsars. Moreover, it can be seen that the posterior of S16 is more heavily distributed in regions where $\alpha +\delta \gt 90^\circ$. The coverage of the radio cone increases with α and δ, but decreases when $\alpha +\delta \gt 90^\circ$ because the projections of the bipolar cones on opposite hemispheres overlap. This again implies a higher fraction of RQ pulsars than RL pulsars.

Using Bayesian inference, assuming pure geometrical effects, we showed that the distribution of the magnetic inclination angle α is skewed toward larger values, which does not exhibit the reported unexpected skew of α toward smaller values from Rookyard et al. (2015a), and is consistent with the conclusion of Rookyard et al. (2015b) that α should be skewed toward larger values. On the other hand, the values of the radio-cone half-angle δ are constrained to be within 10° and 35°. Because we used Equation (8) and considered pulsars with rotation period within 0.03 and 0.2 s, we should obtain δ to be within 10° and 35°. However, using Bayesian statistics instead of conventional frequentist statistics gives us not only a range estimation for δ, but also one for α and their joint posterior, which gives their uncertainties and correlations for free. From the gamma-ray emission point of view, larger inclination angles are preferentially detected by the observations (c.f., Takata et al. 2011; Watters & Romani 2011).

In the future, the Bayesian inference results shown here could be updated when more RQ and RL gamma-ray pulsars are detected. The uncertainties (expressed in HPDIs) will be reduced because the observed distributions in Figures 2 and 3 will be less patchy and the posterior distributions can be constrained more precisely using more data. To avoid incorrect and/or ambiguous results using conventional frequentist fitting methods and tests, Bayesian inference could be and should be used in the observational and statistical studies of other compact objects and in general every fields of high-energy astrophysics, e.g., RL versus millisecond pulsars, and Redback versus Black Widow pulsars, etc.

We thank K. S. Cheng for valuable theoretical discussions. H.F.Y. thanks J. Michael Burgess, Felix Ryde, and Yan-Ting Lam for insightful discussions on Bayesian statistics. H.F.Y. acknowledges support from the Swedish National Space Board and the Swedish Research Council (Vetenskapsrådet). H.F.Y. is supported by the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. C.Y.H. is supported by the National Research Foundation of Korea through grant 2016R1A5A1013277. A.K.H.K. is supported by the Ministry of Science and Technology of the Republic of China (Taiwan) through grants 105-2112-M-007-033-MY2 and 105-2119-M-007-028-MY3. J.T. is supported by NSFC grants of Chinese Government through grants 11573010, U1631103, and 11661161010.

Facility: Fermi Gamma-ray Space Telescope/Large Area Telescope - (Atwood et al. 2009).

Software: PyMC3 (Salvatier et al. 2016).