Supernova remnants as cosmic ray accelerators. SNR IC 443

We examine the hypothesis that some supernova remnants (SNRs) may be responsible for some unidentified gamma-ray sources detected by EGRET instrument aboard the Compton Gamma Ray Observatory. If this is the case, gamma-rays are produced via pion production and decay from direct inelastic collisions of accelerated by SNR shock wave ultrarelativistic protons with target protons of the interstellar medium. We develop a 3-D hydrodynamical model of SNR IC 443 as a possible cosmic gamma-ray source 2EG J0618+2234. The derived parameters of IC 443: the explosion energy E_o=2.7*10^{50} erg, the initial hydrogen number density n(0)=0.21 cm^{-3}, the mean radius R=9.6 pc and the age t=4500 yr result in too low gamma-ray flux, mainly because of the low explosion energy. Therefore, we investigate in detail the hydrodynamics of IC 443 interaction with a nearby massive molecular cloud and show that the reverse shock wave considerably increases the cosmic ray density in the interaction region. Meantime, the Rayleigh-Taylor instability of contact discontinuity between the SNR and the cloud provides an effective mixing of the containing cosmic ray plasma and the cloud material. We show that the resulting gamma-ray flux is consistent with the observational data.


Introduction
Supernova remnants (SNRs) are believed to be the most promising accelerators of ultrarelativistic particles (electrons, protons and nuclei) in our Galaxy [ 1,2]. They are responsible for the majority of cosmic rays (CRs) with energies up to 10 15 eV (Lorentz factor of proton γ p ∼ 10 6 ). Classical observations of synchrotron radiation from SNRs, mainly in the radio band, give experimental confirmation of the presence of relativistic electrons in these cosmic objects [ 3]. Recently, the presence of ultrarelativistic electrons in SNRs has been confirmed by optical, X-ray and soft γ-ray observations [ 3,4].
There are expected two main CR accelerators in SNR. One of them is a young pulsar (if there is any), i.e., a spinning neutron star with a strong magnetic field. The second one is a strong shock wave as the result of an interaction of supernova ejecta with the interstellar medium (ISM). It is the main object of our investigation.
The theory of CR acceleration at the shock front via the first order Fermi mechanism is quite well elaborated [ 1,2,5]. Its main prediction is a considerably higher efficiency of the proton acceleration in comparison with an electron case, i.e., up to 10% of kinetic energy of the incoming flow may be transformed into the energy of newly accelerated CRs. Experimental evidences of the proton acceleration at shock fronts are limited to the solar system region (a bow shock at the boundary of the Earth's magnetosphere, interplanetary shocks and solar flares). Therefore, it is of great importance to find additional evidences of the proton and nucleus component acceleration at the fronts of cosmic shock waves.
Such a new possibility is offered by the recent observations of considerable γ-ray fluxes in directions on some radio bright SNRs [ 6,7,8]. Among the 32 unidentified Galactic plane γ-ray sources detected by the Energetic Gamma-Ray Experimental Telescope (EGRET) aboard the Compton Gamma-Ray Observatory (CGRO) five sources are coincident with SNRs without known pulsars. Four of these SNRs (W28, W44, γ Cygni and IC 443) exhibit evidences of a strong interaction with interstellar clouds. Therefore, it is believed that the observed γ-ray emission is the result of an interaction of the accelerated by SNR shock wave protons with the SNR matter itself [ 6] or with the matter of the adjacent molecular clouds [ 7].
In this paper we check this hypothesis for the case of well-studied SNR IC 443 and show that according to the atypically low initial energy of a supernova explosion such simple explanations do not hold here. In order to explain the γ-ray emission of IC 443 we elaborate a hydrodynamical model of IC 443 as a result of a supernova explosion in the medium with a large scale density gradient, what allows us to explain observational data concerning X-ray surface brightness distribution and obtain the parameters of remnants, especially shock wave characteristics (section 2). Using these parameters, we consider in detail the SNR shock interaction with a nearby molecular cloud (section 3) and show that the reverse shock enhances the energy and number density of newly accelerated CRs (section 4). Then we show that caused by the Rayleigh-Taylor instability, effective mixing of SNR and molecular cloud matter leads to the generation of an observable flux of γ-rays (section 5). Conclusions are given in section 6.

Hydrodynamical model of SNR IC 443
SNR IC 443 is well studied at different energy bands from radio to γ-rays [ 9,10,11,12,13,14] (figure 1 and 2). It is believed to be the result of a supernova explosion which happened about t = 3000 ÷ 5000 yr ago near a dense molecular cloud at the distance about d ≈ 1.5 kpc from us. Its unusual nonspherical morphology should reflect the inhomogeneity of the surrounding ISM, but until now a self-consistent model of IC 443 is absent. We suppose that the main reason of this is the absence of a detailed numerical modelling of IC 443. So, despite the evident influence of inhomogeneity of the surrounding ISM, IC 443 models are mainly based on the use of the one-dimensional (1D) self-similar Sedov solution [ 15] for SNR evolved in a uniform medium. The role of the density gradient in the surrounding ISM was taken into account only for the explanation of shape asymmetry within the framework of the approximate Kompaneetz [ 16] approach [ 3,13].
Therefore, we carry out here a detailed 2D calculation of the IC 443 evolution in a nonuniform medium using the proposed in [ 17,18,19] new approximate method of a complete hydrodynamical description of a point-like explosion in an arbitrary nonuniform medium. Besides geometrical sizes, we took the total flux and the surface brightness distribution in an X-ray band as free characteristics which give the best possibility for the estimation of plasma characteristics inside the SNR.  [ 14] and 95% confidence contour (solid line) of γ-ray source 2EG J0618+2234 [ 7].
Taking into account the generally accepted distance to IC 443 d = 1.5 kpc [ 3,13,14], we have reproduced the following picture of the supernova explosion. At t = 4500 yr ago, a SN explosion with the total energy E o = 2.7 · 10 50 erg took place in the transition zone between a warm ISM with the hydrogen number density n o and an interstellar HI cloud with the number density n i and the high scale h, so the initial density distribution in this transition zone was n(r) = n o +n i exp(−r/h), with n o = 0.16 cm −3 , n i = 16 cm −3 , h = 2.4 pc, the SN position moved onr = r o = 13.8 pc from the cloud centre and the initial hydrogen number density in the point of explosion n(0) = 0.21 cm −3 . Calculations show the following best fitting parameters of the IC 443 model. Depending on the direction from the explosion centre, the shock radius R, its velocity D, the postshock temperature T and the initial hydrogen number density at shock front position vary in the ranges: 7.5 R (pc) 10.1, 480 D (km/s) 880, 3.1 · 10 6 T (K) 1.1 · 10 7 , 0.16 n(R) (cm −3 ) 1.31. The total swept up mass of the ISM is M SN R = 28.0 masses of the Sun (M ⊙ ), the volume of the disturbed by a shock wave region is V = 1.1 · 10 59 cm 3 , the mean SNR radius is R = (3V /4π) 1/3 = 9.6 pc. In figure 3 the X-ray surface brightness of the IC 443 model is presented. The total equilibrium flux in an X-ray band is calculated for X-ray emissivity data from [ 20] and is obtained to be L x (> 2.4 keV)=1.4·10 34 erg/s (no absorption correction is made) which corresponds to the experimental result L x,obs (> 2 keV)=1.5·10 34 erg/s [ 13]. In figure 3(b) we show the effect of the absorption of interstellar gas with the column density N H, ISM = 3 · 10 21 cm −2 and an elongated cylindrical-like molecular cloud in front of the SNR with the radius R cl = 2 pc and the cloud column density N H, cl = (0÷5)·10 21 cm −2 . The data for the absorption cross-section from [ 21] were used. This absorption is mainly responsible for the unusual nonspherical X-ray shape of IC 443.

Hydrodynamics of the SNR IC 443 -molecular cloud interaction
Our calculations show that the visible X-ray shape of IC 443 is close to spherical, although with different radii in the south-west and the north-east, while the ROSAT X-ray image provides a strong evidence of two half-spheres in the IC 443 structure [ 14]. This difference may be naturally explained by the enhanced absorption of Xrays in the adjacent to IC 443 elongated molecular cloud which is visible in a radio band [ 10,11,12,22].
The same molecular cloud (or, more correctly, the part of the cloud disturbed by the shock wave of SNR) should be responsible for the observable γ-rays. Observations in CO and HCO + emission radiolines [ 22] reveal that the total mass disturbed by the shock gas in the molecular cloud reaches 500÷2000 M ⊙ . Although this molecular cloud covers practically the entire remnant, the interaction region is localized in the southeastern part and coincides with the γ-ray source position (figure 2). It means that the interaction region is somewhat displaced from the SNR symmetry axis (on the angle 0 ÷ 30 o , as viewed from the Supernova progenitor position), i.e., IC 443 has essentially a 3D shape.
From the hydrodynamical point of view, the SNR -molecular cloud interaction is a complicated event accompanied by the arising of two shocks: a forward shock in the cloud matter and a reverse shock in the SNR matter. The boundary between the ISM and the molecular cloud becomes a contact discontinuity between the shocked plasmas of the SNR and the cloud [ 3]. According to the complexity of the resulting flow we cannot calculate it within the framework of the used above method for one shock flow calculations. But simple analytical estimates are possible. So, the velocity of the forward shock wave in the molecular cloud D m is determined by the SNR shock velocity D i at the cloud position and the ratio of hydrogen number densities of the molecular cloud n m and the ISM n i [ 3]: D m = D i βn i /n m , where β depends on the number density contrast and β = 6 for n m /n i → ∞. If a SNR shock wave penetrates into the molecular cloud at the time t p ∼ 3500 yr after the SN explosion, when the shock radius is R p ∼ 7 pc, D i ∼ 600 km/s and, according to the Sedov solution [ 15], the effective thickness of the postshock gas shell is ∆R p ∼ 0.1R p , the typical time scale of the reverse shock wave existence is The forward shock penetrates into the molecular cloud up to the distance L p = D m · ∆t rev ∼ 0.1R p βn i /n m , and the shock-excited mass of the cloud is of order where ρ m = 1.4n m m p is the molecular cloud density (the typical number density of helium atoms n He = 0.1n is taken), m p is the proton mass, S m = πR 2 m is the effective surface of the SNR -cloud interaction. We note that the obtained mass of the shocked matter of the molecular cloud is close to the observed one, if we take typical values n m = 10 4 ÷ 10 5 cm −3 .
To summarise, the proposed here 3D model of IC 443 explains in a consistent way the main characteristics of the remnant, including its X-ray properties and interaction with the adjacent molecular cloud. We use this model for the estimation of the efficiency of cosmic ray acceleration and for the explanation of γ-ray production in IC 443.

Acceleration of cosmic rays in SNR IC 443
CRs, i.e., high energy relativistic particles (electrons, protons, nuclei) with the power-law energy spectrum N ∝ ε −α , α ∼ 2.6 ÷ 3.2 and the mean energy density ω cr ∼ 1 eV/cm 3 ∼ 10 −12 erg/cm 3 are believed to be accelerated by different mechanisms in our Galaxy (CRs with energies ε 10 17 eV) and in powerfull extragalactic sources (CRs with energies ε 10 17 eV, up to the maximal observable energies ε max = 3 · 10 20 eV), such as active galactic nuclei, quasars, radiogalaxies etc. By now the most promising mechanism for the galactic CR generation is the first order Fermi mechanism at shock wave fronts [ 1,2,5]. In this mechanism fast particles gain energy during their diffusive motion in the shock region as a result of scattering on different types of magnetic field fluctuations (magnetohydrodynamical waves, especially Alfvén waves, turbulent pulsation etc.) ahead and behind the shock front. The main reason for effective energy gain is the convergent character of hydrodynamical flows in the shock front frame.
The main characteristics of a shock wave acceleration are the following [ 1,2]. The acceleration time is where u i (i = 1 for upstream and 2 for downstream) is the plasma velocity in the shock frame, k i = l i c/3 is the diffusion coefficient, l i = ηr L,i is the mean free path, r L,i = ε/eH i is the Larmor radius of a particle with the energy ε and the charge e in the magnetic field H i , c is the light velocity, D = −u 1 is the shock velocity in the observer's frame. Parameter η is believed to be in the range 1 η 10 for SNR shocks.
The theoretically predicted spectrum of accelerated particles has the power-law type so α = 2 for plasma with the adiabatic index γ ad = 5/3. Both electrons and protons are involved in the acceleration process but the number of accelerated electrons is considerably less than the number of protons N e (ε)/N p (ε) = (m e /m p ) (α−1)/2 . Further we shall consider only a proton component of SNR CRs.
There are different estimates of the efficiency of energy transformation from the kinetic energy of a hydrodynamical flow into the energy of accelerated particles [ 1,2,5]. The reasonable result is that about 10% of the kinetic energy of a flow may be transformed at the shock front into the energy of CRs.The maximal energy of CRs is restricted by the processes of energy loss, particle escape from the region of acceleration, finite time of the shock existence etc. In the considered here SNR IC 443 case, taking into account the calculated IC 443 shock characteristics, we should expect the following characteristics of CR acceleration. The total energy of cosmic rays in IC 443 is (in order to provide the necessary amount of CRs in our Galaxy due to SNRs we need ν = 0.03). So, for our case E o = 2.7 · 10 50 erg we obtain W cr = 8.1 · 10 48 . For the spectral index of protons we take α = 2.1 which is close to theoretical predictions (1) and corresponds to the observable γ-ray spectrum (see below). Then the energy spectrum of IC 443 CRs (mainly protons) is The process of the pion creation in subrelativistic proton -proton at rest collisions is effective beginning from the proton kinetic energy ε = ε o ≈ 600 MeV for which the pion decay dominates over other mechanisms of γ-ray production [ 8]. The same proton energy is needed to generate γ-photons with the energy ε γ 100 MeV. The total number of such protons in IC 443 is The maximal energy of CRs is restricted in the IC 443 case by the SNR age t = 4500 yr. Then, from t = t acc we obtain The mean energy density of CRs inside IC 443 with volume V is Exactly speaking, CRs should not be uniformly distributed inside IC 443. Rather they will be concentrated in a thin shell of thickness ∆R ∼ 0.1R (∆V ≃ 0.3V) near the shock front, where the majority of SNR mass and about half of its thermal energy are concentrated [ 15], since the necessary time for diffusion on the distance of the order of a shock radius considerably exceeds the IC 443 age. It means that at least half of the total energy of CRs W cr is concentrated in this shell, so, the CR energy density near the shock front is ω sh ≈ 0.5W cr /0.3V ≈ 1.7ω cr . As we shall see below, this value is too low for the generation of the necessary gamma-flux from IC 443. Therefore, we propose and consider in detail the promising mechanism of the CR energy density enhancement in IC 443, connected with the reverse shock action.
In the reverse shock the plasma of SNR is subjected to another shock and its final values of the number density n f = γ ad /(γ ad − 1)n 2 and magnetic field H f = (n f /n 2 )H 2 are 2.5 times enhanced for γ ad = 5/3 plasma [ 23], which results in the increase of CR number density n cr,f /n cr,2 = n f /n 2 = 2.5. Meantime the energy of individual ultrarelativistic particles increases according to adiabatic invariant conservation (P 2 ⊥ /H = const where P ⊥ = m p v ⊥ is the normal to the magnetic field orientation component of particle's momentum) ε f /ε 2 = H f /H 2 . Therefore, we expect that behind the reverse shock front the energy density of CRs will be ω rev = ω sh n f n 2 3/2 = 6.6 ω cr (R p ) = 101 ν 0.01 For our model we will have ω rev = 818 eV/cm 3 . The maximal energy of CRs in the reverse shock region will be of the order of maximal energy in the main shock, i.e. about 10 12 eV, while acceleration conditions in the both shocks are similar.

γ-ray emission from IC 443
The γ-ray flux from source 2EG J0618+2234 which coincides with the SNR IC 443 position is well approximated by the power-law dependence [ 24] S γ = (5.4 ± 0.4) · 10 −10 ε γ 293 MeV The total γ-ray luminosity in ε γ 100 MeV band is L γ , obs = 2.2 · 10 35 erg/s and the rate of photon creation isṄ γ,tot = 7.6 · 10 38 phot/s. In [ 6] it was shown that in the case of a pulsar absence, the most promising mechanism for ε γ > 100 MeV γ-ray production is inelastic interaction of relativistic protons with protons at rest resulting in the creation of pions and their consequent decay into γ-rays. The γ-ray luminosity in this case is equal to the rate of energy transformation from relativistic protons to neutral pions [ 1] L γ = c n N ∞ ε min N(ε)σ pp (ε)ε π o (ε)dε (5) where n N = 1.4n is the mean number density of target nuclei in the region of interaction, ε min ≈ 600 MeV is the minimal proton kinetic energy of the effective pion creation (with the cross-section σ pp (ε) close to the mean value σ pp = 3 · 10 −26 cm 2 ), ε π o (ε) = ε/6 is the mean energy transformed into the pion. Substituting in (5) the cosmic ray spectrum (2) we obtain: So, as it was mentioned earlier in [ 6], it is necessary to have the initial number density in the IC 443 region n ∼ 10 cm −3 and the total energy of CRs of W cr ∼ 10 50 ergs in order to explain the observable γ-ray flux in such a simple model. But, as we showed in section 2 in the IC 443 case supernova exploded in a low density region with n ∼ 0.3 cm −3 and the total energy of CRs is expected to be only of the order of W cr ∼ 10 48 erg. In [ 7], the γ-flux from IC 443 is explained as the result of the interaction of accelerated by SNR CRs with the matter of the adjacent molecular cloud. The role of the molecular cloud consists in the increase of the number density of target nuclei n N in the interaction region. Namely, we can rewrite (6) in the form where M gas is the total mass of the interacting with CRs gas, m p is the proton mass. In [ 7], this total mass is taken to be M gas = 5 · 10 3 M ⊙ and the deduced mean CR energy density inside IC 443 is ω cr = 96 eV/cm 3 . But, as we can see from (3), this value of ω cr is considerably higher than the one expected for IC 443. Another problem here is a high total mass of clouds M gas ∼ 5 · 10 3 M ⊙ in the region of interaction, while the cloud mass outside SNR is included (as we mentioned above, only 500 ÷ 2000M ⊙ of gas is inside the SNR). The energy density of CRs outside the SNR should be considerably lower than that inside it [ 1,5], and, therefore, the pion production should be depressed.
We suppose that the observable γ-ray flux from IC 443 is caused by enhancement of the CR energy density in the reverse shock from the SNR -molecular cloud interaction. Indeed, using (4) we obtain from (7) the following γ-ray luminosity in the case of enhancement of the CR energy density: Therefore, we obtain the theoretical γ-ray luminosity L γ = 2.3 · 10 35 erg/s for IC 443 despite the low explosion energy E o = 2.7 · 10 50 erg and, accordingly, the low total energy of cosmic rays W cr = νE o = 8.1 · 10 48 erg.
Additional support for this model follows from the existence of an effective mixing mechanism of a hot containing CR plasma and dense shocked plasma of the molecular cloud -the well known Rayleigh-Taylor instability of a contact discontinuity between the both media which should lead to a decay of contact discontinuity and fragmentation of the shocked cloud material onto a set of cloudlets [ 12].
Another possibility of the effective CR -cloud material interaction takes place in the case of a strongly nonuniform cloud, in which dense molecular cores of mass ∼ 10 M ⊙ and number density ∼ 10 5 cm −3 are embedded in a relatively less dense intercloud medium with the number density ∼ 10 ÷ 10 2 cm −3 . As a result of displacement of a contact discontinuity in the molecular cloud direction, such separate cores penetrate into the shocked by a reverse shock plasma region and should be the sites of an intense γ-ray production. Radio observations of the shocked molecular gas in the IC 443 region reveal the existence of few such cloudlets [ 22].

Conclusions
If the association of some unidentified EGRET sources with SNRs is confirmed, it will be a new argument for the CR acceleration in SNRs. It is of great importance that contrary to a radio band, γ-rays give unique information about the acceleration of a proton component of CRs in SNRs. In this work we have proposed a self-consistent model of IC 443 which explains a total set of its main parameters, including shape, size, X-ray and γ-ray fluxes. We have shown that, despite the low supernova explosion energy, account for the SNR -molecular cloud interaction, especially for the role of the reverse shock wave, is the main reason of an intense γ-ray flux from IC 443. Further observational and theoretical efforts are needed for checking the theory of CRs acceleration by SNRs and IC 443 is a prominent candidate for such investigations.