About cosmic ray transport

The modelling of cosmic-ray (CR) diffusion in the Galaxy includes the solution of the transport equation with a given source distribution and boundary conditions for all CR species. The transport equation describes diffusion, convection by the hypothetical galactic wind, energy losses, and possible distributed acceleration (energy gain). In addition, nuclear collisions with interstellar gas atoms resulting in the production of secondary energetic particles and the decay of radioactive isotopes should be taken into account when considering the cosmic-ray nuclei.

The most often used propagation model, the flat halo diffusion model, has a simple geometry, but reflects the most essential features of the real system (Ginzburg & Ptuskin 1976). In this model, it is assumed that the Galaxy has the shape of a cylinder with a radius R (~20 kpc) and total height 2H (H > 1 kpc). The cosmic ray sources are distributed within a thin disk about the Galactic mid-plane having characteristic thickness ~200 pc. The Sun is at a distance ~8.5 kpc from the center of the Galaxy. The diffusion of cosmic rays averaged over the scale of a few hundred parsecs is isotropic. The particles escape freely through the halo boundaries into intergalactic space where the density of cosmic rays is assumed to be negligible.

The wealth of information contained in the cosmic-ray isotopic abundances makes it possible to study various aspects of CR acceleration and propagation in the interstellar medium as well as the source composition. Stable secondary nuclei tell us about the diffusion coefficient and Galactic winds (convection) and/or re-acceleration in the interstellar medium. Long-lived radioactive secondaries provide constraints on global Galactic properties such as, e.g., the Galactic halo size. Abundances of K-capture isotopes, which decay via electron K-capture after attaching an electron from the interstellar medium, can be used to probe the gas density and acceleration time scale. Details of the Galactic structure, such as the non-uniform gas distribution (spiral arms, the Local Bubble), radiation field and magnetic field (regular and random) distributions, and close SNRs may also affect the local fluxes of cosmic-ray particles.

It has been recently realized that direct information about the fluxes and spectra of cosmic-ray in distant locations is provided by the Galactic diffuse gamma rays, therefore, complementing the local cosmic-ray studies, but this connection requires extensive modeling and is yet to be fully explored in detail.

(For more details see Publications)

ISOTOPES IN COSMIC RAYS

The abundances of stable secondary nuclei (Li, Be, B, Sc, Ti, V) in CR allows one to determine the ratio (halo size)/(diffusion coefficient) and the incorporation of radioactive secondaries (⁴Be, ²⁶Al, ³⁶Cl, ⁵⁴Mn) is used to find the diffusion coefficient and the halo size separately. The derived source abundances of CR may provide some clues to mechanisms and sites of CR acceleration. However, the interpretation of CR data, e.g., the sharp peak in the secondary/primary nuclei ratio, is model dependent. The leaky-box model fits the secondary/primary ratio by allowing the path-length distribution vs. rigidity to vary. The diffusion models are more physical and explain the shape of the secondary/primary ratio in terms of diffusive reacceleration (distributed energy gain) in the ISM, convection by the Galactic wind, or by the damping of the interstellar turbulence by CR on a small scale.

The secondary/primary nuclei ratio is strongly affected by the value of the diffusion coefficient and its energy dependence. The larger value of the diffusion coefficient leads to the lower ratio because of faster escape of primary nuclei from the Galaxy results in smaller amount of secondaries produced and vise versa. The relative abundance of radioactive secondaries (e.g., ¹⁰Be/⁹Be) is used to derive the Galactic halo size, the Galctic volume filled with CR. The larger the halo the longer it takes for radioactives to reach us. The ratio ¹⁰Be/⁹Be thus decreasing with the halo size increase.

K-capture isotopes in CR (e.g., ⁴⁹V, ⁵¹Cr) can serve as important energy markers and can be used to study the energy-dependent effects such as diffusive reacceleration in the interstellar medium and heliospheric modulation. Such nuclei usually decay via electron-capture and have a short lifetime in the medium. In CR they are stable or live essentially longer as they are created bare by fragmentation of heavier nuclei while their β⁺ decay mode is suppressed. At low energies, their lifetime depends on the balance between the processes of the electron attachment from the interstellar gas and stripping. The probability of attachment is strongly energy-dependent increasing toward low energies, while the probability of stripping is flat. This makes the abundances of K-capture isotopes in CR energy-dependent. Without the electron K-capture, the ratio ⁵¹V/⁵¹Cr in CR would be flat because both nuclei are secondary. The electron K-capture transforming ⁵¹V/⁵¹Cr makes the ratio increasing at low energies. The reacceleration (energy gain) increases it at low energies even further making more consistent with the CR data. The solar modulation acts in the opposite direction. In the interstellar medium the ratio ⁵¹V/⁵¹Cr is very steep, while the solar modulation flattens it to some degree dependently on the modulation level.

Study of the light stable nuclei in CR (Li--O) allows us to determine propagation parameters averaged over the large Galactic region, but the local interstellar medium is not necessary the same and the local propagation parameters may significantly differ. A way to study the local medium is to look at radioactive isotopes with shorter lifetimes (e.g., ¹⁴C) and the abundances of heavy nuclei Z>28. Large fragmentation cross sections of heavy nuclei result in relatively small "collection area".

The CR source composition is derived from direct CR measurements by correcting for the effects of propagation, spallation, and solar modulation. The elements with low first-ionization potential (FIP) appear to be considerably more abundant in CR sources relative to the high-FIP elements, when compared with the solar system material. This implies that the source material for CR is the atmosphere of stars, where typical temperatures are ~10⁴ K. A strong correlation between FIP and "volatility", most of low-FIP elements are refractory while high-FIP elements are volatile, suggests that the material of CR may originate in the interstellar dust pre-accelerated by shock waves. Na, Ga, Rb and some other elements Z>28 break this correlation. CR data tend to prefer volatility over FIP, but uncertainties in the derived source abundances (cross sections!) are preventing an unambiguous solution.

Isotopic pecularities of CR composition are also important. The ²²Ne/²⁰Ne enrichment might tell us that CR are produced in cores of superbubbles created by multiple correlated SNe.

For more details and references see Moskalenko & Strong (2004).

DIFFUSE GAMMA RAYS

The diffuse gamma-ray continuum emission is the dominant feature of the gamma-ray sky. This emission in the range 50 keV - 50 GeV has been systematically studied in the experiments Oriented Scintillation Spectrometer Experiment (OSSE), Imaging Compton Telescope (COMPTEL), and Energetic Gamma-Ray Experiment Telescope (EGRET) on the CGRO as well as in earlier experiments. The diffuse gamma-ray continuum emission arises in interactions of energetic CR particles, mainly protons, alpha, electrons, and positrons, with interstellar gas through the decay of neutral pions and bremsstrahlung, and with radiation field through the inverse Compton (IC) scattering. The observation of diffuse gamma rays provides the direct information necessary to access the CR spectra and densities in the distant locations. However, what we see is the sum of, at least, three independent components (pion decay, IC, bremsstrahlung) which are produced in hadronic and leptonic interactions, besides the gamma-ray intensity in certain direction in the sky is the integral along the line of sight while the CR particle spectra, gas density, and interstellar radiation field, all vary with the position. A fit of the gas (H2, H I, H II) templates, dust, and generic IC templates though possible does not produce meaningful information. A proper modeling of the diffuse emission including disentangling different components of the diffuse emission requires extensive modeling and developed propagation models. The launch in 2007 of the Gamma-ray Large Area Space Telescope (GLAST) will tremendously increase the quality and accuracy of the gamma-ray data.

GALPROP computes the gas-related γ - ray intensities from the emissivities as a function of (R,z,Eγ) using the column densities of H I and H2 for Galactocentric annuli based on 21-cm and CO surveys. Neutral pion production is calculated using a formalism by Dermer (1986a,b) as described in Moskalenko & Strong (1998); bremsstrahlung is calculated using a formalism by Koch & Motz (1959) as described in Strong et al. (2000). The inverse Compton scattering is treated using the formalism for an anisotropic radiation field developed by Moskalenko & Strong (2000a).

For more details and references see Moskalenko, Strong, & Reimer (2004), McEnery, Moskalenko, & Ormes (2004).

TRANSPORT EQUATION

The transport equation includes a Galactic wind (convection), diffusive reacceleration in the interstellar medium, energy losses, nuclear fragmentation, and decay:

The source term includes primary CR particles and secondary particles and isotopes produced in the course of CR propagation.

GALPROP solves this transport equation with a given source distribution and boundary conditions for all cosmic-ray species. The numerical solution of the transport equation is based on a Crank-Nicholson (Press et al. 2002) implicit second-order scheme. The spatial boundary conditions assume free particle escape. Since we have a 3-dimensional (R,z,p) or 4-dimentional (x,y,z,p) problem (spatial variables plus momentum) we use "operator splitting" to handle the implicit solution. The reaction network is solved starting at the heaviest nucleus (i.e., ⁶⁴Ni). The propagation equation is solved, computing all the resulting secondary source functions, and then proceeds to the nuclei with A-1. The procedure is repeated down to A=1. To account for some special β^--decay cases (e.g., ¹⁰Be->¹⁰B) the whole loop is repeated twice.

For a given halo size the diffusion coefficient as a function of momentum and the reacceleration or convection parameters is determined by boron-to-carbon ratio data. The spatial diffusion coefficient is taken as D_xx = β D₀(ρ/ρ₀)^δ if necessary with a break (δ=δ_1,2 below/above rigidity ρ₀), where the factor β(=v/c) is a consequence of a random-walk process. For the case of reacceleration the momentum-space diffusion coefficient D_pp is related to the spatial coefficient D_xx (Berezinskii et al., 1990; Seo & Ptuskin, 1994), where δ=1/3 for a Kolmogorov spectrum of interstellar turbulences. The convection velocity (in z-direction only) V(z) is assumed to increase linearly with distance from the plane (dV/dz>0 for all z); this implies a constant adiabatic energy loss. The linear form for V(z) is consistent with cosmic-ray driven MHD wind models (Zirakashvili et al. 1996).

SOURCE DISTRIBUTION

FULL SIZE

The distribution of cosmic-ray sources (Strong & Moskalenko 1998) is chosen to reproduce the cosmic-ray distribution determined by analysis of EGRET γ-ray data (histogram, Strong & Mattox 1996). The injection spectrum of nucleons is assumed to be a power law in momentum. The distribution derived from CR is different from that of SNR  presumably CR sources. Even with 20 kpc halo, SNR produce CR distribution distinctly different from the measurements.

A solution to the apparent discrepancy between the radial gradient in the diffuse Galactic gamma-ray emissivity and the distribution of SNRs (Case & Bhattacharya 1998), believed to be the sources of cosmic rays, was proposed by Strong et al. (2004). The problem is plausibly solved if the X-factor (≡ n_H2 / W_CO) increases by a factor of 5-10 from the inner to the outer Galaxy, as expected from the Galactic metallicity gradient. The source distribution required to match the radial gradient of gamma rays can be reconciled with the distribution of SNRs as traced by current studies of pulsars (Lorimer 2004).

INTERSTELLAR HYDROGEN DISTRIBUTION

FULL SIZE

The gas distribution is important in production of secondary CR species and energy losses. It consists of 3 different components with distinctly different distributions. Molecular gas H₂ concentrated mostly in the plane. The curves show the distributions at different altitudes above the galactic plane (0, 0.1, 0.2 kpc). The atomic hydrogen H I distribution is broader. Ionized hydrogen H II has a very broad distribution but it is only a small fraction of the gas mass. Still it is important for diffuse gamma-ray emission since it contributes to the flux at high galactic latitudes.

The interstellar hydrogen distribution uses H I and C0 (tracer of H₂) surveys and information on the ionized component (Moskalenko et al., 2002); the helium fraction of the gas is taken as 0.11 by number. The H₂ gas number density is defined in the form of table (Bronfman et al., 1988), which is interpolated linearly, and the conversion factor is taken as X ≡ n_H2 /  W_C0 = 1.9 × 10²⁰ mols·cm^-2/(K km s^-1) (Strong & Mattox, 1996). The H I gas number density in the Galactic plane is defined by a table (Gordon & Burton,1976) which is renormalized to agree with the total integral perpendicular to the plane by Dickey & Lockman (1990). The z-dependence is calculated using the approximation by Dickey & Lockman (1990) for R< 8 kpc, using the approximation by Cox et al., (1986) for R>10 kpc, and interpolated in between. The ionized component H II (atom cm^-3) is calculated using a cylindrically symmetrical model (Cordes et al., 1991). See Moskalenko et al. (2002) for more details.

ENERGY LOSSES

Energy loss time-scales (Strong & Moskalenko, 1998) of nucleons (left) and electrons and positrons (right) in neutral and ionized hydrogen. The curves are computed for gas densities n_H =n_H II =0.01 cm^-3, and equal energy densities of photons and magnetic field U=U_B =1 eV cm^-3 (in the Thomson limit). In the left panel solid lines correspond to ionization losses and dashed lines to Coulomb losses (the dotted line is that for protons).

FULL SIZE

INTERSTELLAR RADIATION FIELD

Calculations of the spectrum of γ-rays generated by inverse Compton scattering and electron energy losses require the full ISRF as function of (R,z,ν) is required, which is not available in the literature. Our ISRF calculation uses emissivities based on stellar populations and dust emission. The infrared emissivities per atom of Hi, H2 are based on COBE/DIRBE data from Sodroski et al. (1997), combined with the distribution of Hi, H2. The spectral shape is based on the silicate, graphite and PAH synthetic spectrum using COBE data from Dwek et al. (1997). For the distribution of the old stellar disk component we use the model of Freudenreich (1998) based on the COBE/DIRBE few micron survey. The stellar luminosity function is taken from Wainscoat et al. (1992). For each stellar class the local density and absolute magnitude in standard optical and near-infrared bands is given, and these are used to compute the local stellar emissivity by interpolation in wavelength. The z-scale height for each class and spatial functions (disk, halo, rings, arms) given by Wainscoat et al. (1992) then give the volume emissivity as a functiolength. All their main-sequence and AGB types were explicitly included.

NUCLEAR REACTION NETWORK AND CROSS SECTIONS

FULL SIZE

The study of transport of CR nuclear component requires consideration of nuclear spallation. Hundreds of isotopes are included in the calculation of nuclear fragmentation and production of secondary isotopes. The most of the reactions can be devised into direct channels where we get the stable nucleus (shown in black) and indirect channels where there are radioactive intermediate states.

Most of decays are simple,  they are β^+/−-decays which move a nucleus in diagonal direction, and emission of a proton and neutron, which move a nucleus in vertical or horizontal directions. There are however several dozens of more complicated decays which require special handling.The radioactive isotopes, namely ¹⁰Be, ²⁶Al, ³⁶Cl, ⁵⁴Mn are long lived isotopes with the half-life time of about 1 Myr and they are very important in astrophysics of CR and used as radioactive clocks.

The code includes cross-section measurements and energy dependent fitting functions (Strong & Moskalenko 2001). The nuclear reaction network is built using the Nuclear Data Sheets. The isotopic cross section database is built using an extensive T16 Los Alamos compilation of the cross sections (Mashnik et al. 1998) and modern nuclear codes CEM2k and LAQGSM (Mashnik et al., 2004). The most important isotopic production cross sections (  ²H, ³H, ³He, Li, Be, B, Al, Cl, Sc, Ti, V, Mn) are calculated using the authors' fits to major production channels (e.g., Moskalenko et al. 2001, Moskalenko et al. 2003, Moskalenko & Mashnik 2003). Other cross sections are calculated using phenomenological approximations by Webber et al. (1990) (code WNEWTR.FOR version of 1993 and updated version of 2003) and/or Silberberg and Tsao (code YIELDX_011000.FOR version of 2000) renormalized to the data where it exists. The cross sections on the He target are calculated using a parametrization by Ferrando et al. (1988).

For pp and pA total inelastic cross sections we adapted parametrizations by Tan & Ng (1983a) and Letaw et al. (1983), or Barashenkov-Polanski code CROSEC. Total fragmentation cross section of  ⁴He is calculated using the authors' fit to the data. For a summary of inelastic cross sections employed, see (Mashnik et al., 2002). The reaction network is solved starting at the heaviest nuclei (i.e., ⁶⁴Ni). The propagation equation is solved, computing all the resulting secondary source functions, and then proceeds to the nuclei with A - 1. The procedure is repeated down to A = 1. Our preliminary results for all CR species Z ≤ 28 are given in Strong & Moskalenko (2001) and re-evaluation of the radioactive isotopes of Be, Al, Cl, Mn is given in Moskalenko et al. (2001), and new calculation of isotopes of H, He, Li, Be, B in CR is presented in (Moskalenko et al., 2003).

SECONDARY ANTIPROTONS

FULL SIZE

The code calculates production and propagation of secondary antiprotons as described in Moskalenko et al. (1998) and Moskalenko et al.(2001b). Antiproton production in pp - collisions has been calculated using the parametrization of the invariant p - production cross section given by Tan & Ng(1983b). The antiproton production by nuclei with Z ≥ 2 is calculated using effective nuclear factors; the latter are computed using the Monte Carlo event generator DTUNUC (Roesler et al., 1998; Simon et al., 1998) or scaling factors similar to Gaisser & Schaefer (1992). For the inelastic proton cross sections we adapted parametrizations by Tan & Ng (1983a) and Letaw et al. (1983). The total pp inelastic cross section has been calculated using a fit from Tan & Ng (1983a) and parametrization by Groom et al. (2000). The antiproton absorption cross section on nuclear targets is calculated following to Moiseev & Ormes (1997).

SECONDARY POSITRONS AND ELECTRONS

Secondary positrons and electrons in CR are the final product of decay of charged pions and kaons which in turn created in collisions of cosmic-ray particles with gas. Pion production in pp-collisions is considered following a method developed by Dermer (1986a,b), which combines isobaric (Stecker, 1970) and scaling (Badhwar et al., 1977; Stephens & Badhwar, 1981) models of the reaction. Secondary positron and electron production is computed as described in Moskalenko & Strong (1998). The electron and positron distributions from the muon decay have been corrected according to Kelner et al. (2006).

HELIOSPHERIC MODULATION

As solar activity changes with a period of ~11 years, so does the intensity of Galactic cosmic rays, but in the opposite direction: with an inverse correlation. The solar modulation is very important since it changes the interstellar spectra of cosmic ray particles below ~20 GeV/nucleon during their propagation inside the heliosphere. These modulated spectra are the actual ones measured by balloon-borne and spacecraft instruments. The "solar modulation" is a combination of the effects of convection by the solar wind, diffusion, adiabatic cooling, drifts, and diffusive acceleration. The theory of solar modulation is far from being complete (Fichtner 2005).

Current heliospheric modulation models are based on the numerical solution of Parker's transport equation (Parker 1965). In the past 35 years many analytical and semi-analytical solutions have been found. The most often used are spherically symmetric solutions such as the "force-field" (Gleeson & Axford, 1968), and Fisk (1969) approximations. The Pioneer, Voyager, and Ulysses missions contributed significantly to understanding the global aspects of modulation and limiting the number of free parameters in models. Recent progress in the modeling of heliospheric cosmic ray transport has been stimulated by the new data provided by the Voyager 1 spacecraft which is approaching the outer boundary of the solar system,the solar wind termination shock and the heliosheath (Stone et al. 2005). This has resulted in the development of much more realistic time-dependent 3-dimensional hydrodynamic models of the heliosphere (Ndiitwani et al. 2005, Langner et al. 2006).

NUMERIC SOLUTION OF PROPAGATION EQUATION

The diffusion, reacceleration, convection and loss terms in eq. (1) can all be finite-differenced for each dimension (R,z,p) in the form

∂ψi	=	ψit+Δt − ψit	=	α1ψi-1t+Δt − α2ψit+Δt + α3ψi+1t+Δt	+	qi
∂t		Δt		Δt

where all terms are functions of (R,z,p). In the Crank-Nicholson implicit method (Press et al. 2002) the updating scheme is

ψit+Δt = ψit + α1ψi-1t+Δt − α2ψit+Δt + α3ψi+1t+Δt + qiΔt

Thus the tridiagonal system of equations,

− α1ψi-1t+Δt + (1 + α2)ψit+Δt - α3ψi+1t+Δt = ψit + qiΔt

Table 1: Coefficients for the Crank-Nicholson method.

Process	Coordinate	α1 /Δt	α2 /Δt	α3 /Δt
Diffusion	R	Dxx 2Ri – ΔR 2Ri (ΔR)2	Dxx 2Ri 2Ri (ΔR)2	Dxx 2Ri + ΔR 2Ri (ΔR)2
	z	Dxx /(Δz)2	2Dxx /(Δz)2	Dxx /(Δz)2
Convectiona	z>0(V>0)	V(zi-1 )/Δz	V(zi )/Δz	0
	z<0(V<0)	0	− V(zi )/Δz	− V(zi+1 )/Δz
	p( dV >0) dz	0	1 dV pi 3 dz Pii+1	1 dV pi+1 3 dz Pii+1
Diffusive reaccelerationa	p	− Dpp,i− Dpp,i-1 + Pii−12 + 2 ( Dpp,i + Dpp,i−1 ) Pii−1 Pi−1i+1 pi–1	− Dpp,i − Dpp,i−1 + Pii−12 + 2Dpp,i ( 1 + 1 ) + Pi−1i+1 Pii+1 Pii−1 + 2Dpp,i Pii−1 pi	2Dpp,i+1 Pi−1i+1Pii+1
Energy loss	p	0	− 1 dpi Pii+1 dt	− 1 dpi+1 Pii+1 dt
Fragmentation	R,z,p	0	1/3 τf	0
Radioactive decay	R,z,p	0	1/3 τr	0
a Pji ≡ pi – pj

DIFFUSION IN R

As an example, the coefficients for the radial term are derived here.

1

∂

(

RDxx

∂ψ

)

=

2

Dxx

{

Ri+1

ψi+1 – ψi

– Ri-1

ψi – ψi-1

}

R

∂R

∂R

Ri

Ri+1 – Ri-1

Ri+1 – Ri

Ri – Ri-1

Setting R_i+1 − R_i = R_i − R_i-1 = ΔR, one can obtain the following expressions in terms of our standard form (eq. [1])

α1/Δt = Dxx 2Ri – ΔR 2Ri (ΔR)2

α2/Δt = Dxx 2Ri 2Ri (ΔR)2

α3/Δt = Dxx 2Ri + ΔR 2Ri (ΔR)2

DIFFUSIVE REACCELERATION

In terms of 3-D momentum phase-space densityƒ(p), the diffusive reacceleration equation is

∂ƒ(p)	= divp •	[	Dppgradp ƒ(p)	]	=	1	∂	[	p2Dpp	∂ƒ(p)	]
∂t						p2	∂p			∂p

The distribution is asumed isotropic so ƒ(p) = ƒ(p), where p = |p|. First we rewrite the equation in terms of ψ(p) = 4πp²ƒ(p) instead of ƒ(p) and expand the inner differential:

∂ψ

=

∂

[

p2Dpp

∂

ψ

]

=

∂

Dpp

[

∂ψ

−

2ψ

]

∂t

∂p

∂p

p2

∂p

∂p

p

The differencing scheme is then

2

[

Dpp,i+1

(

ψi+1 – ψi

–

2ψi+1

)

–

Dpp,i-1

(

ψi – ψi-1

–

2 ψi-1

)]

pi+1 – pi-1

pi+1 – pi

pi+1

pi – pi-1

pi-1

In terms of our standard form (eq.[1]) the coefficients for reacceleration are

α1/Δt =	2 Dpp,i-1	(	1	+	2	)
	pi+1 – pi-1		pi – pi-1		pi-1

α2/Δt =	2	(	Dpp,i+1	+	Dpp,i-1	)
	pi+1 – pi-1		pi+1 – pi		pi – pi-1

α3/Δt =	2 Dpp,i+1	(	1	+	2	)
	pi+1 – pi-1		pi+1 – pi		pi+1

One more scheme (#2) comes from further detalization

dψ	=	∂ Dpp	∂ψ	+	Dpp	∂2ψ	– 2	∂	Dppψ
dt		∂p	∂p			∂p2		∂p	p

Here it is

α1/Δt = –	Dpp,i – Dpp,i-1	+	2 Dpp,i	+	2 Dpp,i-1
	( pi – pi-1) 2		( pi+1 – pi-1) ( pi – pi-1)		( pi – pi-1) pi-1

α2/Δt = –	Dpp,i – Dpp,i-1	+	2 Dpp,i	(	1	+	1	)	+	2 Dpp,i
	( pi – pi-1) 2		pi+1 – pi-1		pi+1 – pi		pi – pi-1			( pi – pi-1) pi

α3/Δt =	2 Dpp,i+1
	( pi+1 – pi-1) ( pi+1 – pi)