The Gamma-Ray Emission from the Supernova Remnant RX J1713.7-3946 Interacting with Two-phase Medium

We study the origin of gamma rays from the supernova remnant (SNR) RX J1713.7-3946. Using an analytical model, we calculate the distribution of cosmic rays (CRs) around the SNRs. Motivated by the results of previous studies, we assume that the SNR is interacting with two-phase interstellar medium (ISM), where dense clumps are surrounded by tenuous interclump medium. We also assume that only higher-energy protons (~>TeV) can penetrate the dense clumps. We find that pi^0-decay gamma rays produced by protons reproduce the observed gamma-ray spectrum peaked at ~TeV. On the other hand, it has recently been indicated that the observed ISM column density (N_p), the X-ray surface brightness (I_X), and the gamma-ray surface brightness (I_g) at grid points across the SNR form a plane in the three-dimensional (3D) space of (N_p, I_X, I_g). We find that the planar configuration is naturally reproduced if the ISM or the CR electron-to-proton ratio is not spherically uniform. We show that the shift of the observed data in the 3D space could be used to identify which of the quantities, the ISM density, the CR electron-to-proton ratio, or the magnetic field, varies in the azimuthal direction of the SNR.


INTRODUCTION
Supernova remnants (SNRs) are the most promising accelerators of cosmic rays (CRs) below the knee (∼ 3 × 10 15 eV). CRs are believed to be accelerated through a diffusive shock-acceleration mechanism (DSA; Blandford & Ostriker 1978;Bell 1978;Blandford & Eichler 1987). Nonthermal emissions from the accelerated CRs have been observed. The detections of synchrotron X-rays in some SNRs have been regarded as evidence for the acceleration of electrons to ultrarelativistic energies at SNR shocks (Koyama et al. 1995). On the other hand, it is still uncertain whether the gamma rays from SNRs are generated through either leptonic (inverse Compton, IC, scattering of low-energy photons by high-energy electrons) or hadronic (π 0 -decay photons generated through pp-interaction) processes. The sharp cutoff at low energies (∼ 100 MeV) in the gamma-ray spectra of the middle-aged SNRs IC 443 and W44 can be regarded as the so-called π 0 -bump, which is direct proof for the hadronic origin of the gamma rays (Giuliani et al. 2011;Ackermann et al. 2013). On the other hand, for the young SNR RX J1713.7-3946 (RX J1713 hereafter), Ellison et al. (2010) indicated that the hadronic model of gamma-ray emission cannot reproduce the observed X-ray emission because of an overproduction of thermal X-ray line emission, using 1D hydrodynamic simulations of a supernova blast wave. However, it has been claimed that the thermal X-ray line emission can be faint if the SNR is interacting with inhomogeneous interstellar medium (ISM; Zirakashvili & Aharonian 2010;Inoue et al. 2012;Gabici & Aharonian 2014). These studies also indicated that the spectrum of the π 0 -decay gamma rays generated by CR protons can mimic that of leptonic gamma rays if the ISM is inhomogeneous.
In general, angular resolutions of gamma-ray detectors are worse than those of radio and X-ray telescopes. Thus, gamma-ray observations with high angular resolutions could change the situation. Recently, Fukui et al. (2021) analyzed H.E.S.S gamma-ray data of RX J1713 with improved angular resolutions. They assumed that the gamma-ray counts are given as a linear combination of two terms: one is proportional to the ISM column density, and the other proportional to the X-ray counts. By fitting the expression to the data pixels, that is to say, the data at grid points across the SNR, they discovered that the gamma-ray counts are well represented by a plane in 3D space formed by the ISM column density, the X-ray counts, and the gamma-ray counts (see Figure 4 in Fukui et al. 2021 or Figure 7 in this paper). They indicated that the plane angle suggests that the hadronic and leptonic components constitute (67 ± 8)% and (33 ± 8)% of the total gamma rays, respectively.
In this paper, we study gamma ray emissions from an SNR interacting with inhomogeneous ISM taking RX J1713 as an example. We also discuss the implications of the plane discovered by Fukui et al. (2021). In particular, we take account of the fact that an SNR is rather spherically symmetric, which was not explicitly considered by Fukui et al. (2021). This paper is organized as follows. In Section 2, we describe our models for the CR distributions and the ISM. In Section 3, we show the results of our model and discuss the spectral energy distribution and the plane. The conclusion of this paper is presented in Section 4. We derive the distribution function of CR protons outside the SNR based on the model by Ohira et al. (2011). It can be obtained by solving a diffusion equation, where r is the position, p is the CR momentum, n p is the distribution function, D ISM (p) is the diffusion coefficient for CRs, and q s (t, r, p) is the source term of CRs. We assume that the SNR is spherically symmetric and r is the distance from the SNR center. It is also assumed that CRs with a momentum p escape from the SNR at t = t esc (p) (Ptuskin & Zirakashvili 2005;Ohira et al. 2010). For a point source, the source term is written as , and the solution is where and N esc (p) = dt d 3 r q s (t, r, p) , which represents the spectrum of the whole CRs.
In the case of a spherical SNR, CRs escape from a surface, R esc (p). Thus, the source term is given by The solution of equation (1) can be obtained using equation (2) as the Green function, We need to specify R esc (p), t esc (p), and N esc (p). In this study, we consider the case where the SNR is in the Sedov phase at least until recently. Thus, the shock radius is represented by and the escaping radius is given by where we adopt κ = 0.04 (Ptuskin & Zirakashvili 2005;Ohira et al. 2010). As long as κ 1, the difference between R esc and R sh can be ignored in the following discussion. We assume that CRs with p = p esc escape from the surface (r = R esc ) at t = t esc . We expect that the escape momentum p esc is a decreasing function of the shock radius. Thus, we adopt a phenomenological power-law relation, where p max and R Sedov are the escape momentum and the shock radius at the beginning of the Sedov phase (t = t Sedov ), respectively. We set R Sedov = 2.1 pc, t Sedov = 210 yr, and α = 6.5 following Ohira et al. (2011), although these values are not particularly adjusted to RX J1713. The maximum momentum p max is a parameter, and we fix this at p max c = 10 15 eV to reproduce the observed gamma-ray spectrum. Eliminating R sh from equations (7) and (9), and replacing p esc and t with p and t esc , respectively, we obtain We assume that the CR spectrum at the shock front is given by where A is the normalization. In the following, we adopt s = 2 and β = 3(3 − s)/2; the latter is the relation obtained for a thermal leakage model for CR injection (Ohira et al. , 2011. The spectrum of the escaped CRs (p > p esc ) can be written as ). The normalization is determined from the total energy of the CR protons with pc > 1 TeV (E tot,CR ), which is treated as a parameter.
Outside the SNR (r > R esc ∼ R sh ), the diffusion coefficient is assumed to be D ISM (p) = 10 28 χ pc 10 GeV δ cm 2 s −1 (13) (Ohira et al. 2011). We adopt Kolmogorov-type turbulence (δ = 1/3), which is theoretically motivated (see also Tang & Liu 2021) and close to those derived from observations (δ ∼ 0.4; Evoli et al. 2015;Genolini et al. 2015). The constant χ(≤ 1) is introduced because the coefficient around SNRs can be significantly reduced by waves created through the streaming of escaping CRs (e.g. Skilling 1975;Fujita et al. 2010Fujita et al. , 2011. In this study, we fix it at χ = 0.01 (Fujita et al. 2009). In this way, equation (6) is determined. We note that the equation is applied to escaped CRs (p > p esc ). It is assumed that the functional form of CR electron density n e (t, r, p) is also represented by equation (6) if radiative cooling is inefficient. The ratio to n p is assumed to be K ep = 0.05 in the fiducial model, which is consistent with the observed spectral energy distribution (SED; see Figure 3). In contrast with protons, however, electrons suffer from radiative cooling such as synchrotron emission, IC scattering, and nonthermal bremsstrahlung. Thus, if the cooling time of the electrons t cool,e is longer than t − t esc (p), the density is n e (t, r, p) = 0.

Inside the Supernova Remnant
For the sake of simplicity, we assume that the diffusion coefficient for CRs is almost zero at r < R sh because of efficient scattering of CR protons by turbulence developed down the shock. In this case, the distribution of CRs depends on advection and adiabatic loss and can be obtained by solving the following equation (Ptuskin & Zirakashvili 2005;Ohira et al. 2010): where u is the gas velocity, and we simply set where σ is the compression ratio of the shock and u sh is the shock velocity (Ptuskin & Zirakashvili 2005). In the following, we adopt σ = 4. Given equation (11), the solution of equation (14) is written as In equation (16), the normalization A can be determined from N esc (p) (equation (12)), because it is written as where first term on the right-hand side describes the particles that run away from the shock front, and the second term describes the particles escaping from the shock interior (Ptuskin & Zirakashvili 2005). With the help of equations (9) and (11), they are represented as and respectively . In equation (18), ξ cr is the ratio of the CR pressure to the shock ram pressure, and we adopt ξ cr = 0.5 (Ptuskin & Zirakashvili 2005). We note that Equation (19) is realized when β > (s − 1)(1/σ − 1), which is consistent with our chosen parameters. Thus, once the normalization of equation (12) is given, normalization A in equation (16) is derived from equations (17)-(19). We note that equation (16) is applied for CR protons that have not escaped from the SNR (p < p esc ). Thus, n p = 0 for p > p esc at r < R sh because they have escaped. Moreover, considering that the advection time of CRs, is finite, we also assume that n p = 0 at the radius r that satisfies t < t adv . The distribution of CR electrons is given by n e = K ep n p , where K ep is the electron-toproton ratio, if cooling is ignored.

Interclump Medium (n ic )
Clump (n cl ) Cavity Figure 1. Schematic of the SNR. The ISM surrounding the cavity consists of dense clumps with density n cl and interclump medium with density nic. The inner and outer boundaries are represented by r = R1 and R2, respectively. The SNR only recently hit the inner boundary, and the shock radius (R sh ) and the escape radius (Resc) are located inside the ISM

Interstellar Medium
In our model, the supernova explosion occurred inside a low-density cavity produced by stellar winds driven by the progenitor massive star. We also assume that the SNR only recently hit the dense ISM surrounding the cavity, which we simply refer to as "the ISM" from now on. Motivated by the results of numerical simulations (Inoue et al. 2012), we consider a two-phase medium, in which dense clumps are immersed in the tenuous interclump medium ( Figure 1). Thus, we assume that the density of the interclump medium (n ic ) is much lower than that of the clumps (n cl ), and the volume filling factor of the clumps is very small (x 1). This extreme density contrast is required to avoid strong thermal emission from the interclump medium (see Section 3.1), and it can be realized if stellar winds from the progenitor star have blown out the low-density ISM before the explosion of the supernova (see Figure 9 in Inoue et al. 2012). The average density of the ISM is represented by The inner and outer boundaries of the ISM are located at R 1 = 6 pc and R 2 = 9 pc, respectively ( Figure 1). We assume that the stellar winds had blow away clumps at r < R 1 ; this is required so that the profile of N p is consistent with observations (see Figure 4). For R sh < r < R 2 (upstream), the average density of the ISM is n av = 250 cm −3 and the density of the interclump medium is n ic = 0.025 cm −3 in the fiducial model. Assuming that the filling factor is x = 0.01, the density of the clumps is n cl ≈ 2.5 × 10 4 cm −3 (equation (21)). We choose these values to approximately reproduce observational results for the ISM (e.g. Sano et al. 2015;Fukui et al. 2021). The density contrast between the clumps and the interclump medium is fairly high (n cl /n ic ≈ 10 6 ). The magnetic fields that are compressed during the formation of the clumps are likely to support the clumps. It is assumed that n ic at R 1 < r < R sh (downstream) is σ (= 4) times larger than that of the upstream value, which means n ic = 0.1 cm −3 . 1 This is consistent with the value obtained by Katsuda et al. (2015), and it does not affect the results as long as n ic n av . On the other hand, we assume that n cl at the downstream side is the same as that at the upstream side because numerical simulations showed that the density of the clumps is not much affected by the passage of the shock (Inoue et al. 2012). Thus, the average density n av at the downstream side is almost same as that at the upstream side as long as x 1 and x ∼ const (equation (21)). The size of individual clumps does not appear in our formulation. Numerical simulations showed that when a shock front passes a clump, the clump maintains its density contrast and the front is not much distorted (Celli et al. 2019, see also Slavin et al. 2017;Wang et al. 2018). This suggests that the clump size does not significantly influence the clump sustainability and the shock propagation, although clumps that are too small may be destroyed at the shock passage.
The magnetic fields in the ISM (B) affect synchrotron emissions from CR electrons. We assume that the effective magnetic field strength in the ISM is B u = 12 µG at the upstream side (r > R sh ) and B d = σB u = 48 µG at the downstream side (r < R sh ). We do not mean that the magnetic fields in the interclump space are uniform; they can be intensified locally. In fact, numerical simulations showed that they are amplified around the clumps due to turbulence that develops through the passage of the shock (Inoue et al. 2012, see also Uchiyama et al. 2007). Synchrotron emissions from inside the clumps can be ignored because of the small filling factor.

Spectral Energy Distribution
In this section, radiations from the model SNR are calculated and they are compared with observations. We emphasize that we do not aim to reproduce the observations perfectly, and instead rather attempt to identify physics behind various observed relations. The shock and escape radii at present are set at R sh = 6.7 pc and R esc = 7.0 pc, respectively. The current age of the SNR is estimated to be t age = 1460 yr (equation (7)), and the energy of CRs currently escaping  is E esc ≈ p esc c = 34 TeV (equation (9)). The total energy of CR protons above 1 TeV is E tot,CR = 5.80 × 10 49 (250 cm −3 /1 cm −3 ) −1 erg (H. E. S. S. Collaboration et al. 2018).
Using the models in Section 2.1, we calculate the CR distributions and show them in Figure 2. CRs with energies of E < E esc are distributed at the downstream side of the shock. They have not been swept far downstream and are confined in the ISM (R 1 r < R sh ) because the SNR is relatively young. On the other hand, CRs with E > E esc have escaped from the SNR and are distributed at r > R esc . However, CR electrons with E E esc have cooled down mainly because of synchrotron radiation. Figure 3 shows the SED of the model SNR compared with observations of RX J1713. Here, we set a threshold energy at E th = 0.1 TeV, and assume that only CR protons with energies of E > E th can enter the clumps. This is because numerical simulations showed that magnetic fields and/or plasma waves develop around clumps, and they prevent lower-energy CRs from intruding into the clumps (Inoue et al. 2012;Inoue 2019). This means that at a given radius r, CR protons with E > E th interact with the ISM with an average density n av , while those with E < E th remain in the interclump medium with a density n ic . We ignore the volume of the clumps because the filling factor is x 1. On the other hand, CR electrons may behave differently because amplified magnetic fields around clumps may cool them down and prevent them from entering the clumps. In fact, observations showed that synchrotron emissions originated from electrons decline in large clumps in contrast with gamma-ray emissions (Sano et al. 2013;Sano et al. in preparation). Thus, we assume that while electrons with E > E th can enter the clumps, 90% of them have cooled down and cannot radiate. The following results do not much change as long as the fraction is large. Most electrons with E > E th remain in and interact with the interclump medium. Electrons with E < E th cannot enter the clumps and stay in the interclump space. Here, we implicitly assume that only a ignorable fraction of electrons have cooled down in the magnetic fields covering the clumps because of their small volume.
Although we chose E th = 0.1 TeV to reproduce the observed SED, it is consistent with the results of numerical simulations (Inoue 2019). This suggests that the Fe I Kα line emission at 6.40 keV is unlikely to be observed from the clumps inside the SNR, because while the line emission is produced through the interaction of MeV CR protons with the clump gas (Nobukawa et al. 2018;Makino et al. 2019), these low-energy CRs are not allowed to enter the clumps. 2 We calculate radiative processes for electrons using the models by Fang & Zhang (2008), and we derive gammaray spectra using the models by Kamae et al. (2006), Kelner et al. (2006), and Karlsson & Kamae (2008). The cosmic microwave background radiation and a farinfrared component with a temperature T = 26.5 K and a density of 0.415 eV cm −3 are the seed photon fields for the IC emission. The latter values are obtained from GALPROP by Shibata et al. (2011) at a distance of 1 kpc. The emissions from secondary electrons can be ignored. Figure 3 shows that the fiducial model reproduces the overall trend of the observed SED. While most of the gamma-ray flux is generated through the hadronic process (π 0 decay), a small fraction is attributed to the radiation through the IC scattering. The contribution of the nonthermal bremsstrahlung to the gamma-ray emission is negligible. The hadronic gamma-ray flux at ∼GeV is lower than that at ∼TeV, which is consistent with the results of the Fermi observations (Abdo et al. 2011). The suppression is attributed to our assumption that only CR protons with E > E th (= 0. Although thermal X-ray emissions from RX J1713 are very dim (Katsuda et al. 2015), they are not inconsistent with our model. This is because the density of the hot interclump medium is low and because the temperature of the clumps does not increase strongly because the shocks that propagate in them are weak (Inoue et al. 2012).
In Figure 3, the synchrotron flux we predict is lower than that of the Suzaku observations at E 10 keV. This may be caused by the oversimplification of our model. For example, all CRs with energies of > E esc = 34 TeV have escaped from the SNR, while those with < 34 TeV are confined within the SNR in our model. If we make a smoother transition, that is, if we allow some of CRs with > 34 TeV to be inside the SNR, the discrepancy disappears. NuSTAR observations may be useful to discuss this issue (e.g. Tsuji et al. 2019). In order to reproduce both radio and X-ray synchrotron flux observations, the energy spectrum of electrons rather than the magnetic field strength needs to be modified. However, the elections responsible for the radio emission do not contribute to X-rays and gamma rays that are our focus in this paper. Thus, we neither tune the spectrum nor add another component for these electrons. Figure 4(a) shows the profile of hydrogen column density N p derived from n av . To estimate the projected radius r proj , the distance to the SNR (RX J1713) is taken to be 1 kpc (Fukui et al. 2003;Cassam-Chenaï et al. 2004;Moriguchi et al. 2005). The profile of N p roughly reproduces observations (Fukui et al. 2021). In Figure 4(b), we present the profile of the gamma-ray surface brightness I g at 2 TeV and that of the X-ray surface brightness I X at 2 keV. While they are generally similar, the gamma-ray emission is noticed even at r proj 0.4 deg, which is not observed for the X-rays. The former is produced by CR protons escaped from the SNR (r > R sh and E > E esc ). Since the magnetic field strength at the upstream side is weaker than that at the downside (B u < B d ), CR electrons cannot effectively create synchrotron X-ray emission at r > R sh (Ohira & Yamazaki 2017). This is qualitatively consistent with the results of the H.E.S.S. observations that showed gamma-ray emissions extending beyond the Xray-emitting shell (H. E. S. S. Collaboration et al. 2018).

A Plane Formed in N p -I X -I g Space
Recently, Fukui et al. (2021) indicated that the pixel data of RX J1713 form a plane in the 3D space of the hydrogen column density (N p ), the X-ray surface bright-  Figure 6. Cross-section of the SNR. Unescaped CRs are permeated at R1 r < R sh , while the ISM is distributed at R1 < r < R2. For the line of sight L1,Ñp is the column density along the segments BC+C'B' (= ), while Np is that along AC+C'A'. For the line of sight L2,Ñp is the column density along the segment BB' (= ), while Np is that along AA'. ness (I X ), and the gamma-ray surface brightness (I g ). 3 They indicated that the plane is formed if both hadronic and leptonic IC components contribute to the observed gamma rays. In this section, we discuss the plane in detail using our model SNR.
In Figures 5 we plot (N p , I X , I g ) along the radial direction for the fiducial-model SNR; I X and I g are the values at 2 keV and 2 TeV, respectively. The projected radius (r proj ) increases in the directions of the arrow, and r proj of the individual points, which we call the "data", corresponds to those of the black dots in Figure 4(b). We use the data only for r < R sh because I X and I g are very small outside the shock (Figure 4(b)) and would be discarded in actual observations. However, this does not mean that our calculations for CRs at r > R sh are useless because we have explicitly confirmed that their contribution can be ignored, even to the projected emissions from the inside of the SNR (r proj < R sh ). Because Figure 5(a) shows the results for the fiducial model, the data are the same as those in Figure 4. Following Fukui et al. (2021), we fit the following plane to the data: where a and b are parameters and the best-fitting plane is shown in Figure 5(a). The best-fitting parameters are a = 0.0341 and b = 0.439. Figure 5(b) is a different view of Figure 5(a) and shows that the plane nicely fits to the data. The standard deviation in the direction of I g from the plane is 1.1 × 10 −16 erg cm −2 s −1 arcmin −2 . However, the plane normal is almost on the I x -I g plane in Figure 5, which contradicts Figure 4 in Fukui et al. (2021). This can be explained as follows. The hadronic gamma-ray brightness is written I gp ∝Ñ p n p , whereÑ p is the column density of the ISM hydrogens interacting with CRs. This is because while the ISM is distributed at R 1 < r < R 2 , unescaped CRs are permeated only at R 1 r < R sh (Figures 2 and 6). We emphasize that N p < N p and their relation is nonlinear (notÑ p ∝ N p ); the difference betweenÑ p and N p was not considered in Fukui et al. (2021). The leptonic gamma-ray brightness of bremsstrahlung origin and that of IC origin are represented by I ge,br ∝Ñ p n e and I ge,IC ∝ n ph n e , respectively, where n ph is the seed photon field and is the depth of the region containing CRs along the line of sight (Figure 6). On the other hand, because the brightness of X-ray synchrotron emission is given by we obtain I ge,br ∝Ñ p I X /(B 2 ) and I ge,IC ∝ n ph I X /B 2 . Thus, the gamma-ray brightness is represented as where f , g, and h are constants. The first, second, and third terms on the right-hand side represent the hadronic, IC, and bremsstrahlung components, respectively. For the fiducial model, Figure 3 indicates that the hadronic component is the main contributor to the gamma-ray emission. Thus, the gamma-ray brightness is approximated by I g ∼ I gp ∝ n pÑp . Moreover, because we have assumed that n p = K −1 ep n e (Section 2.1.2), the gamma-ray brightness should be where we used relation (23) and n av =Ñ p / . We emphasize that this relation represents the hadronic component, although the last one is represented by I X and apparently does not includeÑ p . Because n av = const and B = const in the region in which most CRs exist (R 1 r R sh ), we finally obtain the relation of I g ∝ I X , which explains the similarity between I g and I X in Figure 4(b). In Figure 5(a), the linear sequence L1 corresponds to the data around the line of sight L1 in Figure 6 and r proj ∼ 0. • 1-0. • 3 in Figure 4. Along this sequence, all of N p , I X and I g gradually increase as r proj increases. On the other hand, the linear sequence L2 in Figure 5(a) corresponds to the data around the line of sight L2 in Figure 6 and r proj ∼ 0. • 35-0. • 37 in Figure 4. Along this sequence, while N p changes little, I X and I g rapidly decreases as r proj increases. These two sequences form the plane in Figures 5(a) and (b).
Relation (25) suggests that the I X -I g relation depends on f Xg ≡ K −1 ep n av B −2 . This means that if f Xg varies, the plane shown in Figures 5(a) and (b) (I g ∝ f Xg I X ) should rotate around the N p -axis. This happens when the SNR and f Xg are not azimuthally uniform (i.e. not uniform around the SNR center), for example. First, we vary n av while keeping K ep and B unchanged. This modifies N p because it depends on n av . Specifically, we calculate the emissions from the SNR when n av = 150 cm −3 (low-density model) and n av = 350 cm −3 (high-density model); the other parameters are the same as those in the fiducial model (n av = 250 cm −3 ). In Figure 7, we show the data for the fiducial model and the high-and low-density models. As n av changes, the sequence shifts to the direction of N p and I g (red and green sequences in Figures 7(a) and (c)). We fit equation (22) to all the data, and the resultant plane is shown in Figures 7(a) and (b). The standard deviation of the data in the direction of I g from the plane is 2.4 × 10 −15 erg cm −2 s −1 arcmin −2 . The best-fitting parameters are a = 0.147 and b = 0.192. The plane normal has a substantial N p -component, which is consistent with Figure 4 in Fukui et al. (2021). Figure 7(d) clearly shows that the I g -I X relation rotates (I g ∝ f Xg I X ) as n av (i.e. f Xg ) changes. The shift and rotation of the sequence is the reason of the plane inclination, and they have nothing to do with CR electrons. Thus, the plane configuration discovered by Fukui et al. (2021) may simply reflect that the ISM around RX J1713 is not uniform in the azimuthal direction. In fact, Fukui et al. (2012) found that the ISM has both atomic and molecular hydrogens and that N p in the northwest region is larger than in the southeast region. Therefore, the observed plane could be reproduced even if the gamma-ray emission from RX J1713 is purely hadronic and if n av changes in the azimuthal direction (see also Hirai et al. 2020). We note that we do not perform quantitative comparison of the plane angle with that obtained by Fukui et al. (2021). This is because the plane fitting to the data is not necessarily good (Figure 7(b)), and thus the angle depends on the bias of the chosen data pointa. The SED of the whole SNR can be reproduced if the average of n av across the SNR is ∼ 250 cm −3 . Second, we vary K ep while keeping n av and B unchanged. Specifically, we calculate the emissions from the SNR when K ep = 0.01 (low-K ep model) and K ep = 0.1 (high-K ep model); the other parameters are the same as those in the fiducial model (K ep = 0.05). The results are shown in Figure 8. This time, the sequence shifts mainly in the direction of I X (Figures 8(a) and (d)) because the CR electron density n e depends on K ep . The sequence also shifts in the direction of I g (Figure 8(c)) because the contribution of the leptonic IC component changes in the gamma-ray band. We fit equation (22) to all the data, and the resultant plane is shown in Figures 8(a) and (b). The standard deviation of the data in the direction of I g from the plane is 2.5 × 10 −15 erg cm −2 s −1 arcmin −2 . The best-fitting parameters are a = 0.177 and b = 0.108. The plane normal has a substantial N p -component. Figure 8(d) shows that the I g -I X relation rotates as K ep changes. Again, the shift and rotation of the sequence is the reason for the plane inclination. If the plane angle for an SNR is determined by the variation of K ep as it is in Figure 8, it could be used to estimate the average respective contributions of the hadronic and leptonic components to the total gamma rays, as was done by Fukui et al. (2021).
Third, we vary B while keeping n av and K ep unchanged. Specifically, we calculate the emissions from the SNR when the upstream magnetic fields are B u = 6 µG (low-B model) and B u = 16 µG (high-B model); the other parameters are the same as those in the fiducial model (B u = 12 µG). The results are shown in Figure 9. This time, the sequence shifts only in the direction of I X (Figures 9(a) and (d)) because the synchrotron radiation depends on B, while N p and I g do not depend on B. We fit equation (22) to all the data, and the resultant plane is shown in Figures 9(a) and (b). The standard deviation of the data in the direction of I g from the plane is 3.0×10 −15 ergcm −2 s −1 arcmin −2 . The bestfitting parameters are a = 0.206 and b = 0.0423. The plane normal has a much smaller I X -component than in Figure 4 in Fukui et al. (2021). This may indicate that the magnetic field variation is not the cause of the plane inclination for RX J1713. Figure 9(d) shows that the I g -I X relation for the fiducial model rotates as B (i.e. f Xg ) changes. Figures 7-9 indicate that the rotation and the direction of the shift of the radial sequence (red, blue, and green) could be used to identify the parameter that azimuthally varies in the SNR, if individual radial sequences are similar. For example, the sequence shifts both in the N p and I g directions when n av changes (Figures 7(c)), while it does not shift in the N p direction when K ep or B change (Figures 8(c) and 9(c)). If individual radial sequences are not similar, it means that the SNR has a complex structure. Although it may not be easy, we could still explore the structure by solving an inversion problem. For example, we could constrain the structure from the observed pixel data of (N p , I X , I g ) by comparing them with a number of simulation data using deep learning.
Finally, we note that if the gamma-ray emission from the SNR is genuinely produced by IC scattering, it should be I g ≈ I ge,IC ∝ n ph I X /B 2 , which does not depend on N p explicitly. If B = const and I X is irrelevant to N p , we obtain a plane I g ∝ I X , which is inconsistent with the direction of the observed plane (Fukui et al. 2021). However, it could be interesting future work to investigate the possibility that B and/or n e have a special dependence on N p or n av so that the observed plane is reproduced.

CONCLUSIONS
We have studied nonthermal emissions from an SNR interacting with inhomogeneous ISM taking RX J1713 as an example. Based on an analytical model, we obtained the distribution of CRs around the SNR. We assumed that the ISM is composed of high-density clumps immersed in low-density interclump medium, as is suggested by numerical simulations. Moreover, we assumed that only high-energy (> 0.1 TeV) CRs can enter the clumps because magnetic fields are amplified around the clumps. While high-energy protons can radiate gamma rays through pp-interaction, most of the high-energy electrons have cooled down because of the amplified magnetic fields around the clumps, and they do not radiate gamma rays. Lower-energy CRs cannot penetrate the clumps and remain in the interclump medium. Our findings are summarized as follows: 1. Our model can broadly reproduce the SED of RX J1713. In particular, the observed gamma-ray peak at ∼TeV can be explained by the exclusion of the lower-energy protons from the clumps.
2. If the SNR is spherically symmetric and if the hadronic component is dominant in the gamma rays, the radial sequence of pixel data forms a plane in the 3D space formed by the ISM column density, the X-ray surface brightness, and the gamma-ray surface brightness. However, the plane angle is inconsistent with that observationally obtained by Fukui et al. (2021).
3. If the ISM density, the electron-to-proton ratio, or the magnetic field strength is not spherically uniform in the SNR, the plane angle significantly changes from the one at which the SNR is spherically uniform. In particular, if the ISM density or the electron-to-proton ratio is not uniform, the plane angle is qualitatively consistent with observed angles. By studying the rotation and the shift of the radial sequence of data in the 3D space, we could identify which parameter is not uniform.
In the era of the Cherenkov Telescope Array (CTA; Acharya et al. 2013), the intrusion of CRs into the clumps and the structure of the plane will be studied in more detail for even more SNRs because gamma-ray observations with higher angular resolutions become available (Acero et al. 2017).
We would like to thank the anonymous referee for a constructive report. We thank Y.