EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 573 (January 2015) A&A, 573 (2015) A134 Full HTML
Free Access
Issue
A&A
Volume 573, January 2015
Article Number A134
Number of page(s) 5
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201424442
Published online 14 January 2015

© ESO, 2015

1. Introduction

Supernova remnants (SNRs) are believed to be the sources of Galactic cosmic rays (GCRs) from the argument of energetics and the diffusive shock acceleration mechanism of particles expected in the shocks associated with the remnants (e.g., Blandford & Eichler 1987). This view is corroborated by recent TeV gamma-ray observations of SNRs by CANGAROO (e.g., Muraishi et al. 2000), HESS measurement (e.g., Aharonian et al. 2006), and the other experiments. Supernovae (SNe) are believed to occur once per 30 to 100 years in our Galaxy. Hence we must take into account the discreteness of SNe in space and time in the investigation of the propagation process of GCRs from parent SNRs to the solar system. However, the numerical codes widely used to calculate the propagation of the GCRs in the Galaxy such as GALPROP (e.g., Strong & Moskalenko 1998) and DRAGON (e.g., Evoli et al. 2008) may not necessarily be suitable to fully include this discreteness in the calculation of propagation of GCRs. The important effects of this discreteness imprinted in the nature of GCRs observed in the solar system have been pointed out by various authors (Higdon & Lingenfelter 2003; Taillet et al. 2004; Mertsch 2011; Blasi & Amato 2012; Bernard et al. 2012). However, this is the first time that simulations of the propagation of GCRs that fully take into account the stochasticity coming from the spatial and temporal discreteness of occurrences of SN explosion have been attempted.

In this paper, we propose a new and fully three-dimensional stochastic numerical method to calculate age distribution and path length distribution (PLD) of GCRs reflecting the discreteness of the occurrence of supernova explosions. Our stochastic method is based on the equivalence of a coupled set of stochastic differential equations (SDEs) to the Parker convection-diffusion equation describing the propagation of GCRs. The power of this SDE method has been demonstrated in the study of the solar modulation phenomena where the method revealed some important physical properties, which are difficult – if not impossible – to obtain by other numerical methods such as the distribution of the arrival time from the heliospheric boundary to the earth and the distribution of energy loss during this period (Yamada et al. 1998; Zhang 1999). Similarly, it is expected that our method based on SDEs will reveal new aspects of the propagation process of GCRs. In this study, we assume that three-dimensional diffusive propagations of GCR are isotropic and we do not impose free escape boundaries. We also assume axisymmetric configurations for the density distributions of the interstellar matter and for the surface density of supernovae. We present various results of simulations with SDEs including the energy dependence of the B/C ratio to be compared with recent observations by CREAM (Ahn et al. 2008), AMS-02 (Oliva & AMS Collaboration 2013), and PAMELA (Adriani et al. 2014).

2. Models and numerical simulation

In this study we numerically investigate the propagation of GCR protons which originated in SNRs in the Galaxy. We assume that the acceleration of GCRs in SNRs continues uniformly in time up to 105 years after the explosion of the parent supernova, although it has been suggested in recent studies that the escape time of GCRs from SNRs depends on their energy (e.g., Ohira et al. 2010; Caprioli et al. 2010; Drury 2011). The radius of SNRs at the age of 105 years is assumed to be 30 pc by following the Sedov model where the values of 1051 erg, 109 cm s-1, and 1 proton cm-3 for the total explosion energy, the velocity of the shock wave, and the ambient matter density, respectively, are adopted. We also assume that the supernova (SN) occurs randomly both in time and in space within the radius of 20 kpc from the center of the Galaxy. The frequency of occurrences of SN is assumed to be 3 times in 100 years in the Galaxy. The surface density of the SN rate has a galactocentric radial dependence scaled to the molecular gas given by Williams & McKee (1997), and has a Gaussian height distribution where αSN = 0.070 kpc (Boulares & Cox 1990).

The transport and distribution of GCRs is described by the Fokker-Planck equation (FPE) also called the Parker convection-diffusion equation, (1)where f(r,p,t) is the cosmic-ray distribution function, r and p indicate the coordinate values of a generic point in the phase space, t is the time, V is the velocity of the Galactic wind, κ is the spatial diffusion coefficient tensor, Dp is the diffusion coefficient in the momentum space, and ṗc ≡ dpc/ dt is the momentum gain or loss rate. The full set of SDEs equivalent to the FPE (1) is written as (2)where r and p indicate the position and the momentum of pseudo-particle respectively (hereafter we refer to “pseudo-particle” as just “particle” if needed), , and dWs and dWp are the Wiener processes given by the Gaussian distribution. In this study, we assume an isotropic three-dimensional diffusive propagation without imposing free escape boundaries and neglect any energy changes and the existence of the Galactic wind for simplicity. Then, the full set of SDEs (2) is reduced to a simple SDE (3)We solved this SDE backward in time (see Yamada et al. 1998; Zhang 1999). In other word, we follow the sample path of the particles (protons) with various fixed energies at the earth backward in time until the particle arrives at some active SNR.

If we assume an isotropic diffusion, the diffusion tensor κμν is reduced to a scalar κ represented as , where lm.f.p. and v are the diffusion mean free path and the velocity. The mean free path lm.f.p. is represented as a { p/ (1 GeV / c) } 2−δ with two parameters a and δ. We adopted the value (a [cm],δ) as , , and (2.8 × 1017,1) for the Kolmogorov-type, the Kraichnan-type, and the Bohm-type models, respectively. The value of a for each model is determined so that our result of the B/C ratio at 1 GeV coincides with the observed value by AMS-02 as shown later. This normalization may be allowed because the effect of the solar modulation on the ratio of B/C is estimated less than 6% with reasonable values for the modulation parameter according to the force field theory.

The path length X is calculated simultaneously with the SDE as dX = ρvdt where ρ is the density of the interstellar matter. The density ρ is calculated with the model of interstellar medium adopted by Higdon & Lingenfelter (2003). In this model, the interstellar mediums consist of five major components, the molecular gas, cold neutral gas, warm neutral gas, warm ionized gas, and hot tenuous medium. The density of the molecular gas and the other components are given by Williams & McKee (1997) and Ferrière (1998) respectively, as a function of the Galactocentric radius and the height from the Galactic plane.

In our algorithm, first we generate a list of SNe which might have occurred during the last 5 × 108 yr in our Galaxy so that their spatial position rSN and occurrence time tSN are randomly distributed to be consistent with the conditions as stated before. Next we start to follow the sample path r(t) with certain energy by the SDE with an appropriate time step Δt. The path length for each time step can be calculated by vρΔt. Sample paths start at the solar system which is located at (8.5 kpc,0 kpc,0 kpc) for a galactocentric Cartesian coordinate system (x,y,z), and then run backward in time until they hit some active SNR centered at the parent SN on the prepared list where the acceleration of CRs is still going on. If the sample path r(t) at some time t get into some SNR, in other words r(t) satisfies | r(t) − rSN | ≤ 30 pc, tSN − 105 yr ≤ t ≤ tSN, then we stop to calculate the sample path and we record the position, the arrival time, and the total path length, and then we restart the same procedure again for another particle. We follow the sample path up to 5 × 108 yr while the particle fails to hit some active SNR. We adopt this figure of 5 × 108 yr so that it should be long enough compared with the age of GCR ~ 15 Myr estimated in the frame of the leaky box model based on direct observations of radioactive secondary nuclei in the energy range of a few 100 MeV/n by ACE (Yanasak et al. 2001). The adequacy of this figure as the time limit is shown later by our simulation results. By this numerical experiment, we can simultaneously obtain the age distribution and the PLD for particles observed with various energies at the earth. We want to emphasize that the age of GCRs obtained by this method is not the escape time from the Galaxy but the time spent by the particles during their journey to the earth from the SNR where they were born1. Accordingly, the value of this age should depend on the location of observers in our Galaxy.

3. Results and discussions

The age distribution and the PLD of protons are shown in Figs. 1 and 2, respectively, for some fixed energies in the solar system overlaid with the expected distributions by the leaky box model. It can be seen that the resultant age and the path length are distributed over a significantly wide range depending on the observed energy even for particles with the same observed energy. These distributions for each energy are certainly not obtained directly by observations; however, they play a crucial role in determining the model of propagation process of GCRs by the weighted slab method. It may be difficult, if not impossible, to estimate these distributions by other methods than ours. As we stated before, it may be recognized that our choice of 5 × 108 yr as the time limit to follow sample paths is justified by the resultant distributions shown in Fig. 1 where almost all of sample paths are distributed in the range lower than 108 yr even for the lowest energy particles of 100 MeV, which have been revealed to have the longest mean age within a few representative energies shown in Fig. 1.

thumbnail Fig. 1

Calculated age distribution of GCR protons with some fixed energies in the solar system for the case of . The thick line, the long-dashed line, and the thin line indicate the age distribution of protons which arrive at the solar system with a fixed energy of 100 MeV, 10 GeV, and 1 TeV, respectively. The dotted line indicates the age distribution expected by the leaky box model with the escape lifetime of 20 Myr.

Open with DEXTER

thumbnail Fig. 2

Calculated PLD of GCR protons with some fixed energies in the solar system for the case of . For an explanation see the caption to Fig. 1. The dotted line indicates the PLD expected by the leaky box model with the escape path length of 10 g cm-2.

Open with DEXTER

It is remarkable that a rapid decrease or a cut-off can be seen in the distributions shown in Figs. 1 and 2 below a certain age or path length depending on the energy of particles. Obviously no such cuts can be seen for the leaky box model. The origin of the cuts is traced back to the discreteness of the occurrence of SN as the source of GCRs. The position of the cut-off shifts to a lower range when the energy of the particle becomes higher. These tendencies arise because the diffusion time of GCRs from their birth place to the solar system becomes shorter when the energy of the particle becomes high. The effects of the discreteness of the GCR sources on the age distribution and the PLD have been investigated by the myriad-source model by Higdon & Lingenfelter (2003) and Taillet et al. (2004); however, they failed to discover the existence of the cuts in the distributions because they averaged out the spatial position of SNe. The existence of cuts in the PLD was pointed out by Cowsik & Wilson (1975) from quite a different argument of the nested leaky box model. Garcia-Munoz et al. (1987) suggested a lack in short path lengths in the PLD from a comparison of the B/C ratio to the sub-Fe/Fe ratio. Our model succeeded in predicting the missing short path lengths quite naturally and quantitatively.

thumbnail Fig. 3

Mean values of the age distribution of GCR protons. The thin line, the thick line, and the dashed line indicate the mean values for the case of , , and 1, respectively. The error bars shown for some points indicate typical statistical errors of ± 3σ in our simulation.

Open with DEXTER

Both distributions shown in Figs. 1 and 2 for 100 MeV particles exhibit clear peaks near the cut-off regions. We have checked that the peaks are due to particles that originated in several relatively young nearby SNRs. As mentioned before, it is unfortunate that these remarkable features cannot be checked directly by observations. However, the effects may be reflected indirectly in observed ratios of secondary to primary nuclei especially for unstable secondary nuclei. Clear peaks are not seen in the distributions for 10 GeV and 1 TeV particles. This is understandable from the fact that the fraction of the particles from the same nearby SNRs that make the peaks for 100 MeV particles becomes smaller for the higher energy particles because the number of SNRs where these particles were born becomes higher. This is caused because the value of the diffusion coefficient increases with the particle energy.

We note that the ages obtained by observations for GCRs with certain energies may be the values averaged over distributions similar to those shown in Fig. 1. Calculated average ages to be compared with observations are shown in Fig. 3 for three cases with different forms of the diffusion mean free path. As anticipated, the ages decrease monotonically with energies above about 1 GeV. The mean value of ~5 × 107 yr around at 1 GeV agrees well with the age obtained by the observation of radioactive secondary nuclei (e.g., Yanasak et al. 2001).

thumbnail Fig. 4

Mean values of the PLD of GCR protons. The dotted line indicates the mean values of the PLD assumed by Jones et al. (2001). The other line symbols are as in Fig. 3.

Open with DEXTER

Similarly, mean values of PLD as usual notion are calculated from such distributions as shown in Fig. 2 for each energy. The results are shown in Fig. 4 for the same three cases as are shown in Fig. 3. One of the notable features in the PLD is that a peak value is attained at a few 100 MeV. Many authors have assumed that the PLD is piecewise continuous as a function of energy with a peak at around a few GeV to explain phenomenologically the observed energy dependence of secondary to primary ratios in GCR (e.g., Jones et al. 2001; Ptuskin 2012). As an example we overlaid the PLD by Jones et al. in Fig. 4. Our model reproduces reliably, though qualitatively, the PLD with such a form without introducing a break point by hand. The obtained PLD decreases monotonically with energies around 1 GeV depending on the form of the diffusion mean free path, as was anticipated.

Once we get the PLDs for protons with arbitrary energies, we can calculate abundance ratios of the secondary to primary elements at corresponding energies by the weighted slab method if we are allowed to assume the PLDs for the heavy elements are the same as those of protons. In the calculations, we should use the PLDs shown in Fig. 2 prepared in advance for each energy and should not use the averaged PLD shown in Fig. 4 because, as mentioned earlier, the PL is distributed in a wide range even for particles with the same energy. If we use the averaged PLD as is usually done, there are subtle differences in the ratios even though we do not show the details here. We present the calculated energy dependence of B/C ratios in Fig. 5 for the same three different forms of the mean free path as shown in Fig. 3 and Fig. 4 together with recent observational results. Here we calculated the B/C ratio in the energy range where the effect of solar modulation is estimated to be small by the force field theory. In our calculation, only C, N, and O are assumed to be the parent nuclei of B. The source abundance of each parent nuclei was taken from Strong & Moskalenko (2001). Relevant nuclear data were referred to Garcia-Munoz et al. (1987), Webber et al. (2003), and Ramaty et al. (1997).

The numerical values multiplied to each form of the diffusion mean free paths presented in Fig. 5 were chosen so as to reproduce best the value of B/C at 1 GeV/n obtained by AMS-02 as mentioned earlier. As clearly seen, the observational results are most reliably reproduced in the whole range of energy when the power law index of momentum dependence is 1/2. This result indicates that the value of δ is 3/2 and supports the Kraichnan model for the interstellar turbulence.

We did not address here discussions on the B/C in the energy range below ~1 GeV / n where the ratio may be influenced to somewhat high degrees by the convection from the Galactic wind (e.g., Strong & Moskalenko 1998), the reacceleration in the interstellar medium (e.g., Simon et al. 1986; Giler et al. 1989; Seo & Ptuskin 1994), and the various plasma processes reflected in the diffusion coefficient (e.g., Shalchi et al. 2004; Evoli & Yan 2014).

So far we have only discussed nuclear components in the GCR; however, we should not overlook the effects of the discreteness of SNe occurrences on electron components in the GCR. Kobayashi et al. (2004) pointed out the possible peculiar energy spectrum of electrons in the TeV energy range where electrons from only a few nearby SNRs could contribute because of severe energy losses suffered by electrons from other distant SNRs. These interesting features will be revealed by scheduled experiments such as CALET (Akaike et al. 2010) in the very near future.

thumbnail Fig. 5

Derived B/C ratios for three cases with different forms of the diffusion mean free path overlaid with the observational data by HEAO3, CREAM (compilation by Maurin et al. 2014, and the references cited therein), AMS-02 (Oliva & AMS Collaboration 2013), and PAMELA (Adriani et al. 2014). The line symbols are as in Fig. 3.

Open with DEXTER

4. Summary

We demonstrated a new numerical method to investigate the propagation process of GCRs in our Galaxy which is based on the equivalence of the diffusion convection equation to a coupled stochastic differential equations. Assuming SNRs are the birthplace of the GCRs, we examined the effects on the nature of the GCRs arriving at the solar system caused by spatial and temporal discreteness of the SN occurrences in a fully 3D scheme though for a simplified model of an isotropic diffusion.

The resultant PLD for particles arriving at the solar system with the same energy exhibits a distribution spread over a wide range of path length values with a sharp cut below some lower path length depending on the energy. The existence of PLDs with a cut has been anticipated by studies of the abundance ratios of secondary to primary elements. The cut is clearly brought about by the discrete occurrences of SN explosions and accordingly may not be attained by models with continuous source distributions as usually assumed. The age distribution obtained simultaneously with PLDs exhibits a characteristic distribution with a cut also caused by the discreteness of SN explosions. The calculated average age agrees well with the value obtained by observations of radioactive secondary nuclei.

The energy dependence of B/C ratios from the obtained PLDs was calculated by the weighted slab method for a few forms of the diffusion mean free path as a function of momentum. It was found that the model with the mean free path proportional to the square root of particle momentum reproduces best the recent observational results. This finding supports the idea of the Kraichnan model for the interstellar turbulence.

Results by simulations which take into account important factors not considered in this work such as the Galactic wind, various energy changes including reaccelerations, and the more detailed global structure of the Galaxy with the magnetic fields will be presented in another paper in the near future.

1

After we submitted our paper, Lipari (2014) presented an interesting article in arXiv where he develops the concept of cosmic ray ages that is similar to ours.

Acknowledgments

We thank J. Nishimura, T. Kamae, T. Shibata and S. Torii for their invaluable comments for this work. We are grateful to the referee for critical comments and useful suggestions. This work is supported by JSPS KAKENHI Grant Number 26800145 (S.M.) and 22540264 (S.Y.).

References

All Figures

thumbnail Fig. 1

Calculated age distribution of GCR protons with some fixed energies in the solar system for the case of . The thick line, the long-dashed line, and the thin line indicate the age distribution of protons which arrive at the solar system with a fixed energy of 100 MeV, 10 GeV, and 1 TeV, respectively. The dotted line indicates the age distribution expected by the leaky box model with the escape lifetime of 20 Myr.

Open with DEXTER
In the text
thumbnail Fig. 2

Calculated PLD of GCR protons with some fixed energies in the solar system for the case of . For an explanation see the caption to Fig. 1. The dotted line indicates the PLD expected by the leaky box model with the escape path length of 10 g cm-2.

Open with DEXTER
In the text
thumbnail Fig. 3

Mean values of the age distribution of GCR protons. The thin line, the thick line, and the dashed line indicate the mean values for the case of , , and 1, respectively. The error bars shown for some points indicate typical statistical errors of ± 3σ in our simulation.

Open with DEXTER
In the text
thumbnail Fig. 4

Mean values of the PLD of GCR protons. The dotted line indicates the mean values of the PLD assumed by Jones et al. (2001). The other line symbols are as in Fig. 3.

Open with DEXTER
In the text
thumbnail Fig. 5

Derived B/C ratios for three cases with different forms of the diffusion mean free path overlaid with the observational data by HEAO3, CREAM (compilation by Maurin et al. 2014, and the references cited therein), AMS-02 (Oliva & AMS Collaboration 2013), and PAMELA (Adriani et al. 2014). The line symbols are as in Fig. 3.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.