LETTER TO THE EDITOR
Lepton models for TeV emission from SNR RX J1713.73946
Z. H. Fan^{1}  S. A. Liu^{2}  Q. Yuan^{3}  L. Fletcher^{2}
1  Department of Physics, Yunnan University, Kunming 650091, Yunnan, PR China
2  Department of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK
3  Key Laboratory of Particle Astrophysics, Institute of High
Energy Physics, Chinese Academy of Sciences, Beijing 100049,
PR China
Received 8 June 2010 / Accepted 4 July 2010
Abstract
Aims. SNR RX J1713.73946 is perhaps one of the best
observed shelltype supernova remnants with emissions dominated by
energetic particles accelerated near the shock front. The nature of the
TeV emission, however, is an issue still open to investigation.
Methods. We carried out a systematic study of four lepton models for the TeV emission with the Markov chain Monte Carlo method.
Results. It is shown that current data already give good
constraints on the model parameters. Two commonly used parametric
models do not appear to fit the observed radio, Xray, and ray
spectra. Models motivated by diffusive shock acceleration and by
stochastic acceleration by compressive waves in the shock downstream
give comparably good fits. The former has a sharper spectral cutoff in
the hard Xray band than the latter. Future observations with the HXMT and NuSTAR may distinguish these two models.
Key words: acceleration of particles  plasmas  shock waves  turbulence
1 Introduction
Energetic particles produce radiations over a broad energy range under typical astrophysical circumstances. To study the energetic particle population in highenergy astrophysical sources, observations over a broad spectral range are needed. Multiwavelength observations of individual sources usually give sparse spectral data points, which can be fitted reasonably well with a simple powerlaw or broken powerlaw function. A featureless powerlaw distribution implies the lack of characteristic scales or distinct processes in the system. The nature of particle acceleration processes are still a matter of debate nearly a century after the discovery of highenergy particles from the outer space (Butt 2009).
Recent progress in observations of highenergy astrophysical sources has resulted in more detailed spectral coverage, and spectral features start to emerge. These features are produced by the underlying physical processes with welldefined characteristic scales in space, time, and/or energy, and they play essential roles in advancing our understanding of the acceleration mechanism (Berezhko & Krymsky 1988). Diffusive astrophysical shocks are considered as one of the most important particle accelerators (Berezhko 1996; Drury 1983; Toptygin 1980). The acceleration of particles by compressive motions of astrophysical flows is also a very generic process (Ptuskin 1988). Both processes may be important in the acceleration of particles dominating the radiative characteristics of the shock associated with SNR RX J1713.7394 (Fan et al. 2010). Detailed Xray and ray observations of this source have revealed several features: 1) the Xray and ray emissions are well correlated in space (Plaga 2008; Aharonian et al. 2006); 2) the emission spectrum shows a clear highenergy cutoff in both the Xray and TeV bands (Aharonian et al. 2007), and the cutoff in hard Xrays appears to be sharp (Tanaka et al. 2008), implying a sharp cutoff in the electron distribution producing the observed radio to Xray spectrum through the synchrotron process; 3) the ray spectrum has a broad convex shape (Funk 2009); 4) there are bright Xray filaments with a width of 0.1 light year varying on a timescale of about a year (Uchiyama et al. 2007). Liu et al. (2008) showed that these emissions may be attributed to a single population of electrons with the rays produced via the inverse Compton scattering of the background photons by relativistic electrons.
In light of recent highresolution observations with Suzaku and Fermi, we use the Markov chain Monte Carlo (MCMC) method to constrain parameters of four possible lepton models. It is shown that simple models with a powerlaw electron distribution, which cuts off at TeV energies, give a poor fit to either the Xray or ray spectrum. Motivated by the mechanism of diffusive shock acceleration (DSA) and stochastic acceleration (SA) by compressive motions, we introduce two alternative models, both of which give acceptable fits to the observed spectra (Sect. 2). In Sect. 3, we discuss the implications of these results and future work needed to improve these models. Conclusions are drawn in Sect. 4.
2 Spectral fits with lepton models
In the lepton scenario, the radio to Xray emissions are produced through the synchrotron process of relativistic electrons, and the ray emission is produced through the inverse Compton scattering of the background radiation fields by the same electron population. If the background magnetic fields are relatively uniform, similar to the background radiation fields, the model can naturally explain the spatial correlation between the Xrays and rays (Plaga 2008). As pointed out by Liu et al. (2008), the magnetic field required to produce the Xray flux agrees with what is required to produce an Xray spectral cutoff by an electron population cutting off in the TeV energy range (Ellision et al. 2010; Aharonian et al. 2007).
There are at least four parameters in a specific model. In this work we employ an MCMC technique suitable for multiparameter determination to search all the model parameters. The MetropolisHastings sampling algorithm is adopted when determining the jump probability from one point to the next in the parameter space (MacKay 2003). The MCMC approach, which is based on the Bayesian statistics, is superior to the grid approach with a more efficient sampling of the parameter space of interest, especially for high dimensions. The algorithm ensures that the probability density functions (PDF) of model parameters can be asymptotically approached with the number density of samples. Starting with an initial parameter set , which should lie in a reasonable parameter space, one obtains the likelihood function of : . Then another parameter set is randomly selected with the corresponding likelihood function . This parameter set is accepted with the probability . If the new point is accepted, it becomes the starting point for the next step, otherwise, we reset . This procedure is repeated to derive a parameter chain that determines the PDF of model parameters. More details of this procedure can be found in Neil (1993), Gamerman (1997), and MacKay (2003).
I: Without detailed modeling of the particle acceleration
processes, a powerlaw electron distribution with an exponential
highenergy cutoff is often assumed to fit the observed spectrum:
(1) 
where p is the spectral index, the Lorentz factor of the highenergy cutoff , m_{e} and c are the electron mass and the speed of light, respectively. The integration of over the Lorentz factor gives the particle number density. It has been shown that this model gives a sufficient fit to previous observations of SNR RX J1713.73946 when the interstellar radiation field is considered (Porter et al. 2006); however, more detailed studies suggest that such a model systematically underestimates the photon flux at the lowenergy end of HESS observations (Tanaka et al. 2008; Zirakashvili & Aharonian 2010). In this paper, we adopt the interstellar radiation field spectrum of Porter et al. (2006). There are then four model parameters. Besides the magnetic field B, p, and , we need a normalization for the relativistic electrons , which is the total energy of electrons with . Figures 1 and 2 give the best fit of this model to the data and the probability density of the model parameters. All four model parameters are well constrained. However, detailed examination of residuals of the spectral fit shows that there are several data points below 1 and near 10 TeV with normalized residuals significantly greater than 3.
II: Liu et al. (2008) suggested that the fit to the ray
spectrum can be improved by considering a more gradual cutoff in
the electron distribution:
(2) 
Several mechanisms can cause variations in the shape of the highenergy cutoff (Blasi 2010; Zirakashvili & Aharonian 2007; Becker et al. 2006; Berezhko & Krymsky 1988). However, given the sharp cutoff of the Xray spectrum, Tanaka et al. (2008) claim that the electron distribution must have a very steep cutoff. Indeed, by using their data and performing the MCMC fit, we obtain radio fluxes nearly one order of magnitude below the observed values. The model also significantly overestimates the hard Xray fluxes. To increase the weight of the radio and ray data in the model fit, we artificially increase the errors of the Xray data from Suzaku by a factor of 2^{} (Figs. 1 and 2). The normalized residuals are indeed improved significantly with a reduced of 1.2, in agreement with Liu et al. (2008), where cruder Xray data from ASCA observations were used. However, as emphasized by Tanaka et al. (2008), the sharp cutoff in the Xray spectrum is significant and the error estimate of the Suzaku observation appears to be robust (Takahashi et al. 2008). This simple model therefore cannot fit the Xray spectrum.
III: The particle distribution is determined by the injection
process and the spatial diffusion coefficient
in the DSA model.
The timedependent electron distribution with a constant injection
rate may be approximated as (Drury 1991; Berezhko & Krymsky 1988; Forman & Drury 1983)
where is the supernova lifetime since the onset of particle acceleration and U the shock speed. A very sharp highenergy cutoff can be produced when gyroradii of highenergy particles exceed the characteristic length of the magnetic field in the background plasma l_{d}. Particles at such high energies do not perform gyromotions around the chaotic magnetic field and instead interact randomly with the background fields. The spatial diffusion coefficient of these particles increases as the square of the relativistic Lorentz factor in the absence of waves on scales greater than l_{d}. Consequently, they can readily cross the acceleration site without interactions. We assume a spatial diffusion coefficient
where e is the elementary charge, depends on the turbulence intensity, and 1/3 for particle interactions with a spectrum of plasma waves following the Kraichnan (1965) and Kolmogorov (1941) phenomenology for the turbulence spectral evolution, respectively. The characteristic length of the magnetic field is then given by the gyroradius of electrons at : .
Compared with the previous two models, an extra parameter is introduced here: , which is proportional to the electron acceleration timescale. Figures 1 and 2 show the results with . The probability density of peaks near a typical value of 12. Since has to be greater than 1, we have lyr. The gyroradius of electrons at the cutoff Lorentz factor is 0.01 pc comparable to the width of the observed highly variable Xray filaments. The Xray variability may then be attributed to spatial diffusion of Xray emitting electrons from a highdensity region formed presumably via an intermittent process. It is interesting to note that the best fit spectral index p=1.92 is less than 2, implying acceleration by multiple shocks or strong nonlinear effects (Zirakashvili & Aharonian 2010). The model with gives a similar fit with a slightly higher value of and sharper hard Xray cutoff.
IV: Fan et al. (2010) have carried out
detailed modeling of electron acceleration by compressive motions in
the shock downstream and argued that it might also explain these
observations.
Here we follow the same treatment and also consider the effect of
incompressive motions, which can dominate
the spatial diffusion and therefore the escape of energetic electrons
from
the acceleration region. The acceleration rate is given by (Bykov & Toptygin 1993)
where , k, v_{F}, and S are the spatial diffusion coefficient, the wavenumber, speed, and intensity of compressive waves, respectively. The overall wave intensity with for subsonic turbulence^{}. The turbulence is generated on a scale L near the shock front, and l_{d} corresponds to the scale where the turbulence speed is comparable to the Alfvén speed . Relativistic electrons with a gyroradius less than l_{d} have a meanfree path comparable to l_{d} since waves on smaller scales are subject to strong transit time damping by thermal background ions. We then have the diffusion coefficient
where x indicates the spatial location in the downstream, , and we have assumed Bohm diffusion for . Due to gradual decay of the turbulence intensity, related quantities evolve with x in the downstream. Equation (6) is quite different from Eq. (4), where particle scattering by the background turbulence is parameterized without considering details of the turbulence structure. Considering the presence of strong turbulence with between the scales of l_{d} and L, Bohm diffusion may be appropriate, and as we show below, faster diffusion may not be able to accelerate electrons efficiently to account for the observations. The intensity of the incompressive motions T(k) can be determined from the energy conservation: , where 3u^{2}/2 is the total turbulence energy density per unit mass. Further downstream, where v_{F} is high (6^{1/2} u), the incompressive motions vanish completely. To derive the escape time of electrons from the acceleration region, we need to consider the enhancement of the spatial diffusion coefficient by the turbulence motions. According to Eq. (4.30) in Bykov & Toptygin (1993), the effective diffusion coefficient can be derived from
=  (7)  
where we have ignored the decay of incompressive motions. Then we have
(8) 
With the above treatment, we can study particle acceleration by compressive motions all the way to the shock front, which is different from the case studied by Fan et al. (2010), where only the subsonic phase with v_{F}>u is considered. Near the shock front, we assume so that Eq. (5) is applicable. Since the compressive wave intensity is lower than in Fan et al. (2010) and the acceleration timescale of TeV electrons becomes comparable to the supernova lifetime, a timedependent treatment is necessary^{}. Both the acceleration and escape timescale explicitly depend on the diffusion coefficient , which is independent of for . Assuming that electrons are injected at the Lorentz factor with a constant rate, the timedependent results of Becker et al. (2006; see also Park & Petrosion 1995; Cowsik & Sarkar 1984) can be used to describe the evolution of the electron distribution function at energies below . Near the cutoff energy , the nontrivial dependence of the acceleration and escape timescales on the Lorentz factor demand numerical solutions. We assume the timedependent distribution function has a highenergy cutoff given by (For details, see Fan et al. 2010). Figures 1 and 2 show the results for U=4000 km s^{1} with the Kraichnan phenomenology. The quality of this fit is comparable to the DSA model. The deduced model parameters also agree with those obtained in Fan et al. (2010).
Figure 1: Model fit to the spectrum of SNR RX J1713.73946. The TeV data is from Aharonian et al. (2007), the Fermi data is from Funk (2009), and the Xray data is from Tanaka et al. (2008) (Thanks to Jun Fang). The thick solid, dashed line and the thin solid, dashed line in the upper panel correspond to the diffusive shock, stochastic acceleration model, and the model with an exponential, more gradual cutoff in the electron distribution, respectively. The low and highenergy spectral hump are produced by relativistic electrons through the synchrotron and inverse Comptonization process, respectively. The lower panels show the normalized residuals. From top to bottom, they are for the model with an exponential and more gradual cutoff, the DSA, and SA model. For the model with a gradual cutoff, the Xray errors have been artificially increased by a factor of 2 to improve the fit. The reduced of these fits are indicated in the topleft corner. See the text for details. The instrumental sensitivity of NuSTAR and HXMT are indicated by the dotted and dotdashed lines, respectively. 

Open with DEXTER 
Figure 2: Probability density function of model parameters normalized at the peak value. Different line types correspond to different models as explained in Fig. 1 and the text. 

Open with DEXTER 
3 Discussion
The SA model (Fig. 1) predicts a higher hard Xray flux than the DSA model since the diffusion coefficient of the DSA model has a stronger dependence on the electron energies leading to a sharper highenergy cutoff. The electron distribution in the SA model includes contributions from broad regions in the downstream, which also makes the overall electron distribution broader. Future observations with NuSTAR and HXMT may be able to distinguish these two models.
In actual turbulent astrophysical shocks, both the DSA and SA mechanisms may contribute to the electron acceleration (Bykov & Toptygin 1993; Achterberg 1990). The distinction between the two resides in the dissipation structure (Berezhko & Krymsky 1988). Supersonic shock fronts dominate the acceleration in the DSA model. The SA model dominates in the subsonic phase. It is possible that the shock downstream of SNR RX 1713.73946 has both a supersonic and subsonic phase turbulence with the former closer to the shock front. The electron acceleration is therefore a continuous process in the turbulent downstream. In the SA model of this paper, we intentionally remove the DSA by requiring that the compressive waves be subsonic in the downstream. A selfconsistent treatment of the turbulence evolution in the shock downstream is needed for particle acceleration near astrophysical shock fronts.
All the models appear to systematically underproduce ray flux in the Fermi energy band and overproduce flux near 10 TeV. The latter could be caused by inhomogeneities of the background magnetic field, which implies lower cutoff energy of the electron distribution and less highenergy ray flux (Tanaka et al. 2008). The former requires a broader electron distribution from the GeV to TeV energy range. A combination of DSA and SA model with the former dominating the higher energy particle acceleration and the latter enhancing the lower energy ones may address this issue. Such a scenario is possible if both supersonic and subsonic turbulence are produced by the supernova shock (Bykov & Toptygin 1993; Berezhko & Krymsky 1988). Contributions to GeVTeV emission from decay of neutral pions produced by inelastic collisions of relativistic protons with the background ions may also explain these residuals. However, as shown by Ellision et al. (2010) and Katz & Waxman (2008), relativistic protons can not be the dominant TeV emission component (Zirakashvili & Aharonian 2010). Including this component is beyond the scope of the current investigation but may give a good constraint on the relative acceleration efficiency of relativistic electrons and protons (Katz & Waxman 2008).
4 Conclusions
With recent high spectralresolution observations of SNR RX 1713.73946 in Xrays and rays, we show that, in general, models with simple parametric descriptions of energetic particle distribution give a poor fit to the broadband spectrum. Details of relativistic electron acceleration may be probed with advanced models. Although both the DSA and SA model give acceptable fits to the spectrum, the SA model predicts a higher hard Xray flux than the DSA model, which may be tested by future observations. The spatial distribution of energetic electrons in the SA and DSA models is also different, and the source structure can be used to distinguish these two models as well. The results agree with the scenario where energetic electrons are accelerated by interacting diffusively with a turbulent electromagnetic field and produce nearly all of the observed emissions from the shell of SNR RX 1713.73946. Turbulence evolution at astrophysical shocks appears to play a key role in advancing our understanding of particle acceleration processes. AcknowledgementsThis work is supported in part by the SOLAIRE research and training network at the University of Glasgow (MTRNCT2006035484), the National Science Foundation of China (grants 10963004 and 10778702), Yunnan Provincial Science Foundation of China (grant 2008CD061) and SRFDP of China (grant 20095301120006). S.L. thanks Eduard Kontar for helpful discussion. Q.Y. thanks Jie Liu for helping to develop the MCMC code adapted from the COSMOMC code of Lewis & Bridle (2002).
References
 Achterberg, A. 1990, A&A, 231, 251 [NASA ADS] [Google Scholar]
 Aharonian, F. A., Akhperjanian, A. G., BazerBachi, A. R., et al. 2006, A&A, 449, 223 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2007, A&A, 464, 235 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Becker, P. A., Le, T., & Demer, C. D. 2006, ApJ, 647, 539B [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Berezhko, E. G. 1996, Astropart. Phys., 5, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhko, E. G., & Krymsky, G. F. 1988, Sov. Phys. Usp., 31, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Blasi, P. 2010, MNRAS, 402, 2807 [NASA ADS] [CrossRef] [Google Scholar]
 Butt, Y. M. 2009, Nature, 460, 701 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bykov, A. M., & Toptygin, I. N. 1993, Phys. Usp., 36, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 CassamChenaï, G., Decourchelle, A., Ballet, J., et al. 2004, A&A, 427, 199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cowsik, R., & Sarkar, S. 1984, MNRAS, 207, 745 [NASA ADS] [Google Scholar]
 Drury, L. 1983, SSRv, 36, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Drury, L. 1991, MNRAS, 251, 340 [NASA ADS] [Google Scholar]
 Ellison, D. C., Patnaude, D. J., Slane, P., & Raymond, J. 2010, ApJ, 712, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Z. H., Liu, S. M., & Fryer, C. L. 2010, MNRAS, 406, 1337 [NASA ADS] [Google Scholar]
 Forman, M. A., & Drury, L. O. 1983, ICRC, 2, 267 [Google Scholar]
 Funk, S. 2009, Fermi Symp. [Google Scholar]
 Gamerman, D. 1997, Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference (London: Chapman and Hall) [Google Scholar]
 Katz, B., & Waxman, E. 2008, JCAP, 01, 018 [NASA ADS] [CrossRef] [Google Scholar]
 Kolmogorov A. 1941, Dokl. Akad. Nauk SSSR, 30, 301 [NASA ADS] [Google Scholar]
 Kraichnan 1965, Phys. Fluids, 8, 1385 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., & Bridle, S. 2002, Phys. Rev. D 66, 103511 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, S. M., Fan, Z. H., Fryer, C. L., Wang, J. M., & Li, H. 2008, ApJ, 683, L163 [NASA ADS] [CrossRef] [Google Scholar]
 MacKay, D. 2003 Information Theory, Inference, and Learning Algorithms, (Cambridge University Press) [Google Scholar]
 Neal, R. M. 1993 Probabilistic Inference Using Markov Chain Monte Carlo Methods, Technical Report CRGTR931, Department of Computer Science, University of Toronto [Google Scholar]
 Park, B. T., & Petrosian, V. 1995, ApJ, 446, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Plaga, R. 2008, NewA, 13, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Porter, T. A., Moskalenko, I. V., & Strong, A. W. 2006, ApJ, 648, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Ptuskin, V. S. 1988, Soviet Astron. Lett., 14, 255 [NASA ADS] [Google Scholar]
 Takahashi, T., Tanaka, T., Uchiyama, Y., et al. 2008, PASJ, 60, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, T., Uchiyama, Y., Aharonian, F. A., et al. 2008, ApJ, 685, 988 [NASA ADS] [CrossRef] [Google Scholar]
 Toptygin, I. N. 1980, SSRv, 26, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Uchiyama, Y., Aharonian, F. A., Tanaka, T., et al. 2007, Nature, 449, 576 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Zirakashvili, V. N., & Aharonian, F. A. 2007, A&A, 465, 695 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zirakashvili, V. N., & Aharonian, F. A. 2010, ApJ, 708, 965 [NASA ADS] [CrossRef] [Google Scholar]
Footnotes
 ... 2^{}
 Given the complexity of the instrumental background and calibration, and unknowns, such as the estimated cosmic Xray background, such an artificial increase may well be reasonable.
 ... turbulence^{}
 Here we adopt the formula given by Bykov & Toptygin (1993), which gives an acceleration rate a factor of 2 lower than that of Ptuskin (1988).
 ... necessary^{}
 At the point, where v_{F}=u, the compressive wave intensity used here is 3 times lower than in Fan et al. (2010).
All Figures
Figure 1: Model fit to the spectrum of SNR RX J1713.73946. The TeV data is from Aharonian et al. (2007), the Fermi data is from Funk (2009), and the Xray data is from Tanaka et al. (2008) (Thanks to Jun Fang). The thick solid, dashed line and the thin solid, dashed line in the upper panel correspond to the diffusive shock, stochastic acceleration model, and the model with an exponential, more gradual cutoff in the electron distribution, respectively. The low and highenergy spectral hump are produced by relativistic electrons through the synchrotron and inverse Comptonization process, respectively. The lower panels show the normalized residuals. From top to bottom, they are for the model with an exponential and more gradual cutoff, the DSA, and SA model. For the model with a gradual cutoff, the Xray errors have been artificially increased by a factor of 2 to improve the fit. The reduced of these fits are indicated in the topleft corner. See the text for details. The instrumental sensitivity of NuSTAR and HXMT are indicated by the dotted and dotdashed lines, respectively. 

Open with DEXTER  
In the text 
Figure 2: Probability density function of model parameters normalized at the peak value. Different line types correspond to different models as explained in Fig. 1 and the text. 

Open with DEXTER  
In the text 
Copyright ESO 2010