PeVatron Candidate SNR G106.3+2.7 in a Low-density Cavity: a Multiwavelength Test

In this paper, we constrain the density of the interstellar medium (ISM) around the hadronic PeVatron candidate, supernova remnant (SNR) G106.3+2.7, based on X-ray and $\gamma$-ray observations. The purpose of this investigation is to understand the influence of the gaseous environment on this SNR as a proton PeVatron candidate. By modelling the self-regulated propagation of the CRs injected from the SNR, we calculate the $\gamma$-ray emission of CRs via the hadronuclear interactions with the molecular cloud and the ISM, and use the measured $\gamma$-ray flux to constrain the ISM density around the SNR. Our results support the picture that the SNR is expanding into a low-density ($n<0.05 cm^{-3}$) cavity, enabling the SNR to be a potential proton PeVatron despite that it presently is not in the very early phase.


INTRODUCTION
Cosmic Rays (CRs) are high-energy charged particles (predominantly composed of protons and atomic nuclei).While CRs below 100 TeV are believed to be accelerated mainly at the blast wave of supernova remnants (SNRs) via diffusive shock acceleration (Bell 1978;Blasi 2013), the origin of PeV CRs remains under debate (Aharonian et al. 2019).It has been suggested that protons can be accelerated to PeV energies only at the very early phase ( 0.1 kyr) of an SNR in a dense environment (Schure & Bell 2013;Cardillo et al. 2015;Cristofari et al. 2020).SNR G106.3+2.7 has a cometary shape heading towards the northeast in the radio band, where the so-called Boomerang Nebula is in spatial coincidence with the head (Kothes et al. 2001).The tail region of the SNR seems to be expanding in a low-density wind cavity (Kothes et al. 2001).Observation on gas structure put the SNR at a distance of = 800 pc (Kothes et al. 2001).The size of the SNR is ∼0.3 • (∼4 pc at a distance of 800 pc) wide and ∼ 0.9 • (∼12 pc long.The Boomerang Nebula is powered by an energetic pulsar PSR J2229+6114 with a characteristic age of 10 kyr (Halpern et al. 2001).While the SNR and the pulsar may be generated in the same SN explosion, the characteristic age of the pulsar is not necessary the true age of the system.Based on SNR expansion model (Truelove & McKee 1999) and the current size of SNR, Bao & Chen (2021) and Yang et al. (2022) suggest that the age of SNR is only ∼ 1-2 kyr.On the other hand, an age of ∼ 4 kyr is suggested based on the spectral break in the radio spectrum of the PWN (Kothes et al. 2006).
Nevertheless, as an SNR aged at least a few kyr, SNR G106.3+2.7 ★ E-mail: ryliu@nju.edu.cn† E-mail: ygchen@nju.edu.cn is not supposed to be an efficient particle accelerator.However, intriguingly, this SNR has been measured in -ray band (Xin et al. 2019;Liu et al. 2020;MAGIC Collaboration et al. 2023;Albert et al. 2020;Tibet AS Collaboration et al. 2021;Cao et al. 2021;Fang et al. 2022) with the spectrum continuing up to several hundred TeV (Cao et al. 2021).The observation of MAGIC has revealed two distinct -ray emitting regions in the SNR, one in the head region while the other in the tail region (MAGIC Collaboration et al. 2023).In particular, the one in the tail region is in spatial association with a molecular cloud (MC), suggesting a possible hadronic origin of the gamma rays.While the GeV emission in the tail region may be explained as the inverse Compton (IC) scattering of electrons, the emission above TeV energy is difficult to be ascribed to the same process because of the spectral softening caused by the Klein-Nishina effect, especially when taking into account the simultaneous modelling of the X-ray spectrum (Ge et al. 2021).The MC is probably not directly interacting with the SNR, but may be illuminated by the protons escaped from the SNR via the proton-proton ( ) collisions (Liu et al. 2022), implying the SNR to be a proton PeVatron.This possibility is backed up by the X-ray observation.The X-ray emission of the entire SNR region is dominated by the nonthermal component (Ge et al. 2021;Fujita et al. 2021).The X-ray spectrum in the tail region continues to 7 keV without an obvious cutoff, implying a high speed of the SNR shock of ≥ 3000 km/s at the present time (Ge et al. 2021), and thus the CR acceleration is expected to be efficient.Indeed, from the radio morphology of the SNR, we may estimate the average speed of the shock to be ≃ 6000( /800 pc)( SNR /1 kyr) km/s, by simply dividing half of the length by the SNR's age.Note that this is a conservative estimation, because if we assume the SN explosion center is at the pulsar's location or considering the projection effect, the inferred shock speed would be even higher.
Apparently, the UHE -ray observation and the X-ray observation contradict with the theory for SNR acceleration, and raise an intriguing question that what makes SNR G106.3+2.7 so special.Ge et al. (2021) speculated that it is the low-density environment that makes the SNR shock in the tail region not significantly decelerated (although it is not as fast as that at the very early phase), and the protons can be gradually accelerated to PeV energy given a sufficiently long lifetime of the SNR.The low-density environment in the tail region is revealed by the observation of 21 cm HI emission, which may have been excavated by the stellar wind of a massive star or the blast wave of SN of previous generation (Kothes et al. 2001).If this scenario is true, the ISM density around the SNR tail is the key parameter.However, this scenario has not been seriously examined yet.
In this study, we aim to test this scenario with observations in other wavelengths and obtain a self-consistent picture for SNR G106.3+2.7 as a proton PeVatron.As will be shown in the following sections, we will constrain the density of ISM around the SNR in three channels.First, the downstream shocked ISM emits thermal X-rays, the emissivity of which is related to the ISM density.We will re-analyse the X-ray data in the SNR tail and extract the information of the density.Second, in addition to illuminating the MC via the collision when the escaped CRs reach there, these CRs will also radiate -rays when colliding with the ISM permeating the entire region.In this case the ISM density cannot be higher than a critical value, otherwise the hadronic emission will overshoot the -ray data especially at the GeV band.Finally, the ISM density will affect the Alfvén speed which determines the growth rate of the turbulent magnetic field (as will be detailed later).The latter controls the propagation of CRs in the ISM and hence influences the hadronic -ray emission of the MC.
The rest of the paper is structured as follows: in Section 2, we reanalyse the X-ray data in the tail region and derive the ISM density; in Section 3, we model the self-regulated CR propagation taking into account the streaming instability; we discuss and summarise our study in Section 4 and 5, respectively.Ge et al. (2021) show that the surface brightness profile of the X-ray emission in the head region declines monotonically with increasing distance from the pulsar in the head region, and the spectrum presents a gradual softening as well, implying that the X-ray emission is powered by the pulsar wind.Such a tendency continues until the boundary between the head and tail.Upon reaching this boundary, the surface brightness and spectral index become almost constant in the tail region, as revealed by the XMM-Newton and the Suzaku observations.This abrupt transition indicates a different origin of the X-ray emission in the tail region from that in the head region.The X-ray emission in the SNR head is dominated by the pulsar and PWN.We constrain the ISM density from the spectra extracted in the SNR tail, which is mainly covered by the XMM-Newton observations and listed in Table 1.We reduce the XMM-Newton data with the Extended Source Analysis Software (ESAS) that is integrated into the XMM-Newton Science Analysis System (SAS; version 17.0.0).

CONSTRAINTS FROM THE X-RAY OBSERVATION
The XMM-Newton data reanalysis essentially follows the procedure detailed in Ge et al. (2021).Because the XMM-Newton observations were originally proposed for the Solar Wind Charge Exchange (SWCX) X-ray emission, we ignore the X-rays below 1 keV that may be highly affected by the SWCX emission.After subtracting the point sources and background components, the spectra from the SNR tail are best fitted by a power-law model, which suggests that the X-rays are mainly of nonthermal origin.In order to constrain the ISM density, we add an additional non-equilibrium ionization (NEI) collisional plasma component to fit the underlying thermal emission from the ISM.We use solar abundance and zero redshift for the NEI thermal component, and let the temperature and ionization timescale be free parameters.The spectral fit is insensitive to changes in , which can thus actually not be constrained, and therefore we only use the volume emission measure of the thermal component derived from the NEl normalization to constrain the gas density.After setting free the best-fitting parameters of POWERLAW model, we refit the model.The fitting statistics are 2 /dof=1122/892 for POWERLAW and 2 /dof=1122/888 for POWERLAW+NET model.The 3 upper limit for the NEI normalization is also calculated.We then estimate the ISM density from the NEI normalization: We assume a uniform ISM density distribution in the tail region, then the upper limit of electron density is = 3.2 × 10 −2 cm −3 .We can also give a more conservative constraint on the ISM density if we assume all the X-ray emission in the tail is from the ISM.Though the fitting gets worse when using a single NEI model ( 2 /dof=1181/891 for NEI v.s. 2 /dof=1122/892 for POWERLAW), we find a conservative upper limit for the ISM density as = 7.4 × 10 −2 cm −3 .

CONSTRAINTS FROM -RAY OBSERVATIONS
Based on the -ray fluxes from the nearby associated MCs, the diffusion coefficients in the upstream of SNR shocks are usually believed to be suppressed by a large factor (see e.g., Gabici et al. 2007Gabici et al. , 2009;;Ohira et al. 2011;Li & Chen 2012).Since the energy density gradient near SNRs is large, the propagation of CRs in this region is strongly affected by their self-generated non-linear turbulence driven via the streaming instability.The streaming instability has two branches: the resonant Alfvén mode (see e.g., Ptuskin et al. 2008;D'Angelo et al. 2016;Nava et al. 2016;D'Angelo et al. 2018;Nava et al. 2019) and the non-resonant Bell mode (Bell 2004;Amato & Blasi 2009).Here we discuss the resonant branch, which is the most frequently discussed self-regulated mechanism to suppress diffusion coefficient in the neighbourhood of CR sources.
The growth rate of the resonant branch is proportional to the gradient of CR differential number density, and therefore it has a dramatic effect on low-energy (< 10 TeV) CRs whose differential number density is high.On the other hand, the resonant turbulence driven by high-energy CRs is expected to be weaker since the CR differential number density is small unless the age of the source is smaller than the escape timescale of high-energy CRs.

Self-regulated Transport
Following D'Angelo et al. ( 2016), we consider an SNR lying inside the Galactic disk.The coherence length c (over which the background magnetic field tube roughly keeps the same direction) of the background magnetic field in the Galactic disk is ∼ 100 pc (Beck et al. 2016).For simplicity, we consider one-dimensional diffusion.The 1D diffusion-advection equation reads where ( , , ) represents the CR distribution function on CR momentum , the distance from the release site of CRs, ( , , ) the diffusion coefficient, the advection velocity, and ( , , ) the injection term.
The magnetic field is assumed to be weakly perturbed at beginning, and thus is + A at > 0 and − A at < 0, where A = 0 / 4 i i is the Alfvén speed, d /d = 2 A ( ), i the number density of ions, and i the mean ion mass ( i = 1.4 p , with p the proton mass).The diffusion coefficient ( , , ) is defined as where L denotes the Larmor radius of the particles.The dependence of the spectral power of the turbulence F ( , , ) on the resonant wave number = 1/ L ( ) is defined by 2 ( , ) = 2 0 ∫ F ( , , ) d ln .The equation of turbulence evolution writes The growth rate of the self-generated turbulence is (Skilling 1971) For the damping rate Γ D , we consider the damping term which is equivalent to cascade Γ cas = (2 ) −3/2 A F ( , , ) (where ≈ 3.6, Ptuskin & Zirakashvili 2003).As can be seen, both Γ CR and Γ D is proportional to A and thus is proportional to −1/2 i , which means the ISM density can affect the growth and damping (mainly growth here, i.e., Γ CR > Γ D , because we are discussing the turbulence near the source) of the Alfvénic turbulence.For the boundary conditions, we set F | = c /2 = F Gal and | = c /2 = 0, where F Gal = L /[3 Gal ( )], and Gal ( ) = 3.6 × 10 28 1/3 GeV cm 2 s −1 (Ptuskin et al. 2009).
Bao & Chen (2021) considered the situation in which only the CRs with the highest energy can escape due to the turbulence generated via non-resonant branch of streaming instability.Here we consider another situation in which the accelerated CRs are released into the ISM all at once.Then, the injection is described as where cut represents the cutoff energy (∼ 100 TeV, Bell et al. 2013), the spectral index (=4 for canonical diffusive shock acceleration), and 1 the injection constant which can be determined by the total energy of the CRs CR,tot , i.e., with inj the injection scale of the CRs (comparable to the radius of the SNR).Prior to commencing complicated numerical calculations, it is worthwhile to get some hints from the analytical methods.As is shown in Ptuskin et al. (2008), in fully ionized plasma, if we drop the advection and adiabatic terms in Equation 2, the analytical solution of the diffusion coefficient reads (Ptuskin et al. 2008) where = ( g ) 1/3 4/3 0 2 5/3 3 1/3 7/3 K 8/3 (8) (with g being the Lamour radius of the protons) and = ∫ ∞ −∞ d , and the inverse of reads (Ptuskin et al. 2008) Obviously, when 2 / is large enough (≫ (2 4 /3 6 )Γ 8 (1/4)Γ −4 (1/2) 6 5 / 2 ; in other words, the CR energy is low enough), becomes energy-independent.For example, for canonical values CR,tot = 10 50 erg and = 4, Then, we obtain from Equations 9 and 8 That is, independent of the power-law index of injected CR spectrum , at the location in which 2 / is large enough, the energy spectrum of CR will be as hard as d /d ∝ −1.5 .In Figure 1, we compare the analytical solution and the numerical solution of and .As can be seen, the numerical solution roughly goes as ( ) ∝ −3.5 as is expected, and the diffusion coefficient is almost energy-independent before the escape dominates.In Figure 2, we also plot the spatial and energy dependence of particle distribution function for reference.As can be seen from the left panel, most of the particles are confined at ∼ 2 pc from the injection position.In the right panel, it can be seen that 1 PeV particles have partially escaped from the system, while the relatively low-energy particles are still confined well.

Application to G106.3+2.7
We use the code "AAfragpy" (Koldobskiy et al. 2021) to calculate the pionic −ray flux.The absorption of -rays during propagation through the -pair production on CMB photons is considered.The MC is modeled as a cylinder starting at = 1 and ending at = 2 saturating the magnetic flux tube, with a total mass of MC .We also consider the ISM within the magnetic flux tube as another target for the collision.In addition to the ISM density, there are other parameters that can affect the predicted -ray flux.To study the effect of the ISM density, we first model the -ray emission in a reference scenario, and then examine the influence of some parameters on the allowable range of the ISM density.The reference scenario uses intermediate-level parameters to establish a baseline, enabling the subsequent investigation of the influence of individual parameters.It is worth noting that the MAGIC data and the LHAASO/Tibet AS data are somewhat inconsistent with each other.As is discussed by MAGIC Collaboration et al. ( 2023), the MAGIC data may miss some -ray emission in the tail region.On the other hand, HAWC, Tibet AS and LHAASO data may contain the emissions from both the head region and the tail region.As suggested by some previous studies (e.g.Liang et al. 2022), the Boomerang Nebula might partially contribute to -ray flux above several hundreds of TeV.Therefore, we will fit the GeV flux, which is concentrated at the tail region of the SNR, observed by Fermi-LAT, and the flux measured by LHAASO up to 100 TeV.However, we consider the flux above 100 TeV measured by LHAASO as an upper limit for the SNR's emission.
In Figure 3, we show the reference scenario and the parameters used to fit the data are listed in Table 2.We explore the dependence of the -ray spectrum upon the ISM density as is shown in Figure 4.The ISM density can affect the GeV -ray spectrum, because the bulk of low-energy GeV protons have not arrived at the MC located at 20-30 pc away within the age of the SNR, and most of GeV emissions are produced by interactions with the intervening ISM between the source and the MC.Besides, the ISM density can also affect the resulting high-energy cutoff energy, since the growth rate (Equation 5, and the effective growth rate Γ CR − Γ D ) is proportional to the Alfvén speed.The increase in the ISM density will reduce the Alfvén speed and thus the growth of the turbulence will decrease correspondingly.As a consequence, it would aggravate the escape of the high-energy particles from the system before the balance between the growth and the damping of the turbulence can be built.Note that we do not aim to fit all the GeV-PeV -ray data with the model, but rather use them as the upper limits to constrain the ISM density.This is because the GeV flux as measured by Fermi-LAT may be partially contributed by the IC radiation of electrons.If so, the hadronicray flux need be adjusted to a lower level in order not to overshoot the measured flux.On the other hand, the -ray flux above 100 TeV may partially originate from the Boomerang PWN (Liang et al. 2022;Liu et al. 2020).As shown in Figure 4, we can see that an ISM density in the range of i = 0.01-0.05cm −3 are consistent with the data in the reference scenario.A higher ISM density will lead to an overshoot in the GeV -ray flux, while a lower ISM density would make a too high theoretical flux compared to the data around 100 TeV and above.However, the excess in the high-energy end may be avoided if the maximum proton acceleration energy is assumed to be smaller.So the constraint mainly poses an upper limit of 0.05 cm −3 for the ISM density, supporting the picture that the SNR breaks into a low-density cavity, as revealed by the X-ray observation.
The constraints on the ISM density from the -ray data are also relevant with some other model parameters, which are not necessarily equal to the values in the reference scenario.In Figure 5, we show the -ray spectra obtained by setting different values of the coherence length of the interstellar magnetic field (the top-left panel), the age of the SNR (the top-right panel), the position of the MC (the middle-left panel), the spectral index of injected protons (the middle-right panel), and the radius of the tube (the bottom-left panel), other than those in the reference scenario.Under each of the variation scenario, we then increase the ISM density until the resulting -ray flux reaches the 95% confidence level upper bound of any spectral data point of Fermi-LAT.As such, we obtain an upper limit of the ISM density in each scenario.As can be seen from Figure 5 and will be discussed below, our results do not vary significantly with these parameters, as long as the choice of parameters is reasonable.Note that, effects of changing different parameters may also cancel each other.So uncertainties of other parameters basically do not affect our conclusion that the tail of the SNR is expanding in a low-density cavity.
The coherence length determines the escape timescale of the particles, and therefore affects the cutoff energy of the -ray spectra.Meanwhile, the SNR age will also affect the cutoff energy of the spectra as the high-energy particles will escape gradually over time.The influence of the increased/decreased SNR age can be compensated by the increased/decreased coherence length.Also, the excess in the high-energy spectral end can be always cancelled by intro-ducing a smaller cutoff energy at the injection spectrum.These two parameters almost do not have influence on the low-energy -ray flux and hence the constraint on the ISM density from the GeV data is unaffected.
The position of the MC will also affect the -ray spectra.The larger distance from the SNR/release site will lead to a harder -ray spectra, because a smaller fraction of low-energy particles can reach the MC within the SNR age.As shown in Figure 2, the spatial distribution of particles are highly energy-dependent.As the distance increases, the proton spectra becomes harder and the low-energy gamma-ray flux decreases.An upper limit of = 0.1cm −3 for the ISM density can be obtained in the case of 1 = 36 pc.As can be seen, the SNR-MC distance will not have much effect on the density.
The injection power-law index determines the fraction of energy occupied by high-energy particles.While the index is not supposed to deviate from 2 too much according to the canonical shock acceleration theory, we here vary it slightly to see the influence.The resonant branch of the streaming instability is determined by the number density and number density gradient of high-energy particles.Therefore, a softer injection spectrum will lead to fewer high-energy particles and, consequently, weaker generated turbulence resonating with high-energy particles.The escape timescale of high-energy particles would become shorter remarkably for a softer injection spectrum.On the low-energy side, although the total number of particles (i.e., ) would increase significantly with a softer injection spectrum, selfregulation of particles can be established more quickly than in the case of a harder injection spectrum.The second term in the right hand side of Equation 10 dominates and most of the low-energy particles are confined within a small distance (∼ 1 pc) from the injection position, due to self-similarity of the system.Therefore, the increase in low-energy particles is not substantial at large distance.In the case of = 3.8 and 4.2, the upper limits of ISM density are found to be 0.02cm −3 and 0.1cm −3 , respectively.
The radius of the magnetic tube has a systematic influence on the energy density and the energy density gradient of high-energy particles.A larger tube radius would increase the volume of the tube and hence dilute the particle density, and vice versa.As a result, for the smaller tube size the cutoff energy is higher, because a smaller tube size will lead to a higher particle density and density gradient (in order to have the same total particle number).As is shown in Figure 1, the numerical solution will deviate from the analytical solution as the energy increases.This deviation occurs because the particle number density and density gradient are not large enough to reach selfregulation at high energies, and these particles diffuse more quickly than that predicted by the analytical solution (which assumes the self-regulation has been built).Hence, a larger tube size will relax the constraint on the ISM density.We find that with a tube radius of 4 pc, the upper limit of the ISM density is increased to 0.3 cm −3 .Although this upper limit is 6 times higher than the one obtained in the reference scenario, it is still consistent with the picture of a low-density cavity around the tail region of the SNR.In fact, the ISM density derived from the X-ray observation may be also used as a constraint on tube radius in this case.
We also explore the difference between an impulsive injection at = 0 and a constant injection of CRs.If the SNR has not been significantly decelerated, it is probably still processing as a proton PeVatron, and hence we may expect the particle injection is continuous.For simplicity, we assume a constant particle injection rate.In this case, more high-energy particles remain in the tube compared to the impulsive injection case, as more particles are injected recently and do not have time to escape the tube.Therefore, the resulting -ray spectra at the MC show a larger cutoff energy.On the other hand, the density of low-energy particles at the MC is not significantly affected as self-regulation is established at early time.Thus, the obtained constraint on the ISM density remains unchanged even if the particle injection rate is constant.Finally, we discuss the total amount of the injected proton energy.In the above calculation, we fix the total CR energy to 4 × 10 49 ergs which is relatively small compared with the typical proton energy budget of an SN explosion, i.e., a few times 10 50 erg.In general, it is believed that the SNR blast wave can convert ∼ 10% of its kinetic energy to the accelerated nonthermal particles, so the total injected proton energy may have an uncertainty of a factor of a few.From the perspective of -ray production, the amount of total proton energy is in degeneracy with the target gas density.However, the influence is not linear when considering the self-regulated transport of injected protons.A smaller injection normalisation will lead to less particle density gradient in weaker turbulence at high energies.It would then lead to an insufficient amount of high-energy protons at the MC to explain the multi-TeV data, and this deficiency cannot be compensated by employing a higher MC mass which is estimated in the range of 300 -800 ⊙ (Liu et al. 2022).On the contrary, if we adopt a larger total proton energy, it would only lead to a more stringent constraint on the ISM density (i.e., a lower upper limit of ).

Farmer-Goldreich damping and ion-neutral damping
Here we discuss Farmer-Goldreich damping (Farmer & Goldreich 2004;Lazarian 2016) and ion-neutral damping.Subsequent analysis will demonstrate that the Farmer-Goldreich damping is negligible.
Similarly, we will show that ion-neutral damping is also negligible.
The damping arise from interactions with pre-existing turbulence can be written as where A is the Alfvén Mach number, and MHD the injection scale.
The damping rate Γ IND is derived by solving the equation outlined in (Zweibel & Shull 1982).
where is the complex frequency of the wave, = 0 / 4 i i , ion-neutral momentum transfer ratio and = i / n is the ratio between the densities of ions and neutrals.With = R + Γ IND /2, Equation 13 can be written as X-ray observations have uncovered that the cavity density is remarkably low.Correspondingly, the environment's temperature is expected to be exceptionally high, leading to near-complete ionization.Due to this low density, the gas's cooling timescale exceeds 1 Myr (exceeding 10 6 K, Spitzer 1998).Additionally, the pressure balance also requires a high temperature, as the nearby denser gas should has density of ∼ 1 cm −3 and temperature of 10 4 K.The ionization rate can be estimated via Saha equation.The result is that, under 10 6 K, the hydrogen almost fully ionized, and about 90% of the helium are also ionized, i.e., the neutral fraction is ∼ 1%.
In Figure 6, we show the difference in whether considering Γ FG (with an unity A and MHD = c ) and Γ IND .As can be seen, both damping can be neglected.In considering IND, n = 0.0003, i = 0.0297, = 10 6 K.As can be seen, the difference is negligible.

DISCUSSION
The low ISM density indicates that the swept-up mass by the SNR blast wave is ∼ 0.1 ⊙ .With such a low swept mass, the SNR can hardly enter the Sedov-Taylor phase at its current age and the shock is not severely decelerated.In this case, the shock velocity can simply be obtained via dividing the SNR size by the SNR age.If we consider the current position of the pulsar as the explosion position, the shock velocity is about ∼ 6000( /800 pc)( SNR /1 kyr) km/s, which is consistent with the X-ray observation and the fact that PeV particles have been accelerated.In this paper, we constrain the ISM density in the tail region of the SNR considering an asymmetric expanding shock.In Bao & Chen (2021) and Yang et al. (2022), larger densities (0.2 cm −3 and 0.5 cm −3 ) are used to explain the size of the SNR with an isotropic expanding model.The larger density arises from the fact that there is a denser HI region in the northeast "head" region as is suggested by Kothes et al. (2001), and the density is averaged over 4 solid angle.While the tail region is still in the free expansion phase, the head region should have already entered the Sedov phase, and can no longer accelerate PeV protons.Although in the first hundred years the shock in head region could accelerate protons to PeV, these protons may have escaped from the shock already and are diffusing in the ambient ISM.It may suggest that the 100 TeV -rays (Cao et al. 2021) seen in the head region stems from the pulsar.Young SNRs such as SN 1006, RXJ 1713, and Vela Jr., characterised by fast shock velocities, typically exhibit bright synchrotron X-ray filaments amplified by strong magnetic fields (see e.g., Ressler et al. 2014).However, filaments are absent in SNR G106.3+2.7.The explanation is that, for perpendicular shocks, electron injection is hard to trigger (Bohdan 2023), and a quasiperpendicular SNR shock may form for SNRs expanding in stellar wind cavities (Voelk & Biermann 1988;Zirakashvili & Ptuskin 2018).So even if the magnetic field is enhanced, it may also be hard to see the filaments for perpendicular shock because there are not enough accelerated electrons.

CONCLUSION
In this paper, we constrain the ISM density around the PeVatron candidate, SNR G106.3+2.7, with the X-ray and -ray observations.We develop a model for the self-regulated propagation of CR protons injected from the SNR and calculate the -ray emission of the CRs via hadronuclear interactions with the MC and the ISM.We found that non-detection of the thermal X-ray emission from the tail region may impose an upper limit of 0.02 − 0.07 cm −3 for the ISM density.The measurement of GeV-PeV -ray flux may also constrain the density of the ISM where the SNR is expanding.
Our calculations demonstrate that, apart from the influence on the X-ray emission, the ISM density also plays a crucial role in shaping the GeV -ray flux and the high-energy cutoff of the resulting -ray spectrum.In order to achieve an acceptable fit to the GeV-PeV -ray data and a reconciliation with the X-ray observations, the ISM density should lie ∼ 0.01 − 0.05 cm −3 for the reference scenario.Taking the observed -ray data as upper limits of the hadronic emission from the CR interactions with the ambient medium, we find that most of the other parameters have minor impact on the constraint of the ISM density up to a factor of a few.Our results is in agreement with the scenario that the SNR is expanding into a low-density cavity, which may enable it to be a proton PeVatron.