TF SPIM StoppingPower

Stopping Power

Bethe Equation

Classical Energy Loss

Consider the energy lost when a particle of charge ([math]ze[/math]) traveling at speed [math]v[/math] is scattered by a target of charge ([math]Ze[/math]). Assume only the coulomb force causes the particle to scatter from the target as shown below.

Notice

as [math]ze[/math] is scattered the horizontal component of the coulomb force ([math]F[/math]) flips direction; ie net horizontal force for the scattering

[math]F_{vertical} = k \frac{zZe^2}{r^2} \sin(\theta) = k \frac{zZe^2}{r^2} \frac{b}{r}[/math]

where

k =[math]\frac{1}{4 \pi \epsilon_0}[/math]

r = distance between incident projectile and target atom

b= impact parameter of collision

Using the definition of Impulse one can determine the momentum change of [math]ze[/math] as

[math]\Delta p = \int F dt[/math]

Let's assume that the energy lost by the incident particle [math]ze[/math] is absorbed by an electron in the target atom. This energy may be cast in terms of the incident particles momentum change as

[math]\frac{(\Delta p)^2}{2m_e}[/math]

By calculating the change in momentum ([math]\Delta p[/math]) of the incident particle we can infer that the energy lost by the incident particle is absorbed by one of the target material's atomic electrons.

[math]\Delta P = \int F dt = \int k \frac{zZe^2b}{r^3} dt[/math]

using [math]dt = \frac{dx}{v} = \frac{d x}{\beta c}[/math] we have

[math]= k \frac{zZe^2b}{\beta c} \int_{-\infty}^{+\infty} \frac{ dx}{(x^2+b^2)^{3/2}}[/math]

[math]=\frac{kzZe^2b}{\beta c b^2} \int_{-\infty}^{+\infty} \frac{ dx/b}{(1+\frac{x^2}{b^2})^{3/2}}[/math]

[math]\int_{-\infty}^{+\infty} \frac{ dx/b}{(1+\frac{x}{b^2})^{3/2}}=2[/math]

[math] \Delta p = \frac{2kzZe^2b}{\beta c b^2}[/math]

casting this in terms of the classical atomic electron radius [math]r_e[/math]

[math]r_e = \frac{k e^2}{m_e v^2} \sim \frac{k e^2}{m_e c^2}[/math] just equate [math]F = \frac{ke^2}{r_e^2} = m \frac{v^2}{r_e}[/math]

Then

[math] \Delta p = \frac{2zZr_e m_e c}{\beta b}[/math]

and

[math]\Delta E = \frac{(\Delta p)^2}{2m_e} = 2 \left ( \frac{r_e m_e}{\beta b}\right )^2 \frac {z^2 Z^2 c^2}{m_e}[/math] : [math]Z[/math] = 1 here because I shall assume the energy is lost to just the electron and the Atom is a spectator

Now let's calculate an expression representing the AVERAGE energy lost for an incident particle traversing a material of some thickness.

Let

[math]P(\Delta E)[/math] = Probability of an interaction taking place which results in an energy loss [math]\Delta E[/math]

If we let

Z = Atomic Number = # electrons in target Atom = number of protons in an Atom

N = Avagadros number = [math]6.022 \times 10^{23} \frac{Atoms}{mol}[/math]

A = Atomic mass = [math]\frac{g}{mole}[/math]

[math]dP(\Delta E)[/math] = probability of hitting an atomic electron in the area of an annulus of radius ([math]b + db[/math]) with an energy transfer between [math]\Delta E[/math] and [math]\Delta E + d(\Delta E)[/math]

Then

[math]\frac{-dE }{dx}= \int_0^{\infty} dP(\Delta E) \Delta E[/math] = energy lost by the incident particle per distance traversed through the material

I am just adding up all the energy losses weighted by the probability of the energy loss to find the average (total) energy loss.

What is [math]dP(\Delta E)[/math]

[math]dP(\Delta E)[/math] = probability of an energy transfer taking place = probability of an interaction = [math]\frac{N}{A} d \sigma[/math] [ Atoms [math]cm^2[/math]/g]

[math]dP(\Delta E) = \frac{N}{A} d \sigma =\frac{N}{A} (2 \pi b db) Z[/math]

In practice [math] \sigma[/math] is a measured cross-section which is a function of energy.

classically [math]\sigma = \pi b^2 ; d \sigma = 2\pi b db[/math] so let's use this as a first approximation

[math]\Rightarrow \frac{-dE}{dx} = \int_0^{\infty} \frac{N}{A} (2 \pi b db) Z \Delta E = \frac{2 \pi N Z}{A} \int_0^{\infty} \Delta E b db[/math]

= [math]\frac{2 \pi N Z}{A} \int_0^{\infty} \left [ \frac{2 r_e^2 m_e c^2 z^2}{\beta^2 b^2}\right ] b db[/math]

= [math]4 \pi N r_e^2 m_e c^2 \frac{z^2 Z}{A \beta^2} \int_0^{\infty} \frac{db}{b}[/math]

=[math]\frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \int_0^{\infty} \frac{db}{b}[/math]

where [math]\frac{\mathcal{K}}{A} = \frac{4 \pi N r_e^2 m_e c^2}{A} = 0.307 \frac{MeV cm^2}{g}[/math] if A=1

The limits of the above integral should be more physical in order to reflect the limits of the physics interaction. Let [math]b_{min}[/math] and [math]b_{max}[/math] represent the minimum and maximum possible impact parameter where the physics is described, as shown above, by the coulomb force.

What is [math]b_{min}[/math]?

if [math]b \rightarrow 0[/math] then [math]\frac{d E}{dx}[/math] diverges and the energy transfer [math]\rightarrow \infty : \Delta E \sim \frac{1}{b}[/math]. Physically there is a maximum energy that may be transferred before the physics of the problem changes (ie: nuclear excitation, jet production, ...). The de Borglie wavelength of the atom is used to estimate a value for [math]b_{min}[/math] such that

[math]b_{min} \sim \frac{1}{2} \lambda_{de Broglie} = \frac{h}{2p} = \frac{h}{2 \gamma m_e \beta c}[/math]

What is [math]b_{max}[/math]?

As [math]b[/math] gets bigger the interaction is "softer" and longer. If the interaction time ([math]\tau_i[/math]) is so long that it is equivalent to an electron orbit ([math]\tau_R[/math]) then the atom looks more like it is neutrally charged. You move from an interaction in which the electron orbit is perturbed adiabatically such that there is no orbit change and the minimum amount of energy is transferred to no interaction taking place because the atom is neutral.

Let

[math]\tau_i = \frac{b_{max}}{v} (\sqrt{1-\beta^2})[/math] : fields at high velocities get Lorentz contracted

[math]\tau_R \equiv \frac{h}{I}[/math] : I [math]\equiv[/math] mean excitation energy of target material ( [math]E = h \nu = h/ \tau[/math])

Condition for [math]b_{max}[/math] :

[math]\tau_i = \tau_R[/math]

[math]\Rightarrow b_{max} = \frac{h \gamma \beta c}{I}[/math]

[math]-\frac{dE}{dx} = \frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \int_0^{\infty} \frac{db}{b}[/math]

[math]= \frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \ln \frac{b_{max}}{b_{min}}[/math]

[math]= \frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \ln \frac{2 \gamma^2 m_e \beta^2 c^2}{I}[/math]

Example 5: Find [math]\frac{dE}{dx}[/math] for a 10 MeV proton hitting a liquid hydrogen ([math]LH_2[/math]) target

A = Z=z=1
[math]m_e c^2[/math] = 0.511 MeV
I = 21.8 eV : see gray data point for Liquid [math]H_2[/math] From Figure 27.5 on pg 6 of PDG below.

Just need to know [math]\gamma[/math] and [math]\beta[/math]

"a 10 MeV proton" [math]\Rightarrow[/math] Kinetic Energy (K.E.) = 10 MeV = [math](\gamma - 1) mc^2[/math]

[math]\Rightarrow \gamma = \frac{K.E.}{mc^2} + 1 = \frac{10 MeV}{938 MeV} + 1 \sim 1 = \frac{1}{\sqrt{1-\beta^2}}[/math]

Proton is not relativistic

[math]v^2 = \frac{2 K.E.}{m} = \frac{2 \cdot 10 MeV}{938 MeV/c^2} = 2 \times 10^{-2} c^2 \Rightarrow \beta^2 = \frac{v^2}{c^2} = 2\times 10^{-2}[/math]

Plugging in the numbers:

[math]\frac{dE}{dx} = \left ( 0.307 \frac{MeV \cdot cm^2}{g}\right ) (1)^2 (1) \frac{1}{2 \times10^{-2}} \ln \left( \frac{2 (1) (0.511 MeV) (2 \times10^{-2})}{21.8 eV} \frac{10^6 eV}{MeV}\right)[/math]

[math]= 105 \frac{MeV cm^2}{g}[/math]

How much energy is lost after 0.3 cm?

Notice that the units for energy loss are normalized by the density of the material [math]\rho_{LH_2}[/math] = 0.07 [math]\frac{g}{cm^3}[/math]

To get the actual energy lost I need to multiply by the density. So for any given atom the energy loss will depend on the state (solid, gas, liqid) of the atom as this effects the density of the material.

[math]\Delta E = (105 \frac{MeV cm^2}{g}) (0.07 \frac{g}{cm^3}) (0.3 cm)[/math] = 2.2 MeV

File:SPIM HydrogenStoppingPower.pdf Compare with Triumf Kinematics Handbook, 2nd edition, September 1987, L.G. Greeniaus

Bethe-Bloch Equation

While the classical equation above works in a limited kinematic regime, the Bethe-Bloch equation includes the corrections needed to cover most kinematic regimes for heavy particle energy loss.

[math]-\frac{dE}{dx} = \mathcal{K} z^2 \frac{Z}{A} \frac{1}{\beta^2} \left [ \frac{1}{2} \ln \left (\frac{2 m_e c^2 \beta^2 \gamma^2 }{I} \frac{ T_{max}}{I} \right) - \beta^2 - \frac{\delta}{2}\right ][/math]PDG reference Eq 27.1 pg 1

where

[math]T_{max} = \frac{2 m_e c^2 \beta^2 \gamma^2}{1+ \frac{2 \gamma m_e}{M} + \frac{m_e}{M}}[/math]

= Max K.E. transferable to the Target of mass [math]M[/math] in a single collision.

[math]-\beta^2[/math]

= correction for electron spin and very distant collisions which deforms the electron atomic orbits each process reducing dE/dx by [math]\frac{\beta}{2}[/math]

[math]\frac{\delta}{2}[/math]

= density correction term: in the classical derivation the material is treated as just a system of [math]N[/math] atoms uniformly distributed in space. These Atoms, however, give the material polarizability which can reduce the electric field (dielectric).

GEANT 4 implementation

The GEANT4 file (version 4.8.p01)

source/processes/electromagnetic/standard/src/G4BetheBlockModel.cc

is used to calculate hadron energy loss.

line 132 (line 257 in version 4.9.5) [math]\Rightarrow[/math]

[math]-\frac{dE}{dx} = \log \left ( \frac{2 m_e c^2 \tau (\tau +2) E_{min}}{I^2}\right) - \left (1 - \frac{E_{min}}{E_{max}} \right ) \beta^2[/math]

where

[math] \tau = \frac{K.E.}{M}[/math]

line 143 (line 267 in version 4.9.5) [math]\Rightarrow[/math]

[math]\frac{dE}{dx} -= \log ( \tau (\tau + 2) ) -cden[/math] = density corection = [math]\frac{\delta}{2}[/math]

line 148 (line 270 in vers. 4.9.5) [math]\Rightarrow[/math]

[math]\frac{dE}{dx} -= \frac{2c}{Z_{target}}[/math] = shell correction, corrects for the classical asumption that the atomic electron's velocity is initially zero; or the the incident particles velocity is far greater than the atomic electron's velocity.

line 154 (line 273 in version 4.9.5) [math]\Rightarrow[/math]

[math]\frac{dE}{dx} *= \frac{2 \pi m_e c^2 r_e^2 z^2}{\beta^2} \rho_e \;\;\;\; \rho_e \propto \frac{NZ}{A}[/math]

Energy Dependence

The above curve shows the energy loss per distance traveled ([math]\frac{dE}{dx}[/math]) as a function of the incident particles energy. There are three basic regions. At low incident energies ( < 10^5 eV) the incident particle tends to excite or even ionize the atoms in the material it is penetrating. The maximum amount of energy loss per distance traveled is defined as the Bragg peak. The region after the Bragg peak, in which the energy loss per distance traveled reaches its smallest value, is refered to as the point of minimum ionizing. Minimimum ionizing particles will have incident energies corresponding to this value or larger. The characteristic of the minimum ionizing particles is that their energy loss per distance traveled is essentially constant making simulations easier until the particle's energy drops below the minimum ionizing energy level as it passes through the material.

In general the Bethe-Bloch equation breaks down at low energies (below the Bragg peak) and is a good description (to within 10%) for

[math]10 \frac{MeV}{a.m.u.} \lt E \lt 2 \frac{GeV}{a.m.u.}[/math] and [math]Z[/math] < 26 (Iron) : a.m.u = Atomic Mass Unit

the [math]\frac{1}{\beta^2}[/math] term in the Bethe-Bloch equation dominates between the Bragg peak and the minimum ionization region.

the [math]\ln[/math] term and its corrections influence the dependence of [math]\frac{dE}{dx}[/math] as you move up in energy beyond the minimum ionization point.

Energy Straggling

While the Bethe-Bloch formula gives you a way to quantify the amount of energy a heavy charged particle loses as a function of the distance traveled, you should realize that when you calculate the total energy lost via

[math] \Delta E = \int_{E_i}^{E_f} \left ( \frac{dE}{dx} \right ) dx[/math]

you are only determining the AVERAGE energy loss. In other words, Bethe-Bloch is the Astochastic process describing energy loss.

In reality the energy loss process is a stochastic process because of the statistical fluctuations which occur in the actual number of collisions which take place.

Thick Absorber

A thick absorber is one in which a large number of collisions takes place. In this situation the central limit theorem from statistics tells you that the larger the number of random variable samples , [math]N[/math], involved the more the observable will follow a Gaussian distribution. The Gaussian distribution is a good approximation to the binomial distribution when the number of trials is large.

Forest_ErrAna_StatDist#Binomial_with_Large_N_becomes_Gaussian

, and to a Poisson distribution when the mean is a lot larger than 1.

Forest_ErrAna_StatDist#Gaussian_approximation_to_Poisson_when

The gaussian probability function is defined as

[math]P(x,\Delta) \propto e^{\frac{(\Delta - \bar{\Delta})^2}{ \sigma^2}}[/math]

where the Full Width at Half Max (FWHM) of the distribution = [math]\left ( 2 \sqrt{2 \ln 2} \right ) \sigma[/math]

In the case of energy loss, the variance using the Bethe-Bloch equation should be

[math]\sigma_0^2 = 4 \pi N r_e^2 (m_e c^2)^2 \rho \frac{Z}{A} x[/math]

the realitivistic variance is

[math]\sigma^2 = [\frac{1-\beta^2/2}{1-\beta^2} ]\sigma_0^2[/math]

for very thick absorbers see

C. Tschaler, NIM 64, (1968) 237 ; ibid, 61, (1968) 141

When simulating energy loss of heavy charged particles the Bethe-Bloch equation may be used to calculate a [math]\frac{dE}{dx}[/math] which can determine the average energy loss at the given kinetic energy of the particle. This average is then smeared according to a gaussian distribution of variance

[math]\sigma^2 =4 \pi N r_e^2 (m_e c^2)^2 \rho \frac{Z}{A} x [\frac{1-\beta^2/2}{1-\beta^2} ][/math]

Thin Absorbers

In thin absorbers the number of collisions is small preventing the use of the central limit theorem to describe the stochastic process of energy loss in terms of a Gaussian distribution. The large energy transfers that are possible cause the energy loss distribution to look like a Gaussian with a high energy tail (or foot).

The skewness of the resulting energy loss distribution is quantified as

[math]\kappa = \frac{\bar{\Delta}}{W_{max}}[/math]

[math]\Delta \equiv 2 \pi N_a r_e^2 m_e c^2 \rho \frac{Z}{A} \left ( \frac{z}{\beta}\right)^2 x [/math] = lead term in Bethe Bloch equation

[math]\rho[/math] = density of absorbing material.

[math]W_{max} = \frac{(pc)^2}{\frac{1}{2} \left [ m_e c^2 + \left ( \frac{M^2 c^2}{m_e} \right ) \right ] + \sqrt{(pc)^2 + (Mc^2)^2}}[/math] = max energy transfered in 1 collision (headon / knock out collision)

This comes from the relativistic kinematics of an Elastic Collision.

[math]\gamma = \frac{E_{tot}}{Mc^2} = \frac{ \sqrt{(pc)^2 + (Mc^2)^2}}{Mc^2}[/math]

[math]\beta= \frac{pc}{\gamma Mc^2} = \frac{pc}{E_{tot}}[/math]

[math]E_k = E_{tot} - Mc^2 = \gamma Mc^2 - Mc^2 = (\gamma - 1 ) Mc^2[/math]

[math]E_k = \sqrt{(pc)^2 + (Mc^2)^2} - Mc^2 [/math]

[math] (p^{\prime}c)^2 = E_k^2 + 2E_km_ec^2[/math]

Conservation of Momentum [math]\Rightarrow[/math] :

[math]\vec{p} = \vec{p}^{\; \prime \prime} + \vec{p}^{\; \prime}[/math]

Conservation of Energy [math]\Rightarrow[/math] :

[math]E_{tot} + m_ec^2 = E_{tot}^{\prime \prime} + E_{tot}^{\prime}[/math]

[math]\sqrt{(pc)^2 + (Mc^2)^2} + m_ec^2 = \sqrt{(p^{\; \prime \prime} c)^2 + (Mc^2)^2} + E_k + m_e c^2[/math]

using conservation of E & P as well as substituting for [math]p^{\prime}[/math] you can show

[math](p^{\; \prime \prime}c)^2 = (pc)^2 - 2E_k\sqrt{(pc)^2 +(Mc^2)^2} + E_k^2[/math] : cons of E

[math]= (pc)^2 + E_k^2 + 2E_km_ec^2 -2pc\sqrt{E_k^2+2E_km_ec^2} \cos(\theta)[/math] : cons of P

[math]\Rightarrow[/math]

[math]pc \cos(\theta) \sqrt{1+\frac{2m_ec^2}{E_k}} = \sqrt{(pc)^2+(Mc^2)^2} + m_ec^2[/math]

solving for [math]E_k[/math]

[math]E_k = \frac{2m_ec^2(pc)^2\cos^2 (\theta)}{[\sqrt{(pc)^2 + (Mc^2)^2} +m_ec^2]^2 - (pc)^2 \cos^2 (\theta)}[/math]

(Landau Theory)

[math]\kappa \leq 0.01[/math]

Landau assumed

1. [math]W_{max} = \infty[/math] is max energy transfer
2. electrons are free (energy transfer is so large you can neglect binding)
3. incident particle maintains velocity (large momentum transfer from big mass to small mass) (bowling ball hits ping pong ball)

L. Landau, "On the Energy Loss of Fast Particles by Ionization", J. Phys., vol 8 (1944), pg 201

instead of a gaussian distribution Landau used

[math]P(x,\Delta) \propto \frac{1}{\bar{\Delta}\pi} \int_0^{\infty} e^{-u \ln u - u \lambda} \sin(\pi u) du[/math]

where

[math]\lambda = \frac{1}{\bar{\Delta}} \left [ \Delta - \bar{\Delta} \ln \bar{\Delta} - \ln \epsilon + 1 -C \right ][/math]

[math]\bar{\Delta} = 2\pi N_a r_e^2 m_e c^2 \rho \frac{Zz^2}{A \beta^2}x[/math]

[math]\ln \epsilon = \ln \left [ \frac{(1-\beta^2)I^2}{2m_ec^2 \beta^2} \right ][/math]

[math]C = 0.577[/math]

(Vavilou's Theory)

Vavilous paper

P.V. Vavilou, "Ionization losses of High Energy Heavy Particles", Soviet Physics JETP, vol 5 (1950? )pg 749

describe the physics for the case

[math]0.01 \lt \kappa \lt \infty [/math]

The distribution function derived is shown below as well as a conceptual overlay of Vavilou's and Landau's distributions. (The [math]\zeta f(x,\Delta)[/math] in the picture should be a [math]\bar{\Delta}P(x,\Delta)[/math] )

[math]P(x,\Delta) = \frac{1}{\bar{\Delta}\pi} x e^{x(1+\beta^2C)} \int_0^{\infty} e^{xf_1} \cos(y \lambda_1 + xf_2) dy[/math]

where

[math]f_1 = \beta^2 \left [ \ln(y) - C_i(y)\right ] - \cos(y) - y S_i(y)[/math]

[math]f_2 = y\left [ \ln(y) - C_i(y)\right ] + \sin(y) + \beta^2 S_i(y)[/math]

[math]C_i(y) \equiv - \int_y^{\infty} \frac{\cos(t)}{t} dt[/math]

[math]S_i(y) \equiv \int_0^{y} \frac{\sin(t)}{t} dt[/math]

[math]C = 0.577[/math]

GEANT4's implementation

GEANT 4 uses the skewness parameter [math]\kappa[/math] to determine if it will use a "fluctuations model" to calculate energy straggling or the gaussian model described in section 3.2.1.

kappa > 10

If

[math]\kappa \equiv \frac{ \bar{\Delta}}{W_{max}}[/math] > 10

and we have a thick absorber ( large step size) then the Gausian function in 3.2.1 is used to calculate energy straggling.

What happens is [math]\Delta E[/math] is calculated via [math]\int_{E_i}^{E_f} \frac{dE}{dx} dx[/math] then the actual energy loss predicted by the simulation is chosen from a Gaussian distribution to account for energy straggling such that the [math]\sigma[/math] of this Gaussian distribution is given by:

[math]\sigma^2 = 2 \pi r_e^2m_ec^2N_{el} \frac{Z_h}{\beta^2} T_C s (1 - \frac{\beta^2}{2})[/math]

where

[math]N_{el}[/math] = electron density of the medium

[math]Z_h[/math] = charge of the incident particle

[math]s[/math] = step size

[math] T_C[/math] = cutoff kinetic energy for [math]\delta [/math]-electrons

[math]T_C[/math] tells GEANT where to put the cutoff for using the Gaussian distribution for energy straggling. This tells the simulation the low energy cutoff where Bethe-Bloch starts to fail due to ionization.

Delta-electrons

What is a [math]\delta[/math] - electron?

[math]\delta[/math] - electrons are also known as "knock -on" electrons and delta rays.

As heavy particles traverse a medium they can ionize electrons from atoms. The ejected electrons can be given enough energy to ionize as well.

In a cloud chamber (a supercooled volume of super saturated water vapor which ionizes as charged particles pass through) such and event would look like:

The blue spiral in the above gas chamber picture is a high energy electron ejected from a collision that spirals in the B-field ejecting low energy electrons at the end. The B-field is directed out of the picture.

The physics of ionization is different from the physics used to calculate Bethe-Bloch energy loss. Remember Bethe-Bloch starts to break down at low energies below the Bragg peak.

Because of this GEANT 4 sets the cutoff for this process to be

[math]T_{cut}[/math] > 1 keV

Note: The BE energies of an electron in Hydrogen is 13.6 ev and the electrons in Argon have binding energies between 15.7 eV and 3.2 keV.

Fluctuations Model: kappa < 10

If [math]\kappa \equiv \frac{ \bar{\Delta}}{W_{max}} \lt \frac{\Delta E}{T_C}[/math]

Then GEANT 4 uses a "Fluctuations Model" to determine energy loss instead of Bethe-Bloch.

Fluctuations Model

1. the atom is assumed to have on 2 energy levels [math]E_1[/math] and [math]E_2[/math]
2. you can excite the atom and lose either [math]E_1[/math] or [math]E_2[/math] or you can ionize the atom and lose energy according to a [math]\frac{1}{E^2}[/math] function [math]u_j[/math].

The total energy loss in a step will be

[math]\Delta E = \Delta E_{exc} + \Delta E_{ion}[/math]

where

[math]\Delta E_{exc} = \eta_1 E_1 + \eta_2 E_2[/math]

[math]\Delta E_{ion} = \sum_{j=1}^{\eta_3} \frac{I}{1 - u_j \frac{T_{up}-I}{T_{up}}}[/math]

[math]\eta_1[/math], [math]\eta_2[/math], and [math]\eta_3[/math] are the number of collisions which are sampled from a poison distribution

[math]u_j = \int_{I}^{E_j} \frac{I T_{up}}{T_{up} - I} \frac{dx}{x^2}[/math]

[math]E_j = \frac{I}{1- rand \frac{T_{up}-1}{T_{up}}}[/math] : rand = random number between 0 and 1

[math]T_{up} = \left \{ {~ 1 keV \; threshold \;energy \;for \; \delta- ray \; production \atop T_{max} \; \;\;\; if \; T_{max} \lt 1 keV} \right .[/math]

[math]I[/math] = mean ionization energy

[math]E_2 \approx (10 eV) Z^2[/math]

[math]\ln E_1 = \frac{\ln (I) - f_2 \ln (E_2)}{f_1}[/math]

[math]f_1 + f_2 =1[/math]

[math]f_2 =\left \{ {0 \; z=1 \atop \frac{2}{z} \; z \ge 2} \right .[/math]

The fluctuation model was comparted with data in

K. Lassila-Perini and L. Urban, NIM, A362 (1995) pg 416

The cross sections used for excitation and ionization may be found in

H. Bichel, Rev. Mod. Phys., vol 60 (1988) pg 663

Range Straggling

Def of Range (R)

The distance traveled before all the particles energy is lost.

[math]R \equiv \int_0^T \frac{dE}{\frac{dE}{dx}}[/math]

= theoretical calculation of the path length traveled by a particle of incident energy [math]T[/math]

Note units: [math]\left [ R \right ] = \frac{g}{cm^2} ; \left [ \frac{dE}{dx} \right ] = \frac{MeV \cdot cm^2}{g}[/math]

the Energy Straggling introduced in the previous section can explain why identical particles penetrate material to different depths. The energy straggling results in Range straggling.

If we do a shielding experiment where we have a source of incident particles of energy E and we count how many "punch" through a material of thickness (x) we would see a transmission coefficient [math]\left ( \frac{N_{out}}{N_{in}} \right) [/math] which would look like

Fractional Range Straggling

[math]\frac{\sigma_R}{R} \equiv[/math] fractional range straggling

Assuming the energy loss of a non-relativistic heavy ion through matter follows a Gaussian (thick absorber)

Then it can be shown that

[math]\frac{\sigma_R}{R} \approx \frac{1}{2} \sqrt{\frac{M}{A}}[/math]

where

[math]M[/math] = mass of the target electrons

[math]A[/math] = atomic mass of the Projectile

since

[math]m_e = 9.11 \times 10^{-31}[/math] kg

and

1 a.m.u. = [math]1.66 \times 10^{-27}[/math] kg

then

[math]\frac{\sigma_R}{R} \approx \frac{1}{2} \sqrt{\frac{9.11 \times 10^{-31}}{1.66 \times 10^{-27}A}}[/math]

= 1.17 % if using a proton (A=1)

The above is a "back of the envelope" estimate. The experimentally measured values for Cu, Al, and Be target using a proton projectile are

If the incident projectile is an electron then [math]\frac{\sigma_R}{R} \approx \frac{1}{2}[/math] making electron range straggling a vague concept.

There are several definitions of electron range

1.) Maximum Range ([math]R_0[/math])

This range is defined using the continuous slowing down approximation (CSDA) where electrons are assumed to have many collisions over very small distances making it appear to be continuous energy loss instead of discrete. The range is then calculated by integrating over these average energy losses [math]\frac{dE}{dx} \cdot s[/math].

2.) Practical Range ([math]R_P[/math])

This stopping distance is defined by extrapolating the electron transmission curve to zero (see below).

Electron Capture and Loss

Bohr Criterion

"A rapidly moving nucleus is fully ionized if its velocity exceeds that of its most tightly bound electron"

The Bohr Model:

[math]\Rightarrow E = \frac{mz^2e^4}{8 \epsilon_0^2 h^2 n^2}[/math]

for the inner most electron ([math]n=1[/math])

Electron K.E. = [math]\frac{1}{2} mv^2 = \frac{mz^2e^4}{2(4\pi \epsilon_0)^2 \hbar^2} \Rightarrow v = \frac{z e^2}{4 \pi \epsilon_0 \hbar}[/math]

the fine structure constant [math]\alpha \equiv \frac{e^2}{4 \pi \epsilon_0 \hbar c} = \frac{1}{137}[/math]

[math] v = zc \alpha[/math]

If [math]v \gt zc \alpha[/math] the nucleus is fully ionized

or

if [math]\frac{z}{v/c} = \frac{z}{\beta} \lt \frac{1}{\alpha} = 137[/math]

alternatively if the ion is moving through a material with a speed such that

[math]\frac{z}{\beta} \gt \frac{1}{\alpha} =137[/math]

Then electrons may be captured by the projectile and lost by the target.

Z-effective

Describing the charge state of your heavy ion traveling through matter at a velocity below the Bohr criterion is very complicated. There is a competition between electron capture and loss. Accurate cross sections are needed to simulate the process reliably.

Some insight into this process can be found using the Thomas-Fermi model

[math]V \propto \frac{Ze^{-r/a}}{r}[/math]

to describe an atom moving slow enough so it has captured many electrons but fast enough so its not neutral. In the Thomas-Fermi model the distribution of electrons in an atom is described as being uniformly distributed such that there are 2 electrons in each discrete volume of phase space( the space in which all possible states of a system are represented) defined using planks constant as [math]h^3[/math].

For the purpose of simulations you would like a relationship for [math]Z_{eff}[/math] in terms of [math]\beta[/math] and [math]Z[/math].

It is usually adequate to use fits for empirical data as long as we know that we are in the kinematic range in which those fits are valid.

when [math]E \lt 10[/math] MeV the data indicates that

[math]Z_{eff} = Z(1 - e^{-\beta\frac{B}{Z^{2/3}}})[/math]

where

[math]B = 130 \pm 5[/math]

[math]Z_{eff} \equiv[/math] effective charge f the projectile = [math]Z - \bar{q}_c[/math]

[math]Z[/math] = number of protons

[math]\bar{q}_c[/math] = average number of captured electrons

When calculating stopping power for E < 10 MeV you use [math]Z_{eff}[/math] in the Bethe-Bloch equation.

Note: As the ions charge state fluctuates while it slows down (or if accelerated through materials) you will need to recalculate the energy loss, and as a result you will get larger energy loss fluctuations in this energy range.

For thin absorber you will look for stripping and loss cross sections.

Here a thin absorber is one whose thickness is less than the charge equilibrium distance defined as the distance traveled until the projectile's velocity is [math] v \ll zc\alpha[/math]

A rule of thumb is that a thin absorber for low energy ions has a thickness [math]\le \frac{5 \frac{\mu g}{cm^2}}{\rho}[/math]

For thick absorbers: The experimentally determined expression for the change in [math]Z_{eff}[/math] from [math]Z[/math] is

[math]\Delta Z_{eff} = \frac{1}{2} \sqrt{ \left [ Z_{eff} \left (1 - \frac{Z_{eff}}{Z} \right )^{1.67}\right ] }[/math]

Multiple Scattering

The Bethe-Bloch equation tells us how much energy is lost and GEANT4s calculation of this energy is described above.

Now we need to know which direction the scattered particle goes after it has lost this energy.

The work of Moliere describes the angular deflection of the particle which lost the energy thereby leading to a prediction of the Cross-section. GEANT4 though uses the more complete Lewis theory to describe Multiple Coulomb Scattering (MCS) sometimes generically referred to as multiple scattering.

There are 3 regions in which coulomb scattering is calculated

1.) Single Scattering

For thin materials.

If the probability of more than 1 coulomb scattering is small

Then use the Rutherford formula for [math]\frac{d \sigma}{d \Omega}[/math]

2.)Multiple Scattering

In this case the number of independent scatterings is large (N > 20) and the energy loss is small such that the problem can be treated statisticaly to obtain a probability distribution for the net deflection angle [math] [P(\theta)][/math] as a function of the material thickness that is traversed.

3.) Plural Scattering

If 1< N [math]\le[/math] 20 then you can't use Rutherford to describe the scattering nor use a normal random statistical description.

see E. Keil, Z. Naturforsch, vol 15 (1960), pg 1031

Reviews of rigorous multiple scattering calculations may be found in

P.C. Hemmer, et. al., Phys. Rev, vol 168 (1968), pg 294

GEANT4's implementation of MSC (N>20)

GEANT4 models MSC when N>20 using model functions to determine the angular and spatial distributions chosen to give the same moments of these distributions as the Lewis theory.

H.W. Lewis, Phys. Rev., vol 78 (1950), pg 526

modern versions of the above are at

J.M. Fernandez-Varea, et. al., NIM, B73 (1993), pg 447

I. Kawrakow, et. al., NIM, B142 (1998) pg 253

When N>20 multiple scattering can be described as a statistical process using a modified version of the Boltzman transport equation from statistical mechanics.

Note

The simulation step size is chosen such that (N>20), If you have materials so thin that N < 20 then GEANT4 will likely skip the material. (one way around this is to increase the thickness and change the density). If the material thickness can't be increased because its sandwhiched between two other materials then you will need to write a special step algorithm for the volume and have GEANT4 use it for the step.

Let [math]f(s,\vec{x},\hat{v}) \equiv[/math] the distribution function for a system of incident particles traveling through a material.

where

[math]s =[/math] arc length of the particle's path through the material

[math]\vec{x} =[/math] position of a charged particle

[math]\hat{v} =[/math] direction of motion of the particle [math]\frac{\vec{v}}{|\vec{v}|}[/math]

The multiple scattering experienced by a single charged particle traveling through the material is then simulated by sampling from the distribution [math]f(s,\vec{x},\hat{v} )[/math]

The governing transport/diffusion equation is based on the continuity equation but with a "sink" term representing the possibility of collisions ejecting particles out of the volume.

[math]\frac{\partial f(s,\vec{x},\hat{v} ) }{\partial s} + \hat{v} \bullet \vec{\nabla}f(s,\vec{x},\hat{v} ) = N \int \sigma(\hat{v} \bullet\hat{v}^{\prime} )\left [ f(s,\vec{x},\hat{v}^{\prime} ) - f(s,\vec{x},\hat{v} ) \right ] d \hat{v}^{\prime}[/math]

where

[math]N[/math] = number of atoms per volume

[math]\sigma(\hat{v} \bullet\hat{v}^{\prime} )[/math] = cross sections for elastic scattering per Solid angle [math]\left ( \frac{d \sigma}{d \Omega} \right )[/math]

To solve the above diffusion equation the distribution function, [math]f(s,\vec{x},\hat{v} )[/math] is expanded in Spherical Harmonics ( [math]Y_{\ell}^m(\theta,\phi)[/math] ) and [math]\sigma[/math] expand in Legendre Polynomials ([math]P_N(cos \theta)[/math]) since it has no [math]\phi[/math] angle dependence.

Note

For Coulomb Scattering in polar coordinates you can write the potential in terms of Legendre Polynomials such that:

[math]U=k \frac{q}{r}[/math]

= [math]k\frac{q}{\sqrt{r^2-a^2-2ar \cos \theta}}[/math] in polar coordinates

= [math]k\frac{q}{r} \sum_{n=0}^{\infty} P_n(\cos \theta) \left ( \frac{a}{r}\right )^n[/math] (the sqrt term above is expanded using binomial series

[math]f(s,\vec{x},\hat{v} ) = \sum_{\ell,m} f_{\ell,m}(\vec{x},s) Y_{\ell}^m(\hat{v})[/math]

after substituting into the diffusion equation and doing the integral on the righ hand side you get

[math]\frac{\partial f_{\ell,m}(\vec{x},s) }{\partial s} + \frac{f_{\ell,m}(s,\vec{x},\hat{v} }{\lambda_{\ell}} = - \sum_{\lambda, i, j} \vec{\nabla} f_{i,j}(\vec{x},s ) \bullet \int Y_{\ell,m}^{\star} \hat{v} Y_{i,j} d \hat{v} \; \; \; \; \; \; \; \;\hat{v} = f(\theta.\phi)[/math]

where

[math]\frac{1}{\lambda_{\ell}} = 2 \pi N \int_0^{\pi} \left [ 1-P_{\ell}(\cos \theta)\right ] \sigma(\theta) \sin(\theta) d \theta[/math] = [math]\ell^{th}[/math] transport mean free path for the [math]f_{\ell}[/math] distribution function ( [math]\phi[/math] symmetry is assumed making it [math]m[/math] independent)

From the above one can find the average distances traveled and the average deflection angle of the distribution. Again, see :

J.M. Fernandez-Varea, et. al., NIM, B73 (1993), pg 447

The "moments" of [math]f(s,\vec{x},\hat{v}) [/math] are defined as

[math]\lt z\gt = 2 \pi \int z f(s,\vec{x},\hat{v}) \sin(\theta) d \theta d |\vec{x}| = \lambda_1 \left [ 1-e^{-s/\lambda_1}\right ][/math] = mean geometrical path length

[math]\lt \cos(\theta)\gt = 2 \pi \int_{-1}^1 \sum_{\ell} P_{\ell}(\cos \theta) \int f(s,\vec{x},\hat{v}) \sin(\theta) d \theta d |\vec{x}| = e^{-s/\lambda_1} [/math]

[math]\frac{1}{\lambda_1} = 2 \pi N \int_0^{\pi} \left [ 1-P_1(\cos \theta)\right ] \sigma(\theta) \sin(\theta) d \theta[/math]

Notice there are 3 lengths

[math]s[/math] = geometrical path length between endpoints of the step =[math] \left \{ {line \; if \; \vec{B} = 0 \atop arc \; if \; \vec{B} \ne 0 } \right .[/math]

[math]t[/math] = true path length = actual length of the path taken by particle

[math]\lt z\gt [/math] - mean geometrical path length along the z-axis

In GEANT4 the [math]\lambda_{\ell}[/math]'s are taken from

If 100 eV < K.E. of electron or positron < 10 MeV

D. Liljequist, J. Applied Phys, vol 62 (1987), 342

J. Applied Phys, vol 68 (1990), 3061

If K.E. > 10 MeV

R. Mogol, Atomic Data, Nucl, Data tables, vol 65 (1997) pg 55

with <z> now known GEANT will try to determine "[math]t[/math]" for the energy loss and scattering calculations.

A model is used for this where

[math]t=\frac{1}{\alpha} \left [ 1 - (1- \alpha \omega z)^{\frac{1}{\omega}})\right ][/math]

where

[math]\omega = 1 + \frac{1}{\alpha \lambda_{10}}[/math]

[math]\alpha =\left \{ {\frac{\lambda_{10} - \lambda_{11}}{s \lambda_{10}}\;\;\;\; K.E. \ge M_{particle} \atop \frac{1}{R}\;\;\;\; K.E. \lt M_{particle}} \right .[/math]

[math]s[/math] = stepsize

[math]\lambda_{10} - \frac{\lambda_1}{1-\alpha s}[/math]

[math]\lambda_{11} = \lambda_1[/math] at end of strep

while [math]\lt cos \theta \gt [/math] is calculable, GEANT4 evaluates [math]\cos (\theta)[/math] from a probability distribution whose general form is

[math]g[\cos(\theta)] = p \left ( qg_1[\cos(\theta)] + (1-q)g_3[\cos(\theta)] \right ) + (1-p)g_2[\cos(\theta)][/math]

where

[math]g_1(x) = C1e^{-a(1-x)}[/math]

[math]g_2(x) = \frac{C_2}{(b-x)^d}[/math]

[math]g_3(x) = C_3[/math]

[math]C_1, C_2, C_3[/math] are normalization constants

[math]p,q,a,b,d[/math] are parameters which follow the work reported in

V.L. Highland, NIM, vol 219 (1975) pg497

The GEANT4 files in version 4.8 were located in

/source/processes/electromagnetic/utils/src/G4VMultipleScattering.cc

and

/source/processes/electromagnetic/standard/src/G4MscModel.cc

/source/processes/electromagnetic/standard/src/G4MultipleScattering.cc

Simulations_of_Particle_Interactions_with_Matter