partmatter.asciidoc

{baseurl}/Manual.asciidoc[Back to Main Page]

Table of Contents

• 1. Physics Models Used in the Particle Matter Interaction Model
  • 1.1. The Energy Loss
  • 1.2. The Coulomb Scattering
    • 1.2.1. Multiple Coulomb Scattering
    • 1.2.2. Large Angle Rutherford Scattering
  • 1.3. Beam Stripping physics
    • 1.3.1. Residual gas interaction
    • 1.3.2. Electromagnetic stripping
  • 1.4. The ScatteringPhysics Substeps
    • 1.4.1. The Flow Diagram of ScatteringPhysics Class in OPAL
  • 1.5. Available Materials in OPAL
  • 1.6. Example of an Input File
  • 1.7. A Simple Test
  • 1.8. References

1. Physics Models Used in the Particle Matter Interaction Model

The command to define particle-matter interactons is PARTICLEMATTERINTERACTION.

TYPE

Specifies the particle-matter interaction handler. Currently, there are the two following implemented methods: SCATTERING for physical processes of beam scattering and energy loss by heavy charged particles and BEAMSTRIPPING for interactions with residual gas and Lorentz stripping.

MATERIAL

The material of the surface (see Available Materials in OPAL).

ENABLERUTHERFORD

Switch to disable Rutherford scattering (default: true).

LOWENERGYTHR

Low-energy threshold [MeV] for energy loss calculation. Particles with lower energy will be removed (default: 0.01 MeV).

The so defined instance has then to be added to an element using the attribute.

1.1. The Energy Loss

The energy loss is simulated using the Bethe-Bloch equation.

-\frac{\mathrm{d} E}{\mathrm{d} x}=\frac{K z^2 Z}{A \beta^2}\left[\frac{1}{2} \ln{\frac{2 m_e c^2\beta^2 \gamma^2 Tmax}{I^2}}-\beta^2 \right],

where Z is the atomic number of absorber, A is the atomic mass of absorber, m_e is the electron mass, z is the charge number of the incident particle, K=4\pi N_Ar_e^2m_ec^2, r_e is the classical electron radius, N_A is the Avogadro’s number, I is the mean excitation energy. \beta and \gamma are kinematic variables. T_{max} is the maximum kinetic energy which can be imparted to a free electron in a single collision.

T_{max}=\frac{2m_ec^2\beta^2\gamma^2}{1+2\gamma m_e/M+(m_e/M)^2},

where M is the incident particle mass.

This expression is valid for energies from 600 keV to 10 GeV for incident beams (see PARTICLE definition) of PROTON, DEUTERON, MUON, HMINUS and H2P; and also for ALPHA particles from 10 MeV to 1 GeV.

The stopping power is compared with PSTAR program of NIST in Figure 1.

Figure 1. The comparison of stopping power with PSTAR.

The energy loss in the low-energy region is calculated using semi-empirical fitting formulas of Andersen and Ziegler [9].

-\frac{\mathrm{d} E}{\mathrm{d} x}=10^{-21}\frac{N_A}{A} \cdot \epsilon,

where the energy loss is in MeV cm^2/g and \epsilon is a fitted function of the stopping cross section.

\epsilon = \frac{\epsilon_{low}\cdot\epsilon_{high}}{\epsilon_{low}+\epsilon_{high}}

In case of incident protons in the material for energies from 10 keV to 600 keV, the fitting functions are given by:

\epsilon_{low} = A2 \cdot T_s^{0.45}

\epsilon_{high} = \frac{A3}{T_s} \ln{\left(1 + \frac{A4}{T_s} + A5 \cdot T_s\right)}

where T_s is the kinetic energy (in keV) divided by the proton mass (in amu). For T_s between 1 and 10 keV, the fitted function is given by:

\epsilon = A1 \cdot T_s^{0.5}

In case of incident alpha particles for energies from 1 keV to 10 MeV, the stopping power functions is expressed as:

\epsilon_{low} = B1\cdot(1000T)^{B2}

\epsilon_{high} = \frac{B3}{T} \ln{\left(1 + \frac{B4}{T} + B5 \cdot T\right)}

where T is the kinetic energy in MeV. The numerical values of coefficients of the empirical formulas are showed in Table 2.

The particles lost due to excessive energy loss (final energy null or lower than LOWENERGYTHR) are recorded in a HDF5 file (or ASCII if ASCIIDUMP is true). The loss file name is compound by the associated element name and the PARTICLEMATTERINTERACTION name.

Energy straggling: For relatively thick absorbers such that the number of collisions is large, the energy loss distribution is shown to be Gaussian in form. For non-relativistic heavy particles the spread \sigma_0 of the Gaussian distribution is calculated by:

\sigma_0^2=4\pi N_Ar_e^2(m_ec^2)^2\rho\frac{Z}{A}\Delta s,

where \rho is the density, \Delta s is the thickness.

1.2. The Coulomb Scattering

The Coulomb scattering is treated as two independent events: the multiple Coulomb scattering and the large angle Rutherford scattering. Using the distribution given in Classical Electrodynamics, by J. D. Jackson [1], the multiple- and single-scattering distributions can be written:

P_M(\alpha) \;\mathrm{d} \alpha=\frac{1}{\sqrt{\pi}}e^{-\alpha^2}\;\mathrm{d}\alpha,

P_S(\alpha) \;\mathrm{d} \alpha=\frac{1}{8 \ln(204 Z^{-1/3})} \frac{1}{\alpha^3}\;\mathrm{d}\alpha,

where \alpha=\frac{\theta}{<\Theta^2>^{1/2}}=\frac{\theta}{\sqrt 2 \theta_0}.

The transition point is \theta=2.5 \sqrt 2 \theta_0\approx3.5 \theta_0,

\theta_0=\frac{{13.6}{MeV}}{\beta c p} z \sqrt{\Delta s/X_0} [1+0.038 \ln(\Delta s/X_0)],

where p is the momentum, \Delta s is the step size, and X_0 is the radiation length.

1.2.1. Multiple Coulomb Scattering

Generate two independent Gaussian random variables with mean zero and variance one: z_1 and z_2. If z_2 \theta_0>3.5 \theta_0, start over. Otherwise,

x=x+\Delta s p_x+z_1 \Delta s \theta_0/\sqrt{12}+z_2 \Delta s \theta_0/2,

p_x=p_x+z_2 \theta_0.

Generate two independent Gaussian random variables with mean zero and variance one: z_3 and z_4. If z_4 \theta_0>3.5 \theta_0, start over. Otherwise,

y=y+\Delta s p_y+z_3 \Delta s \theta_0/\sqrt{12}+z_4 \Delta s \theta_0/2,

p_y=p_y+z_4 \theta_0.

1.2.2. Large Angle Rutherford Scattering

Generate a random number \xi_1, if \xi_1 < \frac{\int_{2.5}^\infty P_S(\alpha)d\alpha}{\int_{0}^{2.5} P_M(\alpha)\;\mathrm{d}\alpha+\int_{2.5}^\infty P_S(\alpha)\;\mathrm{d}\alpha}=0.0047, sampling the large angle Rutherford scattering.

The cumulative distribution function of the large angle Rutherford scattering is

F(\alpha)=\frac{\int_{2.5}^\alpha P_S(\alpha) \;\mathrm{d} \alpha}{\int_{2.5}^\infty P_S(\alpha) \;\mathrm{d} \alpha}=\xi,

where \xi is a random variable. So

\alpha=\pm 2.5 \sqrt{\frac{1}{1-\xi}}=\pm 2.5 \sqrt{\frac{1}{\xi}}.

Generate a random variable P_3,
if P_3>0.5

\theta_{Ru}=2.5 \sqrt{\frac{1}{\xi}} \sqrt{2}\theta_0,

else

\theta_{Ru}=-2.5 \sqrt{\frac{1}{\xi}} \sqrt{2}\theta_0.

The angle distribution after Coulomb scattering is shown in Figure 2. The line is from Jackson’s formula, and the points are simulations with Matlab. For a thickness of \Delta s=1e-4 m, \theta=0.5349 \alpha (in degree).

Figure 2. The comparison of Coulomb scattering with Jackson’s book.

1.3. Beam Stripping physics

Beam stripping physics takes into account two different physical processes: interaction with residual gas and electromagnetic stripping (also called Lorentz stripping). Given the stochastic nature of such interactions, the processes are modeled as a Monte Carlo method.

Both processes are described in the same way: Assuming that particles are normally incident on a homogeneous medium and that they are subjected to a process with a mean free path \lambda between interactions, the probability of suffering an interaction before reaching a path length x is:

P(x) = 1 - e^{-x/\lambda}

where P(x) is the statistic cumulative interaction probability of the process.

Figure 3 summarizes the iterative steps evaluated by the algorithm.

Figure 3. The diagram of BeamStrippingPhysics in OPAL.

1.3.1. Residual gas interaction

The mean free path of the interaction is related with the density of interaction centers and the cross section: \lambda=1/N\sigma. Assuming a beam flux incident in an ideal gas with density N (number of gas molecules per unit volume under the vacuum condition), the number of particles interacting depends on the gas composition and the different reactions to be considered. Thus, in agreement with Dalton’s law of partial pressures:

\frac{1}{\lambda_{total}}=\sum \frac{1}{\lambda_k}=N_{total}\cdot\sigma_{total}=\sum_jN_j\,\sigma^{j}_{total}=\sum_j\left(\sum_iN_j\,\sigma^{j}_{i}\right)

where the first summation is over all gas components and the second summation is over all processes for each component.

The fraction loss of the beam for unit of travelled length is: f_g=1-e^{-\delta _s /\lambda_{total}}. For a individual particle, f_g represents the interaction probability and it is evaluated through an independent random variable each step \delta _s.

Gas stripping could be applied for four different types of incoming particles: negative hydrogen ions (HMINUS), protons (PROTON), neutral hydrogen atoms (HYDROGEN), hydrogen molecule ions (H2P) and deuterons (DEUTERON). Single / Double - electron detachment or capture reactions are considered for each of them.

Cross sections are calculated according to energy of the particle employing analytic expressions fitted to cross section experimental data. The suitable function is selected in each case according to the type of incident particle and the residual gas under consideration. There are different fitting expressions:

• Nakai function [2]

\sigma_{qq'} = \sigma_0 \left[ f(E_1) + a_7\!\cdot\!f(E_1/a_8) \right]

f(E) = \frac{ a_1\!\cdot\!\left(\!\displaystyle\frac{E}{E_R}\!\right)^{\!a_2} }{ 1+\left(\!\displaystyle\frac{E}{a_3}\!\right)^{\!a_2+a_4}\!\!+\left(\!\displaystyle\frac{E}{a_5}\!\right)^{\!a_2+a_6} }

E_R = hcR_{\infty}\!\cdot\!\frac{m_H}{m_e}

E_1 = E_0 - E_{th}

where \sigma_0 is a convenient cross section unit (\sigma_0 = 1\cdot 10^{-16}\,\text{cm}^2), E_0 is the incident projectile energy in keV, E_{th} is the threshold energy of reaction in keV, and the symbols a_i denote adjustable parameters.

• Tabata function [3]: A linear combination of the Nakai function, f(E), improved and fitted with a greater number of experimental data and considering more setting parameters. The enhancement of the function makes it possible to extrapolate the cross section data to some extent.

• Barnett function [4]:

\ln{[\sigma(E)]}=\frac{1}{2}a_0 + \sum_{i=1}^{k} a_i \cdot T_i(X)

X=\frac{(\ln{E}-\ln{E_{min}})-(\ln{E_{max}}-\ln{E})}{\ln{E_{max}}-\ln{E_{min}}}

where T_i are the Chebyshev orthogonal polynomials.

• Bohr function [5]:

\sigma =
\left\{\begin{array}{l}
4\pi a_0^2 \displaystyle\frac{z_t+z_t^2}{z_i}\left(\frac{v_0}{v}\right)^2 \hspace{8mm} z_t < 15\\
\\
\pi a_0^2 \displaystyle\frac{z_t^{2/3}}{z_i}\left(\frac{v_0}{v}\right) \hspace{14mm} z_t > 15
\end{array}\right.

where z_i and z_t are the charge of the incident particle and the charge of the target nuclei, a_0 is the Bohr radius, v_0=e^2/4\pi\varepsilon_0\hslash is the characteristic Bohr velocity and v is the incident particle velocity.

1.3.2. Electromagnetic stripping

In case of HMINUS particles, the second electron is slightly bounded, so it is relevant to consider the detachement by the magnetic field. The orthogonal component of the magnetic field to the median plane (read from CYCLOTRON element), produces an electric field according to Lorentz transformation, E\!=\!\gamma\beta cB. The fraction of particles dissociated by the electromagnetic field during a time step \delta _t is:

f_{em}=1-e^{\,-\,\delta _t/\gamma\tau}

where \tau is the particle lifetime in the rest frame, determined theoretically [6]:

\tau= \frac{4mz_T}{S_0N^2\hslash\,(1+p)^2\left(1-\displaystyle\frac{1}{2k_0z_t}\right)}\cdot\exp\!{\left(\frac{4k_0z_T}{3}\right)}

where z_T\!=\!-\varepsilon_0/eE is the outer classical turning radius, \varepsilon_0 is the binding energy, p is a polarization factor of the ionic wave function, k_0^2\!=\!2m(-\varepsilon_0)/\hslash^2, S_{\!0} is a spectroscopy coefficient, and the normalization constant is N=(2k_0(k_0+\alpha)(2k_0+\alpha))^{1/2} / \alpha, where \alpha is a parameter for the ionic potential function.

The electromagnetic stripping calculation is restricted to OPAL-cycl.

1.4. The ScatteringPhysics Substeps

Small step is needed in the routine of ScatteringPhysics.

If a large step is given in the main input file, in the file ScatteringPhysics.cpp, it is divided by a integer number n to make the step size using for the calculation of scattering physics less than 1.01e-12 s. As shown by Figure 4 and Figure 5 in the following subsection, first we track one step for the particles already in the element and the newcomers, then another (n-1) steps to make sure the particles in the element experience the same time as the ones in the main bunch.

Now, if the particle leave the element during the (n-1) steps, we track it as in a drift and put it back to the main bunch when finishing (n-1) steps.

1.4.1. The Flow Diagram of ScatteringPhysics Class in OPAL

Figure 4. The diagram of ScatteringPhysics in OPAL.

Figure 5. The diagram of ScatteringPhysics in OPAL (continued).

1.5. Available Materials in OPAL

Different materials have been implemented in OPAL for scattering interactions and energy loss calculation. The materials that are supported are listed in Table 1, including the atomic number, Z, the atomic weight, A, the mass density, \rho, the radiation lenght, X0, and the mean excitation energy, I. In addition, the coefficients of the Andersen-Ziegler empirical formulas for the stopping power in the low-energy region are illustrated in Table 2.

Table 1. List of materials with their atomic and nuclear properties [7] [8].

Material (OPAL Name)	Z	A	\rho [g/cm^3]	X0 [g/cm^2]	I [eV]
Air	7	14	1.205E-3	36.62	85.7
Aluminum	13	26.9815384	2.699	24.01	166.0
AluminaAl2O3	50	101.9600768	3.97	27.94	145.2
Beryllium	4	9.0121831	1.848	65.19	63.7
BoronCarbide	26	55.251	2.52	50.13	84.7
Copper	29	63.546	8.96	12.86	322.0
Gold	79	196.966570	19.32	6.46	790.0
Graphite	6	12.0172	2.210	42.7	78.0
GraphiteR6710	6	12.0172	1.88	42.7	78.0
Kapton	6	12	1.420	40.58	79.6
Molybdenum	42	95.95	10.22	9.80	424.0
Mylar	6.702	12.88	1.400	39.95	78.7
Titanium	22	47.867	4.540	16.16	233.0
Water	10	18.0152	1	36.08	75.0

Table 2. List of materials with the coefficients of the Andersen-Ziegler empirical formulas for the stopping power in the low-energy region [9].

Material (OPAL Name)	A1	A2	A3	A4	A5	B1	B2	B3	B4	B5
Air	2.954	3.350	1.683e3	1.900e3	2.513e-2	1.9259	0.5550	27.15125	26.0665	6.2768
Aluminum	4.154	4.739	2.766e3	1.645e2	2.023e-2	2.5	0.625	45.7	0.1	4.359
AluminaAl2O3	1.187e1	1.343e1	1.069e4	7.723e2	2.153e-2	5.4	0.53	103.1	3.931	7.767
Beryllium	2.248	2.590	9.660e2	1.538e2	3.475e-2	2.1895	0.47183	7.2362	134.30	197.96
BoronCarbide	3.519	3.963	6065.0	1243.0	7.782e-3	5.013	0.4707	85.8	16.55	3.211
Copper	3.969	4.194	4.649e3	8.113e1	2.242e-2	3.114	0.5236	76.67	7.62	6.385
Gold	4.844	5.458	7.852e3	9.758e2	2.077e-2	3.223	0.5883	232.7	2.954	1.05
Graphite	0.0	2.601	1.701e3	1.279e3	1.638e-2	3.80133	0.41590	12.9966	117.83	242.28
GraphiteR6710	0.0	2.601	1.701e3	1.279e3	1.638e-2	3.80133	0.41590	12.9966	117.83	242.28
Kapton	0.0	2.601	1.701e3	1.279e3	1.638e-2	3.83523	0.42993	12.6125	227.41	188.97
Molybdenum	6.424	7.248	9.545e3	4.802e2	5.376e-3	9.276	0.418	157.1	8.038	1.29
Mylar	2.954	3.350	1683	1900	2.513e-02	1.9259	0.5550	27.15125	26.0665	6.2768
Titanium	4.858	5.489	5.260e3	6.511e2	8.930e-3	4.71	0.5087	65.28	8.806	5.948
Water	4.015	4.542	3.955e3	4.847e2	7.904e-3	2.9590	0.53255	34.247	60.655	15.153

1.6. Example of an Input File

particlematterinteraction.in

FX5 is a slit in x direction, the APERTURE is POSITIVE, the first value in APERTURE is the left part, the second value is the right part. FX16 is a slit in y direction, the APERTURE is NEGATIVE, the first value in APERTURE is the down part, the second value is the up part.

1.7. A Simple Test

A cold Gaussian beam with \sigma_x=\sigma_y=5 mm. The position of the collimator is from 0.01 m to 0.1 m, the half aperture in y direction is 3 mm. Figure 6 shows the trajectory of particles which are either absorbed or deflected by a copper slit. As a benchmark of the collimator model in OPAL, Figure 7 shows the energy spectrum and angle deviation at z=0.1 m after an elliptic collimator.

Figure 6. The passage of protons through the collimator.

Figure 7. The energy spectrum and scattering angle at z=0.1 m

1.8. References

[1] J. D. Jackson, Classical Electrodynamics, John Wiley & Sons, 3rd ed. (1999).

[2] Y. Nakai et al., Cross sections for charge transfer of hydrogen atoms and ions colliding with gaseous atoms and molecules, At. Dat. Nucl. Dat. Tabl., 37, 69 (1987).

[3] T. Tabata and T. Shirai, Analytic Cross Section for Collisions of H+, H2+, H3+, H, H2, and H- with Hydrogen Molecules, At. Dat. Nucl. Dat. Tabl., 76, 1 (2000).

[4] C. F. Barnett, Atomic Data for Fusion Vol. I: Collisions of H, H2, He and Li atoms and ions atoms and molecules, Tech. Rep. ORNL-6086/V1, Oak Ridge National Laboratory (1990).

[5] H.-D. Betz, Charge States and Charge-Changing Cross Sections of Fast Heavy Ions Penetrating Through Gaseous and Solid Media, Rev. Mod. Phys. 44, 465 (1972).

[6] L. R. Scherk, An improved value for the electron affinity of the negative hydrogen ion, Can. J. Phys., 57, 558 (1979).

[7] Atomic Weights of the Elements 2019, International Union of Pure and Applied Chemistry (IUPAC).

[8] Particle Data Group (PDG), Atomic and Nuclear Properties of Materials for more than 350 materials.

[9] Stopping Powers and Ranges for Protons and Alpha Particles, Tech. Rep. ICRU-49, International Commission on Radiation Units and Measurements (1993).