National Aeronautics and

Space Administration

IS20

George C. Marshall Space Flight Center

Huntsville, Alabama 35812
The NASA STI Program...in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

Program Office plays a key part in helping NASA 

maintain this important role.

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for 

NASA’s scientific and technical information. The 

NASA STI Program Office provides access to 

the NASA STI Database, the largest collection of 

aeronautical and space science STI in the world. 

The Program Office is also NASA’s institutional 

mechanism for disseminating the results of its 

research and development activities. These results 

are published by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA programs and include extensive data 

or theoretical analysis. Includes compilations 

of significant scientific and technical data 

and information deemed to be of continuing 

reference value. NASA’s counterpart of peerreviewed formal professional papers but has less 

stringent limitations on manuscript length and 

extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical conferences, 

symposia, seminars, or other meetings sponsored 

or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, technical, 

or historical information from NASA programs, 

projects, and mission, often concerned with 

subjects having substantial public interest.

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results...even 

providing videos.

For more information about the NASA STI Program 

Office, see the following:

• Access the NASA STI program home page at 

<http://www.sti.nasa.gov>

• E-mail your question via the Internet to 

<help@sti.nasa.gov>

• Fax your question to the NASA STI Help Desk 

at 443 –757–5803

• Phone the NASA STI Help Desk at 

443 –757–5802

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076–1320
i

NASA/TP—2012–217464

A Review of the Ginzburg-Syrovatskii’s 

Galactic Cosmic-Ray Propagation 

Model and Its Leaky-Box Limit

A.F. Barghouty

Marshall Space Flight Center, Huntsville, Alabama

July 2012

National Aeronautics and

Space Administration

Marshall Space Flight Center • Huntsville, Alabama 35812
ii

Available from:

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076 –1320

443–757–5802

This report is also available in electronic form at

<https://www2.sti.nasa.gov/login/wt/>
iii

TABLE OF CONTENTS

 1. INTRODUCTION ........................................................................................................... 1

 2. THE GINZBURG-SYROVATSKII EQUATION ............................................................ 3

 3. THE DIFFUSION APPROXIMATION ......................................................................... 5

 4. PARTICLE-WAVE INTERACTIONS ............................................................................. 7

 5. GREEN’S FUNCTIONS SOLUTION TO THE GINZBURG-SYROVATSKII 

EQUATION ..................................................................................................................... 8

 6. TRYING OUT DIFFERENT G’S ................................................................................... 10

6.1 The Slab Model .......................................................................................................... 10

6.2 The Venerable Leaky-Box Model ............................................................................... 11

6.3 Nonstandard Leaky-Box G ......................................................................................... 14

 7. THE CONTACT POINT: xesc⇔ Dˆ ................................................................................ 16

 8. PUTTING THE BOX TOGETHER ................................................................................ 17

 9. SAMPLE NUMERICAL IMPLEMENTATION ............................................................ 20

10. REMARKS ...................................................................................................................... 29

REFERENCES ....................................................................................................................... 31
iv
v

LIST OF FIGURES

 1. A schematic (not to scale) showing the long-winded path of a high-energy 

cosmic-ray particle in the Milky Way as it makes its way out of the ‘box’ ................ 12

 2. The source function, eq. (54), for protons and alphas as used 

in the sample illustration .......................................................................................... 22

 3. The escape function, eq. (55), for protons and alphas as used in the illustration ...... 22

 4. Energy loss by protons traversing a hydrogen gas—Bethe-Block relation, 

eq. (56). Discrete points are as simulated by the code SRIM/TRIM ........................ 23

 5. The total reaction cross sections, eq. (58), for protons and alphas as used 

in the illustration ...................................................................................................... 24

 6. The spallation cross section, eq. (59), for alphas on ISM protons to produce 

secondary CR protons as used in this illustration .................................................... 24

 7. The arriving number densities for protons and alphas, eq. (63), as predicted 

by this sample illustration using xesc=7 g/cm2 and Φ = 400 MV .............................. 25

 8. Ratio of the secondary protons (produced by the spallation of the primary 

CR alphas in the ISM) to the primary CR protons .................................................. 26

 9. Ratio of arriving alphas to protons (including secondary protons) .......................... 26

10. Observed proton and alpha flux by BESS ................................................................ 27

11. Ratio of alphas to protons as observed by BESS ..................................................... 27
vi
vii

LIST OF ACRONYMS AND SYMBOLS

BESS balloon-borne with superconducting spectrometer

CR cosmic ray

GCR galactic cosmic ray

G-S Ginzburg-Syrovatskii

He helium

IMF interplanetary magnetic field

ISM interstellar medium

OB spectral type O or early type B (stars)

PDE partial differential equation

PDF probability density function

PLD path length distribution

QLT quasi-linear theory

SDE stochastic differential equation

SEP solar energetic particles
viii

NOMENCLATURE

A mass number

A area-like scale

a endpoint 

Bo	 uniform component of magnetic field

b endpoint 

C(x,t) drift term

c speed of light

Dˆ diffusion coefficient

D(x,t) diffusive term

dW Wiener process

E kinetic energy

e electron charge

G characteristic solution

G(x) Green’s function

hg characteristic height of the galactic disc

hh characteristic height of the galactic halo

I mean excitation energy

Ii source intensity of CR nuclide i 

J flux

j differential particle flux

k nuclide index

M species transformation matrix

m mass

mec2 mass of the electron (eV)

N macroscopic distribution function

n ISM number density

nˆ unit vector

n average density

ng disc’s matter density
ix

NOMENCLATURE (Continued)

ni constant that carries the relative abundance factor

O differential operator

p parameter; nonrelativistic momentum of the particle; probability

Q;q source term

R rigidity parameter 

r position vector

resc escape distance

S surface

Tesc escape time

Tmax maximum kinetic energy imparted

t time

V homogeneous volume

v velocity

W(E) energy loss function; stopping power

x amount of ISM matter; particle’s random position; independent variable

xesc escape path length

Z atomic number

z charge of the transversing particle

a exponent associated with turbulence spectrum

b Lorentz factor

s cross section

DE energy change

Dt time step

d density factor (Bethe-Block)

d



B	 fluctuating irregular component of magnetic field

g Lorentz factor

k Kolmogorov turbulence decay exponent

l nuclear interaction mean free path (g/cm2)

m pitch angle

r branching ratio
x

NOMENCLATURE (Continued)

t half-life of species

Φ potenial-like function; modulation potential (MV)

Ω	 solid angle

Subscripts

1 alpha particles

2 protons

esc escape

i coefficient for species; nuclide

j nuclide
1

TECHNICAL PUBLICATION

A REVIEW OF THE GINZBURG-SYROVATSKII’S GALACTIC COSMIC-RAY 

PROPAGATION MODEL AND ITS LEAKY-BOX LIMIT

1. INTRODUCTION

Galactic cosmic rays (GCRs), while mostly protons, are energetic nuclei representing perhaps all of the nuclides of the periodic table. Their origin remains unclear even as the centennial of 

their discovery by Victor Hess is celebrated in 2012. In his pioneering balloon-borne experiments, 

Hess discovered that the level of ionizing radiation was actually higher up in the atmosphere than 

at sea level (e.g., by a factor of 2 at a height of 5 km). He concluded that the increase was due 

to ‘radiation’ penetrating the atmosphere. Robert Millikan confirmed Hess’s discovery in 1925 and 

coined the term ‘cosmic rays’ (CRs). Hess was awarded the Noble prize in physics in 1936 for his 

discovery of ‘cosmic radiation.’ (Hess shared the 1936 prize with Carl Anderson for Anderson’s 

discovery of the positron.) The histories and disciplines of nuclear and particle physics have been 

interwined with CR physics ever since. The 1947 discovery of the pi-meson in ‘atmospheric CRs’ is 

a telling example of this common heritage.

Today, a host of balloon-borne, space-borne, as well as ground-based observations of (primary) CRs reveal they are nuclei of regular (i.e., baryonic) matter from the lightest, hydrogen, through 

the very heavy elements like uranium. In addition, CRs also appear to include electrons, positrons, 

and antiprotons. Their arrival, for the most part, i.e., better than 1 in 104, is isotropic. The chemical 

abundances—elemental as well as isotopic—of these primary CRs appear to reflect the chemical 

composition of massive stars’ winds (as opposed to solar system-like abundances, for example).

Their arrival spectrum near Earth’s orbit reveals an almost structureless power law spanning 

many decades of energy, down to ≈1018 eV. Primary CRs arriving at the boundaries of the heliosphere with energy below a few hundred MeV/nucleon are unable to penetrate the outward-flowing 

solar wind plasma. Those that do lose a fraction of their energy as they drift and diffuse into the 

solar wind—an effect known as solar modulation. Near the inner heliosphere and Earth’s orbit, the 

primary CR nuclei arriving at the top of the atmosphere interact with its atomic constituents (mostly 

carbon, nitrogen, and oxygen), losing (as they cascade through), in most cases, all of their kinetic 

energy while producing a host of secondary radiation and particles (like gamma rays, neutrons, and 

pions). Very few primary CRs make it to sea level. Their secondary products, however, like muons, 

are able to penetrate the hard surface of the Earth down to depths of tens to hundreds of meters.
2

From an astrophysical perspective, some of the outstanding fundamental questions regarding CRs pertain to the nature and location of their source(s) as well as to the precise mechanism(s) 

of their acceleration. Their journey after synthesis and acceleration—their propagation in the interstellar medium (ISM)—appears to be reasonably well understood. Currently, the so-called diffusive shock acceleration theory appears to be the accepted physical theory behind their acceleration. 

The ‘standard’ source is usually taken to be a supernova explosion. However, galactic superbubbles, 

where a collection of massive OB stars tend to congregate in close proximity in space and time, are 

increasingly becoming a more likely source.

The propagation stage of that cosmic journey will be the focus in this Technical Publication. 

The medium of propagation is the ISM and, for purposes here, this is idealized as a gas of mostly 

hydrogen and helium (He) with a number density on the order of 1 per cm3 and a kinetic temperature of ≈104 K, permeated by a weak (~μG) but turbulent magnetic field. The CRs themselves are 

expected to contribute to the energy density of the ISM (which is ≈2 eV/cm3), but these and similar 

effects will be ignored here. Finally, the particles’ transversal of this medium is describe only statistically, and classically, taking into account only basic energy gain and loss processes as well as regular 

(baryonic) particles’ sinks and sources.
3

2. THE GINZBURG-SYROVATSKII EQUATION

In transport and acceleration of energetic (suprathermal) charged particles in space and 

astrophysical plasmas, the relevant length scales are typically much larger than the particle’s Larmor 

radius. As a result, the motion of these particles is effectively determined by the structure of the electromagnetic field via both its regular and stochastic components. In the absence of particle-particle 

collisions (the low density limit), the electromagnetic field is also responsible for accelerating the particles up to relativistic energies. The fluctuating component of the field can scatter particles, especially 

when the particle’s Larmor radius is comparable to the wavelength of the scattering hydromagnetic 

wave, ‘resonant scattering.’ If scattering is frequent and strong, it can isotropize the particle’s density 

function and, in the diffusive limit, the motion of those particles can then be described statistically.

As CR nuclei traverse the ISM and interact with its constituents and fields, they can lose or 

gain energy as well as suffer nuclear interactions. The Ginzburg-Syrovatskii (G-S) equation1 follows 

the evolution of a CR particle density distribution function, N(E,



r,t), that evolves in space and time 

once created at t=0 and by some source somewhere in the galaxy. N is assumed to depend explicitly 

on the kinetic energy, E, (or momentum) of the CR particle. For a particular CR nuclide i, the G-S 

equation for Ni(E,



r,t) is a continuity equation in phase space and reads (eq. (14.8) in ref. 1): 

∂Ni

∂t − ∇⋅ Dˆ ( i

∇N) +

∂

∂E

b( iNi) − 1

2

∂2

∂E2 d( i

Ni) = Qi − piNi + pi

k

k<i

∑ Nk . (1)

On the left-hand side of the G-S equation, Dˆ

i is the diffusion coefficient (or tensor) for species 

i in the ISM, which, in general, can depend on E, t, and 

r. The coefficient or function bi=dE/dt

describes the systematic energy change, ∆E, while di describes the fluctuations in this change, i.e., di(E)

=d/dt(ΔE)

2. On the right-hand side of the G-S equation, Q (E,r,t) i

 is the source intensity, i I , of the 

CR nuclide i, where, for an isotropic distribution of particles, Ni(>E) = 4π Ii

(E) / v dE, E

∞

∫ where v is 

the particle’s speed. pi is the probability (per unit time) for a CR nuclide i to suffer a collision in the 

ISM with a cross section σi , i.e., pi=nvσi, where n is the ISM number density. In the last term on the 

right-hand side, pi

k is the probability per unit time for the appearance of a nuclide i due to collisions 

of nuclide k in the ISM, and hence, pi

k can be written as =nvσik, where σik, is, for example, the production cross section for nuclide k from spallation reactions of nuclide i. In this notation, index i=1 

refers to the nuclide with the highest atomic number. Note that the G-S equation also applies to CR 

electrons with some modifications.1

In eq. (1), terms related to further addition and subtraction of species can be added as well 

as terms related to further energy losses and gains. But note the clear distinction of space-time 

evolution, (∂N/∂t) and (∇ i(Dˆ ∇N), from chemical transformation in which energy appears only as 

a parameter rather than as an independent variable, if the energy loss or gain terms are ignored. 

Actually, for CR particles with kinetic energy greater than a few GeV/nucleon, energy loss (or gain) 
4

mechanisms can, for propagation purposes, be neglected to a good approximation. As will be seen 

shortly, the transformation part can be separated and rewritten in terms of the amount of ISM matter, x, traversed by the CR nucleus as 

dNi

S

dx = −λi

−1

Ni

S + λik−1

Nk

S

k<i

∑ + qi δ(x) , (2)

where, under the homogeneity assumption, x=nvt which is typically expressed in units of g/cm2, and 

n being the (average, in this case) matter density of the ISM. In equation (2), λi is the nuclear interaction mean free path for species i in the ISM and λik in the interaction mean free path to produce 

species k upon the interactions of species i with the ISM; λi and λik are also typically expressed in 

g/cm2. Note that λ ∝ 1/σ. Also implicity assumed in equation (2) is that Ni

s

(x=0)=qi, qi being the 

fraction of nuclide i in the CR source and that Ni

s is only a function of x.

Note that eq. (2) is a special case of eq. (1) for Dˆ

i =bi=di=0 and for Qi=Qi(t)=qiδ (t). Essentially then, eq. (1) follows the evolution of Ni

 as it evolves in time (or x) and diffuses in space while 

suffering, while not necessarily in parallel; however, nuclear transformations due to interactions in 

the ISM which can alter the particular species Ni describes (for a thorough treatment of nuclear 

spallation processes of CRs, see ref. 2.) The above space-time versus chemical transformation separation (which will allow for a general solution to the G-S equation to be written down) is made possible 

by a number of assumptions, most important of which is the assumption that the space-time and 

transport parameters Dˆ

i , bi

, di do not depend on nuclear species while the transformation parameters pi

 and k pi do not depend on either time or space coordinates. Before the general solution and 

various simplifying assumptions of the G-S equation are examined, it is worthwhile to appreciate, at 

least qualitatively, the justification as to what makes a diffusion-like equation appropriate to model 

the propagation of CRs in the ISM.
5

3. THE DIFFUSION APPROXIMATION

Phenomenologically, the propagation of GCRs in the ISM is thought of as an electrodynamic 

interaction of a relativistic charged particle in a turbulent, magnetized cosmic plasma. The subject, 

of course, is rich in history and accomplishments, as well as in complexity. It is typically formulated within the context of a wave-particle interaction using the tools of weak turbulence theory. On 

a lower level still, i.e., more formal, one can think of it as a problem (and a solution) belonging to 

time-dependent perturbation theory. None of these tools or theorems will be recited here. Rather, the 

aim is to gain insight into how one can abridge a notion about the motion of a fast charged particle 

in a turbulent field with the diffusion approximation so as to arrive at the G-S differential equation, 

i.e., making the passage from the basic notion of motion in a turbulent magnetic field to a governing 

equation (like the G-S equation) with reductionist and predictive powers.

Diffusive motion due to some randomness in the problem, e.g., Brownian motion, is a well 

studied notion.3 At some level, one envisages a particle making a small and discrete number of random walks, due to irregular scattering, for example. A phenomenological description amounts to 

writing down a probability density function that can give a measure of where the particle would be 

and with what probability. A Bernoulli distribution function, for example, for the particle’s motion 

in one-dimension and for a small number of discrete steps can deliver that. (Keep in mind that the 

Bernoulli probability density function (PDF) has an asymptotic Gaussian form in the limit of a 

large number of steps. This is helpful because the Gaussian form also happens to be a solution to 

the simple diffusion equation.) When one proceeds to formulate the problem for a large number of 

random walks, the process (solution) begins to take on a simple form. Further, instead of addressing 

the motion of a single particle and its PDF, one thinks of a distribution function of a large number 

of particles that share the same initial conditions. To make the passage then to a differential equation 

governing the motion (with the appropriate initial and/or boundary conditions), one assumes some 

time step ∆t that is large enough so that the particle will have time to make a large number of steps 

but short enough so that the rms of the displacement is small compared with large-scale motion of 

the particle. Now, in the limit ∆t→ 0 and by taking a Taylor-series expansion of our macroscopic 

distribution function N around ∆t, one can arrive3 at a differential equation for N that resembles 

a ̨diffusion equation, i.e., a partial differential equation (PDE) first order in time and second order in 

space.

One can also arrive at a similar PDE by invoking the macroscopic Fick’s law along with the 

continuity equation. Fick’s law reads: 



J = −Dˆ∇N , (3)

where 

J is flux and Dˆ is the diffusion coefficient (or tensor). The continuity equation for N reads 

∂N

∂t

+ ∇⋅



J = 0. (4)
6

Now, combining eq. (3) with eq. (4) for N yields:

∂N

∂t − ∇⋅(Dˆ∇N) = 0. (5)

One can add a convection term, a source term, energy loss/gain terms, particle loss/gain terms, etc. to 

eq. (5) to arrive at a PDE more or less resembling the G-S equation. Note that eq. (5) is a linear PDE 

for N. (The same linear PDE (for N above), however, can be nonlinear if formulated in the ‘wrong’ 

variable, e.g., for 

v where 

J = 

v N,4 as in eq. (28).)

A Langevin5,6 approach, in contrast, would specify instead, a stochastic differential equation 

(SDE) for the particle’s own random position, x, for example, as:5,6

dx = C(x,t) dt + D(x,t) dW , (6)

where for the scalar random variable x(t), C(x,t) is its drift (or mean) term and D(x,t) is its diffusive (or dispersive) term, and dW is a Wiener process (i.e., a Gaussian white noise). While the two 

approaches are numerically different, they are nonetheless considered equivalent.5
7

4. PARTICLE-WAVE INTERACTIONS

Being very qualitative here, a few words about the so-called quasi-linear approximation is 

in order. When a group of particles interact with a wave (magnetic field) and they remain in phase 

with the wave over many cycles, the interaction is termed a resonant wave-particle process. In this 

process, the particles do have time to exchange energy with the wave while the distribution function 

will only suffer slow changes. Such changes are attributed to a diffusion process termed quasi-linear 

diffusion. The process is kinetic in nature, with the characteristic plateau and tail formations in the 

distribution function. Quasi-linear theory (QLT) has been most successfully—but still not without 

some difficulties—applied to the transport of solar energetic particles (SEPs) as well as CR particles 

in the interplanetary magnetic field (IMF).7,8

In the case of the IMF, the magnetic field is modeled as having a uniform component 

Bo (an 

Archimedian spiral emanating from the Sun) and a fluctuating irregular component δ



B. The motion 

of a typical SEP or CR particle then in this field will comprise of a smooth, regular component along 



Bo and a chaotic component due to δ



B. The latter component results in the scattering of the particle 

in pitch-angle μ (μ is defined as the cosine of the angle between the particle velocity vector and the 

magnetic line of force). If the scattering is assumed to be ‘sufficiently strong,’ then a diffusion process 

is a good approximation for the evolution of N. To quantitatively describe this diffusion process, 

a Fokker-Planck equation can be constructed that follows the evolution of N. 

Determination of the diffusion coefficient from knowledge of the magnetic fluctuations is the 

essence of any statistical description of the transport process. In QLT formalism(s), the assumption 

is that those fluctuations affect the motion of the particle on timescales slower than the characteristic 

timescale of the motion in the ambient field in general. In other words, a perturbative approach is 

assumed to hold true and the motion can be looked at as being uniform to the zeroth approximation plus an irregular one to the first approximation. In the zeroth approximation, the particle ‘sees’ 

a stationary magnetic field, so a magnetostatic approximation can be used to derive the diffusion 

coefficient in terms of the observable statistical properties of the IMF like its mean value and power 

spectrum. This is an example of a Kolmogorov formalism wherein an incomplete knowledge of the 

nature and geometry of the magnetic fluctuations is replaced by a hierarchy of correlation functions 

and their Fourier transforms.

Other processes that can occur, but typically not included in (early) standard QLT of CR 

transport, are nonlinear wave-wave interactions and wave-particle-wave interactions.9 The former, 

also a resonant process but no particle resonances are involved, can give rise to the so-called parametric instabilities, where two waves beat together to interact with a third wave with the interaction 

becoming unstable above a certain threshold. The latter is both a resonant and a kinetic process 

where the particles keep in constant phase relative to the beat frequency of two waves, giving rise to 

the so-called plasma wave echoes.10
8

5. GREEN’S FUNCTIONS SOLUTION TO THE GINZBURG-SYROVATSKII EQUATION

The Green’s functions approach to a general solution of the G-S equation is appropriate here 

by virtue of the fact that the G-S equation, eq. (1), can be shown to satisfy the nonhomogeneous 

Sturm-Liouville equation which, in one independent variable x, reads: 

O[z(x)] + f (x) = 0 , (7)

where O is a differential operator of the form 

O[] = d

dx p(x) ddx

⎡

⎣

⎢ ⎤

⎦

⎥[] + q(x) , (8)

and p(x) and q(x) are some functions of x. Differential operators satisfying eq. (8) are called selfadjoint operators. One then seeks a solution, z(x), satisfying eq. (7) in addition to certain boundary 

conditions at endpoints a and b over the interval [a,b], with the help of another function (Green’s 

function). The sought-after Green’s function, G (x), is required to satisfy the homogeneous SturmLiouville equation: 

O[G(x)] = 0 (9)

in addition to other conditions.11,12 Most notable is that G has to satisfy the same boundary conditions that z(x) satisfies and that G be a continuous function of x while its first derivative with respect 

to x be discontinuous. In terms of G, the general (Green’s functions) solution to eq. (7) is now written

z(x) = dx′G(x,x′) f (x′) , a

b

∫ (10)

where the integrated-over variable x′ is not exactly a dummy variable; rather, it is the variable that 

G′ needs to be discontinuous in when p(x), eq. (8), is evaluated at.

Back to the G-S equation, a general, Green’s function solution13 for N satisfying eq. (1) can 

now be written as 

Ni(t,



r) = dxG(t,



r,x)Ni

s

(x) , 0

∞

∫ (11)

where x is an independent variable. Note, however, that x and t will no longer be independent under 

the homogeneity assumption and the Ni

s (x) are functions of x only that satisfy 

dNi

s

dx = −σi

Ni

s + σijN js + ci δ(x) .

j<i

∑ (12)
9

The source term is now a Dirac delta function in x with strength ci (to simulate relative abundance of 

the different species i at synthesis). Collecting the interaction terms together and ignoring the energyloss term, G is now required to satisfy 

∂

∂t

− ∇⋅(Dˆ∇) + nv

∂

∂x

⎡

⎣

⎢ ⎤

⎦

⎥G(t,



r,x) = cδ(x) , (13)

in addition to the boundary conditions that N satisfies. Solving for N amounts to constructing G,

then solving eq. (11). Intuitively, one can think of G as that key part of the diffusion model that 

embodies all space-time as well as x information since eq. (12) carries no space-time information, i.e., 

the standard notion of a Green’s function as an ‘influence’ function is clear here as well.

Constructing a general (or a class of) G, however, is not an easy task, and the carrying of 

space-time information will, to a large degree, be model dependent. As a prime example, in the literature,13,14 a steady-state model, i.e., ∂N/∂t=0, with cylindrical symmetry has been put forward, 

meaning a general G is written down with space information as well as other astrophysical input, e.g., 

relative thickness and density of the galactic disc to halo, etc. The model assumes that CR particles 

can freely escape the galactic volume (disc+halo) by crossing the surface S that encloses this volume, 

i.e., the assumed boundary condition at the surface is that Ni |

S=0 (as opposed to a reflecting boundary, ∂N/ ∂z |

S=0, where z is some space coordinate). Also, it is further assumed that Ni

 is continuous 

across the boundary of the disc and halo. While the complete formulation of this construction of G

is rather involved, for purposes here, the limit of the formulation in one-dimension is what is more 

instructive. It is helpful to think of that limit as when the diffusion process takes place mainly along 

a tube of magnetic field lines, and the only relevant space coordinate in this case will be the one along 

the field lines of the galactic magnetic field. Also, for purposes here, z (say z=0) can be pinned down 

and the essential information will still be intact. In this case, and for a constant diffusion coefficient, 

Dˆ , G takes the simple form: 

Gi(z = 0)  A

Dˆ 1+ ngvAσi / ⎡ Dˆ ⎣ ⎤

⎦

−1

. (14)

Above, A is an ‘area’ of the order of hg hh, disc (gas) and halo heights, si is the total reaction 

cross section for nuclide i, and ng is the disc’s matter density. The above approximation is valid for 

σi  Dˆ / (ngvhh

2 ).
10

6. TRYING OUT DIFFERENT G’S

One can construct a global model for which one derives the appropriate G, as in the 

above model of Ginzburg and Ptuskin,14 or one can try out various forms for G and pick the 

one that seems to fit the modeling purpose. It is instructive to go through some examples, 

however, as they will help us connect with the various notions, models, etc. that are widely used—

especially in data analyses.

6.1 The Slab Model

The simplest form one can pick is that of a Dirac delta function: 

G → δ x − x0 ( ). (15)

This particular choice of G carries no explicit space or time information. The only information it carries is that of a path length x=x0, i.e., species i is propagated through ISM matter equivalent to path 

length x0. For this G, eq. (11) becomes 

Ni t,

 ( r)→ Ni

s

x0 ( ) . (16)

Ni

s

x0 ( ), in turn, is the solution to eq. (12) evaluated at x=x0. Incidentally, eq. (12) has a general 

solution of the form:15

Ni

s

(x) = exp +Mx ( ) ⎡

⎣ ⎤⎦ij

Nj

s

(x = 0) ,

j

∑ (17)

where Mx is a species transformation matrix, i.e., nuclear chemistry information. Thus, with the 

choice (limit) of a Dirac delta function for G, the solution is reduced to the so-called ‘slab’ model of 

GCR propagation. Also obvious is that this limit suppresses all spatial information about the diffusion process, i.e., nothing is being diffused, only species transformed. Or, diffusion was so rapid, suggesting space homogeneity throughout the transit process. The limit, however, does serve the purpose 

of pointing out the complete separability of terms carrying the space-time information (diffusion), 

G, and those that are responsible for species transformation, Ni

s

(x) , in eq. (11). This feature, which 

comes about from the form of the G-S equation, is useful for studies of the GCR source composition as it allows one to ‘propagate back’ to the source, insofar as species transformation is concerned, 

simply by solving for Ni

s

(x = 0) given Ni

s

x0 ( ): 

Ni

s

(x = 0) = exp(−Mx ) ⎡

⎣ ⎤⎦ ij

Nj

s

x( 0 ) .

j

∑ (18)
11

But since the diffusion process is irreversible (initial conditions about space-time are lost in the process), propagating back cannot refer to this part of eq. (11). Fortunately, the separability feature 

alluded to above and the fact that the various nuclear information (cross sections) are not expected 

to be space-time dependent, one can, with some confidence, speak of propagating back to the GCR 

source.

6.2 The Venerable Leaky-Box Model

Starting here with the results, then the steps will be justified. The G in the leaky-box limit is 

that of a simple exponential: 

G(x) = exp −x / x ( esc ) , (19)

where xesc is the so-called escape path length related to the escape time (from the galaxy) of GCR by 

xesc = nvTesc , where n is average density. Back to eq. (11), with the above G, one has a solution for 

N in the so-called leaky-box limit. One recognizes that, instead of a simple slab model, i.e., a singular 

x, here is a distribution of x (path length distribution or PLD). This is termed the ‘weighted slab’ 

model due to the appearance of G in eq. (11) as a weighing function in x. Thus, the leaky-box limit 

(or model now) and the weighted slab model are mathematically equivalent in the absence of any 

energy losses, as assumed from the outset.15,16 In the presence of such losses, one can still treat them 

equally (at least for numerical purposes) as long as energy appears as a parameter in the equation to 

be solved.

Let us begin our justification by making the leaky-box ansatz that GCRs, for the most part, 

are confined to the galactic volume except for those GCRs that are endowed with enough kinetic 

energy to escape the confining effect of the galactic magnetic field (see fig. 1). Note also that the characteristic size of our galaxy is of the order of 100 kly (kilo light years). Since these particles travel 

at almost the speed of light, the characteristic transit time should then be of the order of thousands 

of years. Cosmic-ray data, however, as inferred mostly from the secondary to primary ratios which 

are directly affected by the amount of ISM matter traversed (of the order of 10 g/cm2) suggest that 

this timescale is actually of the order of 10 My. Hence, CRs travel in long-winded (diffused) paths 

between their source(s) and detection. Also note that CRs with energies above 1019 eV have a gyroradius in the ISM’s magnetic field (strength on the order of μG) larger than the galactic diameter 

of ≈100 kly; they can escape. Those that do escape are assumed to be able to do so at any point on 

surface S enclosing the galactic volume.

There is a subtle point here with regard to boundary conditions. In the diffusion model, the 

boundary condition is such that N drops to zero at the boundary point (Dirichlet boundary condition), any point on S. Now, in this leaky-box ansatz, a strong reflection, i.e., the gradient of N drops 

to zero at any point on S (Neumann boundary condition), is assumed for most of the CR particles 

except for those that can escape, of course. As such, i.e., different boundary conditions, the diffusion model and the leaky-box model cannot, strictly speaking, be considered related (even though a 

fundamental relationship in what follows will be established). What this means is that the leaky-box 

limit can be considered a truly mathematical limit (however physical or unphysical that limit may be) 

of the diffusion model at any point inside S; but not on S itself. (It is not clear what implications, if 
12

Figure 1. A schematic (not to scale) showing the long-winded path of a high-energy

cosmic-ray particle in the Milky Way as it makes its way out of the ‘box.’

any, this ‘technical’ point may have on conclusions about the size of the confining volume reached 

via the leaky-box model.) Recall that in addition to this, in the leaky-box model, one is also assuming 

spatial homogeneity.

Now, the first link that ought to be established between the two models is the simple exponential form of G (in the diffusion model), or equivalently, the simple exponential PLD characteristic 

of the (standard) leaky-box model. One arrives at the leaky-box equation by replacing the diffusion 

term with an escape term (the homogeneity assumption), 

−∇⋅(Dˆ∇N) → +

NT

Tesc

, (20)

and with no explicit time dependence, ∂N / ∂t = 0 , eq. (13) for N will read,

N

Tesc

+ nv

∂N

∂x = q , (21)

where both the density and source terms now denote average density and average source so as to 

be consistent with the homogeneity assumption. Equation (21) is the leaky-box equation without 

energy loss/gain or decay terms. Now, G must satisfy a similar equation: 

G

Tesc

+ nv

∂G

∂x = δ(x) , (22)

where, for simplicity, the strength of the delta function source term has been set to unity. Solving 

eq. (22) for G gives, 
13

G(x) = exp −x / x ( esc ) , (23)

which, after inserting it back into eq. (11), yields the solution of the G-S equation in the leaky-box 

limit. The escape path length xesc will be related to the diffusion coefficient as will be shown when the 

connection between the two models is established more firmly.

But, before that is done, the two underlying assumptions—steady state and homogeneity—

need some attention. Steady state is readily argued for on the basis of relevant GCR propagation timescales; i.e., relevant times are much longer than it takes a typical GCR particle to traverse the confining 

volume once, implying that a steady-state situation can be safely assumed. Homogeneity deserves 

a little algebra, and will be looked at from two different angles.

Species transformation aside, recall the separability feature of the G-S equation and examine 

the following term:

−∇⋅ Dˆ∇N(

 ( r)) = q , (24)

where now N is a function of space only due to the steady-state assumption. Now, for a homogeneous volume V ′ (enclosed by S′ rather than S because of the different boundary conditions), and 

after averaging over that volume, gives for the left-hand side of eq. (24): 

−∇⋅ Dˆ∇N(

 ( r)) = d3

r′ ∫ −∇⋅( ⎡ Dˆ∇N) ⎣ ⎤

⎦ /V ′

= ∫ dS′(−Dˆ∇N) ⋅ nˆ /V ′

= dS′ Jesc ( ) /V ′ , ∫ (25)

where nˆ is a unit vector ⊥ S′. Integrating over S′ the J’s vanish except for Jesc, which is the leakybox ansatz. Now, if Jesc = vescN and also assuming that where the particle escapes (on S′) is independent of S′, then

−∇⋅ Dˆ∇N(

 ( r))  vescNS′

S′r′

esc

= N

Tesc

. (26)

In eq. (26), the volume V ′ ∼ S′ r′

esc , where r′

esc is the characteristic escape distance and Tesc is 

the characteristic escape time. Averaging over V ′ on the right-hand side of eq. (24) gives q . Thus, 

eq. (21), after inserting back the species transformation term, is derived under the homogeneity 

assumption.
14

Also, for a constant diffusion coefficient, Dˆ , the diffusion term in the G-S equation under the 

leaky-box ansatz reduces to N/const where the constant (with the dimension of time) can be identified with Tesc. Start by writing 

−∇⋅(Dˆ∇N) |

esc= ∇⋅ N

vesc , (27)

which is valid for any point either on or inside the enclosing surface, be it S or S′ . The right-hand 

side of eq. (27) can be expanded to

∇⋅ Nvesc = N∇⋅



vesc + vesc ⋅ ∇N

= N

Dˆ ∇⋅



vesc − vesc

2 ( )

= N

Dˆ Φ , (28)

where ∇⋅



vesc − vesc

2 ( ) = Φ, which has the look of a ‘potential.’ Rewriting the diffusion equation in 

terms of 

v and at v

esc gives

∂



vesc

∂t = ∇ ∇⋅



vesc − vesc

2 ( )

= ∇Φ .

(29)

Now, if escape is not an explicit function of time, i.e., ∂



vesc / ∂t = 0, then Φ =const, which for a constant Dˆ gives, 

−∇i (Dˆ∇N) |

esc = N / const , (30)

which is again the leaky-box ansatz.

The difference between the two derivations, using S′ and Φ (even though they both lead to 

the same leaky-box limit) is that the question of boundary conditions does not come up when using 

Φ. What this means, though, is that the identification of the constant in eq. (30) as that related to 

some (constant) characteristic escape time seems valid, while not unique. This points again to the 

question of the different boundary conditions in the two models and why one should be somewhat 

careful when using conclusions of one model to infer or support assumptions about another.

6.3 Nonstandard Leaky-Box G

As seen above, the Green’s function for the standard leaky-box limit is that of a simple exponential, eq. (19). Variants17–19 of this simple function have also been used and found to fit some of 

the observed CR data. These include the double exponential, G=G1–G2, where both G1 and G2 have 
15

the form of eq. (19), but with different xesc. This is the so-called nested leaky-box model and the 

physics behind it is that of a ‘shrouded’ CR source. Also, there is the truncated exponential PLD: 

G(x) = 0,

exp −(x − δ / xesc ⎡

⎣ ⎤⎦

⎧

⎨

⎪

⎩

⎪

x < δ

x ≥ δ , (31)

where one assumes a significant near-source grammage is encountered by CRs, or, that of a CR 

source that is farther away than in the standard leaky-box limit.

The point to appreciate here is that, mathematically, if the simple exponential appears to be 

reasonable, data-wise that is, there is every reason to suspect that small departures and/or generalizations in the form of G from that of a simple exponential should still be reasonable. Or, one could even 

ask if there is a general form for G (of which the simple exponential, and variants, are but special 

cases) that is consistent with the leaky-box ansatz. One is tempted to answer this question by going 

back to the simple one-dimensional diffusion equation and to begin relaxing the main assumptions 

that lead to the leaky-box limit as a general solution is rederived.
16

7. THE CONTACT POINT: xesc ⇔ Dˆ

To recite the main conclusion of section 6 (with the question of boundary conditions aside): 

The leaky-box ansatz is the mathematical limit of the one-dimensional diffusion model for a homogeneous confining volume (with the ‘right’ initial and boundary conditions) and implicit time dependence. Naturally, one needs to identify any points of contact between the two models here. The 

diffusion model, in one or more dimensions, is characterized by a diffusion coefficient (or tensor) 

Dˆ. The leaky-box limit is characterized by the escape time Tesc or the escape path length xesc. Clearly 

then, a connection between the two models is manifested in a connection between Dˆ and xesc . Examining eq. (13), one-dimensional diffusion, and equation (21), leaky-box limit, one can easily show14

that 

xesc 

ngvhhhg

Dˆ , (32)

which is valid as long as hg/hh  1, i.e., tube-like diffusion along the galactic magnetic field line of 

force (recall that hg is a characteristic height scale of the disc and hh for the halo).

Equation (32) is also somewhat obvious from the second derivation of the leaky-box limit, 

i.e., the one using the potential-like function Φ. One can see from eqs. (29) and (30) that the constant 

must somehow be related to (a constant) Dˆ. Even though it appears absent in the literature, one 

should be able to arrive at a relation like eq. (32) using this derivation and invoking (chaotic) adiabatic invariance, in some simple model of (spatially) trapped charged particles.19

The significance of the contact point is that one can inject an energy dependence in xesc, which 

appears needed by the data, via any energy dependence in Dˆ. From a Kolmogorov-type description (statistical description) of a fast charged particle’s motion in turbulent fields,7,8,20,21 the dependence of Dˆ on energy has the form 

Dˆ(E) ∝ E−κ , (33)

where κ is a typical Kolmogorov decay exponent, characteristic of diffusive transport processes. In 

terms of rigidity R (total momentum per charge), 

Dˆ(R) ∝ βRα , (34)

with β=v/c and the exponent a=1/3. With this dependence of Dˆ on R, xesc is written 

xesc = x0β R0 ( / R)

α′

, (35)

with R0 being some rigidity parameter and x0 a constant that appears to depend on the assumed 

composition of the ISM.22 The exponent α′ seems to depend on the particular species propagated 

and whether it happens to be a primary or a secondary CR nucleus.23 
17

8. PUTTING THE BOX TOGETHER

In this section and section 9, how one goes about implementing the leaky-box model is briefly 

sketched. From the outset, the evolution of N as a particle density distribution function has been 

discussed. What typically one measures though is the differential particle flux j (number of particles 

observed per unit area (dA), per solid angle (dΩ), per energy interval (dE), and per unit time (dt)). 

However, the measurement of j is equivalent to a measurement of N. What makes this true is the 

simple fact that the number of particles is a conserved quantity, i.e., 

j dE dAdΩdt = N p

2 dpv dAdΩdt , (36)

where p is the nonrelativistic momentum of the particle with velocity v and mass m. Since pdp=mdE, 

eq. (36) simplifies to 

j=p2N . (37)

The same can be shown for relativistic particles.

First, reinject any energy loss/gain terms in the model for realistic numerical implementation 

of the model. Energy losses here primarily refer to those due to the ionizing effects of the charged 

CR particles as they pass through the ISM matter. Such losses are calculated using the well-known 

Bethe-Bloch expression for the stopping power W(E). Qualitatively, and at low energies, 

W(E) = − dE

dx ∝ 1

v

2 . (38)

At high energies, W(E) scales as lng, where g is the Lorentz factor. For ultrarelativistic particles in 

dense media, the situation is slightly more complicated, where a relativistic correction has to be 

carefully injected in the expression of W(E).24 While W(E) is quite sensitive to the particle’s kinetic 

energy, it has no explicit dependence on the particle’s (rest) mass. This implies that particles with different mass but same momentum will have similar rates of energy loss. Also, W(E) has a minimum, 

below this minimum, and because of the dependence depicted by eq. (38), the energy loss increases 

(quite rapidly) as v decreases. Above this minimum, however, the loss rate increases quite slowly. To 

a good approximation and up to two decades of E above the minimum, the loss is more or less the 

same and is of the order of a few MeV per g/cm2 (for 1 GeV protons, for example) of matter, independent of the exact characteristics of the matter.

Energy losses due to nuclear interactions—those that are responsible for species transformation—are typically neglected; the so-called straight-ahead approximation. This is a rather good 

approximation for energies above a few GeV/nucleon. Relaxing the straight-ahead approximation 

suggests25 that the approximation tends to introduce an error on the order of only 5%.
18

The energy loss term enters in eq. (21), the leaky-box equation, looking like 

+

∂

∂E

Wi(E)Ni ⎡ (E) ⎣ ⎤

⎦ = qi . (39)

Energy gain terms, e.g., diffusive acceleration,16 stochastic acceleration,18 or reacceleration,26 enter 

eq. (21) in a similar manner. 

Finally, decay terms, loss as well as gain of species due to radioactive decay, also need to be 

included in eq. (21) as they tend to alter the balance of Ni given qi. Even though such terms are 

functions of time, and there is no explicit time dependence in eq. (21), they are taken as functions of 

energy only and a cutoff time is typically assumed, i.e., the characteristic age of CRs, on the order 

of few million years. Meaning, radioactive species that decay with half-lives longer that this time are 

considered stable here. Now, one can inject this sort of information in eq. (21) via two ways. The first 

is to write, in eq. (21), the decay term explicitly as a function of path length x as: 

dNi

dx decay∝ − Niτ i

+

ρijNj

τ j j

∑ ⎛

⎝

⎜

⎜

⎞

⎠

⎟

⎟ , (40)

where ti is the half-life of species i and rij is the branching ratio (fraction of decays of species j yielding species i). Another way is to inject the same information about half-lives and branching ratios 

in the cross-section terms, i.e., in si and sij appearing in eq. (2). If no explicit time dependence is 

assumed, the two ways are equivalent.

Equation (21), in terms of the various escape, interaction, and decay path lengths, along with 

ionization energy losses, is now written as 

 qi(E) − Ni(E)

xi

esc (E)

− Ni(E)

xi

int(E)

+

Nj (E)

xji

int(E) j>i

∑ − Ni(E)

xi

dec (E)

+

N j(E)

xji

dec (E) j>i

∑ − ∂

∂E

Wi(E)Ni ⎡ (E) ⎣ ⎤

⎦ = 0 , (41)

where the plus or minus sign of each term is obvious as the process (described by that term) either 

enhances the source term qi (+) or depletes it (–). Equation (41) is the standard leaky-box propagation equation.

Solving eq. (41)—for Ni at a given energy, i.e., energy is treated as a parameter—proceeds in 

two equivalent ways. One can either solve for Ni

s

(x,E) in 

∂Ni

s

(x,E)

∂x = left-handside of eq. (41) , (42)
19

with the source term understood as Ni (x=0,E), and then integrates 

Ni(E) = dxNi

s

0

x

∫ (x,E) . (43)

Or, one solves eq. (41) for Ni (x,E) without the escape loss term, then integrates 

Ni(E) = dx Ni

s

0

∞

∫ (x,E)P(x,E) , (44)

withP(x,E) being the normalized PLD function, or G earlier, e.g., eq. (19). The source term, 

Ni(x=0,E) or qi(E), is typically taken to be a power law in rigidity for each species either due to 

a single shock27,28 or multiple shocks29,30 but with an added flattening term (power law in velocity) 

for low energies.
20

9. SAMPLE NUMERICAL IMPLEMENTATION

In this simple—but hopefully illustrative—example, the leaky-box model will be implemented for 

a two-species CR source made up of only protons (90%) and alpha particles (10%). At the ‘point source,’ 

these are synthesized and energized with a ‘prepropagation and spallation’ spectrum determined solely 

by the particle’s rigidity. In the ISM (made up of protons only), both the protons and the alpha particles 

suffer some energy losses but only the alphas suffer spallation. Neither suffers any decay. A simple formula will be used to account for modulation before comparing—qualtitatively only—the leaky-boxgenerated spectra to observed spectra at 1 AU.

The two-species system of equations to be solved is then, 

∂N1

∂x = q1(E) − N1(E)

x1

esc (E)

− N1(E)

x1

int(E)

+

N2 (E)

x2→1(E)

− ∂

∂E

W1(E)N1 ⎡

⎣ ⎤⎦ (45)

and

∂N2

∂x = q2 (E) − N2 (E)

x2

esc (E)

− N2(E)

x2

int

(E)

+

N1(E)

x1→2 (E)

− ∂

∂E

W2 (E)N2 ⎡

⎣ ⎤⎦ , (46)

where ‘1’ refers to alpha particles and ‘2’ to protons. If the energy loss function W does not depend 

on x explicitly, which it usually does not, then the structure of this two-species system suggests the 

method of characteristics for solution. Multiply the equation for N1 by W1 so that it can be rewritten 

as: 

∂G1

∂x

+W1

∂G1

∂E = W1 q1 − λ1

−1 ( G1) , (47)

where G1(x,E) is defined as W1(E)N1(x,E) and where, for the alpha particles, 

1

λ1

= 1

x1

esc

+

1

x1

int , (48)

and x2→1→∞. Using the standard method-of-characteristics solution for this linear, first-order PDE, 

gives

N1(x,E) = !1(E)

W1(E)

W1(E)q1(E) " exp "x / !1 ( (E)) W1(E#)q1(E#) 1" n1!1

"1 ( (E#)) $

% & { '} , (49)
21

where E′ = E + xW1(E) . Note that the assumption is used that limx→0 Ni

(x,E) = niqi(E),i=1,2, 

where ni is a constant that carries the relative abundance factor (i.e., 0.9 for protons and 0.1 for 

alphas) as well as the right dimensional conversion, without any loss of generality, the latter can be 

taken to be unity. For N2(x,E), if its source term is rewritten as 

( , ) = ( ) ( , ) , 1

1 Q2 x E q2 E x1 2 N x E − + → (50)

its PDE will resemble that for N1(x,E). With G2 defined, as before, as W2N2,

∂G2

∂x

+W2

∂G2

∂E = W2Q2 − λ2

−1 ( G2 ) , (51)

with a solution, 

N x E

E

W E

W E Q x E x E 2

2

2

2 2 2 ( , ) ( )

( ) = − ( ) ( , ) exp( ) − / ( ) 





λ λ

× W2 (E′)q2(E′) 1− λ2

−1

(E′) + λ1

−1

(E′)

q1(E′)

q2(E′)

⎛

⎝

⎜

⎞

⎠

⎟

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

⎫

⎬

⎪

⎭⎪

,

(52)

where, for the protons, E′ = E + xW2 (E), and 

1

λ2

= 1

x2

esc

+

1

x2

int .

(53)

Next, the various source, escape, energy loss, and interaction terms appearing in the above

solutions for N1(x,E) and N2(x,E) are addressed. The source function is typically parameterized31

using a power law in rigidity as (see fig. 2): 

( ) = , 2.5

0

− q R q R (54)

with R = E + 2m0 ( )E A/Z, m0 being the nucleon rest mass and A and Z are the mass and charge

number of the particles, respectively.

The escape function is also typically parameterized32 using a power law in rigidity as well. If 

R is expressed in units of GV, the escape path length in units of g/cm2 can be well approximated as 

(see fig. 3): 

xesc (R) = 26.7β

(βR)

0.5 + (βR / 1.4)−1.4 ⎡

⎣ ⎤⎦

. (55)
22

10–6

10–8

10–10

10–12

10–14

102 103 104 105

Kinetic Energy (MeV/Nucleon)

Number Density (arb. units)

Protons

Alphas

Figure 2. The source function, eq. (54), for protons and alphas as used 

in the sample illustration.

14

12

10

8

6

4

2

0

102

 103 104 105

Protons

Alphas

Kinetic Energy (MeV/Nucleon)

esc (g/cm2

) χ

Figure 3. The escape function, eq. (55), for protons and alphas as used

in the illustration.

In units of MeV per g/cm2, the energy loss function W(E), i.e., the Bethe-Block relation, can 

be written33 as 

W(E) = Kz

2 Z

A

1

β2

⎛

⎝

⎜

⎞

⎠

⎟

1

2

ln

2mec

2

β2γ 2Tmax

I 2 − β2 − δ

2

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

, (56)

where K=0.307075 g/cm2 for A=1 g/mol, z is the charge of the traversing particle, Z the atomic 

number of the medium (1, in this case, for an ISM of protons only), A the atomic mass of the 
23

medium (also 1), mec2 is the mass of the electron, b and g are the standard Lorentz factors for the 

traversing particle, Tmax is the maximum kinetic energy that can be imparted to a free electron in a 

single collision, I is the mean excitation energy, and d is a density effect term. For a traversing particle 

with mass M, 

Tmax = 2mec

2

β2γ 2

1+ 2γ me / M + me

2 / M2 . (57)

The density effect term, δ, can be ignored for such low density media like the ISM (see fig. 4).

30

25

20

15

10

5

0

101

 102 103 104 105

Kinetic Energy (MeV)

Energy Loss Rate (MeV/(g/cm2))

Figure 4. Energy loss by protons traversing a hydrogen gas—Bethe-Block relation,

eq. (56). Discrete points are as simulated by the code SRIM/TRIM.34

For the protons-and-alphas-only leaky-box model, there are three nuclear interaction terms 

that need to be estimated. The first two are for the total reaction cross sections for CR alphas and 

protons impacting the ISM protons as a function of energy; these are needed to infer the interaction mean free paths x1

int

(E) and x2

int

(E) . These cross sections, for kinetic energies in the range 

102 <E<105 MeV/nucleon, are easily parameterized2 using a simple function like: 

σ total(E) = a ⎡1+ b exp(c lnE)cos(d lnE + e) ⎣ ⎤⎦ , (58)

where {a,b,c,d,e} are constant parameters. For alphas, a=100, b=16, c=–1, d=–2.5, and e=–2, for 

σ(E) in units of mb and E in MeV/nucleon. For protons, a=34, b=122, c=–1, d=1.5, and e=1 (see 

fig. 5). The third term is the cross section for production of protons from spallation collisions of CR 

alphas with the ISM protons, to estimate x1→2. This is also relatively easy to parameterize in the 

same energy range as follows (see fig. 6): 

σ prod (E) = a[1+ b exp(c lnE)] , (59)
24

using a=93, b=–37, and c=–1. These approximated cross sections (in mb) (typically as a function 

of energy) are then converted to mean free paths (in g/cm2) using 

x or λ [in g/cm

2 ] = molar mass of the ISM[in g per mole]

Avogadro's number [in No. per mole] × σ [in cm

2 ( ]) ≈ 104

(3σ [in mb]) . (60)

140

120

100

80

60

40

20

0

102

 103 104 105

Kinetic Energy (MeV/Nucleon)

Total Reaction Cross Section (mb)

He + p

p + p

Figure 5. The total reaction cross sections, eq. (58), for protons and alphas

as used in the illustration.

140

120

100

80

60

40

20

0

102

 103 104 105

Kinetic Energy (MeV/Nucleon)

Production Cross Section (mb)

He + p→p + X

Figure 6. The spallation cross section, eq. (59), for alphas on ISM protons 

to produce secondary CR protons as used in this illustration.
25

Next, the solutions need to be averaged, which are now functions of both distance or path 

length (in g/cm2) and energy using an assumed PLD. In the leaky-box approximation, this distribution, eq. (19), when normalized becomes

G(x) = xesc

−1

exp −x / xesc ( ) , (61)

where xesc is the mean escape path length (typically <10 g/cm2). Hence, the arriving solutions (but 

still unmodulated yet by the solar wind) become (i=1,2),

Ni(E) = Ni(x,E)xesc

−1 exp −x / x ( esc ) 0

∞

∫ dx . (62)

Solar modulation, for our purposes, can be illustrated using a potential35 and an energy loss as, 

N

i(E) = E2 + 2m0E

(E + ΔE)

2 + 2m0(E + ΔE)

Ni(E + ΔE) , (63)

where ∆E=ZeΦ /A, Φ being the potential (in V; typically about 400 MV for solar minimum conditions and about 800 or higher for solar maximum conditions) and e the electron charge.

The ‘arriving’ number densities at 1 AU for protons and alpha, N

2 (E) and N

1(E), are 

shown in figure 7, where the averaging over the path length, eq. (62), has been done numerically. 

Here, no atmospheric effects are taken into account. In figure 8, the spallation of alphas’ contribution to the arriving protons are highlighted. This easily is calculated using 100 × N′

2 (E) / N2 ( (E)),

where N′

2 (E) is a solution to the proton equation, eq. (51), but with the spallation term ‘turned off’ 

by letting x1→ 2→ ∞. Figure 9 shows the ratio of the arriving alphas to the arriving protons, i.e., 

N

1(E) and N

2 (E) , which includes the secondary proton component.

10–7

10–9

10–11

10–13

10–15

102 103 104 105

Kinetic Energy (MeV/Nucleon)

Number Density (arb. units)

Protons

Alphas

Figure 7. The arriving number densities for protons and alphas, eq. (63), as predicted 

by this sample illustration using xesc=7 g/cm2 and Φ = 400 MV.
26

1

0.8

0.6

0.4

0.2

0

102

 103 104 105

Kinetic Energy (MeV/Nucleon)

Percent

Figure 8. Ratio of the secondary protons (produced by the spallation of the primary CR

alphas in the ISM) to the primary CR protons.

10

8

6

4

2

0

102

 103 104 105

Kinetic Energy (MeV/Nucleon)

Percent

Figure 9. Ratio of arriving alphas to protons (including secondary protons).

To put these sample calculations in perspective, figure 10 shows actual observed BESS36 data 

for the flux of CR protons and alphas; their ratios are shown in figure 11. BESS balloon-borne spectrometer was flown from Manitoba, Canada, in July 1998 (i.e., under solar minimum conditions) at 

a height of ≈37 km, where residual atmospheric effects are equivalent to <5 g/cm2. For this comparison, 

the measured BESS flux has been multiplied by a constant (=2×10–8) so as to scale it to the calculated 

number density. The results from this simple illustration are not meant to be a model prediction designed 

to be compared directly to data. Only qualitative comparisons are intended and useful here. Already, 

a number of general observations can be made: 
27

10–7

10–9

10–11

10–13

10–15

102 103 104 105

Kinetic Energy (MeV/Nucleon)

Flux×2×10–8(1/(

m2–sr–s–MeV/Nucleon)

BESS-Observed Protons

Alphas

Figure 10. Observed proton and alpha flux by BESS.36

10

8

6

4

2

0

102

 103 104 105

Kinetic Energy (MeV/Nucleon)

Percent

Figure 11. Ratio of alphas to protons as observed by BESS.36

(1) The general shape of the protons’ and alphas’ arriving spectra appear straightforward to 

reproduce assuming simple rigidity dependent escape and source functions. 

(2) Alphas do not appear to be a major contributor to secondary protons.

(3) Spallation of CR nuclei heavier than alphas must contribute to both protons and alphas.

(4) Low energy (less than a few GeV’s per nucleon) behavior of GCRs is more difficult to 

model than at higher energies, where solar modulation effects and nuclear (total reaction and spallation) cross sections appear to subside.
28

(5) Propagation effects (nuclear chemistry of CRs) appear to be well captured by simple models like the leaky-box model mainly because the energy dependence of the relevant nuclear parameters is rather weak as well as being well understood and modeled.
29

10. REMARKS

In ending this brief review focusing on the practical side of modeling the galactic propagation 

of cosmic rays, the following are three general closing remarks:

(1) The intended reader is a newcomer to the field as well as that practitioner who simply 

wants to delve a bit deeper into this immensely interesting area. At various points in this review, I 

suspect both wondered why I had focused and picked areas to elaborate on and paid only scant attention to others. I have no strong general justification to this (accurate) observation except to say the 

following. I had ventured my own study into cosmic rays from a nuclear physics background (and 

predisposition). Learning the eclectic new subject, I had to build my knowledge base and formulate 

my own understanding from what seemed to me like a mosaic of specialities, styles, backgrounds, 

rigors, and aims. It was ‘natural’ for a nuclear physicist to gravitate towards the nuclear aspects of 

cosmic rays. In retrospect, I believe, it was equally ‘natural’ to also look for something ‘different’ to 

do. Galactic propagation aspects of cosmic rays are related (cf. nuclear reactions) yet are somewhat 

different (cf. kinetic and transport formulations).

(2) I was attracted by the Ginzburg-Syrovatskii elegant formulation of the problem from 

early on. I had found their derivation and exposition of the kinetic theory of this problem (‘the 

quantitive theory of the origin of cosmic rays,’ as they had auspiciously termed it) an area I wanted 

to get into. Their theory survives to this day. I strongly recommend that the interested reader pick 

up and study their classic book. The limit of the Ginzburg-Syrovatskii propagation theory, the 

leaky-box model, is now a standard tool for both cosmic-ray data analysis and modeling as well as 

a more or less complete model of the nuclear chemistry of cosmic rays. In other words, the leaky-box 

model, albeit a limit to or a crude approximation of ‘the quantitative theory of the origin of cosmic 

rays’ is essentially a complete framework for the nuclear chemistry of cosmic rays. 

(3) Detailed modeling—with an eye on data comparison and such—of the nuclear chemistry 

of cosmic rays using the leaky-box framework is actually quite straightforward as attempted to be 

demonstrated. ‘Professional’ versions of this leaky-box demonstration (e.g., see ref. 37) are fundamentally similar to this simple one. 

The first 100 years since the discovery of cosmic rays, and especially since the coming of 

the space age, has seen unimaginable progress, refinement, and evolution of understanding of this 

intriguing phenomenon. It is quite likely, however, that the nuclear chemistry aspects of cosmic rays 

may have matured to a point where similar and equally fantastic progress cannot be foretold—in yet 

another sigmoidal display of change in nature.
30
31

REFERENCES

 1. Ginzburg, V.L.; and Syrovatskii, S.I.: The Origin of Cosmic Rays, Pergamon Press, Oxford, 1964.

 2. Silberberg, R.; and Tsao, C.H.: “Spallation Processes and Nuclear Interaction Products of Cosmic 

Rays,” Phys. Rep., Vol. 191, pp. 351–408, 1990.

 3. Chandrasekhar, S.: “Stochastic Problems in Physics and Astronomy,” Rev. Mod. Phys., Vol. 15, 

pp. 1–81, 1943.

 4. Montroll, E.W.; and West, B.J.: Fluctuation Phenomena, Ch. 2, E.W. Montroll and J.L. Lebowitz 

(eds.), North Holland, Amsterdam, 1987.

 5. Gardiner, C.W.: Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences, 

Ch. 4, Springer, Berlin, Germany, 1983.

 6. Lefebvre, M.: Applied Stochastic Processes, Ch. 5, Springer, Berlin, Germany, 2000.

 7. Jokipii, J.R.: “Propagation of Cosmic Rays in the Solar Wind,” Rev. Geophys. Space Phys., Vol. 7, 

pp. 27–87, 1971.

 8. Dröge, W.: “Transfer of Solar Energetic Particles,” Astrophys. J. Supp., Vol. 90, pp. 567–576, 1994.

 9. Lee, M.A.: “Self-consistent Kinetic Equations and the Evolution of a Relativistic Plasma in an 

Ambient Magnetic Field,” Plasma Phys., Vol. 13, pp. 1079–1098, 1971.

10. Swanson, D.G.: Plasma Waves, Ch. 7, Academic Press, Boston, 1989.

11. Fetter, A.L.; and Walecka, J.D.: Theoretical Mechanics of Particles and Continua, Ch. 7, McGraw 

Hill, New York, 1980.

12. Butkov, E.: Mathematical Physics, Ch. 12, Addison Wesley, Reading, Massachusetts, 1968.

13. Ptuskin, V.S.: Astrophysics of Cosmic Rays, Ch. 3, V.L. Ginzburg (ed.), North Holland, Amsterdam, 

1990.

14. Ginzburg, V.L.; and Ptuskin, V.S.: “The Origin of Cosmic Rays,” Sov. Sci. Rev. E Astrophys. Space 

Phys., Vol. 4, pp. 161–256, 1985.

15. Letaw, J.R.; Silberberg, R.; and Tsao, C.H.: “Propagation of Heavy Cosmic-ray Nuclei,” Astrophys. 

J. Supp., Vol. 56, pp. 369–391, 1984.
32

16. Heinbach, U.; and Simon, M.: “Propagation of Galactic Cosmic Rays Under Diffusive Reacceleration,” 

Astrophys. J., Vol. 441, pp. 209–221, 1995.

17. Wefel, J.P.: in Genesis and Propagation of Cosmic Rays, M.M. Shapiro and J.P. Wefel (eds.), Reidel, 

Dordrecht, pp. 1–40, 1988.

18.	 Jones, F.R.; Lukasiak, A.; Ptuskin, V.; and Webber, W.: “The Modified Weighted Slab Technique:

Models and Results,” Astrophys. J., Vol. 547, pp. 264–271, 2001.

19. Brown, R.; Ott, E.; and Grebogi, C.: “The Goodness of Ergodic Adiabatic Invariants,” J. Stat. Phys., 

Vol. 49, pp. 511–550, 1987.

20. Volk, H.J.: “Cosmic Ray Propagation in Interplanetary Space,” Rev. Geophys. Space Phys., Vol. 13, 

pp. 547–566, 1975.

21. Palmer, I.D.: “Transport Coefficients of Low-energy Cosmic Rays in Interplanetary Space,” Rev. 

Geophys. Space Phys., Vol. 20, pp. 335–351, 1982.

22. Ferrendo, P.: “Chemical and Isotopic Composition of Cosmic Rays,” J. Phys. G: Nucl. Part. Phys., 

Vol. 19, pp. S53–S64, 1993. 

23. Gaisser, T.K.: Cosmic Rays and Particle Physics, Ch. 9, Cambridge University Press, Cambridge, 

1990.

24. Jackson, J.D.: Classical Electrodynamics, 2nd Ed., Ch. 13, Wiley, New York, 1975.

25. Tsao, C.H.; Silberberg, R.; Barghouty, A.F.; and Sihver, L.:”Energy Degradation in Cosmic-Ray 

Nuclear Spallation Reactions: Relaxing the Straight-ahead Approximation,” Astrophys. J., Vol. 451, 

pp. 275–283, 1995.

26. Letaw, R.; Silberberg, R.; and Tsao, C.H.: “Comparison of Distributed Reacceleration and Leaky-box 

Models of Cosmic-ray Abundances (Z = 3-28),” Astrophys. J., Vol. 414, pp. 601–611, 1993.

27. Blandford, R.D.; and Ostriker, J.P.: “Supernova Shock Acceleration of Cosmic Rays in the Galaxy,” 

Astrophys. J., Vol. 237, pp. 793–808, 1980.

28. Eichler, D.: “Basic Inconsistencies in Models of Interstellar Cosmic-Ray Acceleration,” Astrophys. J., 

Vol. 237, pp. 809–813, 1980.

29. Silberberg, R.; Tsao, C.H.; Letaw, J.R.; and Shapiro, M.M.: “Distributed Acceleration of Cosmic 

Rays,” Phys. Rev. Lett., Vol. 51, pp. 1217–1220, 1983.

30. Simon, M.; Heinrich, W.; and Mathis, K.D.: “Propagation of Injected Cosmic Rays Under Distributed 

Reacceleration,” Astrophys. J., Vol. 300, pp. 32–40, 1986.
33

31. Ormes, J.F.; and Protheroe, R.J.: “Implications of HEAO-3 Data for the Acceleration and Propagation 

of Galactic Cosmic Rays,”Astrophys. J., Vol. 272 , pp. 756–764, 1983.

32. Yanasak, N.E.; Wiedenbeck, M.E.; Mewaldt, R.A.; et al.: “Measurement of the Secondary 

Radionuclides 10Be, 26Al, 36Cl, 54Mn, and 14C and Implications for the Galactic Cosmic-Ray 

Age,” Astrophys. J., Vol. 563 , pp. 768–792, 2001.

33. Eidelman, S.; Hayes, K.G.; Olive, K.A.; et al.: “Review of Particle Physics (Particle Data Group),” 

Phys. Lett., Vol. B 592, pp. 1–5, 2004.

34. Ziegler, J.F.; Biersackm, J.P.; and Littmark, U.: The Stopping Power and Range of Ions in Solids, 

Pergamon Press, New York, <http://www.srim/org>, 1996.

35. Gleeson, L.I.; and Axford, W.I.: “Solar Modulation of Galactic Cosmic Rays,” Astrophys. J., Vol. 154, 

pp. 1011–1026, 1968.

36. Sanuki, T.; Motoki, M.; Matsumoto, H.; et al.: “Precise Measurement of Cosmic-Ray Proton and 

Helium Spectra with the BESS Spectrometer,” Astrophys. J., Vol. 545, pp. 1135–1142, 2000.

37. “The GALPROP Code for Cosmic-Ray Transport and Diffuse Emission Production,” <http://galprop.

stanford.edu/>, Accessed April 2012.
34

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining 

the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing 

this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operation and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid 

OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S) 5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

 REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITOR’S ACRONYM(S)

11. SPONSORING/MONITORING REPORT NUMBER

12. DISTRIBUTION/AVAILABILITY STATEMENT

13. SUPPLEMENTARY NOTES

14. ABSTRACT

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF:

a. REPORT b. ABSTRACT c. THIS PAGE

17. LIMITATION OF ABSTRACT 18. NUMBER OF

 PAGES

19a. NAME OF RESPONSIBLE PERSON

19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18

A Review of the Ginzburg-Syrovatskii’s Galactic Cosmic-Ray 

Propagation Model and its Leaky-Box Limit

A.F. Barghouty

George C. Marshall Space Flight Center

Huntsville, AL 35812

National Aeronautics and Space Administration

Washington, DC 20546–0001

Unclassified-Unlimited

Subject Category 90

Availability: NASA CASI (443–757–5802)

Prepared for the Astrophysics Office, Science and Mission Directorate

M–1340

Technical Publication

NASA/TP—2012–217464

cosmic rays, interstellar medium, diffusion, nuclear reactions

01–07–2012

UU 48

NASA

U U U

Phenomenological models of galactic cosmic-ray propagation are based on a diffusion equation known as 

the Ginzburg-Syrovatskii’s equation, or variants (or limits) of this equation. Its one-dimensional limit in 

a homogeneous volume, known as the leaky-box limit or model, is sketched here. The justification, utility, 

limitations, and a typical numerical implementation of the leaky-box model are examined in some detail.

STI Help Desk at email: help@sti.nasa.gov

STI Help Desk at: 443–757–5802
The NASA STI Program...in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

Program Office plays a key part in helping NASA 

maintain this important role.

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for 

NASA’s scientific and technical information. The 

NASA STI Program Office provides access to 

the NASA STI Database, the largest collection of 

aeronautical and space science STI in the world. 

The Program Office is also NASA’s institutional 

mechanism for disseminating the results of its 

research and development activities. These results 

are published by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA programs and include extensive data 

or theoretical analysis. Includes compilations 

of significant scientific and technical data 

and information deemed to be of continuing 

reference value. NASA’s counterpart of peerreviewed formal professional papers but has less 

stringent limitations on manuscript length and 

extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical conferences, 

symposia, seminars, or other meetings sponsored 

or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, technical, 

or historical information from NASA programs, 

projects, and mission, often concerned with 

subjects having substantial public interest.

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results...even 

providing videos.

For more information about the NASA STI Program 

Office, see the following:

• Access the NASA STI program home page at 

<http://www.sti.nasa.gov>

• E-mail your question via the Internet to 

<help@sti.nasa.gov>

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Phone the NASA STI Help Desk at 

443–757–5802

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076–1320
National Aeronautics and

Space Administration

IS20

George C. Marshall Space Flight Center

Huntsville, Alabama 35812