Stopping Power Calculation

Introduction
The code
- Translation of the formula
- Physics and code variables
Test with real values
More

Introduction

Here we explain a bit of the physics behind this macro.

According to the Particle Data Group, a particle entering a material will have a mean rate of energy loss (also known as stopping power) defined by the (in)famous Bethe-Bloch formula:

$- \frac{d E}{d x} = K z^{2} \frac{Z}{A} \frac{1}{β^{2}} [\frac{1}{2} \ln \frac{2 m_{e} c^{2} β^{2} γ^{2} T_{m a x}}{I^{2}} - β^{2}]$

where dE represents the amount of energy E that is lost by the particle per unit of distance dx traveled inside a material. Since it is a lost, the decrease is represented by a negative sign (-).

What follows after the equal sign is a series of constants and parameters that define the particle and the material. Not all of them, only the ones that have an influence in the loss of energy. In the following table, these variables are introduced (the numbers in parentheses represent the error of the value):

Symbol	Definition	Units or value
E	Incident particle energy	MeV^[*]
K	4π N_A r_e² m_e c²	0.307 075 MeV cm² / g
z	Charge of the incident particle	unitless
Z	Atomic number of the material	unitless
A	Atomic mass of the material	g / mol
β	velocity in terms of the speed of light c	unitless
γ	Lorentz factor	unitless
m_e c²	Electron mass x c²	0.510 998 918(44) MeV
T	Kinetic energy	MeV
I	Mean excitation energy	eV

[*] Electron-volts (eV) are often used as units of energy and mass in physics, using a system of natural units with c and ℏ set to 1.

This beautiful formula looks something like this when we put real numbers in it. In this case, the values for muons in copper:

The code

Translation of the formula

The macros are written in C++. The function stoppingPower() in simulation1.cpp and simulation2.cpp calculates this formula. To translate it to code, first we calculate gamma (the Lorentz factor), which depends on the energy of the particle:

$γ = 1 + \frac{E}{m}$

and so we have:

g = 1 + E/m;

We can then play with different particles by changing the value of the mass (m), or launch a particle with different values of energy (E), or combine both to make a study on different particles with different energies.

Next, we calculate beta squared, the speed of the particle:

$β^{2} = 1 + \frac{1}{γ^{2}}$

which translates in the code to:

b2 = 1 - (1/(g*g));

The kinetic energy of the particle is defined by:

$T = \frac{2 m_{e} c^{2} β^{2} γ^{2}}{1 + 2 γ \frac{m_{e}}{m} + {(\frac{m_{e}}{m})}^{2}}$

since m_e c² ≈ 0,511 MeV, and so 2m_e c² = 1.022 MeV, we can shorten the expression for T like this:

T = (1.022*b2*g*g) / ( 1 + 2*g*(0.511/m) + pow(0.511/m,2) );

And finally, setting the value of the constant K as double k = 0.307075; (see table), and calling the charge of the particle (z) q, we get to calculate the stopping power:

$- \frac{d E}{d x} = K z^{2} \frac{Z}{A} \frac{1}{β^{2}} [\frac{1}{2} \ln \frac{2 m_{e} c^{2} β^{2} γ^{2} T_{m a x}}{I^{2}} - β^{2}]$

If the mass and energies are given in MeV, the distances in cm, and A in g / mol, the units of the stopping power will be of MeV cm² / g. We can multiply the stopping power by ρ₀, the density of the material in g / cm³, to have it in units of MeV / cm. That way, it will be more straight-forward to calculate the energy loss per unit of travelled distance. Multiplying the stopping power by the travelled distance will automatically give the energy lost by travelling that distance.

The stopping power expression translated into code is:

SP = ((ro*k*q*q*Z)/(A*b2))*(0.5*log( (1.022*b2*g*g*T) / (I*I) ) - b2);

The final subroutine looks like this:

double stoppingPower(double m, double q, double E, double ro, double Z, double A, double I) {

	double k = 0.307075;
	double g, b2, T, SP;

	printf("\n|******** Stopping-Power values for this material ********|\n");

	g = 1 + (E/m);
	printf("\n Gamma  =  %lf", g);

	b2 = 1 - (1/(g*g));
	printf("\n Beta squared  =  %lf", b2);

	T = (1.022*b2*g*g) / (1 + 2*g*(0.511/m) + pow(0.511/m,2) );
	printf("\n Maximum kinetic energy transferred  =  %lf MeV ", T);

	SP = ((ro*k*q*q*Z)/(A*b2))*(0.5*log( (1.022*b2*g*g*T) / (I*I) ) - b2);
	printf("\n Stopping Power  =  %lf MeV/cm \n", SP);

	return SP;
}

Physics and code variables

The simulation is simple: we have a material surrounded by air. At time 0, a particle enters a material of a certain thickness. To track the travelling of the particle in the material, we use steps of a fixed size. Inside a loop, we calculate how much energy the particle loses in each step, and determine if it is stopped inside the material or if it goes out, in which case we repeat the process to calculate the energy loss in air.

Here is a diagram that shows what we have just said:

Diagram of a particle entering a material

The first simulation uses the remaining distance to the exit L:

for(n = 0, L = thickness; L > 0 && E > 0; n++) {
	E = E - SP*step;
	L = thickness - n*step;
}

where n is a counter for the number of steps, and we exit the loop when whether the particle is stopped (E > 0 condition) or it has abandoned the material (L > 0 condition).

The second simulation uses the range:

for(n = 0; range < thickness && E > 0; n++) {
	SP = stoppingPower(m, q, E, ro, Z, A, I, &G, &B, &t);
	E = E - SP*step;
	range += step;
}

Also, the second simulation calculates the stopping power every time the energy changes. The first simulation doesn't because it is intended to be a rough approximation.

Test with real values

Play with the values of energy of the particle and thickness of the material to see what happens. Here is an example of a proton particle entering an aluminium foil. The values can be found in the Wikipedia or in specialized publications, just pay attention to the units and make the necessary conversions:

Parameter	Value
Particle values (proton)
Mass (m)	938.272 MeV / c²
Charge (z or `q`)	1 (in units of e)
Initial kinetic energy (E)	100 MeV^[*]
Material values (Al foil)
Density (ρ₀)	2.7 g / cm³
Atomic number (Z)	13
Mass number (A)	27
Ionization Potential (I)	0.000006 MeV
Material thickness	1 cm
Step size	0.01 cm^[**]
Stopping-Power values for this material
Gamma (γ)	1.106579
Beta squared (β²)	0.183351
Maximum kinetic energy transferred (T_max)	0.229180 MeV
Stopping Power (- dE / dx)	22.572949 MeV/cm
Result	The particle exits the material with an energy of 77.427051 MeV in 100 steps.

[*] You can make a study in energy dependence changing this value systematically.
[**] If more precision is needed, introduce a shorter step size.

If the value of the initial energy is changed, a plot can be made of the dependence of the stopping power with the energy of the incidental particle, or a plot of the energy against the range, like the one below:

Plot of the initial energy against the range