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

The Tibet AS$\gamma$ experiment has measured $\gamma$-ray flux of supernova remnant G106.3+2.7 up to 100 TeV, suggesting it {being} potentially a"PeVatron". Challenge arises when the hadronic scenario requires a hard proton spectrum (with spectral index $\approx 1.8$), while {usual observations and numerical simulations prefer} a soft proton spectrum {(with spectral index $\geq 2$)}. 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 ray {particles} 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 bizarre broadband spectrum of G106.3+2.7 in combination with leptonic emission from the remnant.


INTRODUCTION
Galactic cosmic rays (CRs) are mostly charged particles (mainly protons) with energy up to the so called "knee" (∼ 1-10 PeV). Based on energetic arguments, supernova remnants (SNRs) are usually believed to be the accelerators of Galactic CRs. However, although a large number of SNRs have been detected in γ-rays, none of them 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 a cometary SNR with a tail in the southwest and a compact head (containing PSR J2229+6114) in the northeast, at a distance of 800 pc away from Earth (Kothes et al. 2001). Recently, HAWC and the Tibet ASγ experiments have reported the γ-ray spectrum of SNR G106.3+2.7 above 40 TeV, arguing that the remnant is a promising PeVatron candidate (Albert et al. 2020;Tibet ASγ Collaboration 2021). In particular, Tibet ASγ Collaboration (2021) for the first time measured the γ-ray flux up to 100 TeV, finding that the centroid of γ-ray emissions deviates from the pulsar at a confidence of 3.1σ, and is well correlated with a molecular cloud (MC). The offset of the γ-ray emission centroid from PSR J2229+6114 is measured to be 0.44 • (∼ 6 pc at a distance of 800 pc).
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, Tibet ASγ Collaboration 2021; Liu et al. 2020), or be accelerated by the blast wave directly (Tibet ASγ Collaboration 2021; Ge et al. 2021); in the hadronic scenario, the protons are suggested to be accelerated by the blast wave in earlier ages, or re-accelerated by the PWN adiabatically (Ohira et al. 2018;Tibet ASγ Collaboration 2021). The large offset between PSR J2229+6114 and the γ-ray emission centroid indicates that the γ-rays are more likely to be mainly contributed from the particles accelerated by the blast wave in the southwestern "tail" region (Tibet ASγ Collaboration 2021). In order to figure out the origin of the γ-rays, Ge et al. (2021) separate the PWN-dominated X-ray emitting region in the northeast from the other (the "tail") part of the SNR using XMM-Newton and Suzaku observations. They found that a pure leptonic model can hardly fit the radio, X-ray and γ-ray spectral data of the "tail" region simultaneously. Therefore, there must be a hadronic component in the γ-ray spectrum, and the SED can only be explained with a hadronic (with a proton spectral index ≈ 1.8, Tibet ASγ Collaboration 2021) or a leptonic-hadronic hybrid (with a proton index ≈ 1.5, Ge et al. 2021) model. However, challenge arises when hard proton spectrum is required in both hadronic and leptonic-hadronic hybrid models, while numerical simulations show that diffusive shock acceleration can only give rise to a soft proton spectrum (with spectral index ≥ 2, see e.g., Caprioli et al. 2020).
Tibet ASγ Collaboration (2021) suggest that the very hard proton spectrum can be formed 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 γ-ray spectrum can be explained more naturally. The magnetic field in the upstream of the blast wave 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 γ-ray spectrum well in combination with leptonic emission from the SNR. The model is described in §2 and the SED of the remnant is fitted in §3, discussion is presented in §4 and the conclusion is drawn in §5.

MODEL DESCRIPTION
For simplicity, we follow the approximation that the global maximum energy E max,global is reached at the beginning of the Sedov phase in which R sh ∝ t 2/5 (see e.g., Ohira et al. 2012). The maximum proton energy can be given by (Bell et al. where η is the acceleration efficiency, n e the electron number density of the interstellar medium (ISM, we assume that the ISM near the SNR are fully ionized), v sh the velocity of the shock, and R sh the radius of the shock. Following Cardillo et al. (2015), we calculate the spectrum of escaped CR protons in the Sedov phase. In the upstream, the protons are scattered by the wave generated via escaping of the higherenergy protons, and therefore there is no wave upstream to scatter the protons with the highest energy E esc (t) (Bell 2004;Bell et al. 2013). Hence, the latter escape the system quasi-ballistically at a speed ∼ c, inducing a current j CR = n CR ev sh at the shock (Bell 2004, where n CR represents the number density of the CRs). The differential number of the escaped protons with energy E esc (t) can thus be evaluated via We further assume that the CR pressure at the shock is a fixed fraction ξ of the ram pressure α is the power-law index of the parent CR spectrum at the shock, and E 0 the minimum energy of protons. Finally, N esc can be expressed as (see Cardillo et al. 2015, for derivation) 3. APPLICATION TO SNR G106.3+2.7

γ-ray spectroscopic luminosity
Since kinematic/physical signature of direct MC-SNR contact seems inconclusive yet (Liu et al. 2021), the MC may only be illuminated by the escaped protons. For such a scenario, the MC is sometimes approximated as a truncated cone (see e.g., Li & Chen 2012;Celli et al. 2019). Here the MC is assumed to subtend a solid angle Ω at the SNR center with an inner radius R 1 and an outer radius R 2 . The diffusion coefficient near the SNR is adopted to be D(E) = D 100 TeV (E/100 TeV) δ , where D 100 TeV is the diffusion coefficient for particles at 100 TeV, and δ the energy dependence index.
The diffusion equation for the protons with energy E esc escaping the shock when the shock radius was R sh at time T esc writes where is the injection term, f the proton distribution function, and C Q a free parameter which is proportional to the total proton energy. As is shown in Appendix A, the solution of Equation 6 is (see also Atoyan et al. 1995;Celli et al. 2019) At time T age , the differential number of protons with energy E esc lying inside the conic MC shell can thus be calculated to be Finally, the hadronic γ-ray spectroscopic luminosity Φ γ (E γ ) is evaluated via the cross section presented in Kafexhiu et al. (2014) where dσ/dE γ is the γ-ray differential cross section. In addition to the hadronic component, we add a leptonic component contributed by the electrons accelerated by the blast wave. We approximate the spectrum of electrons in the SNR to be a broken power-law where dN e /dE the differential number of electrons, and E loss the maximum energy determined by energy loss; α 1 and α 2 are the power-law indices in the low energy and high energy band, respectively; E b , the break energy of the electron spectrum, can be constrained well (∼ 9 TeV, assuming a magnetic field of 6 µG, Ge et al. 2021) by the radio (Pineault & Joncas 2000) and X-ray (Ge et al. 2021;Fujita et al. 2021) emissions, which are both dominated by the SNR. The electron spectrum above the break is quite soft, and the leptonic γ-ray emissions are dominated over by the hadronic emissions above 10 TeV. Hence, the total γ-ray spectrum is insensitive to the electron spectrum above E b . We fit the SED as is plotted in Figure 1, with the parameters listed in Table 1. We find that the parameters are insensitive to X-ray flux, since the different sets of X-ray flux data from Ge et al. (2021) and Fujita et al. (2021) can be fitted with very similar parameters.
In Figure 2, we plot the proton spectrum insides the MC. As can been seen, the protons with energy < 75 TeV are still trapped by the turbulence driven by the escaping protons via Bell instability, absent in the MC. Hence, only protons with energies ranging from 75 TeV to 280 TeV can reach the MC within the relatively short lifetime of the SNR and contribute to a novel hadronic γ-ray spectrum. In Figure 3, we explore the dependence of hadronic γ-ray spectrum on α and δ and find that the choice of α (in the range 2.0-2.6) or δ (in the range of 0-1) does not impact the results significantly. The reason is that the protons which can hit the MC within the relatively short lifetime of the SNR at a narrow energy range (75-280 TeV) is almost mono-energetic (as shown in Figure 2). Hence, the energy dependence index δ and the power-law index α can hardly affect the proton spectrum in the MC. Celli et al. (2019) have also proposed a similar phenomenological escaping scenario which can lead to a hard hadronic γ-ray spectrum in middle-aged SNRs, while in our case, we model a younger SNR from the first principle.
Since the hadronic γ-ray luminosity Φ γ ∝ C Q n MC Ω (where n MC is the number density of the gas in the MC, a free parameter), we only list C Q n MC (Ω/4π) as a single parameter in Table 1. Once C Q n MC Ω is obtained from data fitting, the total proton energy E p,tot (including that in the GeV protons trapped near the shock) can be calculated via If we adopt n MC = 100 cm −3 , E p,tot will be 1.6 (4π/Ω) × 10 49 erg, which is quite reasonable. As shown in Figure 1, the leptonic component accounts for the radio (Pineault & Joncas 2000) and X-ray (Ge et al. 2021 for Model A; Fujita et al. 2021 for Model B) emission and dominates the γ-ray flux below 500 GeV. Meanwhile, the hadronic γ-rays (the yellow lines) have a very hard spectrum (dN γ /dE γ ∝ E −1 γ ) below 500 GeV, which stems from the proton cutoff at 75 TeV and the low-energy tail of the pp interaction cross section. The spectrum above 10 TeV is dominated by the hadronic component, and the whole γ-ray spectrum can be explained by the hybrid model naturally.

Age of SNR G106.3+2.7
There remains some uncertainty about the age of SNR G106.3+2.7: although the characteristic age of PSR J2229+6114 (≈ 10 4 yr) is usually adopted to be the age of the remnant, sometimes the remnant is suggested to be very young (∼ 1 kyr, see e.g., Albert et al. 2020). With an adiabatic expansion (R sh ≈ 1.2(E SN /ρ ISM ) 1/5 t 2/5 , where ρ ISM is the density of the ISM), the estimated E SN ≈ 7 × 10 49 erg (Kothes et al. 2001) is unusually low if the age of remnant is ≈ 10 kyr. However, as is estimated in Cardillo et al. (2015), Equation 1 indicates that the global maximum energy E max,global ∝ E SN , and SNRs with a low E SN can hardly accelerate CRs to ∼ 10 2 TeV. In order to reconcile the dilemmas, in this paper we consider a spherically symmetric scenario in which E SN = 10 51 erg (the canonical value) and T age = 1000 yr. Such a remnant age is plausible for a hosted pulsar with a characteristic age of ∼ 10 4 yr. Assuming a canonical braking index n:=νν/ν 2 = 3 (where ν is the spin frequency of the pulsar), the real age of the pulsar writes T age = τ c [1 − (P 0 /P now ) 2 ], where τ c is the characteristic age of the pulsar, P 0 the initial spin period of the pulsar, and P now the spin period of the pulsar at present. For PSR J2229+6114, P now = 50 ms (Halpern et al. 2001), and T age can be as small as ∼ 1 kyr if P 0 is appropriately close to P now . Such value of P 0 is allowed, since in pulsar evolution models P 0 is usually suggested to be in a wide range from ∼ 4 ms to ∼ 400 ms (see e.g., Arzoumanian et al. 2002;Faucher-Giguère & Kaspi 2006;Popov et al. 2010). Instances in which P 0 ≈ P now and T age ≪ τ c can also be found in the catalog: (a) CCO 1E 1207.4−5209 inside SNR G296.5+10.0 has P 0 ≈ P now = 424 ms, leading to τ 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 ≈ P now ≈ 10 2 ms, giving rise to τ c ≈ 2 × 10 8 yr, 4 orders of magnitude larger than T age ∼ 5 × 10 3 yr measured from the observation of the associated SNR (Gotthelf et al. 2005). On the other hand, as can be seen from Equation 1, E esc is determined by the R sh and v sh which can be constrained by observations directly, irrespective of T age . Ge et al. (2021) showed that v sh is at least 3000 km s −1 in the southwestern tail region based on the non-thermal X-ray spectrum up to 7 keV without a clear spectral cutoff, while R sh ≈ 6 pc can be estimated via radio observations (see e.g., Tibet ASγ Collaboration 2021). Although T age does affect R dif , its impacts can be compensated by the free parameter of D 100 TeV . Hence, for simplicity we only consider a spherically symmetric remnant evolving in uniform medium following the well-known Truelove & McKee (1999) model 1 .

Other issues
A sharp cutoff in the novel proton spectrum is an interesting characteristic in our particle escaping scenario with a spherically symmetric morphology applied, which can well fit the Tibet ASγ data in this SNR. Meanwhile, instead of a sharp cutoff, smoother break can be predicted if an asymmetric morphology is considered. In that case, the shock radius R sh and velocity v sh both vary in different directions, and thus E esc ∝ v 2 sh R sh also varies in different direction accordingly, and thus the spatial variation of E esc may turn the sharp cutoff into a softer break. Since the data available at present can not distinguish the asymmetric model, here we adopt a symmetric model to explain the γ-ray spectrum in the zeroth order. Further observation carried out by LHASSO may help to distinguish these models. During the lifetime of SNRs, 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 catalog 2 (Ferrand & Safi-Harb 2012), there are only a few SNRs with γ-ray emissions in 10 TeV which are confirmed to be of hadronic origin, and our novel hadronic spectrum may be expected to be found in more SNRs with the help of LHASSO (see e.g., Aharonian et al. 2021), Tibet ASγ, HAWC, CTA (see e.g., Abdalla et al. 2021) and ASTRI Mini-Array (see e.g., Pintore et al. 2020) experiments in future.

SUMMARY
Although the standard theory of non-linear diffusive shock acceleration can predict hard proton spectrum with α < 2 (Caprioli et al. 2020), it is noted that observations and recent numerical simulations seem to prefer a soft hard proton spectrum (Caprioli et al. 2020). In this paper, we explore an alternative plausible scenario to explain the bizarrely hard γ-ray spectrum of SNR G106.3+2.7 within the current acceleration theory which predicts soft (α ≥ 2) proton spectra. Apart from the common leptonic component which can be constrained by radio and X-ray observations, we invoke a novel hadronic component. This component arises from the escape scenario proposed by Bell (2004) and Cardillo et al. (2015), in which the CRs at the shock are diffused by the turbulence that is generated by the escaping CR particles with higher energies. Consequently, at a given time, only CRs with the highest energy can escape upstream, while CRs with lower energies are confined. Therefore, only protons with energies between 75-280 TeV can reach the MCs at an age 1 kyr. The cutoff of the proton spectrum at 75 TeV then gives rise to a very hard hadronic γ-ray spectrum below 10 TeV because of the low-energy tail of the pp cross section. The hadronic γ-ray spectrum, which dominates above 10 TeV, together with the leptonic component, can explain the bizarrely hard γ-ray spectrum of SNR G106.3+2.7 well.
In our model treatment, the γ-ray spectrum is calculated in zeroth order, assuming that a spherically symmetric SNR expands in a uniform medium. There are seven free parameters in the hadronic component, namely, n ISM , D 100TeV , R 1 , R 2 , M ej , η, and n MC C Q (Ω/4π). Neither of parameters α and δ has strong correlation to the hadronic γ-ray spectrum, because the protons which can hit the MC within the relatively short lifetime of the SNR in the narrow energy range (75-280 TeV) are approximately mono-energetic.
and the boundary condition is ∂f /∂r = 0 and F | r=0 = 0. Defining the differential operator can be further written as LF = Qr. We firstly discuss a simpler equation LG(r, ζ) = δ(r − ζ) with the boundary condition G| r=0 = 0. The solution is apparently Then we multiply Equation A2 by h(ζ) = Q(ζ)ζ on both sides, integrate over ζ, and exchange L and , i.e., Hence, ∞ 0 G(r, ζ)Q(ζ)ζ dζ = F , and we thus have 10 −6 10 −2 10 2 10 6 10 10 10 14   0.3 ... * d is the distance to the SNR, R now the radius of the SNR at present, E SN the explosion energy of the SNR, M ej the ejecta mass, B the magnetic field strength, and W e the total electron energy. T CMB , T FIR , and T NIR are the temperatures of the CMB, far infrared, and near infrared photons, respectively; u CMB , u FIR , and u NIR are energy densities of the CMB, far infrared, and near infrared photons, respectively.