Prospects for Cherenkov Telescope Array Observations of the Young Supernova Remnant RX J1713.7−3946

1.1. Origin of Galactic Cosmic Rays

1.2. Origin of Gamma Rays from SNR RX J1713.7−3946

The identification of the emission mechanism of gamma rays from nonthermal SNRs like RX J1713.7−3946 is of utmost importance for advancing cosmic-ray physics. We will evaluate the capabilities of CTA on serving this purpose as follows. First, we perform morphological studies in order to assess the possibility of using spatial information to pin down the major component of the VHE gamma-ray emission from leptonic and hadronic origins. In the case where the leptonic component is dominant, we can proceed to check the capability of CTA to search for a possible "additional" hadronic component using spectral analyses. We will also evaluate the prospect of detecting a time variation of the spectral cutoff energy over a relatively long timescale, which we believe is a very useful and unique approach for identifying the VHE emission mechanism.

• 1.  
  Imaging with good spatial resolution.CTA is designed to reach angular resolutions of arcminute scale at TeV energies, which is at least a factor of three better than H.E.S.S. (Actis et al. 2011). This improvement is critical in our present study because the resolution of CTA will catch up with that of current CO and H i gas observations in the radio waveband and also get closer to that of X-ray images such as those obtained by XMM-Newton. Synchrotron X-rays are enhanced around the CO and H i clumps on a length scale of a few parsecs,¹⁴¹ but are relatively faint inside the clumps on a subparsec scale. In other words, the peaks in X-ray brightness are offset from those of the CO/H i gas distribution (Sano et al. 2013). Indeed, there exist many gas clumps whose associated X-ray peaks are located more than 0.05 deg away from their centers. If the gamma rays are leptonic in origin and if the energy spectra of the CR protons are spatially uniform, we expect that the gamma rays possess brightness peaks coincident with the CO/H i gas clumps and show similar spatial offsets to the X-rays.
• 2.  
  Broadband gamma-ray spectrum. CTA will cover photon energies from 20 GeV to 100 TeV with a sensitivity about 10 times better than H.E.S.S. (Actis et al. 2011). The spectrum of RX J1713.7−3946 measured by H.E.S.S. is well described by a single smooth component. However, it is possible that there exists another dimmer, and harder, component at the high-energy end of the spectrum that H.E.S.S. cannot discern owing to its insufficient effective area. More specifically, if the gamma-ray spectrum is described by an inverse Compton component with a primary electron index s = 2.0, which has been inferred from the Fermi/LAT observation (Abdo et al. 2011), there can exist an accompanying hadronic component from the CR protons with the same spectral index. This hadronic component should appear flat in νF_ν space, and its expected flux is

  $$\begin{eqnarray}\nu {F}_{\nu } & \approx & 3\times {10}^{-12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\\ & & \times \left(\displaystyle \frac{{E}_{\mathrm{CR}}}{3\times {10}^{49}\,\mathrm{erg}}\right)\left(\displaystyle \frac{n}{1\,{\mathrm{cm}}^{-3}}\right),\end{eqnarray} \tag{ 1 }$$

  where we have assumed the distance to RX J1713.7−3946 to be 1 kpc. The quantities E_CR and n are the total CR proton energy in the emission sites and the average target gas density, respectively. By adopting the values E_CR = 3 × 10⁴⁹ erg and n = 1 cm⁻³ as in Equation (1), this flux is about 10% of the currently observed νF_ν spectral peak at a photon energy of ∼TeV; such a hard component can dominate at the high-energy end of the gamma-ray spectrum and will be discernible as a spectral hardening feature by CTA.Note that a dim hadronic component can be dominant at the highest gamma-ray energies only if the maximum energy E_max of the parent CR proton spectrum is not lower than ∼100 TeV, and that such a high value of E_max is not guaranteed. In particular, if we assume the simplest one-zone model with the DSA around the shock of RX J1713.7−3946, as well as an inverse Compton origin for the gamma rays, then the maximum electron energy, ∼70 TeV, is not limited by energy loss but by the SNR age or CR escape to far upstream because of the weak average magnetic field strength of ∼10 μG required for a leptonic scenario, so that the maximum proton energy is basically similar to that of electrons. However, the emergence of a dim hadronic emission component with a maximum energy of more than 50 TeV is still possible if we consider a multizone model. Multizone models may be motivated by the observational results discussed in Section 1.2, which seem to be inconsistent with each other.¹⁴² For example, a small amount of protons, which were accelerated when the SNR was younger and the magnetic field was stronger, could have escaped the acceleration region and are hitting the surrounding molecular clouds now to produce dim hadronic gamma-ray emission, while most of the gamma rays are produced by either electrons or protons that are currently being accelerated. Assuming a Bohm-like diffusion in magnetic turbulence with an average strength of ∼0.1 mG, the propagation length of PeV protons within 10³ yr is estimated to be ∼1 pc, so that a non-negligible fraction of such high-energy protons can exist inside and around the SNR to interact with dense targets. Spatially resolved spectroscopy with the better angular resolution of CTA will help us investigate such "multizone" structures. Another possible scenario is the following. Both protons and electrons are currently being accelerated, but protons are accelerated in regions with stronger magnetic fields, while electrons are accelerated in regions with weaker magnetic fields. In the strong magnetic field regions, electrons are less efficiently accelerated, because the Alfvén Mach number is small owing to the strong upstream magnetic field, resulting in a small injection rate of electrons (Hoshino & Shimada 2002). Generally speaking, if a young SNR is producing CRs up to the knee energy, we naturally expect a dim hadronic component with gamma-ray energy more than 100 TeV unless the hadrons have escaped into the ISM during an earlier phase. Note that this statement is independent from the detailed acceleration mechanism. Therefore, if such a dim and hard hadronic component were detected, it would provide evidence for proton acceleration up to the knee energy.
• 3.  
  Time variation of maximum CR energies. The CTA observatory will last for more than 10 yr. This motivates us to measure any possible time variation of the gamma-ray spectrum over long timescales. It has been shown that the maximum energy of CRs confined in SNRs depends on the SNR age during the Sedov phase (e.g., Ptuskin & Zirakashvili 2003; Ohira et al. 2012). Since the downstream magnetic field strength around the shock decays with the shock velocity, as the SNR ages it becomes harder to confine high-energy CRs through pitch-angle scatterings near the shock. Accordingly, the maximum energy of CR protons should also decrease with time.¹⁴³ The overall normalization of the gamma-ray spectrum may also change with time, although its behavior is highly uncertain. Assuming that the energy in cosmic-ray protons is proportional to the integration of the kinetic energy flux of upstream matter passing through the shock front, Dwarkadas (2013) derived a simple scaling of hadronic gamma-ray luminosity, which can be written as ${L}_{\gamma }\propto {t}^{5m-2m\beta -2}$, where the SNR radius evolves as R ∝ t^m and the surrounding medium has a power-law profile in density ρ ∝ r^−β. Then, for a uniform surrounding medium (β = 0), we obtain L_γ ∝ t³ for m = 1 (free expansion) and L_γ ∝ t⁰ for m = 2/5 (adiabatic expansion). In the case of β = 2 (an isotropic wind), we find L_γ ∝ t⁻¹ for m = 1 (free expansion) and L_γ ∝ t^−4/3 for m = 2/3 (adiabatic expansion). Another model predicts a secular flux increase at a few hundred GeV at a level around 15% over 10 yr (Federici et al. 2015).In the following, we consider the simplest case in which only the cutoff energy of the gamma-ray spectrum changes with time. The cutoff energy may vary ∼10% in 10–20 yr. To demonstrate this, we consider a power-law behavior for the time evolution of the escape-limited maximum energy, i.e., E_max,p ∝ t^−α. The value of α is important for understanding the average source spectrum of Galactic CRs, which should be steeper than E⁻² (Ohira et al. 2010). Although there are theoretical studies on α (e.g., Ptuskin & Zirakashvili 2003; Yan et al. 2012), its precise value has not yet been determined observationally or theoretically. To reproduce the spectrum of Galactic CRs, we can phenomenologically assume that E_max,p ≈ 10^15.5 eV at the beginning of the Sedov phase (t_Sedov), and that E_max,p ≈ 1 GeV at the end of the Sedov phase (≈10^2.5 t_Sedov), from which we can obtain α ≈ 2.6 (Ohira et al. 2010). The age of RX J1713.7−3946 is around 10³ yr (Section 1.2). Then, the expected evolution rate of the maximum energy at this age is estimated as ${\dot{E}}_{\max ,{\rm{p}}}/{E}_{\max ,{\rm{p}}}\approx -\alpha /t\approx -2.6\times {10}^{-3}\,{\mathrm{yr}}^{-1}$. This implies a decay of the maximum energy of about 5% over a period of 20 yr. If the observed gamma rays are produced by CR protons, the cutoff energy of the gamma rays is proportional to E_max,p. On the other hand, the leptonic gamma-ray scenario predicts the cutoff energy to be proportional to ${E}_{\max ,{\rm{e}}}^{2}$. Therefore, the time variation of the cutoff energy of gamma rays over 20 yr is expected to be about 5% and 10% for CR protons and electrons, respectively. These estimations are based on the assumption of a smooth transition of E_max,p or E_max,e with a power-law scaling. However, RX J1713.7−3946 is thought to be currently interacting with dense regions. If so, a collision between the SNR shock and density bumps causes a sudden deceleration of the shock, which results in a variability of E_max,p and/or E_max,e on shorter timescales.For the CR electrons, if their maximum energy, E_max,e, is also limited by their escape from the SNR, its time evolution should be the same as that of the CR protons. On the other hand, if it is limited by synchrotron cooling, its evolution should differ from that of the CR protons. A cooling-limited E_max,e near the shock can either decrease or increase with time depending on the evolution of the amplified magnetic field (Ohira et al. 2012), for which shock–cloud interactions can play a significant role.

In this section, we evaluate the feasibility of achieving the scientific objectives discussed above using a series of observation simulations for CTA. The simulation software package ctools¹⁴⁴ (Knödolseder et al. 2013, 2016) that we use in this study is developed as an open source project aiming at a versatile analysis tool for the broader gamma-ray astronomy community, and it is very similar to the Science Tools¹⁴⁵ developed for the Fermi/LAT data analysis. ctobssim is an observation simulator that generates event files containing the reconstructed incident photon direction in sky coordinates, reconstructed energy, and arrival time for each VHE photon detected in the field of view. Simulated VHE sky images are produced by ctskymap using the simulated event files. Event selection based on these parameters can be applied by ctselect, and event binning is performed by ctbin. Unbinned or binned likelihood model fitting is executed through ctlike.

3.1. Input Data

3.2. Results and Analysis

3.2.1. Gamma-ray Image

We first look at the simulated images for CTA and study the possibility of using morphological information to determine the major component of the VHE gamma-ray emission. Different gamma-ray images in the energy range of 1–100 TeV are generated by changing A_p/A_e, the ratio between the hadronic and leptonic contributions. Figures 1(a) and (b) show the images for A_p/A_e = 0.01 (lepton dominated; e.g., Abdo et al. 2011) and 100 (hadron dominated; e.g., Fukui et al. 2012), respectively. Each image corresponds to 50 hr of observations with CTA. As per our assumptions described above for the underlying templates, the lepton-dominated model (Figure 1(a)) shows a gamma-ray image that resembles the X-rays, and the hadron-dominated case (Figure 1(b)) delineates the ISM proton distribution including both CO and H i. The difference between them (Figure 1(c)) is significant, as is evident from the subtraction between the two images.

Zoom In Zoom Out Reset image size

To perform a quantitative evaluation, we employ the method of likelihood fitting (e.g., Mattox et al. 1996; Cowan 1998; Feigelson & Babu 2012). We focus here on the ability of CTA to determine the dominant emission component from the prospectively observed morphology. We calculated the maximum log-likelihood, L_e and L_p, by fitting the simulated images with the leptonic and hadronic spatial template, respectively. We note that the templates used during the fits are the same as those used in the simulations, meaning that a rather idealistic case is assumed. To determine the log-likelihoods L_e and L_p individually, we employ either the hadronic or leptonic template one at a time in each fitting model (regardless of the intrinsic A_p/A_e ratio in the simulation data being fitted). We keep the normalization, photon index, and cutoff energy of the power law as free parameters during the fitting process. We also fit the data with the composite model containing both the leptonic and hadronic components in order to obtain the log-likelihood L_ep. Then we calculated the differences 2(L_ep − L_e) and 2(L_ep − L_p) for various A_p/A_e, as summarized in Table 1. These differences in log-likelihood are large when the composite model is a significant improvement over the leptonic or hadronic models, and the difference is intended to be distributed approximately as a χ² variable with 3 degrees of freedom (for the extra parameters in the composite model). However, this is only a rough indication, as the conditions of Wilks's theorem are violated because the simple model is on the physical boundary of allowed parameters of the composite model (Wilks 1938). Table 1 indicates that the composite model is strongly favored over the simple models and that we can easily find the dominant component when the contribution of the second component is small. For example, when A_p/A_e = 100, L_ep is very different from L_e but nearly the same as L_p. On the other hand, in the case where there is a more or less equal contribution from both particle populations (i.e., A_p/A_e = 1), we found that the fitting tends to be biased toward the hadronic component. When A_p/A_e = 1 in our model, the total number of hadronic photons is larger than that of the leptonic photons owing to the differences in E_c, especially at higher energies. Since the angular resolution becomes better at such higher energy ranges, the likelihood result expectedly favors a hadron-dominated scenario. Nevertheless, we note that the case of A_p/A_e > 1 (hadron dominated) shows a more spread-out structure in gamma rays than that of A_p/A_e < 1 (lepton dominated). This trend is consistent with the latest H.E.S.S. results (Abdalla et al. 2016). We conclude that CTA observations will be able to distinguish between hadronic and leptonic gamma rays based on morphological characterizations.

Table 1.  Differences in the Maximum Log-likelihood for Spatial Templates with Various Values of the A_p/A_e Ratio

Ap/Aea	2(Lep − Le)b	2(Lep − Lp)b
0.01	16.8	20085.8
0.1	761.9	15030.1
1.0	12280.2	3494.1
10.0	27083.1	171.9
100	31649.0	0.9

Notes.

^aA_e and A_p are normalization constants for the leptonic and hadronic component, respectively. The details are described in the text. ^bL_e, L_p, and L_ep are the maximum log-likelihoods for the leptonic, hadronic, and composite templates, respectively.

Download table as:  ASCIITypeset image

3.2.2. Spectrum

In the case where the leptonic component is dominant, we can subsequently attempt to search for the "hidden" hard component with a hadronic origin as discussed above. In order to evaluate the capability of CTA to achieve this interesting task, we further perform likelihood analyses over the wider energy band of E > 0.3 TeV using 50 hr of simulation data with various assumed A_p/A_e ratios. The spatial templates and the spectral shapes in the fitting models are the same as those used originally for the simulation. The fitting results are summarized in Table 2. We can significantly detect the hadronic component even for a small A_p/A_e = 0.02. However, the best-fit spectral shape for the hadronic component is generally slightly harder than the true input value of 2.0 for the cases with smaller ratios. This is probably because the dimmer and widely extended hadronic component is more easily confused with the background photons, especially at lower energies. Indeed, when we simulated and analyzed without the GDBG or the cosmic-ray background, the best-fit photon indices were converged consistent with the input value of 2.0. We thus refrain from drawing any strong conclusion on the detectability of the A_p/A_e ratio at this point.

Table 2.  Resulting Parameters of the Likelihood Analyses for the Simulation Data with Various Values of the A_p/A_e Ratio

Ap/Aea	2 (Lep − Le)b	Γec	Eced	Γpe
0.01	9.2	2.03 ± 0.01	17.08 ± 0.56	1.73 ± 0.26
0.02	41.9	2.03 ± 0.01	17.50 ± 0.60	1.86 ± 0.12
0.06	346.5	2.04 ± 0.01	17.62 ± 0.66	1.90 ± 0.04
0.10	998.0	2.03 ± 0.01	17.00 ± 0.67	1.93 ± 0.03

Notes.

^aA_e and A_p are normalization constants for the leptonic and hadronic component, respectively. The details are described in the text. ^bL and L₀ are the maximum log-likelihoods for the model with/without the hadronic component, respectively. ^cΓ_e is a photon index for the gamma-ray spectrum of the leptonic component. ^dE_c^e is a cutoff energy for the gamma-ray spectrum of the leptonic component. ^eΓ_p is a photon index for the gamma-ray spectrum of the hadronic component.

Download table as:  ASCIITypeset image

We proceed to perform maximum likelihood fits for the simulation data (with a ratio A_p/A_e = 0.1) for 18 logarithmically spaced energy bands spanning from 0.3 to 100 TeV. Unlike our analysis above, where we have to specify spectral shapes for the hadronic and leptonic components, in this way we can measure the spectrum independently of any spectral model assumption. Figure 2 shows the resulting spectrum from our "bin-by-bin" analysis of the same 50 hr simulation data. For each spatial template, it is clear that our likelihood fits satisfactorily reproduce the simulated spectrum above 0.3 TeV. The data in each energy bin were fit by power laws with a fixed index of 2.0. As a final check, we compared the total gamma-ray flux points from the likelihood fit over the 18 energy bins with our simulation model in which the emission is purely leptonic (see Figure 3). Although it depends on the energy binning, in the case of A_p/A_e = 0.1 and 50 hr of observation, we can obtain the flux points above 30 TeV deviating from the purely leptonic model with a statistical significance level of >3σ. Thus, we reach the conclusion that 50 hr of CTA observation could detect at above 3σ the hard hadronic component with a ratio A_p/A_e = 0.1, under the idealistic conditions assumed here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.3. Time Variation of Cutoff Energy

Detecting the time variation of the cutoff energy of the gamma-ray spectrum provides a clue to distinguish between different emission scenarios and CR acceleration theories. We start off with the generation of simulated photons for a 100 hr observation of RX J1713.7−3946. The simulations are performed for three sets of intrinsic cutoff energies at 17.9, 19.7, and 16.1 TeV each. The first set with E_c = 17.9 TeV is tagged as our nominal case based on the current best-fit value suggested by H.E.S.S. data. The other two cases represent the possibility of a varying E_c in the coming years within the CTA mission lifetime. As a start, we consider a variation of ΔE_c/E_c = ±10% to obtain an initial impression on how sensitive CTA will be to such fractional changes in the spectral cutoff. Roughly ∼10 yr could be expected for ±10% variation, as discussed in Section 2. Here we consider the purely leptonic scenario as an example (see Section 3.1 for details).

For each case, we explore how the detectability depends on the total exposure time by extracting subsets from the full data set with different exposures shorter than 100 hr. This is achieved by applying a time-based cut using ctselect. We then perform unbinned likelihood fitting for each of these subsets to derive the corresponding best-fit cutoff energies. In the fitting, the normalization, photon index, and ${E}_{{\rm{c}}}^{e}$ are all free parameters and are fitted simultaneously. Figure 4 shows the best-fit ${E}_{{\rm{c}}}^{e}$ and their associated errors from our likelihood analyses. As expected, the fitted ${E}_{{\rm{c}}}^{e}$ are closer to their true values, with smaller uncertainties for longer exposures in all cases.

Zoom In Zoom Out Reset image size

We then proceed to define a significance, s, for any observed ${E}_{{\rm{c}}}^{e}$ variation as

$$\begin{eqnarray}&&{s}_{\pm }(t)=\displaystyle \frac{| {E}_{\pm }(t)-{E}_{0}(t)| }{\sqrt{{\sigma }_{\pm }^{2}(t)+{\sigma }_{0}^{2}(t)}},\end{eqnarray} \tag{ 4 }$$

where E₀ and E_± are the best-fit E_c for the cases with the nominal value and ±10% variation, respectively. σ₀ and σ_± are the corresponding errors. To suppress statistical fluctuations, we repeat all of our simulations for 100 times and take the average of the calculated s(t) from all runs. In the top left panel of Figure 5, our result indicates that a decrease in E_c with time is slightly easier to identify than an increase, which can be understood as follows. In the energy range of E_c (∼10 TeV), the sensitivity of CTA begins to drop with energy. As a result, a lower cutoff energy can actually be measured more easily and precisely for a given exposure. Our result is hence consistent with expectations. We perform a simple fit to the points by a function $\propto \sqrt{t}$, where t is the exposure. If we observe >70 hr in the two epochs, we are able to achieve a 3σ detection for the ΔE_c/E_c = −10% case, whereas ∼100 hr are necessary for the +10% case.

Zoom In Zoom Out Reset image size

We can also estimate s(t) from another point of view by considering the following observation scenario. Supposing that we will observe RX J1713.7−3946 with different exposures during the first and second epochs, how does the detection significance of a variation in E_c depend on the observation time in each epoch? In this case, we can calculate

$$\begin{eqnarray}&&{s}_{\pm }({t}_{1},{t}_{2})=\displaystyle \frac{| {E}_{\pm }({t}_{2})-{E}_{0}({t}_{1})| }{\sqrt{{\sigma }_{\pm }^{2}({t}_{2})+{\sigma }_{0}^{2}({t}_{1})}},\end{eqnarray} \tag{ 5 }$$

where t₁ and t₂ are the exposure times for the first and second epoch, respectively. Our result is presented in the top right and bottom left panels of Figure 5. Here we again arrive at the conclusion that an increase of E_c is harder to detect. A longer exposure in the first year apparently makes it easier for us to detect the variation with a shorter observation in the next epoch. Although several combinations of observation time are available for the 3σ detection, 50 hr is the minimum required for the first observation period.

Since the full operation of CTA will begin roughly 10 yr after the H.E.S.S. observations, it is also of interest to know whether it is possible to detect any variation of E_c by comparing previous H.E.S.S. results and the upcoming CTA observations. In this context, we can calculate s(t) again by Equation (4), fixing E₀ = 17.9 TeV and σ₀ = 3.3 TeV (Aharonian et al. 2007a), and repeat the same procedure described above. The result is plotted in Figure 5 (bottom right panel). Owing to the large uncertainty of the H.E.S.S. results, however, the expected significance remains low and unmeaningful even with long CTA exposures up to 100 hr. This result consolidates the necessity of performing CTA observations of RX J1713.7−3946 in at least two epochs separated by a ∼10 yr time interval.

In this paper, we have studied the feasibility and prospects for achieving a set of key scientific goals through CTA observations of RX J1713.7−3946, based on simulations with exposures of 50 hr or above. We showed that a 50 hr observation is adequate to identify the dominant gamma-ray emission component, namely, leptonic or hadronic, from the morphology of the SNR that is going to be revealed by CTA. And we should be able not only to quantify both the leptonic and hadronic components but also to detect a possible hidden hard component at >3σ level through spectral analyses if they are mixed with a ratio of A_p/A_e ≥ 0.1. Interestingly, we also found that CTA will be able to reveal fractional variations of the spectral cutoff energy over a timescale of ∼10 yr, for the very first time. A variation of ΔE_c/E_c = ±10% is found to be detectable provided that an exposure time longer than 50 hr can be secured for the first epoch. The result may indicate that the detection would be a tough task, but at the same time the achievable science convinces us that it is worth the challenge.

The uncertainties of our results presented above are purely statistical. The energy dispersion and uncertainty in the measured energy scale can affect the fidelity of the spectral analyses especially in the measurement of E_c. We confirmed that the significance of the E_c variations due to the predicted energy dispersion would be negligible, while the expected energy scale uncertainty would cause ∼±1σ deviations, respectively. We also confirmed that the contamination of the very faint thermal X-ray to the leptonic template and variations of the hadronic template due to the uncertainty of the radio observations are both negligible to the results. Katsuda et al. (2015) reported the detection of very faint thermal components, but their spectral analysis region is limited only in the very center of this SNR, with a radius of 53. They also presented the softness map, and the analysis is performed on one of the softest regions where the thermal component is brighter than other areas. Even in such a region, the flux of the thermal emission is far less than 10% of the synchrotron. Concerning the radio map, the uncertainty of the intensity is estimated to be less than 10%. When we modified the hadronic image template by multiplying a random factor conservatively between ±10% to every image pixel, the results are very similar and differences are typically less than 1σ. Estimating other systematic uncertainties for our analysis is not trivial and can be highly model dependent. Possible effects due to a deviation from the actual PSF profile with a tail¹⁴⁸ may exist, since the PSF was assumed simply as a Gaussian in our simulations. Since RX J1713.7−3946 is located at the Galactic plane, the source is expected to be contaminated by the Galactic diffuse gamma-ray emission, which constitutes the majority of the background events (i.e., "noise" for our purpose). The possible existence of currently unknown TeV sources nearby, such as the CCO 1WGA J1713.4−3939, for instance, can contribute to the systematic errors in such a crowded region on the Galactic plane.

On the other hand, CTA will measure gamma rays with energies far below 0.3 TeV, where the Fermi/LAT and H.E.S.S. data points connect. If the gamma rays are mostly of a leptonic origin and the magnetic field strength B is around 10 or a few tens of μG, which has been inferred from theoretical studies (e.g., Yamazaki et al. 2009; Ellison et al. 2010; Lee et al. 2012), we expect that the electron spectrum possesses a cooling break (Longair 1994) at ≈7.5 TeV (B/40 μG)⁻²(t_age/10³ yr)⁻¹, where t_age is the age of RX J1713.7−3946. In this case the gamma-ray spectrum must also exhibit a corresponding break at E_b ≈ 0.3 TeV (B/40 μG)⁻⁴(t_age/10³ yr)⁻² (since ${E}_{\gamma }\propto {E}_{e}^{2}$ for inverse Compton emission in the Thomson regime), if the particle injection is impulsive. Future identification of such a cooling break feature in the gamma-ray spectrum by CTA will equip us with a powerful and independent tool to measure the value of B, which would bring important constraints on the particle acceleration mechanism at high Mach number collisionless shocks.

Our present study is based on a fairly simplified input model, but can be extended to include more physics motivated by theoretical models. Some possibilities are listed in the following, and their incorporation in our model is reserved for future works.

• (1)  
  Hadronic model based on MHD simulations: The shock interaction with dense clouds excites turbulence, which amplifies the magnetic field up to the mG order (Inoue et al. 2012; Sano et al. 2013, 2015). These results may be crucial for our better understanding of the gamma-ray and X-ray images of SNRs. Under these circumstances, the CR electrons will be significantly affected via synchrotron cooling for magnetic fields greater than 100 μG, and the gamma-ray image will follow if the leptonic gamma-ray production is important. In order to fully understand such effects on the gamma-ray images, we need to incorporate the shock–cloud interaction by utilizing the ISM distribution with a reasonably high angular resolution and detailed numerical simulations. Using the full performance of CTA by lowering the energy threshold could be essential in these studies, since more apparent differences between the leptonic and hadronic scenarios are expected.
• (2)  
  Highest-energy CRs from SNRs: Accelerated particles eventually escape from the SNR forward shock to become Galactic CRs. According to recent studies, the highest-energy CRs start to escape from the SNR at the beginning of the Sedov phase (Ptuskin & Zirakashvili 2005; Ohira et al. 2010; Ohira 2012). The escaped CRs make a gamma-ray halo, which is more extended than the SNR shell (Gabici et al. 2009; Fujita et al. 2010; Ohira et al. 2011). Spectral measurements of gamma-ray halos will provide information on the highest-energy CRs accelerated in the SNR in the past. Observed profiles of the gamma-ray halo will also constrain the diffusion coefficient of CRs and density structure in the ambient medium (Fujita et al. 2011; Malkov et al. 2013).
• (3)  
  Hydrodynamical models and spectral modifications via nonlinear effects: More sophisticated models taking into account NLDSA processes predict realistic CR proton electron spectra (e.g., Zirakashvili & Aharonian 2010; Lee et al. 2012), both trapped and escaped, with a full treatment of their time evolution and spatial distribution in a hydrodynamic calculation for RX J1713.7−3946. The comparison of these models with CTA results will provide important new insights on the very complex NLDSA mechanism at young SNRs.

This research made use of ctools, a community-developed analysis package for Imaging Air Cherenkov Telescope data. ctools is based on GammaLib, a community-developed toolbox for the high-level analysis of astronomical gamma-ray data.

We gratefully acknowledge financial support from the following agencies and organizations: Ministerio de Ciencia, Tecnología e Innovación Productiva (MinCyT), Comisión Nacional de Energía Atómica (CNEA), Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Argentina; State Committee of Science of Armenia, Armenia; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Brasil; Croatian Science Foundation, Croatia; Ministry of Education, Youth and Sports, MEYS LE13012, LG14019, 7AMB14AR005, and the Czech Science Foundation grant 14-17501S, Czech Republic; Ministry of Higher Education and Research, CNRS-INSU and CNRS-IN2P3, CEA-Irfu, ANR, Regional Council Ile de France, Labex ENIGMASS, OSUG2020 and OCEVU, France; Max Planck Society, BMBF, DESY, Helmholtz Association, Germany; Department of Atomic Energy, Department of Science and Technology, India; Istituto Nazionale di Astrofisica (INAF), Istituto Nazionale di Fisica Nucleare (INFN), MIUR, Italy; ICRR, University of Tokyo, JSPS, Japan; Netherlands Research School for Astronomy (NOVA), Netherlands Organization for Scientific Research (NWO), Netherlands; The Bergen Research Foundation, Norway; Ministry of Science and Higher Education, the National Centre for Research and Development and the National Science Centre, Poland; Slovenian Research Agency, Slovenia; MINECO support through the National R + D + I, Subprograma Severo Ochoa, CDTI funding plans and the CPAN and MultiDark Consolider-Ingenio 2010 program, Spain; Swedish Research Council, Royal Swedish Academy of Sciences, Sweden; Swiss National Science Foundation (SNSF), Ernest Boninchi Foundation, Switzerland; Durham University, Leverhulme Trust, Liverpool University, University of Leicester, University of Oxford, Royal Society, Science and Technologies Facilities Council, UK; U.S. National Science Foundation, U.S. Department of Energy, Argonne National Laboratory, Barnard College, University of California, University of Chicago, Columbia University, Georgia Institute of Technology, Institute for Nuclear and Particle Astrophysics (INPAC-MRPI program), Iowa State University, Washington University McDonnell Center for the Space Sciences, USA.

The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 262053.

We thank S. Katsuda for providing the original XMM-Newton image.