# On the hard $\gamma$-ray spectrum of the supernova remnant G106.3+2.7

Yiwei Bao  Department of Astronomy, Nanjing University, 163 Xianlin Avenue,
Nanjing 210023, China  Yang Chen  Department of Astronomy, Nanjing University,
163 Xianlin Avenue, Nanjing 210023, China  Key Laboratory of Modern Astronomy
and Astrophysics, Nanjing University, Ministry of Education, Nanjing, China

###### Abstract

Tibet AS$\gamma$ experiment has measured $\gamma$-ray flux of supernova
remnant G106.3+2.7 up to 100 TeV, suggesting it to be potentially a
“PeVatron”. Challenge arises when the hadronic scenario requires a hard proton
spectrum (with spectral index $\approx 1.8$), while the diffusive shock
acceleration can only predict a soft (with spectral index $\geq 2$) proton
spectrum. In this paper, we explore an alternative scenario to explain the
$\gamma$-ray spectrum of G106.3+2.7 within the current understanding of
acceleration and escape processes. We consider that the cosmic rays upstream
of the shock are scattered by the turbulence driven via Bell instability. The
resulting hadronic $\gamma$-ray spectrum is novel, dominating the contribution
to the emission above 10 TeV, and can explain the bizzare broadband spectrum
of G106.3+2.7 in combination with leptonic emission from the supernova
remnant.

ISM: supernova remnants — ISM: individual objects (G106.3+2.7) — diffusion —
(ISM:) cosmic rays

††journal: ApJL

##  1 INTRODUCTION

Galactic cosmic rays (GCRs) are mostly charged particles (mainly protons) with
energy up to the “knee” ($\sim 1$ PeV). Based on energetic arguments,
supernova remnants (SNRs) are usually believed to be the accelerators of GCRs.
However, although a large number of SNRs have been detected in $\gamma$-rays,
none of these SNRs has been confirmed to be a PeV particle accelerator, as
called “PeVatron” (see e.g., Bell et al., 2013). Discovered in the DRAO
Galactic-plane survey (Joncas & Higgs, 1990), G106.3+2.7 is an cometary SNR
with a tail in the southwest and a compact head (containing PSR J2229+6114) in
the northeast. Recently, HAWC and AS$\gamma$ experiments have reported the
$\gamma$-ray spectrum of SNR G106.3+2.7 above 40 TeV, arguing that SNR
G106.3+2.7 (at a distance of 800 pc away from Earth, Kothes et al., 2001) is a
promising PeVatron candidate (Albert et al., 2020; Amenomori et al., 2021). In
particular, Amenomori et al. (2021) for the first time measured the
$\gamma$-ray flux up to 100 TeV, finding that the centroid of $\gamma$-ray
emissions is away from PSR J2229+6114 by 0.44∘ ($\sim 6$ pc at a distance of
800 pc), and is well correlated with a molecular cloud (MC), implying that the
high-energy $\gamma$-rays can be of hadronic origins.

Both leptonic and hadronic models have been proposed to give a plausible
explanation to the spectral energy distribution (SED) of the SNR. In the
leptonic scenario, the electrons are suggested to be transported to its
current position from the pulsar wind nebula (PWN, Amenomori et al., 2021; Liu
et al., 2020), or be accelerated by the SNR shock directly (Amenomori et al.,
2021); in the hadronic scenario, the protons are suggested to be accelerated
by the SNR shock in earlier ages, or re-accelerated by the PWN adiabatically
(Amenomori et al., 2021). In order to figure out the origin of $\gamma$-rays,
Ge et al. (2020) separate the PWN-dominated X-ray emitting region in the
northeast from the other (southwestern) part of the SNR using XMM-Newton and
Suzaku observations. They found that a pure leptonic model can hardly fit the
X-ray and $\gamma$-ray spectral data of the latter part of SNR simultaneously.
Therefore the SED can only be explained with a hadronic (with a proton
spectral index $\approx 1.8$, Amenomori et al., 2021) or a hybrid leptonic-
hadronic (with a proton index $\approx 1.5$, Ge et al., 2020) model. However,
challenge arises when hard proton spectrum is required in both hadronic and
hybrid leptonic-hadronic models, while numerical simulations show that
diffusive shock acceleration can only give rise to a soft proton spectrum
(with spectral index $\geq 2$, see e.g., Caprioli et al., 2020).

Amenomori et al. (2021) suggest that the very hard proton spectrum can be
obtained in very efficient acceleration (which seems to be an extreme case),
and after a very slow particle diffusion. We here explore an alternative self-
consistent scenario in which the $\gamma$-ray spectrum can be explained more
naturally. The magnetic field in the upstream of the SNR shock is believed to
be amplified via the non-resonant Bell instability (Bell, 2004; Amato & Blasi,
2009): protons escaping from the upstream of the shock can drive non-resonant
turbulence whose scale is smaller than the Lamour radius of the escaped
particles, and the relatively-low-energy particles are thus diffused by the
magnetic field amplified by the non-resonant turbulence (for reviews, see
e.g., Blasi, 2013). Based on the numerical simulations (Bell, 2004; Bell et
al., 2013), we calculate the escape process from the first principle instead
of a model with apriori phenomenological assumptions. The resulting proton
spectrum is novel and can explain the hard $\gamma$-ray spectrum well. The
model is described in §2 and the SED of the SNR G106.3+2.7 is fitted in §3,
discussion is presented in §4 and the conclusion is drawn in §5.

##  2 MODEL DESCRIPTION

For simplicity, we follow the approximation that the global maximum energy is
reached at the beginning of the Sedov phase in which $R_{\rm sh}\propto
t^{2/5}$ (see e.g., Ohira et al., 2012). We use the canonical diffusive shock
acceleration scenario in which the spectrum of accelerated particles $f_{\rm
acc}\propto E^{-2}$. The maximum proton energy can be given by (Bell et al.,
2013)

$E_{\rm max}(t)=230\,n_{\rm
e}^{1/2}\left(\frac{\eta}{0.03}\right)\left(\frac{v_{\rm sh}}{10^{4}{\rm
kms^{-1}}}\right)^{2}\left(\frac{R_{\rm sh}}{\rm pc}\right)\,{\rm TeV},$ (1)

where $\eta$ is the acceleration efficiency, $n_{\rm e}$ the electron number
density of the interstellar medium (ISM), ${v_{\rm sh}}$ the velocity of the
shock, and ${R_{\rm sh}}$ the radius of the shock.

Following Cardillo et al. (2015), we calculate the spectrum of escaped protons
in the Sedov phase. In the upstream, the protons are scattered by the wave
generated via escaping of higher-energy protons, and therefore there is no
wave upstream to scatter the protons with the highest-energy $E_{\rm max}(t)$
(Bell, 2004; Bell et al., 2013). Hence, the latter escape the system quasi-
ballistically at a speed $\sim c$, inducing a current $j_{\rm CR}=ev_{\rm
sh}\pi p_{\rm max}^{3}f_{0}(p_{\rm max})$ at the shock, where $f_{0}$ is the
isotropic part of cosmic ray (CR) distribution function in momentum space at
the shock, and $p_{\rm max}$ the maximum momemtum of the non-thermal protons
(Bell et al., 2013). On the other hand, the CR pressure at the shock (Bell et
al., 2013)

$P_{\rm CR}=\frac{4\pi}{3}cp_{\rm max}^{4}f_{0}(p_{\rm max}){\rm
ln}\left(\frac{p_{\rm max}}{p_{0}}\right),$ (2)

where $p_{0}$ is the minimum momentum of non-thermal protons. Following Bell
et al. (2013), we assume that the CR pressure is a fixed ratio of the ram
pressure $\rho v_{\rm sh}^{2}$ (where $\rho$ is the mass density of the
preshock ISM), so $j_{\rm CR}$ can be calculated to be (the coefficient is a
constant which is not important in the calculation)

$j_{\rm CR}\propto\frac{e\rho v_{\rm sh}^{3}}{E_{\rm max}{\rm ln}\left({E_{\rm
max}}/{E_{0}}\right)},$ (3)

where $E_{\rm max}$ and $E_{0}$ are maximum and minimum energies of protons,
respectively.

The total number of the escaped protons with energy $E_{\rm max}(t)$ can thus
be evaluated via

$eN_{\rm esc}(E_{\rm max})\,{\rm d}E_{\rm max}=4\pi R^{2}_{\rm sh}j_{\rm
CR}\,{\rm d}t.$ (4)

So it can be expressed as (see Cardillo et al., 2015, for derivation)

$N_{\rm esc}(E_{\rm max})\propto\rho v^{2}_{\rm sh}R^{3}_{\rm sh}E^{-2}_{\rm
max}\left[\frac{1+{\rm ln}(E_{\rm max}/E_{0})}{{\rm ln}(E_{\rm
max}/E_{0})^{2}}\right].$ (5)

##  3 Application to SNR G106.3+2.7

To characterize the $\gamma$-ray radiation from SNR G106.3+2.7, we consider a
truncated molecular cloud cone subtending solid angle $\Omega$ at the SNR with
inner radius $R_{1}$ and outer radius $R_{2}$. and The diffusion coefficient
near the SNR is adopted to be $D(E)=D_{\rm ISM}/K_{\rm sl}$, where $D_{\rm
ISM}=4.16\times 10^{28}\left(E/4{\ \rm GeV}\right)^{0.5}$ cm2 s-1 (Yuan et
al., 2017, model DR2), and $K_{\rm sl}$ is the modification factor.

The diffusion equation for protons with energy $E$ escaping the shock when the
shock radius was $R_{\rm sh}$ at time $T_{\rm esc}$ writes

$\frac{\partial}{\partial t}f(E,r,t)=\\\
\frac{D(E)}{r^{2}}\frac{\partial}{\partial
r}\left[r^{2}\frac{\partial}{\partial r}f(E,r,t)\right]+Q,$ (6)

where

$Q=\frac{N_{\rm esc}(E)}{4\pi R_{\rm sh}^{2}}\delta(r-R_{\rm
sh})\delta(t-T_{\rm esc})=\frac{C_{\rm Q}\rho v^{2}_{\rm sh}R^{3}_{\rm
sh}E^{-2}_{\rm max}}{4\pi R^{2}_{\rm sh}}\left[\frac{1+{\rm ln}(E_{\rm
max}/E_{0})}{{\rm ln}(E_{\rm max}/E_{0})^{2}}\right]\delta(r-R_{\rm
sh})\delta(t-T_{\rm esc}),$ (7)

with $C_{\rm Q}$ a constant to be considered later.

Defining a new function $F=Rf(E,r,t)$, Equation 6 reads

$\frac{\partial}{\partial t}F=D(E)\frac{\partial^{2}F}{\partial
r^{2}}+Q(E,r,t)r$ (8)

According to Atoyan et al. (1995), the Green function for the Equation 8 is

$G(E,r,t;\xi,T_{\rm esc})=\frac{1}{\sqrt{\pi}R_{\rm
dif}}\left\\{\exp\left[-\left(\frac{r-\xi}{R_{\rm
dif}}\right)^{2}\right]-\exp\left[-\left(\frac{r+\xi}{R_{\rm
dif}}\right)^{2}\right]\right\\},$ (9)

where $R_{\rm dif}=2\sqrt{D(E)(t-T_{\rm esc})}$ is the diffusion length scale.
The solution of Equation 6 is thus

$f(E,r,t)=\frac{1}{r}\int_{0}^{\infty}Q(r-\xi)(r-\xi)G(E,r,t;R_{\rm sh},T_{\rm
esc}){\rm d}\xi.$ (10)

At time $T_{\rm age}$, the differential number of protons with energy $E$
lying inside the concentric spherical molecular cloud shell can thus be
calculated to be

$N_{\rm CR,MC}=\int_{R_{1}}^{R_{2}}f(E,r,t)\Omega r^{2}\,{\rm d}r.$ (11)

Finally the $\gamma$-ray luminosity is evaluated via the cross-section
presented in Kamae et al. (2006). In addition to the hadronic component, we
also add a leptonic component contributed by the electrons accelerated by the
SNR shock. The electron spectrum is constrained by the radio (dominated by the
SNR Pineault & Joncas, 2000) and X-ray (Ge et al., 2020; Fujita et al.,
2021)111The former X-ray flux is extracted from the region essentially
overlapping the molecular cloud, and the latter flux from the regions covering
most of the SNR radiations. We approximate electron spectrum to be a broken
power-law

$\frac{dN_{\rm e}}{dE}\propto\begin{cases}E^{-\alpha_{1}}&{E\leq E_{\rm b}}\\\
E^{-\alpha_{2}}&{E_{\rm b}<E<E_{\rm max,global}}.\end{cases}$ (12)

We fit the SED as is plotted in Figure 1, with the parameters listed in Table
1, where $R_{\rm now}$ is the radius of the SNR at present, $E_{\rm SN}$ the
explosion energy of the SNR, $M_{\rm ej}$ the ejecta mass, $B$ the magnetic
field strength, and $W_{\rm e}$ the total electron energy. $T_{\rm CMB}$,
$T_{\rm FIR}$, and $T_{\rm NIR}$ are the temperatures of the CMB, far
infrared, and near infrared photons, respectively; $u_{\rm CMB}$, $u_{\rm
FIR}$, and $u_{\rm NIR}$ are energy densities of the CMB, far infrared, and
near infrared photons, respectively. We find that the parameters are
insensitive to X-ray flux, since we can fit the X-ray data from Ge et al.
(2020) and Fujita et al. (2021) with very similar parameters.

Since the hadronic $\gamma$-ray luminosity $L_{\rm pp}\propto\Omega n_{\rm
MC}$ (where $n_{\rm MC}$ is the number density of the MC, a free parameter),
and $L_{\rm pp}\propto E_{\rm esc}\propto C_{\rm Q}$, $L_{\rm pp}\propto
C_{\rm Q}n_{\rm MC}\Omega$. Therefore, we only list $C_{\rm Q}n_{\rm
MC}(\Omega/4\pi)$ as a single parameter in Table 1. At each time $t$, the
protons with energy less than $E_{\rm max}(t)$ are trapped by the turbulence
driven by the escaping protons with energy $E_{\rm max}(t)$ via Bell
instability. At present, protons with energy $\lesssim{58}$ TeV are still
trapped by the turbulence. Therefore the non-thermal proton spectrum in the MC
is almost a $\delta$ function, leading to a novel hadronic $\gamma$-ray
spectrum. Dominating the $\gamma$-ray flux in $\lesssim 500$ GeV, the leptonic
component (the cyan line in Figure 1) is consistent with the Fermi
observations (Xin et al., 2019) and X-ray observations (Ge et al., 2020;
Fujita et al., 2021). Meanwhile, the hadronic $\gamma$-ray radiation
contributed by pp interactions (the yellow line) has a very hard spectrum
(d$N/{\rm d}E\propto E^{-1}$) in $\lesssim 500$ GeV. The spectrum above 10 TeV
is dominated by the hadronic component, and the whole $\gamma$-ray spectrum
can be explained by the hybrid model naturally.

##  4 discussion

There remains some uncertainty about the age of the SNR: although the
characteristic age of PSR J2229+6114 ($\approx 10^{4}$ yr) is usually adopted
to be the age of the SNR, sometimes the SNR is suggested to be very young
($\sim 1$ kyr, see e.g., Albert et al., 2020). With an adiabatic model
($R_{\rm sh}\approx 1.2(E_{\rm SN}/\rho_{\rm ISM})^{1/5}t^{2/5}$, where
$\rho_{\rm ISM}$ is the density of the ISM), Kothes et al. (2001) estimated
that $E_{\rm SN}\approx 7\times 10^{49}$ if the age of the SNR were $\approx
10$ kyr. However, as is estimated in Cardillo et al. (2015), Equation 1
indicates that the global maximum energy $E_{\rm max,global}\propto E_{\rm
SN}$, and SNRs with low $E_{\rm SN}$ can hardly accelerate CRs to $\sim
10^{2}$ TeV. Although the dilemma may be solved with asymmetry incoporated, in
this paper we consider a spherically symmetric scenario in which $E_{\rm
SN}=10^{51}$ erg (the canonical value) and $T_{\rm age}=1000$ yr. Assuming a
canonical braking index $n=3$, the real age of the pulsar writes $T_{\rm
age}=\tau_{\rm c}[1-(P_{0}/P_{\rm now})^{2}]$, where $\tau_{c}$ is the
characteristic age of the pulsar, $P_{0}$ the inital spin period of the
pulsar, and $P_{\rm now}$ the spin period of the pulsar at present. For PSR
J2229+6114, $P_{\rm now}=50$ ms (Halpern et al., 2001), and $T_{\rm age}$ can
be as small as $\sim 1$ kyr if $P_{0}$ is close to $P_{\rm now}$. The $P_{0}$
assumed here is possible since in pulsar evolution models, $P_{0}$ is usually
suggested to be in a wide range from $\sim 4$ ms to $\sim 400$ ms (see e.g.,
Arzoumanian et al., 2002; Faucher-Giguère & Kaspi, 2006; Popov et al., 2010).
Instances in which $P_{0}\approx P_{\rm now}$ and $T_{\rm age}\ll\tau_{c}$ can
also be found in the catalog: (a) CCO 1E 1207.4$-$5209 inside SNR G296.5+10.0
has $P_{0}\approx P_{\rm now}=424$ ms, leading to $\tau_{\rm c}>27$ Myr, which
exceeds the age of the SNR by 3 orders of magnitude (Gotthelf & Halpern,
2007), and (b) PSR J1852$-$0040 is suggested to have $P_{0}\approx P_{\rm
now}\approx 10^{2}$ ms, giving rise to $\tau_{\rm c}\approx 2\times 10^{8}$
yr, 4 orders of magnitude larger than $T_{\rm age}\sim 5\times 10^{3}$ yr
measured from the observation on its associated SNR (Gotthelf et al., 2005).

On the other hand, as can be seen from Equation 1, $E_{\rm max}$ is determined
by the $R_{\rm sh}$ and $v_{\rm sh}$ which can be constrained by observations
directly, irrespective of $T_{\rm age}$: Ge et al. (2020) showed that $v_{\rm
sh}$ is at least 3000 km s-1 in the southwest region based on the absence of
thermal X-ray emissions, while $R_{\rm sh}\approx 6$ pc can be estimated via
radio observations (see e.g., Amenomori et al., 2021, for the multi-wave
length morphology). Although $T_{\rm age}$ does affect $R_{\rm dif}$, its
impacts can be compensated by $K_{\rm sl}$ which is a free parameter. Hence,
for simplicity we only consider a spherically symmetric SNR evolving following
the well-known model222See Leahy & Williams (2017) for a fast python
calculator on SNR evolution (Truelove & McKee, 1999).

Our particle escaping scenario has an interesting characteristic: in a
spherically symmetric system, there will be a sharp cutoff in the $\gamma$-ray
spectrum. Meanwhile, smoother break can be predicted if asymmetric morphology
is taken into account: asymmetric shock velocity distribution gives rise to
asymmetric $E_{\rm max,global}$ distribution, and an integration over the
whole SNR consequently smoothen the sharp cutoff, leading to a break in the
$\gamma$-ray spectrum. The escape scenario presented in this paper can
coexists with the former continuous escape scenarios. During the lifetime of
SNR, the shock may break and the CRs can thus escape with continuous power-law
spectra, which is the case as previously adopted (see e.g., Li & Chen, 2010;
Ohira et al., 2011). In the SNR
catalog333http://snrcat.physics.umanitoba.ca/SNRtable.php (Ferrand & Safi-
Harb, 2012), there are only a few SNRs with $\gamma$-ray emissions in
$\lesssim 10$ TeV are confirmed to be of hardronic origin, and our novel
hadronic spectrum may be found in more SNRs with the help of LHASSO,
AS$\gamma$ and HAWC experiments in future.

##  5 conclusion

In this paper, we explore an alternative plausible scenario to explain the
bizzare $\gamma$-ray spectrum of SNR G106.3+2.7 within the current
acceleration theory. Apart from the common leptonic component which can be
constrained by radio and X-ray observation, we add a novel hadronic component.
The hadronic component arises from the escape scenario proposed by Bell (2004)
and Cardillo et al. (2015), in which the CRs at the shock are diffused the by
turbulence generated by higher energy CRs escaping upstream. Consequently, at
time $t$, only CRs with energy $E_{\rm max}(t)$ can escape upstream, while CRs
with lower energy are confined. The resulting hadronic spectrum is quite hard,
and can explain the bizzare $\gamma$-ray spectrum of SED of SNR G106.3+2.7
well.

As can be seen from our calculation, although SNR G106.3+2.7 is sometimes
explained as a “’PeVatron”, current theory can also give rise to a hard
$\gamma$-ray spectrum extending to 100 TeV with $E_{\rm max,global}\approx
4\times 10^{2}$ TeV. Hence detection of 100 TeV $\gamma$-ray emissions is a
necessary condition for “PeVatron” rather than a sufficient condition. More
observations are required to confirm whether G106.3+2.7 is a “PeVatron” or
not.

Figure 1: SED of the emission from SNR G106.3+2.7. The CGPS data are taken from Pineault & Joncas (2000), XMM-Newton data from Ge et al. (2020), Suzaku data from Fujita et al. (2021), Fermi-LAT data from Xin et al. (2019), VERITAS data from Albert et al. (2020), HAWC data from Albert et al. (2020), and AS$\gamma$ data from Amenomori et al. (2021). Table 1: Fitting Parameters Parameter | Model A | Model B
---|---|---
$T_{\rm age}$ (yr) | 1000 | …
$R_{\rm now}$ (pc) | 5.24 | …
$E_{\rm max,global}$ (TeV) | 427 | …
$E_{0}$ (GeV) | 1.0 | …
$R_{1}$ (pc) | 7 | …
$R_{2}$ (pc) | 10 | …
$E_{\rm SN}$ (erg) | $10^{51}$ | …
$d$ (pc) | 800 | …
$M_{\rm ej}$ ($M_{\odot}$) | 1.0 | …
$n_{\rm e}$ (cm-3) | 0.6 | …
$\eta$ | 0.06 | …
$n_{\rm MC}C_{\rm Q}\left(\Omega/4\pi\right)$ (erg cm-3 s-1) | $8\times 10^{-43}$ | $6\times 10^{-43}$
$K_{\rm sl}$ | 100 | …
$E_{\rm b}$ (TeV) | 12 | 16.5
$\alpha_{1}$ | 2.3 | …
$\alpha_{2}$ | 3.8 | 3.45
$B$ ($\mu$G) | 4.6 | …
$W_{\rm e}$ (erg) | $6.8\times 10^{47}$ | …
$T_{\rm CMB}$ (K) | 2.73 | …
$u_{\rm CMB}$ (eV cm-3) | 0.25 | …
$T_{\rm FIR}$ (K) | 25 | …
$u_{\rm FIR}$ (eV cm-3) | 0.2 | …
$T_{\rm OPT}$ (K) | 3000 | …
$u_{\rm OPT}$ (eV cm-3) | 0.3 | …

We are in debt to Pasquale Blasi for fruitful discussion, Yutaka Fujita, Ruo-
Yu Liu, Qiang Yuan and Xiangdong Li for helpful comments. This work is
supported by the National Key R&D Program of China under grants
2017YFA0402600, NSFC under grants 11773014, 11633007, 11851305 and U1931204.

## References

  * Amenomori et al. (2021) Amenomori, M., Bao, Y.W. et al. (Tibet ASg collaboration) 2021, Nat Astron doi:10.1038/s41550-020-01294-9
  * Albert et al. (2020) Albert, A., Alfaro, R., Alvarez, C., et al. 2020, ApJ, 896, L29. doi:10.3847/2041-8213/ab96cc
  * Amato & Blasi (2009) Amato, E. & Blasi, P. 2009, MNRAS, 392, 1591. doi:10.1111/j.1365-2966.2008.14200.x
  * Arzoumanian et al. (2002) Arzoumanian, Z., Chernoff, D. F., & Cordes, J. M. 2002, ApJ, 568, 289. doi:10.1086/338805
  * Atoyan et al. (1995) Atoyan, A. M., Aharonian, F. A., & Völk, H. J. 1995, Phys. Rev. D, 52, 3265. doi:10.1103/PhysRevD.52.3265
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550. doi:10.1111/j.1365-2966.2004.08097.x
  * Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., et al. 2013, MNRAS, 431, 415. doi:10.1093/mnras/stt179
  * Blasi (2013) Blasi, P. 2013, A&A Rev., 21, 70. doi:10.1007/s00159-013-0070-7
  * Caprioli et al. (2020) Caprioli, D., Haggerty, C. C., & Blasi, P. 2020, ApJ, 905, 2. doi:10.3847/1538-4357/abbe05
  * Cardillo et al. (2015) Cardillo, M., Amato, E., & Blasi, P. 2015, Astroparticle Physics, 69, 1. doi:10.1016/j.astropartphys.2015.03.002
  * Faucher-Giguère & Kaspi (2006) Faucher-Giguère, C.-A. & Kaspi, V. M. 2006, ApJ, 643, 332. doi:10.1086/501516
  * Ferrand & Safi-Harb (2012) Ferrand, G. & Safi-Harb, S. 2012, Advances in Space Research, 49, 1313. doi:10.1016/j.asr.2012.02.004
  * Fujita et al. (2021) Fujita, Y., Bamba, A., Nobukawa, K. K., et al. 2021, arXiv:2101.10329
  * Ge et al. (2020) Ge, C., Liu, R.-Y., Niu, S., et al. 2020, arXiv:2012.11531
  * Gotthelf & Halpern (2007) Gotthelf, E. V. & Halpern, J. P. 2007, ApJ, 664, L35. doi:10.1086/520637
  * Gotthelf et al. (2005) Gotthelf, E. V., Halpern, J. P., & Seward, F. D. 2005, ApJ, 627, 390. doi:10.1086/430300
  * Halpern et al. (2001) Halpern, J. P., Camilo, F., Gotthelf, E. V., et al. 2001, ApJ, 552, L125. doi:10.1086/320347
  * Joncas & Higgs (1990) Joncas, G. & Higgs, L. A. 1990, A&AS, 82, 113
  * Kamae et al. (2006) Kamae, T., Karlsson, N., Mizuno, T., et al. 2006, ApJ, 647, 692. doi:10.1086/505189
  * Kothes et al. (2001) Kothes, R., Uyaniker, B., & Pineault, S. 2001, ApJ, 560, 236. doi:10.1086/322511
  * Leahy & Williams (2017) Leahy, D. A. & Williams, J. E. 2017, AJ, 153, 239. doi:10.3847/1538-3881/aa6af6
  * Li & Chen (2010) Li, H. & Chen, Y. 2010, MNRAS, 409, L35. doi:10.1111/j.1745-3933.2010.00944.x
  * Liu et al. (2020) Liu, S., Zeng, H., Xin, Y., et al. 2020, ApJ, 897, L34. doi:10.3847/2041-8213/ab9ff2
  * Ohira et al. (2011) Ohira, Y., Murase, K., & Yamazaki, R. 2011, MNRAS, 410, 1577. doi:10.1111/j.1365-2966.2010.17539.x
  * Ohira et al. (2012) Ohira, Y., Yamazaki, R., Kawanaka, N., et al. 2012, MNRAS, 427, 91. doi:10.1111/j.1365-2966.2012.21908.x
  * Pineault & Joncas (2000) Pineault, S. & Joncas, G. 2000, AJ, 120, 3218. doi:10.1086/316863
  * Popov et al. (2010) Popov, S. B., Pons, J. A., Miralles, J. A., et al. 2010, MNRAS, 401, 2675. doi:10.1111/j.1365-2966.2009.15850.x
  * Truelove & McKee (1999) Truelove, J. K. & McKee, C. F. 1999, ApJS, 120, 299. doi:10.1086/313176
  * Xin et al. (2019) Xin, Y., Zeng, H., Liu, S., et al. 2019, ApJ, 885, 162. doi:10.3847/1538-4357/ab48ee
  * Yuan et al. (2017) Yuan, Q., Lin, S.-J., Fang, K., et al. 2017, Phys. Rev. D, 95, 083007. doi:10.1103/PhysRevD.95.083007
  * Zavlin et al. (2000) Zavlin, V. E., Pavlov, G. G., Sanwal, D., et al. 2000, ApJ, 540, L25. doi:10.1086/312866