Download Theory Manuals Examples Terms of Use F.A.Q. 
About cosmic ray transportThe modelling of cosmicray (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 cosmicray nuclei. COSMIC RAY TRANSPORT

∂ψ  = q(r,p) + div (D_{xx} grad ψ − V ψ) +  ∂  p^{2}D_{pp}  ∂  ψ  −  ∂  [ψ  dp  −  p  (div · V) ψ ] −  ψ  −  ψ 
∂t  ∂p  ∂p  p^{2}  ∂p  dt  3  τ_{f}  τ_{r} 
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 cosmicray species. The numerical solution of the transport equation is based on a CrankNicholson (Press et al. 2002) implicit secondorder scheme. The spatial boundary conditions assume free particle escape. Since we have a 3dimensional (R,z,p) or 4dimentional (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., ^{64}Ni). The propagation equation is solved, computing all the resulting secondary source functions, and then proceeds to the nuclei with A1. The procedure is repeated down to A=1. To account for some special β^{}decay cases (e.g., ^{10}Be>^{10}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 borontocarbon ratio data. The spatial diffusion coefficient is taken as D_{xx} = β D_{0}(ρ/ρ_{0})^{δ} if necessary with a break (δ=δ_{1,2} below/above rigidity ρ_{0}), where the factor β(=v/c) is a consequence of a randomwalk process. For the case of reacceleration the momentumspace 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 zdirection 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 cosmicray driven MHD wind models (Zirakashvili et al. 1996).
The distribution of cosmicray sources (Strong & Moskalenko 1998) is chosen to reproduce the cosmicray 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 gammaray 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 Xfactor (≡ n_{H2} / W_{CO}) increases by a factor of 510 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).
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_{2} 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 gammaray emission since it contributes to the flux at high galactic latitudes.
The interstellar hydrogen distribution uses H I and C0 (tracer of H_{2}) 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_{2} 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^{20} 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 zdependence 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 loss timescales (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).
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 nearinfrared bands is given, and these are used to compute the local stellar emissivity by interpolation in wavelength. The zscale 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 mainsequence and AGB types were explicitly included.
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 ^{10}Be, ^{26}Al, ^{36}Cl, ^{54}Mn are long lived isotopes with the halflife time of about 1 Myr and they are very important in astrophysics of CR and used as radioactive clocks.
The code includes crosssection 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 ( ^{2}H, ^{3}H, ^{3}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 BarashenkovPolanski code CROSEC. Total fragmentation cross section of ^{4}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., ^{64}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 reevaluation 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).
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 in CR are the final product of decay of charged pions and kaons which in turn created in collisions of cosmicray particles with gas. Pion production in ppcollisions 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).
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 balloonborne 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 semianalytical solutions have been found. The most often used are spherically symmetric solutions such as the "forcefield" (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 timedependent 3dimensional hydrodynamic models of the heliosphere (Ndiitwani et al. 2005, Langner et al. 2006).
The diffusion, reacceleration, convection and loss terms in eq. (1) can all be finitedifferenced for each dimension (R,z,p) in the form
∂ψ_{i}  =  ψ_{i}^{t+Δt} − ψ_{i}^{t}  =  α_{1}ψ_{i1}^{t+Δt} − α_{2}ψ_{i}^{t+Δt} + α_{3}ψ_{i+1}^{t+Δt}  +  q_{i} 
∂t  Δt  Δt 
where all terms are functions of (R,z,p). In the CrankNicholson implicit method (Press et al. 2002) the updating scheme is
ψ_{i}^{t+Δt} = ψ_{i}^{t} + α_{1}ψ_{i1}^{t+Δt} − α_{2}ψ_{i}^{t+Δt} + α_{3}ψ_{i+1}^{t+Δt} + q_{i}Δt 
Thus the tridiagonal system of equations,
− α_{1}ψ_{i1}^{t+Δt} + (1 + α_{2})ψ_{i}^{t+Δt}  α_{3}ψ_{i+1}^{t+Δt} = ψ_{i}^{t} + q_{i}Δt 
Process  Coordinate  α_{1} /Δt  α_{2} /Δt  α_{3} /Δt  

Diffusion  R 




z  D_{xx} /(Δz)^{2}  2D_{xx} /(Δz)^{2}  D_{xx} /(Δz)^{2}  
Convection^{a}  z>0(V>0)  V(z_{i1} )/Δz  V(z_{i} )/Δz  0  
z<0(V<0)  0  − V(z_{i} )/Δz  − V(z_{i+1} )/Δz  

0 



Diffusive reacceleration^{a} 
p 




Energy loss  p  0 



Fragmentation  R,z,p  0  1/3 τ_{f}  0  
Radioactive decay 
R,z,p  0  1/3 τ_{r}  0  
^{a} P_{j}^{i} ≡ p_{i} – p_{j} 
As an example, the coefficients for the radial term are derived here.
1  ∂  (  RD_{xx}  ∂ψ  )  =  2  D_{xx}  {  R_{i+1}  ψ_{i+1} – ψ_{i}  – R_{i1}  ψ_{i} – ψ_{i1}  } 
R  ∂R  ∂R  R_{i}  R_{i+1} – R_{i1}  R_{i+1} – R_{i}  R_{i} – R_{i1} 
Setting R_{i+1} − R_{i} = R_{i} − R_{i1} = ΔR, one can obtain the following expressions in terms of our standard form (eq. [1])



In terms of 3D momentum phasespace densityƒ(p), the diffusive reacceleration equation is
∂ƒ(p)  = div_{p} •  [  D_{pp}grad_{p} ƒ(p)  ]  =  1  ∂  [  p^{2}D_{pp}  ∂ƒ(p)  ] 
∂t  p^{2}  ∂p  ∂p 
The distribution is asumed isotropic so ƒ(p) = ƒ(p), where p = p. First we rewrite the equation in terms of ψ(p) = 4πp^{2}ƒ(p) instead of ƒ(p) and expand the inner differential:
∂ψ  =  ∂  [  p^{2}D_{pp}  ∂  ψ  ]  =  ∂  D_{pp}  [  ∂ψ  −  2ψ  ] 
∂t  ∂p  ∂p  p^{2}  ∂p  ∂p  p 
The differencing scheme is then
2  [  D_{pp,i+1}  (  ψ_{i+1} – ψ_{i}  –  2ψ_{i+1}  )  –  D_{pp,i1}  (  ψ_{i} – ψ_{i1}  –  2 ψ_{i1}  )] 
p_{i+1} – p_{i1}  p_{i+1} – p_{i}  p_{i+1}  p_{i} – p_{i1}  p_{i1} 
In terms of our standard form (eq.[1]) the coefficients for reacceleration are
α1/Δt =  2 D_{pp,i1}  (  1  +  2  ) 
p_{i+1} – p_{i1}  p_{i} – p_{i1}  p_{i1} 
α2/Δt =  2  (  D_{pp,i+1}  +  D_{pp,i1}  ) 
p_{i+1} – p_{i1}  p_{i+1} – p_{i}  p_{i} – p_{i1} 
α3/Δt =  2 D_{pp,i+1}  (  1  +  2  ) 
p_{i+1} – p_{i1}  p_{i+1} – p_{i}  p_{i+1} 
One more scheme (#2) comes from further detalization
dψ  =  ∂ D_{pp}  ∂ψ  +  D_{pp}  ∂^{2}ψ  – 2  ∂  D_{pp}ψ 
dt  ∂p  ∂p  ∂p^{2}  ∂p  p 
Here it is
α1/Δt = –  D_{pp,i} – D_{pp,i1}  +  2 D_{pp,i}  +  2 D_{pp,i1}  
( p_{i} – p_{i1}) ^{2}  ( p_{i+1} – p_{i1}) ( p_{i} – p_{i1})  ( p_{i} – p_{i1}) p_{i1} 
α2/Δt = –  D_{pp,i} – D_{pp,i1}  +  2 D_{pp,i}  (  1  +  1  )  +  2 D_{pp,i} 
( p_{i} – p_{i1}) ^{2}  p_{i+1} – p_{i1}  p_{i+1} – p_{i}  p_{i} – p_{i1}  ( p_{i} – p_{i1}) p_{i} 
α3/Δt =  2 D_{pp,i+1} 
( p_{i+1} – p_{i1}) ( p_{i+1} – p_{i}) 
Terms of Use  Page last modified: September 5, 2012 09:00:36 PDT 