Introduction

lindhard is a clean-room, MIT-licensed Rust engine for Monte Carlo transport of ions in matter, in the binary-collision approximation (BCA). For a beam of ions entering a layered amorphous target it computes where the ions stop, how much damage their cascades make, and what is sputtered and backscattered. An electron engine is planned and not yet started.

This book has two parts.

• Physics manual. One page per model: the equations, the assumptions, the range over which the model is meant to be used, and the papers it comes from. Each page names the input keys and the library types that select the model, so you can go from a line in an input file to its physics and back.
• User guide. Installing, a first run from start to finish (5 keV boron into silicon), the reference for the TOML input, and the output formats.

The source of truth for every equation is the code and its doc comments (cargo doc -p lindhard --open); this manual is written from them and checked against them in CI where that can be automated (every model variant must be named in the manual). The design notes, the validation results and the data provenance table stay in the repository’s docs/ directory:

• docs/validation.md: analytic checks, code-to-code comparisons and comparisons with measurement, with the known deviations.
• docs/data-provenance.md: where every number in the tree comes from, and how far it has been verified.
• docs/architecture.md: the crate layout and the milestone plan.

lindhard is early, pre-release software. The project’s rules (the clean-room license tiers in particular) are in CONTRIBUTING.md.

Conventions and the clean-room position

Units

The library works in SI internally. At the input and output boundary the units are in the key names: energies in eV (_ev), lengths in nm (_nm), angles in degrees (_deg), densities in g/cm³. Electronic stopping cross sections are reported per target atom, in the customary eV·10⁻¹⁵ cm² (to_ev_1e15_cm2 converts from J m²).

Throughout the manual \( Z_1, M_1 \) are the atomic number and mass of the moving particle, \( Z_2, M_2 \) those of the target atom, \( E \) the laboratory kinetic energy, \( N \) the atom density, \( a_0 \) the Bohr radius, \( v_0 = \alpha c \) the Bohr velocity, and \( e^2 \) stands for \( e^2 / 4\pi\varepsilon_0 \). Physical constants are CODATA 2022 (Mohr et al. 2025).

Reduced variables

The collision models use the reduced variables of Lindhard, Scharff and Schiøtt (1963). With a screening length \( a \) (see Screening lengths):

\[ x = \frac{r}{a}, \qquad \beta = \frac{b}{a}, \qquad \varepsilon = \frac{a\, E_\mathrm{cm}}{Z_1 Z_2 e^2}, \qquad E_\mathrm{cm} = E \frac{M_2}{M_1 + M_2}, \]

with \( r \) the separation and \( b \) the impact parameter. In these variables the scattering angle depends only on \( (\varepsilon, \beta) \) and on the screening function.

Determinism

Every model is a pure function of its arguments. Random numbers come from a counter-based stream (ChaCha8) keyed on the run seed and the index of the primary history, so a run gives the same bits on one thread or many. A model page that introduces randomness says which draws it makes and in which order.

The clean-room position

lindhard implements published physics from the papers. It does not use the code or the data of closed or copyleft programs (the tiers are set out in CONTRIBUTING.md). In particular:

• No SRIM stopping tables, and nothing interpolated or fitted from one, and no ICRU stopping tables. Electronic stopping is computed from closed-form models; tabulated stopping enters only as a user-supplied table that carries its own provenance, and is never committed with SRIM- or ICRU-derived numbers.
• No ZBL stopping tables. The ZBL universal screening function (eight published coefficients) and the ZBL reduced nuclear stopping fit are used; the electronic stopping coefficient sets of the same book are not.
• Terms whose only published coefficients are tabulated (the Barkas term, shell corrections, multi-oscillator density-effect parameters, Chu and Yang-O’Connor-Wang straggling) are declined, and the declines are listed in Validity ranges and declined terms.

Every coefficient that enters the code as a number has a row in docs/data-provenance.md saying where it comes from and how far it has been verified. Where a value has only been checked against a secondary source, its page says so.

How a model page is laid out

Each page has the same sections: what the model is, the equations, the assumptions, the validity range, how to select it (the TOML key and the Rust type), its verification status, and its references. New models start from the template book/src/models/_template.md in the repository.

References

• J. Lindhard, M. Scharff, H. E. Schiøtt, Mat. Fys. Medd. Dan. Vid. Selsk. 33 (14) (1963).
• P. J. Mohr, D. B. Newell, B. N. Taylor, E. Tiesinga, Rev. Mod. Phys. 97, 025002 (2025), doi:10.1103/RevModPhys.97.025002 (CODATA 2022).
• D. J. Bernstein, ChaCha, a variant of Salsa20 (2008); J. K. Salmon, M. A. Moraes, R. O. Dror, D. E. Shaw, Proc. SC’11 (2011) (counter-based random streams).

Interatomic potentials

Code: lindhard/src/ion/potential.rs (Screening, Potential).

Model

Two atoms at separation \( r \) interact through a screened Coulomb potential

\[ V(r) = \frac{Z_1 Z_2 e^2}{r}\, \phi\!\left(\frac{r}{a}\right), \]

where \( \phi \) is a universal screening function (\( \phi(0) = 1 \), decreasing to 0) and \( a \) a screening length that depends on \( Z_1, Z_2 \) (Screening lengths). The scattering code works in \( x = r/a \), where \( \phi \) is the only input; the screening length only converts to and from SI.

Three of the four screening functions are sums of exponentials,

\[ \phi(x) = \sum_i c_i\, e^{-b_i x}, \]

and the fourth is a polynomial times an exponential.

The screening functions

The default length is the one each function was introduced or fitted with (Screening::default_length); screening_length overrides it.

ZBL universal

Four exponentials fitted by Ziegler, Biersack and Littmark (1985, ch. 2) to Hartree-Fock-Slater solid-state pair potentials:

These are the commonly printed four-digit rounding of the published values.

Kr-C

Three exponentials fitted to the Hartree-Fock Kr-Kr pair interaction (Wilson, Haggmark and Biersack 1977):

Molière

Molière’s three-exponential approximation to the Thomas-Fermi screening function (Molière 1947):

Lenz-Jensen

A polynomial-times-exponential approximation to the Thomas-Fermi-Jensen statistical model (Lenz 1932; Jensen 1932), in the form printed by Möller (2017, p. 11, eq. (28)):

\[ \phi(x) = \left(1 + q + 0.3344\, q^2 + 0.0485\, q^3 + 0.002647\, q^4\right) e^{-q}, \qquad q = \sqrt{9.67\, x}. \]

Assumptions

• Pair potentials: the interaction of two atoms does not depend on any third atom.
• The screening function is universal: \( Z_1 \) and \( Z_2 \) enter only through the screening length.
• The potentials are purely repulsive. There is no attractive well, so they are not meant for energies comparable with chemical binding (a few eV), where the cutoffs and binding energies of the BCA take over.

Validity

Screened Coulomb potentials describe the repulsive part of the interaction, from the close collisions of keV ions down to the tens of eV of cascade atoms. The ZBL universal function was fitted to pair potentials over a broad set of ion-target pairs and is the usual default. It is not exact for any particular pair: the computed ranges of B in amorphous Si run long against the measurement, and the validation attributes the offset below about 5 keV mainly to the nuclear stopping (the ZBL function being too soft for B on Si) (docs/validation.md, level 3).

Verification status

From docs/data-provenance.md: the ZBL, Kr-C and Molière coefficients have been cross-checked against an independent MIT-licensed implementation (ir2-lab/screened_coulomb, commit f84c3c8), a secondary source; the Lenz-Jensen set agrees with a textbook tabulation (Möller 2017) but has not been checked against the primary papers, which were not accessible.

References

• J. F. Ziegler, J. P. Biersack, U. Littmark, The Stopping and Range of Ions in Solids (Pergamon, New York, 1985), ch. 2.
• W. D. Wilson, L. G. Haggmark, J. P. Biersack, Phys. Rev. B 15, 2458 (1977).
• G. Molière, Z. Naturforsch. A 2, 133 (1947).
• W. Lenz, Z. Phys. 77, 713 (1932), doi:10.1007/BF01342150.
• H. Jensen, Z. Phys. 77, 722 (1932), doi:10.1007/BF01342151.
• W. Möller, Fundamentals of Ion-Solid Interaction, HZDR-073 (Helmholtz-Zentrum Dresden-Rossendorf, 2017), p. 11, eqs. (26)-(30), https://www.hzdr.de/publications/PublDoc-10091.pdf.

Screening lengths

Code: lindhard/src/ion/potential.rs (ScreeningLength).

Model

The screening length \( a(Z_1, Z_2) \) scales the separation in the screening function, \( x = r/a \). Three forms are implemented. Two of them use the Thomas-Fermi constant

\[ C_\mathrm{TF} = \left(\frac{9\pi^2}{128}\right)^{1/3} = 0.8853, \]

in the closed form of Firsov (1958, p. 535), printed as 0.8853 by Lindhard, Scharff and Schiøtt (1963, p. 8).

Without screening_length the length paired with the potential is used: universal for ZBL, Firsov for Kr-C and Molière, Lindhard for Lenz-Jensen. The run echoes the length it used in summary.json (input.physics.screening_length).

Assumptions

• The universal length is the one the ZBL screening function was fitted with; the other two come from Thomas-Fermi scaling arguments for the two-atom system.
• Mixing a screening function with a length other than its own is allowed (it is a model choice), but the fitted functions were fitted with their own length.

Validity

As for the potentials: screened Coulomb collisions from cascade energies upward.

Verification status

From docs/data-provenance.md: the Firsov and Lindhard lengths have been verified against scans of the primary papers (page references above). The universal-length prefactor has been seen only in secondary sources, which agree on the exponent 0.23 and print the prefactor as either 0.8854 or 0.8853 (a 10⁻⁴ relative difference).

References

• O. B. Firsov, Sov. Phys. JETP 6, 534 (1958), pp. 535-536.
• J. Lindhard, M. Scharff, H. E. Schiøtt, Mat. Fys. Medd. Dan. Vid. Selsk. 33 (14) (1963), p. 8.
• J. F. Ziegler, J. P. Biersack, U. Littmark, The Stopping and Range of Ions in Solids (Pergamon, New York, 1985), ch. 2.

The scattering integral

Code: lindhard/src/ion/scattering.rs (theta_quadrature, theta_magic, ScatteringTable).

Model

A binary collision in a central potential is classical elastic scattering. In the reduced variables \( x = r/a \), \( \beta = b/a \), \( \varepsilon = a E_\mathrm{cm} / (Z_1 Z_2 e^2) \), the radial motion in the centre-of-mass frame obeys

\[ G(x) = 1 - \frac{\phi(x)}{\varepsilon x} - \frac{\beta^2}{x^2}, \]

with the distance of closest approach \( x_0 \) the root of \( G(x_0) = 0 \). The centre-of-mass scattering angle is (Goldstein 1980; Ziegler, Biersack and Littmark 1985, ch. 2)

\[ \theta = \pi - 2\beta \int_{x_0}^{\infty} \frac{dx}{x^2 \sqrt{G(x)}}. \]

Distance of closest approach

\( G \) is strictly increasing, so the root is unique and lies in \( [\beta,\, (1/\varepsilon + \sqrt{1/\varepsilon^2 + 4\beta^2})/2] \). Newton’s method runs inside that bracket, and a step that leaves it is replaced by bisection.

Gauss-Mehler quadrature

With \( u = x_0/x \) and then \( u = \cos t \) the integrable endpoint singularity is removed, leaving

\[ \int_0^{\pi/2} \frac{dt}{\sqrt{H(\cos t)}}, \qquad H(u) = \frac{\beta^2}{x_0^2} + \frac{\phi(x_0) - u\, \phi(x_0/u)}{\varepsilon x_0 (1 - u^2)}, \]

which the midpoint rule in \( t \) (the Gauss-Mehler, or Gauss-Chebyshev of the first kind, nodes; Mendenhall and Weller 1991) integrates with exponential convergence. The default is 64 nodes. \( H \) is evaluated in this cancellation-free form.

Precomputed angle table (the transport hot path)

The BCA does not integrate at every collision. Once per run and per \( (Z_1, Z_2, \text{potential}) \) pair it builds a table of \( y = \ln\tan(\theta/2) \) on a grid uniform in \( \ln\varepsilon \) and \( \ln\beta \), by quadrature, and interpolates bilinearly in those coordinates. \( y \) is close to linear in \( \ln\beta \) at both ends, which is what makes bilinear interpolation accurate, and an error \( \delta y \) in the table is at most \( \delta y \) radians in \( \theta \). At build time the interpolation error is measured against direct quadrature at cell centres on a deterministic subset of cells; the grid and the largest absolute and relative errors found are written to summary.json under physics.scattering_table. These are empirical bounds from sampling, not proofs. The run grid is \( 10^{-6} \le \varepsilon \le 10^{4} \), \( 10^{-5} \le \beta \le 10^{2} \), 32 points per decade (lindhard::input::TABLE_SPEC). Beyond the largest tabulated \( \beta \) the angle is taken as 0 (the constant-free-path \( p_\mathrm{max} \) is about 20 screening lengths or less at solid densities); below the smallest it is clamped; and at a reduced energy outside the table the engine falls back to direct quadrature (deterministic, only slower).

Magic formula (cross-check only)

The Biersack-Haggmark “magic formula” is kept as a fast approximate cross-check, not used in transport:

\[ \cos\frac{\theta}{2} = \frac{\beta + \rho + \Delta}{x_0 + \rho}, \qquad \Delta = \frac{A (x_0 - \beta)}{1 + G}, \]

\[ A = 2\alpha\varepsilon\beta^{b}, \quad \alpha = 1 + C_1 \varepsilon^{-1/2}, \quad b = \frac{C_2 + \varepsilon^{1/2}}{C_3 + \varepsilon^{1/2}}, \quad G = \frac{\gamma}{\sqrt{1 + A^2} - A}, \quad \gamma = \frac{C_4 + \varepsilon}{C_5 + \varepsilon}, \]

with \( \rho \) the radius of curvature of the trajectory at closest approach. Constant sets exist for the ZBL universal function (Ziegler, Biersack and Littmark 1985) and for Molière (Biersack and Haggmark 1980).

Nuclear stopping

The reduced nuclear stopping cross section follows from the angle,

\[ s_n(\varepsilon) = 2\varepsilon \int_0^\infty \sin^2\frac{\theta}{2}\, \beta\, d\beta, \]

integrated over \( \ln\beta \) by Simpson’s rule. It is used in tests and validation, not in transport (the BCA samples the collisions themselves). Two references are checked against it: the ZBL universal fit \( s_n = \ln(1 + 1.1383\varepsilon) / [2(\varepsilon + 0.01321\varepsilon^{0.21226} + 0.19593\varepsilon^{0.5})] \) for \( \varepsilon \le 30 \) (Ziegler, Biersack and Littmark 1985), and the high-energy limit for a sum of exponentials from the first-order impulse approximation of Lindhard, Nielsen and Scharff (1968, eq. (3.4)).

Assumptions

• Classical, elastic, binary collisions in a central potential; no inelastic energy loss inside the collision itself (electronic loss is handled separately, see Electronic stopping).
• The angle depends only on \( (\varepsilon, \beta) \) and the screening function; \( Z_1, Z_2 \) and the screening length enter through the reduced variables.

Validity

The classical treatment holds when the collision is well localised compared with the screening length, which is the case for heavy particles from cascade energies (a few eV) to well above the MeV range.

Verification status

The quadrature is checked against the small-angle perturbation formula of Lindhard, Nielsen and Scharff (1968) and against the ZBL nuclear stopping fit; the magic formula is checked against the quadrature for both constant sets. The ZBL magic-formula constants were cross-checked against ir2-lab/screened_coulomb (MIT); the Molière set is unverified against the paper and is checked only against the quadrature. See docs/data-provenance.md.

References

• H. Goldstein, Classical Mechanics, 2nd ed. (Addison-Wesley, 1980).
• J. F. Ziegler, J. P. Biersack, U. Littmark, The Stopping and Range of Ions in Solids (Pergamon, New York, 1985), ch. 2.
• M. H. Mendenhall and R. A. Weller, Nucl. Instrum. Methods B 58, 11 (1991).
• J. P. Biersack and L. G. Haggmark, Nucl. Instrum. Methods 174, 257 (1980), doi:10.1016/0029-554X(80)90440-1.
• J. Lindhard, V. Nielsen, M. Scharff, Mat. Fys. Medd. Dan. Vid. Selsk. 36 (10) (1968), eqs. (3.3)-(3.4).

Electronic stopping

Code: lindhard/src/ion/stopping/.

A moving atom loses energy to the target electrons as well as to the nuclei. In the BCA the electronic loss is handled separately from the collisions: continuously along each free flight (nonlocal), at the collisions as a function of the distance of closest approach (local), or as a mix of the two (BCA transport, “Electronic loss”).

Conventions

• Every model returns the stopping cross section per target atom, \( S_e = -\frac{1}{N}\frac{dE}{dx} \), in J m² internally and eV·10⁻¹⁵ cm² at the boundary.
• Energies are the laboratory kinetic energy of the moving atom.
• Every model is a pure function of (ion, target \( Z_2 \), energy): no random numbers, no state.
• Each model exposes an advisory validity range (ElectronicStopping::validity). The CLI warns on stderr when the beam energy lies outside it; it does not refuse the run.
• Compounds and mixtures combine the element cross sections by Bragg additivity.

Choosing a model

A [stopping] table replaces the chosen model for the one (ion, target element) pair it declares (User stopping tables). Energy-loss straggling is described on its own page. It is a library function that the BCA engine does not call at present: the electronic loss along each flight is deterministic, \( N S_e(E)\, s \).

The validity ranges side by side, and the terms that are deliberately not implemented, are on Validity ranges and declined terms.

Lindhard-Scharff

Code: lindhard/src/ion/stopping/lindhard_scharff.rs (LindhardScharff).

Model

At ion velocities below \( v_0 Z_1^{2/3} \) the electronic stopping is proportional to the velocity (Lindhard and Scharff 1961). In the reduced units of Lindhard, Scharff and Schiøtt (1963):

\[ a = 0.8853\, a_0 \left(Z_1^{2/3} + Z_2^{2/3}\right)^{-1/2}, \qquad \varepsilon = \frac{E\, a\, M_2}{Z_1 Z_2 e^2 (M_1 + M_2)}, \qquad \rho = N x\, 4\pi a^2 \frac{M_1 M_2}{(M_1 + M_2)^2}, \]

the reduced electronic stopping is

\[ \left(\frac{d\varepsilon}{d\rho}\right)_e = k_L\, \varepsilon^{1/2}, \qquad k_L = \xi_e\, \frac{0.0793\, Z_1^{1/2} Z_2^{1/2} (A_1 + A_2)^{3/2}} {\left(Z_1^{2/3} + Z_2^{2/3}\right)^{3/4} A_1^{3/2} A_2^{1/2}}, \qquad \xi_e = Z_1^{1/6}, \]

with \( A_1, A_2 \) the masses in u. The equivalent dimensional form is

\[ S_e = \xi_e\, \frac{8\pi e^2 a_0\, Z_1 Z_2}{\left(Z_1^{2/3} + Z_2^{2/3}\right)^{3/2}}\, \frac{v}{v_0}. \]

The factor \( \xi_e \approx Z_1^{1/6} \) (order 1 to 2) is part of the Lindhard-Scharff result. A test checks the reduced form, with the rounded 0.0793, against the dimensional form to 1 % for all tested ion-target pairs.

A per-element multiplicative correction \( f(Z_2) \) can be attached in the library (LindhardScharff::with_correction, default 1; it must be finite and non-negative, and 0 switches the element’s stopping off). The CLI does not expose it; use a user table for measured stopping instead.

Selecting it

stopping = "lindhard-scharff" in [physics] (the default), StoppingChoice::LindhardScharff. All of the loss is applied nonlocally, along the free flights. It is also the nonlocal half of stopping = "equipartition-ls-or" (Oen-Robinson).

Assumptions

• A free-electron-gas picture with Thomas-Fermi scaling; the stopping is smooth in \( Z_1 \) and \( Z_2 \) and has no shell structure (the measured \( Z_1 \) oscillations are not reproduced).
• No dependence on the chemical or physical state of the target.

Validity

\( v < v_0 Z_1^{2/3} \), that is \( E/A_1 \) below about \( 25\ \mathrm{keV} \cdot Z_1^{4/3} \). This is the regime of keV-range implantation and of cascade atoms. Above it the stopping peaks and turns over, which this model does not describe; see Bethe-Bloch for high velocities and Validity ranges and declined terms for why no interpolation joins the two.

The validation pages report a known offset: for B in Si the computed ranges run long against the measurement, and at 10 to 20 keV part of the offset is attributed to LS stopping being too small for that pair (docs/validation.md, level 3).

Verification status

The constants 0.8853, 0.0793 and \( \xi_e = Z_1^{1/6} \) are checked against the dimensional form in a unit test (1 %); they have not yet been checked digit by digit against the papers (docs/data-provenance.md).

References

• J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961).
• J. Lindhard, M. Scharff, H. E. Schiøtt, Mat. Fys. Medd. Dan. Vid. Selsk. 33 (14) (1963).

Oen-Robinson and the equipartition mix

Code: lindhard/src/ion/stopping/oen_robinson.rs (OenRobinson), lindhard/src/ion/stopping/mix.rs (EquipartitionMix).

Model

Oen and Robinson (1976) make the electronic loss local: it is taken at each collision and depends on how close the two atoms came. A collision with distance of closest approach \( r_\mathrm{min} \) loses

\[ \Delta E_e(E, r_\mathrm{min}) = S_\mathrm{LS}(E)\, \frac{c^2}{2\pi a^2}\, \exp\!\left(-\frac{c\, r_\mathrm{min}}{a}\right), \qquad c = 0.3, \]

with \( S_\mathrm{LS} \) the Lindhard-Scharff cross section and \( a \) the Firsov screening length \( 0.8853\, a_0 (Z_1^{1/2} + Z_2^{1/2})^{-2/3} \). The prefactor normalises \( \int \Delta E_e\, 2\pi p\, dp = S_\mathrm{LS} \) when \( r_\mathrm{min} \) is identified with the impact parameter \( p \); that identification is the approximation of the original paper, so the impact-averaged stopping equals the LS value.

Equipartition mix

The equipartition mix splits the Lindhard-Scharff loss in two equal halves: half nonlocal, continuous along the free flight, and half local, at each collision, by the Oen-Robinson formula evaluated at the distance of closest approach. This combination is the one described in the BCA literature that builds on Oen and Robinson (1976). Averaged over impact parameter the mix equals the LS stopping (a unit test checks this to 10⁻¹⁴).

Selecting it

stopping = "equipartition-ls-or" in [physics], which is StoppingChoice::EquipartitionLsOr in the input and ElectronicLoss::EquipartitionLsOr in the engine (see BCA transport). The stopping model passed to the engine is then not used: this mode carries its own LS and Oen-Robinson losses, so it cannot be combined with [stopping] tables. The local losses are reported in the energy budget as electronic_local.

Assumptions

• \( r_\mathrm{min} \approx p \) in the normalisation (Oen and Robinson 1976).
• The local part is only sampled at the collisions the BCA makes, out to \( p_\mathrm{max} \) (to \( p_\mathrm{max}\sqrt{K + 1} \) with \( K \) weak collisions). Where that radius is not large compared with \( a / 0.3 \), the decay length of the local loss, the total electronic stopping comes out somewhat below the LS value.

Validity

The same as Lindhard-Scharff: velocities below \( v_0 Z_1^{2/3} \).

Verification status

The constants \( c = 0.3 \) and the choice of the Firsov length were entered from the contributor’s recollection of the paper and are not verified against it; spot-check them before relying on the local loss (docs/data-provenance.md). The averaging identity with LS is tested.

References

• O. S. Oen and M. T. Robinson, Nucl. Instrum. Methods 132, 647 (1976).
• O. B. Firsov, Sov. Phys. JETP 6, 534 (1958) (screening length).
• J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961).

Bethe-Bloch and the density effect

Code: lindhard/src/ion/stopping/bethe.rs (BetheBloch, EffectiveCharge, density_effect_single_oscillator).

Model

At high velocity the electronic stopping cross section per target atom is (Bethe 1930, 1932; Fano 1963; the same equation as the Particle Data Group review of the passage of particles through matter):

\[ S = \frac{4\pi (e^2)^2 z^2 Z_2}{m c^2 \beta^2} \left[ \frac{1}{2}\ln\frac{2 m c^2 \beta^2\gamma^2 W_\mathrm{max}}{I^2}

• \beta^2 - \frac{\delta}{2} - \frac{C}{Z_2} + L_1 + L_2 \right], \]

\[ W_\mathrm{max} = \frac{2 m c^2 \beta^2\gamma^2}{1 + 2\gamma m/M + (m/M)^2}, \]

with \( m \) the electron mass, \( M \) the ion mass, \( z \) the projectile effective charge, \( I \) the mean excitation energy, \( \delta \) the density-effect correction and \( C \) the shell correction. \( \beta^2 \) is computed as \( \tau(\tau + 2)/(1 + \tau)^2 \) with \( \tau = E/Mc^2 \), which stays accurate at low energy.

• Bloch term (Bloch 1933), on by default: \( L_2 = -y^2 \sum_{n \ge 1} \frac{1}{n(n^2 + y^2)} \), \( y = z\alpha/\beta \).
• Mean excitation energy: by default the Bloch rule \( I = 10\ \mathrm{eV} \cdot Z_2 \) (Bloch 1933), a rough estimate (about 20 % in \( I \), about 2 % in \( S \)). The library takes a measured, cited value through BetheBloch::with_mean_excitation_ev.
• Shell correction \( C/Z_2 \) and density effect \( \delta \): 0 by default; the library accepts constants from a cited source (shell_over_z, density_delta).
• Barkas term \( L_1 \): not implemented (see Validity ranges and declined terms).

Effective charge

The CLI uses the bare charge; the Barkas form is available in the library and is not applied by default.

Single-oscillator density effect (opt-in)

The dielectric formulation of the density effect (Fermi 1940; Sternheimer 1952), in the form given by Fano (1963), is

\[ \delta = \sum_i f_i \ln\!\left(1 + \frac{L^2}{\omega_i^2}\right) - \frac{L^2 (1 - \beta^2)}{\beta^2\omega_p^2}, \qquad \sum_i \frac{f_i\, \omega_p^2}{\omega_i^2 + L^2} = \frac{1}{\beta^2} - 1. \]

Specialised to one oscillator (\( f = 1 \), \( \omega_0 = I/\hbar \)) the constraint solves in closed form, \( L^2 = \omega_p^2(\beta\gamma)^2 - \omega_0^2 \), and

\[ \delta = \ln\!\left[(\beta\gamma)^2 \left(\frac{\hbar\omega_p}{I}\right)^2\right] - 1

• \frac{(I/\hbar\omega_p)^2}{(\beta\gamma)^2} \quad \text{for } \beta\gamma > \frac{I}{\hbar\omega_p}, \]

and 0 below, with \( \hbar\omega_p = \hbar\sqrt{n e^2/\varepsilon_0 m} \) the free-electron plasma energy. The one-oscillator reduction is ours, not quoted from a paper. It is enabled with BetheBloch::with_density_effect (library only).

Selecting it

stopping = "bethe-bloch" in [physics], StoppingChoice::BetheBloch: bare charge, Bloch term on, Bloch-rule \( I \), no shell or density correction, all loss nonlocal.

Use it only for problems that stay above its range. Where the bracket is not positive the model returns NotApplicable, and the transport stops with that error rather than silently switching model. For any ion that slows down in the target (a semi-infinite substrate always) the run fails once the energy falls below the bracket’s zero: for example, 300 keV H into Si with stopping = "bethe-bloch" stops with “bethe-bloch is not applicable at 85306 eV”. No low-to-high energy interpolation is implemented, for the reasons given in Validity ranges and declined terms.

Assumptions

• First Born approximation plus the Bloch correction; the projectile charge is fixed by the effective-charge model.
• Without shell and Barkas corrections, expect errors of a few percent at 1 to 10 MeV/u.

Validity

\( v \ge 3 v_0 Z_1^{2/3} \) up to 1 GeV/u (the upper end because the density effect is off by default). The single-oscillator density effect switches on only at \( \beta\gamma > I/\hbar\omega_p \), about 5 to 6 for Si (a proton of about 5 GeV), so it is zero over the whole advisory range. Real materials switch the effect on much earlier (around \( \beta\gamma \approx 1.5 \) for Si), because their oscillator spectrum is spread: the one-oscillator form underestimates \( \delta \) through the transition and is only the correct asymptote for ultra-relativistic ions.

Verification status

Hand-calculated proton-in-Si values (1, 10 and 100 MeV, \( I = 140 \) eV, CODATA constants) are unit tests; they are our own arithmetic, not taken from any stopping table. The density-effect form is tested through its limits (\( \delta \to 0 \) continuously at threshold, and \( \delta \to 2\ln(\beta\gamma\, \hbar\omega_p / I) - 1 \) at large \( \beta\gamma \)). The Bloch rule, the Bloch series and the Barkas effective charge have not been checked against the original papers (docs/data-provenance.md).

References

• H. Bethe, Ann. Phys. 397, 325 (1930); H. Bethe, Z. Phys. 76, 293 (1932).
• F. Bloch, Ann. Phys. 408, 285 (1933).
• U. Fano, Ann. Rev. Nucl. Sci. 13, 1 (1963).
• Particle Data Group, “Passage of particles through matter”, in the Review of Particle Physics.
• W. H. Barkas, Nuclear Research Emulsions I (Academic Press, 1963).
• L. C. Northcliffe, Ann. Rev. Nucl. Sci. 13, 67 (1963).
• E. Fermi, Phys. Rev. 57, 485 (1940).
• R. M. Sternheimer, Phys. Rev. 88, 851 (1952).

User stopping tables

Code: lindhard/src/ion/stopping/table.rs (StoppingTable, TableOverride).

Model

A user table gives \( S_e(E) \) for one projectile (atomic number and mass) in one target element, as data the user supplies. Between the tabulated points the cross section is interpolated piecewise linearly in \( \ln S \) against \( \ln E \), which keeps the monotonicity of each segment of the data. Nothing is extrapolated or clamped: a query outside the table’s energy range is an error (OutOfTableRange) that stops the run.

provenance = "Author, Journal vol, page (year), Table N"   # required
ion_z = 5
ion_mass_amu = 11.0093            # optional; default: standard atomic weight
target_z = 14
energy_ev = [1.0e2, 1.0e3, 1.0e4]                 # strictly increasing, eV
stopping_ev_1e15_cm2 = [10.0, 30.0, 60.0]         # eV 1e-15 cm^2 per atom

(The numbers are a format illustration, not data.)

Provenance is mandatory

A table without a non-empty provenance string is rejected at load time: data without an origin is not admitted. The run records each table in summary.json (physics.stopping_tables: the path as written, the resolved path, the SHA-256 of the file’s bytes, the provenance string, the pair and the energy range), and lists it under physics.models as user-table with its provenance as the source.

Do not load SRIM- or ICRU-derived tables into anything committed to a repository. A table is the user’s own data and its terms are the user’s concern; this project’s own tree never contains such tables.

Selecting it

Declare the files under [stopping] tables (see TOML input reference). A table replaces the [physics] stopping model for exactly the (ion_z, target_z) pair it declares, including recoils of that species when recoils are followed; every other pair uses the model. A table cannot be combined with stopping = "equipartition-ls-or", which carries its own Lindhard-Scharff and Oen-Robinson loss.

Assumptions

• The table is tied to the projectile mass it was declared for. A query for a mass that differs by more than a small relative tolerance (tight enough to separate neighbouring isotopes) is an error (TableMassMismatch): the same energy at a different mass is a different speed, and no energy-axis conversion is attempted.
• Compounds are built from element tables by Bragg additivity; a compound table is not a supported input.

Validity

Exactly the table’s energy range. The CLI checks up front that the range contains the beam energy, and for a recoil species that it starts at or below recoil_cutoff_ev; it warns if a table starts above primary_cutoff_ev, because the run fails if a projectile slows below it.

Verification status

The loader, the interpolation and the range and mass checks are unit-tested. The tests use tables generated at test time from our own Lindhard-Scharff model; no table is committed (docs/data-provenance.md).

References

The table’s own provenance string is its reference. The format and the rules are this project’s.

Bragg additivity

Code: lindhard/src/ion/stopping/bragg.rs (bragg_cross_section_per_atom, CompoundCorrection).

Model

The stopping cross section of a compound or mixture, per average atom, is the atom-fraction-weighted sum of the element cross sections (Bragg and Kleeman 1905):

\[ S_\mathrm{compound} = f \sum_j x_j S_j, \qquad -\frac{dE}{dx} = N f \sum_j x_j S_j, \]

with \( x_j \) the atom fractions, \( S_j \) the element cross sections from the chosen model (or a user table for that pair), \( N \) the total atom density and \( f \) a per-compound correction factor.

Selecting it

Always on: every layer’s electronic loss goes through the Bragg sum, which for a pure element is just that element’s cross section. summary.json lists it under physics.models as bragg-additivity. The CLI applies no correction (\( f = 1 \)). The library hook CompoundCorrection (a constant or energy-dependent factor, finite and non-negative) is there for chemical or phase corrections, such as cores-and-bonds schemes, supplied by the caller from a cited source.

Assumptions

• Each atom stops the ion independently of its chemical environment.

Validity

Wherever the element models apply. Bragg additivity ignores the chemical and physical state of the target, and no correction for it is applied unless one is supplied.

Verification status

The sum and the rejection of invalid correction factors are unit-tested.

References

• W. H. Bragg and R. Kleeman, Phil. Mag. 10, 318 (1905).

Energy-loss straggling

Code: lindhard/src/ion/stopping/straggling.rs (StragglingModel, variance_per_atom).

Model

The energy loss of an ion along a path fluctuates. For a thick target in the free-electron picture, the variance of the energy loss per unit areal density is (Bohr 1948)

\[ \Omega_B^2 = 4\pi z^2 Z_2 (e^2)^2 \quad \text{(J}^2\,\text{m}^2 \text{ per target atom)}, \]

independent of energy in the non-relativistic regime, with \( z \) the projectile charge from an effective-charge model. The relativistic variant multiplies it by \( (1 - \beta^2/2)/(1 - \beta^2) \) (Bethe and Livingston 1937; Fano 1963). For a compound the variances add by atom fraction, \( \sum_j x_j \Omega_j^2 \).

Selecting it

Library only (variance_per_atom, variance_per_atom_material). The BCA engine does not use it at present and there is no TOML key: the electronic loss along each free flight is deterministic.

Assumptions

• Free, stationary target electrons; every electron of the target contributes, whatever its binding.

Validity

High energy, \( v \gg v_0 Z_1^{2/3} \). Below about 1 MeV/u Bohr’s value is an overestimate, because the bound electrons contribute less than free ones. The corrections that reduce it at intermediate energy (Chu; Yang-O’Connor-Wang) are declined, and the Lindhard-Scharff low-velocity correction is omitted (not verified); see Validity ranges and declined terms.

Verification status

Not verified against the original papers (docs/data-provenance.md).

References

• N. Bohr, Mat. Fys. Medd. Dan. Vid. Selsk. 18 (8) (1948).
• H. A. Bethe and M. S. Livingston, Rev. Mod. Phys. 9, 245 (1937).
• U. Fano, Ann. Rev. Nucl. Sci. 13, 1 (1963).
• W. K. Chu, Phys. Rev. A 13, 2057 (1976) (declined correction).
• Q. Yang, D. J. O’Connor, Z. Wang, Nucl. Instrum. Methods B 61, 149 (1991) (declined correction).

Validity ranges and declined terms

This page is docs/stopping-models.md, included verbatim so there is one copy. It collects the validity ranges of the electronic stopping models side by side, and records the terms that are deliberately not implemented and why (the clean-room rules exclude coefficient sets that are only published as tables). The model pages are Lindhard-Scharff, Oen-Robinson, Bethe-Bloch, user tables, Bragg additivity and straggling; the references cited below are listed in full on those pages and at the end of this one.

Code: lindhard/src/ion/stopping/. All models return the stopping cross section per atom in J m² (to_ev_1e15_cm2 converts to eV·10⁻¹⁵ cm²) and expose the ranges below at run time through ElectronicStopping::validity. Ranges are advisory; the models do not refuse energies outside them (Bethe-Bloch returns NotApplicable where its bracket is not positive, and user tables return OutOfTableRange). v0 = α c is the Bohr velocity, Z1 the ion charge number.

What is deliberately not implemented

The clean-room rules (CONTRIBUTING.md) rule out tabulated coefficient sets from SRIM/ZBL, ICRU reports and similar. As a result:

• Barkas (L1) term: declined (see “Literature findings” below).
• Shell correction C/Z2: declined. Callers may set a constant from a cited source (BetheBloch::shell_over_z); default 0.
• Density effect δ: only the single-oscillator form is implemented (opt-in, BetheBloch::with_density_effect); the tabulated multi-oscillator parameter sets are declined. density_delta can still be set as a constant.
• Mean excitation energy I: defaults to the Bloch rule 10 eV · Z2 (rough); supply a cited measured value with with_mean_excitation_ev.
• Biersack-Varelas interpolation joining the low- and high-energy regimes: documented but not implemented (issue #44, closed; see “Biersack-Varelas interpolation” below). Fitted ZBL/Biersack-Varelas coefficients are not used.
• Straggling: Chu and Yang-O’Connor-Wang corrections declined (see below) and the Lindhard-Scharff low-velocity correction (not verified) is omitted; Bohr plus an optional relativistic factor only.
• Heavy-ion effective charge: the Barkas empirical form is available in bethe and is not applied by default.

Biersack-Varelas interpolation (issue #44, closed)

Status: documented, not implemented.

The joining form

The interpolation joins a low-energy and a high-energy stopping branch harmonically:

1/S = 1/S_low + 1/S_high     (equivalently S = S_low S_high / (S_low + S_high))

Two open-access sources were read (rendered pages) and show this form:

• M. V. Moro, PhD thesis, Universidade de Sao Paulo (2017), doi:10.11606/T.43.2017.tde-18092017-095345, eq. (3.7), printed p. 40. It gives s_low = A1 E^0.45 and s_high = (A2/E) ln(1 + A3/E + A4 E) (E in reduced units as defined there), attributes the form to Varelas and Biersack (1970), refined by Andersen and Ziegler (1977), and says A1 to A4 are fitted to experimental data.
• P. de Vera et al., arXiv:2608.15368 (2026), eq. (9), p. 3: the same s = s_low s_high / (s_low + s_high), citing Varelas and Biersack, and Andersen and Ziegler.

NISTIR 4999 (Berger 1992), section 3.5.2, p. 8, also mentions the “fitting formula of Varelas and Biersack (1970)” used with ICRU 49 coefficients (no equation given).

The original papers were not consulted (closed access; Unpaywall and Semantic Scholar report no open copy, ScienceDirect returned 403):

• C. Varelas and J. P. Biersack, Nucl. Instrum. Methods 79, 213 (1970), doi:10.1016/0029-554X(70)90141-2.
• J. P. Biersack and L. G. Haggmark, Nucl. Instrum. Methods 174, 257 (1980), doi:10.1016/0029-554X(80)90440-1.

Citation correction: “Biersack and Varelas, NIM 194, 93 (1982)” is not a Biersack-Varelas paper. Crossref resolves NIM 194, 93-100 (1982) to Biersack and Ziegler, “Refined universal potentials in atomic collisions” (doi:10.1016/0029-554X(82)90496-7). It is not a source for this join.

Why it is not implemented

The published s_high = (A2/E) ln(1 + A3/E + A4 E) is a fully fitted functional form (de Vera et al. eqs. (10)-(11), following ICRU 49; A1 to A4 are all fitting parameters) built so that it never crosses zero. For positive coefficients the 1 + keeps the logarithm’s argument above 1, so s_high > 0 at every E whatever A3 is. The harmonic join tends to s_low at low energy even with A3 = 0: s_high then tends to the finite value A2 A4 while s_low = A1 E^0.45 goes to 0. The fitted A3/E term only sets how fast s_high grows as E goes to 0. The coefficients come from the Andersen-Ziegler / ICRU 49 fits, which the clean-room policy (CONTRIBUTING.md) does not allow.

lindhard’s only sourced high-energy branch is BetheBloch, which has no such structure: its bracket crosses zero inside the crossover region. As S_high goes to 0+, the join S_low S_high / (S_low + S_high) goes to 0, not to S_low. Evaluated with LindhardScharff and BetheBloch defaults (Bloch term on, I = 10 eV Z2), in units of 1e-15 eV cm^2/atom (computed by a Python mirror of the formulas, not by the Rust models):

Bethe-Bloch is negative for H in Si below roughly 90 keV, for P in Si up to about 400 keV/u, and for He in Au up to about 450 keV/u. Neither workaround is acceptable. Falling back to S_low below the Bethe zero makes S jump from about 0 to the full LS value (about 27 for H in Si near 90 keV), which breaks continuity. Propagating NotApplicable leaves the model undefined over the whole low-energy regime, so ions slowing down through it cannot be transported. Restricting the join to where Bethe-Bloch is positive violates the low-energy limit and gives badly wrong stopping near the zero. A fix would need a high branch that stays positive at low E. The options are the fitted A1 to A4 (Tier C data), or a regularisation of our own, such as grafting the ln(1 + ...) structure onto Bethe by identifying A2 and A4 with Bethe quantities and setting A3 = 0. That identification is our own construction, which neither source makes. Neither option is admissible.

When it could be revisited

• A published, coefficient-free high-energy branch that stays positive below the Bethe zero (for example a Bethe form with a closed-form, non-tabulated shell correction; the review of the omitted correction terms under #41, now closed, is recorded below) is found and can be cited to an equation.
• The A1 to A4 coefficients become available from a source whose licence permits use, or are derived from our own fits to data that may be used.
• The primary papers are read and show a high branch that stays positive without fitted coefficients, or that ties the ln(1 + ...) form to Bethe quantities in a way that can be cited.

Literature findings (issue #41, closed)

Each of the four omitted items was reviewed for a closed form whose coefficients come from a paper itself. The review was from the contributor’s knowledge of the literature; no paper text was available to check digits against, so every statement below about a paper’s content is “as recalled” and the declines are on that basis. Nothing was taken from ICRU/SRIM/NIST tables.

• Barkas L1: declined. Ashley, Ritchie & Brandt (Phys. Rev. B 5, 2393 (1972)) give L1 = F(b / x^(1/2)) / (Z2^(1/2) x^(3/2)) with F a function that is tabulated (and later fitted to it) rather than closed-form. Lindhard’s 1976 treatment of the Barkas effect (NIM 132, 1) gives asymptotic forms from a harmonic-oscillator model, whose numerical coefficients we could not reproduce with confidence from memory. Jackson & McCarthy (Phys. Rev. B 6, 4131 (1972)) and later fits carry fitted coefficients that cannot be verified here. No form with verifiable coefficients exists to us; an invented one is not acceptable.
• Shell correction: declined. Walske and Bichsel shell corrections come from hydrogenic-shell calculations published as tables in η = βγ and I; the compact empirical forms (e.g. the ICRU 49 expression) are fits with coefficients tuned to data and tables, which Tier C of CONTRIBUTING.md excludes. A first-principles computation (hydrogenic shell sums) is a research task of its own, not a closed form.
• Density effect: partly implemented. The Sternheimer parameter sets (x0, x1, a, m, C̄) are tables, which we do not ingest. What is admissible is the dielectric formulation itself, computed from a model; we implement its one-oscillator specialisation (see density_effect_single_oscillator). Its high-energy limit is exact, but it switches on too late to matter in range. A multi-oscillator version needs per-material oscillator strengths, which are tabulated material data. Callers wanting a realistic δ in the range 1 to 100 GeV/u should supply one from a cited source via density_delta.
• Chu / Yang-O’Connor-Wang straggling: declined. Chu’s correction (Phys. Rev. A 13, 2057 (1976)) is built from Hartree-Fock-Slater charge densities and is published as tables; Yang, O’Connor & Wang (NIM B 61, 149 (1991)) is a fit whose coefficients are fitted to data tables. Neither has a closed form with paper-sourced coefficients that we could verify.

Deviations from the issue text

• The Lindhard-Scharff k_L quoted in the issue omits the factor ξ_e = Z1^(1/6). It is included here, because with it the reduced form k_L ε^(1/2) agrees with the dimensional LS expression to better than 1 % for all tested ion/target pairs (and without it Z1 > 1 disagrees by Z1^(1/6)).
• The Oen-Robinson local-loss constants are flagged unverified in data-provenance.md.

References

• J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961).
• O. S. Oen and M. T. Robinson, Nucl. Instrum. Methods 132, 647 (1976).
• H. Bethe, Ann. Phys. 397, 325 (1930); H. Bethe, Z. Phys. 76, 293 (1932).
• F. Bloch, Ann. Phys. 408, 285 (1933).
• U. Fano, Ann. Rev. Nucl. Sci. 13, 1 (1963).
• W. H. Bragg and R. Kleeman, Phil. Mag. 10, 318 (1905).
• E. Fermi, Phys. Rev. 57, 485 (1940).
• R. M. Sternheimer, Phys. Rev. 88, 851 (1952).
• N. Bohr, Mat. Fys. Medd. Dan. Vid. Selsk. 18 (8) (1948).
• M. V. Moro, PhD thesis, Universidade de São Paulo (2017), doi:10.11606/T.43.2017.tde-18092017-095345.
• P. de Vera et al., arXiv:2608.15368 (2026).
• M. J. Berger, NISTIR 4999 (1992).
• C. Varelas and J. P. Biersack, Nucl. Instrum. Methods 79, 213 (1970), doi:10.1016/0029-554X(70)90141-2 (not consulted).
• J. P. Biersack and L. G. Haggmark, Nucl. Instrum. Methods 174, 257 (1980), doi:10.1016/0029-554X(80)90440-1 (not consulted for this join).
• J. P. Biersack and J. F. Ziegler, Nucl. Instrum. Methods 194, 93 (1982), doi:10.1016/0029-554X(82)90496-7.
• J. C. Ashley, R. H. Ritchie, W. Brandt, Phys. Rev. B 5, 2393 (1972).
• J. Lindhard, Nucl. Instrum. Methods 132, 1 (1976).
• J. D. Jackson and R. L. McCarthy, Phys. Rev. B 6, 4131 (1972).
• W. K. Chu, Phys. Rev. A 13, 2057 (1976).
• Q. Yang, D. J. O’Connor, Z. Wang, Nucl. Instrum. Methods B 61, 149 (1991).

BCA transport and free-flight conventions

Code: lindhard/src/ion/bca/ (Bca, BcaConfig, MeanFreePath, ElectronicLoss; kinematics.rs).

Model

The target is a stack of homogeneous, structureless (amorphous) layers, the first starting at depth \( x = 0 \); a semi-infinite substrate may close it. A moving atom alternates straight free flights and binary elastic collisions with target atoms, losing energy to electrons along each flight. This is the amorphous-target BCA of Biersack and Haggmark (1980), within the general BCA framework of Robinson and Torrens (1974) and the treatment of Eckstein (1991). One history:

1. The primary enters at the origin of the front face with the beam’s direction.
2. Free flight of length \( \tau\lambda \) (\( \lambda \) and the distribution of \( \tau \) depend on the free-path convention below), with the nonlocal electronic loss \( N S_e(E)\, s \) taken at the energy at the start of the segment (with Bragg additivity over the layer composition).
3. Collision with a partner drawn by stoichiometry, at impact parameter \( p = p_\mathrm{max}\sqrt{R} \) (uniform over a disc) and a uniform azimuth. The centre-of-mass angle comes from the scattering table.
4. Recoil. If the transfer \( T \) exceeds the partner’s displacement energy \( E_d \), the atom is displaced with energy \( T - E_b \) (the lattice binding \( E_b \) stays in the lattice); otherwise \( T \) stays in the lattice at the site (Biersack and Haggmark 1980; Eckstein 1991). Displaced atoms are followed in turn as full cascades.
5. A particle stops when its energy falls below its cutoff, or leaves through the front or back face if it can overcome the surface barrier; otherwise it is reflected specularly back into the target.

Kinematics

Classical elastic two-body kinematics (Goldstein 1980; Robinson and Torrens 1974), with \( \mu = M_1/M_2 \):

\[ \tan\psi = \frac{\sin\theta}{\cos\theta + \mu}, \qquad \varphi = \frac{\pi - \theta}{2}, \qquad T = \gamma E \sin^2\frac{\theta}{2}, \qquad \gamma = \frac{4 M_1 M_2}{(M_1 + M_2)^2}, \]

where \( \psi \) is the projectile’s laboratory deflection and \( \varphi \) the recoil’s, on the opposite azimuth.

Surface barrier

A planar barrier (Sigmund 1969; Eckstein 1991): a particle of energy \( E \) reaching a face with direction cosine \( c \) to the normal escapes if \( E c^2 > E_s \). Outside, its energy is \( E - E_s \), the parallel momentum is unchanged, and the normal component satisfies \( E’ c’^2 = E c^2 - E_s \) (refraction away from the normal). Target atoms use the surface binding energy \( E_s \) of the material at the face; the beam species uses primary_surface_binding_ev (default 0).

Free-path conventions

Constant

The flight length is the mean interatomic distance \( l = N^{-1/3} \) and impact parameters are uniform over a disc of radius \( p_\mathrm{max} = (\pi N^{2/3})^{-1/2} \), so \( N\pi p_\mathrm{max}^2 l = 1 \): each flight sweeps exactly one atom’s worth of target (Biersack and Haggmark 1980). So that every primary does not make its first collision at the same depth, the primary’s first flight is \( R\, l \) with \( R \) uniform in \( [0, 1) \); recoils start at an atom site and fly a full \( l \). Collisions with \( p > p_\mathrm{max} \) are dropped unless weak collisions add them.

Energy-dependent

Collisions that deflect by less than a minimum centre-of-mass angle \( \theta_\mathrm{min} \) (min_cm_angle_deg, required with this convention and only with it) are neglected. For each element \( i \), \( p_i(E) \) is the impact parameter at which the angle equals \( \theta_\mathrm{min} \), capped at the constant-convention \( p_\mathrm{max} \). The mean free path is

\[ \lambda = \frac{1}{N\pi \sum_i x_i p_i^2}, \]

flight lengths are exponential with that mean, the partner is drawn with probability \( x_i p_i^2 / \sum_j x_j p_j^2 \), and \( p = p_i\sqrt{R} \). This is the standard cross-section cut-off treatment of a Poisson collision process (Eckstein 1991). It reduces to the constant convention, with exponential flight lengths, when every \( p_i \) hits the cap. The nuclear loss of the neglected small-angle collisions is dropped, so choose \( \theta_\mathrm{min} \) small.

Layer boundaries

A flight is drawn as a dimensionless number of mean free paths \( \tau \). When it reaches an interface it is truncated there, the electronic loss for the truncated length is applied, the material is switched, and the flight continues in the new layer with the unused part \( \tau - s/\lambda \) (it is not redrawn). For exponential paths this is exact by memorylessness; for the constant path it means that splitting one layer into two of the same material changes nothing but rounding (a test checks this).

Weak collisions (optional)

With the constant free path every flight ends in one collision with \( p \le p_\mathrm{max} \), and the nuclear loss of collisions beyond \( p_\mathrm{max} \) is dropped while the electronic loss of the flight is charged in full. At cascade-tail energies that dropped part is large, so the electronic share of a cascade comes out too high. weak_collisions = K (0 to 3) adds the weak collisions of Möller and Eckstein (1988): before the hard collision, \( K \) collisions with partners at

\[ p_k = p_\mathrm{max}\sqrt{k + R}, \qquad k = 1, \dots, K, \]

each uniform over an annulus of area \( \pi p_\mathrm{max}^2 \), with its own partner and azimuth. A weak collision deflects the particle and takes its transfer, but never makes a recoil (its \( T \) stays in the lattice). A weak collision whose partner would lie in front of the front surface is skipped. Each collision is evaluated at the energy left after the previous one, which differs from the report, where all collisions of a step use the starting energy (see the ion::bca module docs for the measured effect). Weak collisions are available with the constant free path only. Default 0, because there is no single published convention and the option changes sputter yields substantially (see docs/validation.md).

Electronic loss

The local half is evaluated at the distance of closest approach of every collision, weak ones included; see Oen-Robinson.

Cutoffs and energy bookkeeping

primary_cutoff_ev and recoil_cutoff_ev are required: they have no defensible default. A recoil cutoff above the surface binding energies suppresses sputtering, so keep it below the smallest \( E_s \) when sputtering matters (Eckstein 1991). With follow_recoils = false a displaced atom stops where it was created.

Every energy change is subtracted from the particle and added to exactly one field of the energy budget (results.energy_budget_ev_per_ion), so each history conserves energy up to rounding; the largest per-history relative residual is reported. Electronic losses larger than the particle’s energy are clamped to it.

Randomness and determinism

Each primary has its own random stream, keyed on the run seed and the primary’s index; its recoils draw from the same stream, and cascades are followed from an explicit stack in a fixed order. Histories are grouped in chunks of a fixed size (64, independent of the thread count), which fixes the summation order of the tallies. The output is therefore bit-identical on any number of threads (lindhard/tests/determinism.rs).

Assumptions

• Amorphous, structureless targets: no crystal structure (no channeling).
• The target does not change with fluence (no dynamic composition).
• No refraction of the incident beam at the entrance surface (negligible at keV energies).
• Binary collisions only; many-body effects at low energy are represented only through the cutoffs, \( E_d \), \( E_b \) and \( E_s \).

Validity

The BCA is meant for energies well above the binding energies of the target. It follows cascade atoms down to a few eV, where its results depend on the cutoffs and binding energies chosen. Comparisons with other codes and with measured ranges and sputter yields, with the known deviations, are in docs/validation.md.

References

• J. P. Biersack and L. G. Haggmark, Nucl. Instrum. Methods 174, 257 (1980), doi:10.1016/0029-554X(80)90440-1.
• M. T. Robinson and I. M. Torrens, Phys. Rev. B 9, 5008 (1974).
• W. Eckstein, Computer Simulation of Ion-Solid Interactions (Springer, Berlin, 1991).
• H. Goldstein, Classical Mechanics, 2nd ed. (Addison-Wesley, 1980), sec. 3.11.
• P. Sigmund, Phys. Rev. 184, 383 (1969).
• W. Möller and W. Eckstein, TRIDYN - Binary collision simulation of atomic collisions and dynamic composition changes in solids, report IPP 9/64, Max-Planck-Institut für Plasmaphysik, Garching (1988), p. 14, eq. (26).
• O. S. Oen and M. T. Robinson, Nucl. Instrum. Methods 132, 647 (1976).

Displacement damage

Code: lindhard/src/ion/damage.rs (models), lindhard/src/tally/ (counts).

lindhard reports two different quantities side by side and never adds one to the other or substitutes one for the other:

• Model estimates of the number of displacements made by a primary knock-on atom (PKA) of a given energy: NRT and Kinchin-Pease, from the Lindhard partition of the PKA energy (results.damage.nrt).
• Event counts from the full-cascade simulation: vacancies, interstitials and replacements counted collision by collision as the BCA follows every recoil (results.damage.cascade).

Damage energy (Lindhard partition)

Of a PKA’s energy \( T \), the part eventually given to atomic motion is (Lindhard, Nielsen, Scharff and Thomsen 1963)

\[ T_\mathrm{dam} = \frac{T}{1 + k\, g(\varepsilon)}, \qquad g(\varepsilon) = 3.4008\, \varepsilon^{1/6} + 0.40244\, \varepsilon^{3/4} + \varepsilon, \]

with the analytic fit \( g \) of Robinson (1970) as adopted by the NRT standard, and

\[ \varepsilon = \frac{T a A_2}{Z_1 Z_2 e^2 (A_1 + A_2)}, \qquad a = 0.8853\, a_0 \left(Z_1^{2/3} + Z_2^{2/3}\right)^{-1/2}, \]

\[ k = \frac{0.0793\, Z_1^{2/3} Z_2^{1/2} (A_1 + A_2)^{3/2}} {\left(Z_1^{2/3} + Z_2^{2/3}\right)^{3/4} A_1^{3/2} A_2^{1/2}}, \]

the Lindhard reduced energy and the Lindhard-Scharff electronic stopping coefficient of the PKA (\( Z_1, A_1 \)) in the target (\( Z_2, A_2 \)). For a self-ion these reduce to the forms written in the NRT standard, \( \varepsilon = T / (86.931\, Z^{7/3}) \) (\( T \) in eV) and \( k = 0.1337\, Z^{1/6} (Z/A)^{1/2} \); a test checks this.

NRT and Kinchin-Pease

Norgett, Robinson and Torrens (1975):

\[ N_\mathrm{NRT} = \begin{cases} 0 & T_\mathrm{dam} < E_d, \\ 1 & E_d \le T_\mathrm{dam} < 2E_d/0.8, \\ 0.8\, T_\mathrm{dam} / (2 E_d) & T_\mathrm{dam} \ge 2E_d/0.8, \end{cases} \]

with \( E_d \) the displacement threshold and 0.8 the displacement efficiency the authors took from BCA simulations. The older Kinchin and Pease (1955) count is 0, 1, or \( T/(2E_d) \) with thresholds \( E_d \) and \( 2E_d \); the tallies evaluate it on the damage energy (the usual “modified” form).

NRT is not a cascade count. Stoller et al. (2013) discuss how BCA full-cascade vacancy counts and NRT estimates differ, and recommend that a displacement dose comparable with the NRT standard be computed from the damage energy with the NRT formula, not taken from the full-cascade vacancy count. NRT also overestimates the number of defects that survive in-cascade recombination; it is a standard exposure unit, not a prediction of surviving defects.

Cascade counts

During transport a target atom is displaced when it receives \( T > E_d \) (BCA transport, step 4). The counts are defined in lindhard::tally::CascadeDefects:

• Replacement: a moving atom comes to rest, immediately after the collision in which it displaced an atom of its own element, at that atom’s site, so it fills the site. A particle stops when its energy falls below its cutoff, so this count depends on the cutoffs (a recoil cutoff near \( E_d \) gives the most replacements); with follow_recoils = false it is always 0.
• Vacancies = displacements minus replacements.
• Interstitials: recoils that came to rest in the target and did not fill a site. Implanted beam particles are not counted here (they are in the range tally).

These are kept in total, per layer and as a depth profile (damage_profile.csv).

Selecting it

Always on: both are computed in every run. The inputs are the per-element energies: \( E_d \) (e_d_ev), \( E_b \) (e_b_ev) and \( E_s \) (e_s_ev), set per material or for an element in every layer with [physics.energies.<symbol>]. There are no model variants to choose.

Assumptions

• Compounds (this crate’s own convention). NRT and the Lindhard partition are defined for a monatomic target. For a layer with several elements the tally uses the PKA’s own \( Z_1, A_1 \), the atom-fraction-weighted mean \( Z_2, A_2 \) of the layer the PKA starts in (non-integer in general), and the PKA element’s own \( E_d \) in that layer. No published source is claimed for this averaging; it reduces exactly to standard NRT for a monatomic layer. Treat NRT numbers for compounds as an exposure index comparable within lindhard, not as a value following any compound-target standard.
• The cascade counts depend directly on \( E_d \), \( E_b \) and the recoil cutoff, and on the amorphous-target BCA itself (no recombination, no crystal structure).

Validity

NRT is a convention for damage accounting, defined for monatomic targets. The cascade counts are BCA counts of displacement events, not of surviving defects.

Verification status

The self-ion reductions of \( \varepsilon \) and \( k \) are checked in a test. The default \( E_d \) values of the element table follow the ASTM E521 convention and are not yet verified against the standard (docs/data-provenance.md); Si has no default, so an input must choose one.

References

• J. Lindhard, V. Nielsen, M. Scharff, P. V. Thomsen, “Integral equations governing radiation effects”, Mat. Fys. Medd. Dan. Vid. Selsk. 33 (10) (1963).
• M. T. Robinson, in Nuclear Fusion Reactors (British Nuclear Energy Society, London, 1970), p. 364.
• M. J. Norgett, M. T. Robinson, I. M. Torrens, Nucl. Eng. Des. 33, 50 (1975).
• G. H. Kinchin and R. S. Pease, Rep. Prog. Phys. 18, 1 (1955).
• R. E. Stoller, M. B. Toloczko, G. S. Was, A. G. Certain, S. Dwaraknath, F. A. Garner, Nucl. Instrum. Methods B 310, 75 (2013).

Dynamic composition (bookkeeping only)

Code: lindhard/src/ion/dynamic.rs (CompositionGrid, Relaxation).

Status. This is the start of milestone M3 (targets that change with fluence). Only the bookkeeping exists: a grid of slabs whose composition callers change with explicit inventory deltas. The fluence loop, and the adapter that turns transport tallies into deltas, are not written yet, so no CLI run uses it and there is no TOML key.

Model

The finite layers of a target are held as slabs. Each slab stores an areal inventory \( A_i \) (atoms/m²) per element, kept sorted by \( Z \) so every derived quantity is deterministic; a slab built from a layer of atom density \( n_i \) and thickness \( t \) starts with \( A_i = n_i t \). An optional semi-infinite substrate backs the grid and never changes.

Inventories change only through CompositionGrid::apply. A retained primary adds one atom of its species where it came to rest; a recoil subtracts an atom where it was created and adds it where it stops (an escaping recoil is a loss). The update either succeeds completely or leaves the grid untouched.

Volume relaxation

After every update each slab’s thickness is recomputed from its inventory by one of two stated conventions:

For an element in its own solid, atomic_volume_from_density gives \( v = M / (N_A \rho) \). After relaxation the slab boundaries are rebuilt from the front surface at \( x = 0 \). Swelling moves the interior interfaces and the back face. With sputter erosion on, the sputtered atoms of each element are removed from the front slabs (the loss they would otherwise cause in the slab where they were displaced is cancelled, so nothing is removed twice), the thickness they occupied under the chosen convention is the recession of that step, and the grid is re-anchored so the surface is again \( x = 0 \); a depth plus the cumulative recession is the depth in the original frame. Erosion is off by default and then changes nothing.

Assumptions

• Ideal mixing is additivity of atomic volumes (a Vegard-type rule applied to atoms); it ignores chemistry, voids and amorphisation swelling, and a compound’s real density generally differs from its prediction.
• A fixed number density reproduces a known compound density but cannot respond to composition, and one value serves every slab.
• Neither is an equation of state with a pressure or phase model from the literature; both are this crate’s stated conventions, and the module claims no published source for them.

Validity

Bookkeeping for compositions that stay close to the phase whose volumes or density the caller supplies. The caller must give the volume of every species that occurs (none is inferred, in particular not for gases).

Verification status

Conservation, transactional updates and the relaxation conventions are covered by lindhard/tests/dynamic.rs. The fluence stepping loop (DynamicRun, with the tally-to-delta adapter) is covered by lindhard/tests/dynamic_run.rs: the low-fluence limit equals the static engine, refining the step size converges, steps continue one global random stream, results are bit-identical on 1, 2 and 8 threads, adaptive rejection consumes no indices, and inventory follows the event conventions. There is no comparison with a measured dynamic profile yet.

The fluence loop

DynamicRun delivers the beam’s primaries in steps. A step runs n primaries on the current target, scales the integer atom counts of the events by the fluence one primary stands for, applies them to the grid and relaxes it. The target is held fixed within a step, so the step must be small enough that the composition does not change much inside it: the adaptive policy bounds the largest relative change of a slab per step and retries a too-large step with fewer ions from the same first primary. The front surface stays at x = 0; the time series reports the interface depths and total thickness measured from it, not a receding surface. The step-size policy is this crate’s own design, not a published scheme. In the CLI it is the [dynamic] table (docs/cli.md).

References

No published source is claimed for the two relaxation conventions (see above). The atomic volume uses the Avogadro constant of CODATA 2022 (P. J. Mohr, D. B. Newell, B. N. Taylor, E. Tiesinga, Rev. Mod. Phys. 97, 025002 (2025)) and the element densities of lindhard/src/elements.rs, whose sources are in docs/data-provenance.md.

Installing

lindhard is pre-release and is not yet published as a package; build it from source. It is pure Rust with no C or Fortran dependencies in the default build.

Prerequisites

• A stable Rust toolchain. The repository’s rust-toolchain.toml selects the stable channel, so rustup installs the right one on first use.
• git.

Build

git clone https://github.com/2AMLogic/lindhard.git
cd lindhard
cargo build --release -p lindhard-cli

The binary is target/release/lindhard. Build in release mode for real runs; a debug build is many times slower. You can run it in place, through cargo, or install it into cargo’s binary directory:

target/release/lindhard --version
cargo run --release -p lindhard-cli -- --version
cargo install --path lindhard-cli      # puts `lindhard` on your PATH

--version prints the crate version and the git describe of the source it was built from, for example lindhard 0.0.1 (c2da637). The same two values are written into every summary.json, so a result can be traced to the binary that produced it.

Commands

lindhard check input.toml                  # parse and validate, no transport
lindhard run input.toml --out dir/         # run, write dir/summary.json and CSVs
lindhard run input.toml --out dir/ --ions 200 --seed 7 --threads 4
lindhard --version

--ions and --seed override run.ions and run.seed, and the output records the override. --threads overrides run.threads and never changes the results. The first run walks through both commands.

The library

The physics is in the lindhard library crate; the command line is a thin front end on it. To browse the library API:

cargo doc -p lindhard --open

First run: 5 keV B into Si

This walkthrough runs the example examples/b_5keV_si.toml: 10 000 boron ions of 5 keV into amorphous silicon, 7° off the surface normal. Every block of command output below is taken from an actual run of that file: the blocks are generated by validation/book_walkthrough.py, and CI reruns the example and fails if they no longer match. Run the commands from the root of the repository, with lindhard built as in Installing (or replace lindhard with cargo run --release -p lindhard-cli --).

1. The input

# 5 keV boron into amorphous silicon, 7 degrees off normal.
#
# Run:   lindhard run examples/b_5keV_si.toml --out out/b_5keV_si
# Check: lindhard check examples/b_5keV_si.toml
#
# Schema: docs/cli.md. Energies are in eV, lengths in nm, angles in degrees.

[beam]
ion = "B"
energy_ev = 5000.0
tilt_deg = 7.0

[target]
# An element symbol names the pure element at its tabulated density
# (lindhard/src/elements.rs).
substrate = "Si"

[physics]
potential = "zbl"
stopping = "lindhard-scharff"
free_path = "constant"
primary_cutoff_ev = 5.0
recoil_cutoff_ev = 2.0

# Si has no default displacement energy in the element table, so the input
# must choose one. 15 eV is an illustrative model parameter, not a
# recommendation: choose E_d for your problem.
[physics.energies.Si]
e_d_ev = 15.0

[run]
ions = 10000
seed = 1

[tally]
depth_bin_nm = 0.5
depth_bins = 200

# Fit a dual-Pearson profile to the depth histogram (off by default).
dual_pearson = true

Section by section:

• [beam]: the projectile by element symbol (its mass defaults to the standard atomic weight), its energy in eV and its tilt from the surface normal in degrees.
• [target]: substrate = "Si" names pure silicon at its tabulated density as a semi-infinite substrate. There are no finite layers, so nothing can be transmitted.
• [physics]: the model choices, spelled out even where they are the defaults: the ZBL universal potential, Lindhard-Scharff electronic stopping, and the constant free path. The two cutoffs are required: the beam ion is dropped below 5 eV, recoils below 2 eV (kept below the surface binding energy of Si, so that sputtering is not cut off).
• [physics.energies.Si]: the displacement energy \( E_d \) of Si. The element table has no default for Si, so the input must choose one; 15 eV here is an illustrative value, not a recommendation. Si’s surface binding energy \( E_s \) comes from the element table and its lattice binding \( E_b \) defaults to 0.
• [run]: the number of histories and the seed. Together with the input and the binary, the seed fixes every bit of the output.
• [tally]: a depth grid of 200 bins of 0.5 nm, and a dual-Pearson fit of the depth profile, which is off by default.

The full list of keys is in the TOML input reference.

2. Check the input

lindhard check examples/b_5keV_si.toml

examples/b_5keV_si.toml: OK
  beam: B at 5000 eV, tilt 7 deg, azimuth 0 deg; 10000 ions, seed 1
  layer 0: Si (Si 1.0000), semi-infinite, 4.9940e22 atoms/cm^3
  transport: amorphous-bca
  screening function: zbl-universal
  screening length: universal
  scattering angle: gauss-mehler-quadrature-table
  electronic stopping: lindhard-scharff
  compound stopping: bragg-additivity
  free path: constant
  displacement criterion: e_d-e_b
  surface barrier: planar
  random numbers: chacha8-per-history-stream

check parses and validates the input without transporting anything. It prints the resolved beam and target (the atom density comes from the tabulated mass density of Si) and every model the run will use. An invalid input exits non-zero with a message naming the offending key, for example target.layers[0].thickness_nm: -5 nm must be finite and positive.

3. Run it

lindhard run examples/b_5keV_si.toml --out out/b_5keV_si

10000 ions: 9441 stopped, 559 backscattered, 0 transmitted, 2534 sputtered atoms; wrote out/b_5keV_si

The one-line summary counts the fates of the 10 000 primaries and the target atoms that left through the front face. The run takes a few seconds in a release build. Add --threads N to choose the number of worker threads; the results are bit-identical whatever N is.

The output directory holds:

damage_profile.csv
depth_profile.csv
escape_spectra.csv
ions.csv
lateral_profile.csv
summary.json

4. Read the summary

summary.json starts with the format and software version, then the input exactly as it was run (defaults filled in, so the run can be reproduced from its own header), then the physics and the results. A few parts of it follow; the Output files page describes every key.

The models used. physics.models names every model with its citation (the citations are omitted here):

{
  "models": [
    {
      "role": "transport",
      "name": "amorphous-bca"
    },
    {
      "role": "screening function",
      "name": "zbl-universal"
    },
    {
      "role": "screening length",
      "name": "universal"
    },
    {
      "role": "scattering angle",
      "name": "gauss-mehler-quadrature-table"
    },
    {
      "role": "electronic stopping",
      "name": "lindhard-scharff"
    },
    {
      "role": "compound stopping",
      "name": "bragg-additivity"
    },
    {
      "role": "free path",
      "name": "constant"
    },
    {
      "role": "displacement criterion",
      "name": "e_d-e_b"
    },
    {
      "role": "surface barrier",
      "name": "planar"
    },
    {
      "role": "random numbers",
      "name": "chacha8-per-history-stream"
    }
  ]
}

Where the primaries went. results.primaries:

{
  "primaries": {
    "stopped": 9441,
    "backscattered": 559,
    "transmitted": 0,
    "stopped_depth_mean_nm": 24.389957255951394,
    "stopped_depth_std_nm": 12.728934527954353
  }
}

A few percent of the boron ions backscatter; the rest stop in the silicon.

The range distribution. results.range has the moments of the depth at which the stopped primaries came to rest, with their standard errors. mean_nm is the projected range \( R_p \) and std_dev_nm the straggle \( \Delta R_p \); the kurtosis is 3 for a Gaussian.

{
  "depth": {
    "n": 9441,
    "mean_nm": 24.389957255951394,
    "std_dev_nm": 12.728934527954353,
    "skewness": 0.3718620249214205,
    "kurtosis": 2.713342344381711,
    "mean_std_err_nm": 0.1310035467755083,
    "std_dev_std_err_nm": 0.08573835216188522,
    "skewness_std_err": 0.02024523266795846,
    "kurtosis_std_err": 0.04670779974658438
  },
  "pearson_iv": null,
  "pearson_iv_error": "moments outside the Pearson IV region: skewness = 0.3718620249214205, kurtosis = 2.713342344381711 (needs kurtosis > Some(3.2610741629133053))"
}

When the moments lie outside the region of the Pearson IV family, as they do here, pearson_iv is null and pearson_iv_error gives the reason instead of a silent fallback. The dual-Pearson fit requested in [tally] (results.range.dual_pearson, shortened here) splits the profile into a near-surface head and a tail:

{
  "head_fraction": 0.19614973489188495,
  "head": {
    "mean_nm": 14.200851032315631,
    "std_dev_nm": 10.282737064212109
  },
  "tail": {
    "mean_nm": 26.299554956251598,
    "std_dev_nm": 12.449831675138565
  },
  "chi_square": 142.04173850533618
}

docs/validation.md (level 3) compares the projected range of this problem with a SIMS measurement of B in amorphous Si. With these default models the computed range runs long, and the page discusses why.

Damage. results.damage.per_ion puts the two kinds of damage figure side by side (Displacement damage): the NRT and Kinchin-Pease estimates from the damage energy of the primary knock-on atoms, and the vacancies, interstitials and replacements counted in the simulated cascades. They are different quantities and are not expected to agree.

{
  "per_ion": {
    "pka": 21.3841,
    "pka_energy_ev": 3375.943722511579,
    "damage_energy_ev": 2727.1474995691347,
    "nrt_displacements": 75.1806883174286,
    "kinchin_pease_displacements": 92.1413018627793,
    "displacements": 124.9028,
    "replacements": 8.3302,
    "vacancies": 116.5726,
    "interstitials": 116.3192
  }
}

Sputtering and backscatter.

{
  "sputtering": {
    "yield_per_ion": 0.2534
  },
  "escapes": {
    "backscatter_coefficient": 0.0559,
    "transmission_coefficient": 0.0,
    "energy_reflection_coefficient": 0.009828116927787605
  }
}

Where the energy went. Every history’s energy is accounted for, per ion: electronic loss, energy left in the lattice (sub-threshold transfers; \( E_b = 0 \) here), work against the surface barrier, energy carried out by backscattered ions and sputtered atoms, and the kinetic energy of particles when they fell below their cutoff (rest). The largest relative bookkeeping residual of any history is at the level of rounding.

{
  "energy_budget_ev_per_ion": {
    "incident": 5000.0,
    "electronic_nonlocal": 2172.168730014593,
    "electronic_local": 0.0,
    "lattice": 2628.475261767412,
    "surface_barrier": 1.1732420000000003,
    "backscattered": 49.14058463893802,
    "sputtered": 14.90822822514642,
    "transmitted": 0.0,
    "rest": 134.13395335390996,
    "max_relative_residual": 3.819877747446298e-15
  }
}

5. The profiles

The CSV files hold the profiles on the grids set in [tally]. The first rows of depth_profile.csv, the stopped primaries per 0.5 nm bin and that count per incident ion per nm:

depth_lo_nm,depth_hi_nm,stopped_primaries,fraction_per_nm
0.0,0.5,13,0.0026
0.5,1.0,31,0.0062
1.0,1.5,38,0.0076
1.5,2.0,36,0.0072
2.0,2.5,50,0.01
2.5,3.0,46,0.0092
3.0,3.5,68,0.0136

The last row has an upper edge of inf and collects everything deeper than the grid, so nothing is dropped. damage_profile.csv gives the cascade defects on the same grid, lateral_profile.csv the lateral and radial spread of the stopped primaries, escape_spectra.csv the energy and angle spectra of everything that left the target, and ions.csv the final state of every primary. Their columns are described in Output files.

6. Change something

Some edits to try, each a one-line change to the input:

• potential = "kr-c" and screening_length = "lindhard": another potential and screening length (Interatomic potentials).
• stopping = "equipartition-ls-or": half of the electronic loss taken locally at the collisions (Oen-Robinson); electronic_local in the energy budget becomes non-zero.
• weak_collisions = 3 under [physics]: weak collisions beyond \( p_\mathrm{max} \) (BCA transport).
• A finite layer in front of the substrate, as in examples/as_50keV_si_sio2.toml.

lindhard check shows the effect of each edit on the model list before you spend time on a run.

TOML input reference

An input file is TOML with these tables: [beam], [materials] (optional), [target], [physics], [stopping] (optional), [run] and [tally] (optional). An electron run has an [electron] table instead of [beam], [physics] and [tally] (section “Electron runs” below). The physics behind each [physics] choice is in the physics manual: potentials and screening lengths, electronic stopping, the free-path conventions and weak collisions, and user stopping tables. The first run walks through a complete file, and the repository’s examples/ directory has more.

This page is the input section of docs/cli.md, included here so there is one copy.

Units are in the key names: energies in eV (_ev), lengths in nm (_nm), angles in degrees (_deg), densities in g/cm³. Every table rejects unknown keys. The schema types are lindhard::input (shared with future front ends).

[beam]

[materials.<name>]

Named materials, in the form of lindhard::material::MaterialSpec:

[materials.SiO2]
density_g_cm3 = 2.2                    # required for compounds
elements = [
  { symbol = "Si", atom_fraction = 1.0 },
  { symbol = "O", atom_fraction = 2.0, e_d_ev = 20.0, e_s_ev = 2.0 },
]

Each element takes symbol or z, atom_fraction or mass_fraction (relative weights, normalised), and optional e_d_ev, e_b_ev, e_s_ev. Defaults for the energies come from the element table where one is tabulated; see lindhard/src/material.rs.

[target]

[target]
substrate = "Si"            # optional semi-infinite substrate

[[target.layers]]           # finite layers, front to back
material = "SiO2"
thickness_nm = 10.0

A material (in a layer or as the substrate) is either a key of [materials], an element symbol (the pure element at its tabulated density), or an inline material table with the same keys as [materials.<name>]. A [materials] key wins over an element symbol of the same name. Without a substrate the target has a back face and particles can be transmitted.

[physics]

[physics.energies.<symbol>] sets e_d_ev, e_b_ev and/or e_s_ev for that element in every layer that contains it, after (so overriding) the material’s own values. An element with no tabulated default and no value set is an error naming the layer and the key to set.

tuning (phenomenological calibration, not a published model choice). "none" or omission leaves the physics and the echoed input exactly as without the key. A named set scales the surface binding energy E_s by a per-element factor fitted to measured data. Rules: the factor multiplies the resolved E_s (an explicit [physics.energies] or material value if given, else the elemental default), once per layer, after overrides; the global element table and the collision algorithm are untouched. The pilot supports static ion runs on single-element layers only, with a beam species the set was fitted for and target elements the set lists; compounds, [dynamic] targets, other beams, unlisted elements and unknown set names are rejected with a physics.tuning error. A beam energy outside the set’s fitted range, or a tilted beam, runs with a warning (an extrapolation). summary.json then has physics.tuning with the set, its version and provenance, and per layer the original E_s, the factor and the effective E_s. Tuned results must be reported next to, never in place of, untuned ones.

Shipped sets (fit record and held-out scores: docs/data-provenance.md, “Tuning factor sets”; a new version ships under a new name):

The factors are a calibration of yields under those settings; other settings (potential, E_d, cutoffs, weak collisions) were not part of the fit, and the set does not claim better accuracy for them.

[stopping]

Optional. Supplies user stopping tables for the electronic stopping of particular (ion, target element) pairs.

[stopping]
tables = ["tables/b_in_si.toml", "tables/p_in_si.toml"]

Absent, the input means what it always meant (format.version is unchanged and nothing is echoed).

Table file. The lindhard::ion::stopping::table::StoppingTable format:

provenance = "Author, Journal vol, page (year), Table N"   # required
ion_z = 5
ion_mass_amu = 11.0093            # optional; default: standard atomic weight
target_z = 14
energy_ev = [1.0e2, 1.0e3, 1.0e4]                 # strictly increasing, eV
stopping_ev_1e15_cm2 = [10.0, 30.0, 60.0]         # eV 1e-15 cm^2 per atom

Interpolation is piecewise linear in ln S versus ln E. provenance is mandatory (data without an origin is not admitted). Do not use SRIM- or ICRU-derived tables in anything committed to a repository; a table is the user’s own data and its terms are the user’s concern.

Paths. A relative path resolves against the directory of the input file (not the current directory). The echoed input keeps the path as written.

Composition with [physics] stopping. A table replaces the [physics] stopping model for exactly the pair it declares (ion_z, target_z), including recoils of that species when follow_recoils is on. Every other pair uses the [physics] stopping model. A pair with a table is never silently served by the model: a query outside the table’s energy range, or for a different ion mass, is an error that stops the run. Nothing is extrapolated. Declare each pair once. A table for a pair that cannot occur in the run (including a table for a target element’s recoils when follow_recoils = false) is accepted with a warning that it is unused.

Recoil species. Tables are keyed by (ion_z, target_z) only, so with follow_recoils = true a table whose ion_z is a target element also serves every recoil of that element. Recoils carry the standard atomic weight and are followed down to physics.recoil_cutoff_ev, so such a table is checked up front against both: its ion_mass_amu must be the standard weight, and it must start at or below recoil_cutoff_ev; either failure is an error. A consequence is that an isotopic self-ion beam (e.g. beam.mass_amu = 27.9769 for Si into Si) cannot take a table for its own pair while recoils are followed: drop the table, use the standard weight, or set follow_recoils = false. A recoil-species table that ends below the largest energy the beam can transfer to that element warns. Tables cannot be combined with stopping = "equipartition-ls-or" (that mode carries its own Lindhard-Scharff/Oen-Robinson loss and would ignore them).

Errors name the field (stopping.tables[0]): an unknown key in [stopping], a missing or unreadable file, invalid table contents (including a missing provenance), a duplicate pair, a table whose ion mass differs from the beam ion’s, a table whose range does not contain the beam energy, and the recoil-species checks above. A table that starts above physics.primary_cutoff_ev warns, because the run fails if a projectile slows below it.

Provenance in the output. Tables are user data, so the run records them. summary.json has physics.stopping_tables, one entry per table: path (as written), resolved_path (absolute where possible), sha256 of the file’s bytes, the table’s provenance string, ion_z, ion_mass_amu, target_z and the energy range. physics.models lists each as user-table with the path and provenance as its source. The key is absent without [stopping].

[dynamic] (optional)

Makes the run fluence-dependent: the target composition is updated as the fluence builds up (sputter erosion, build-up of implanted atoms). Without the table nothing changes: the run, its output files and its bytes are those of a static run. The model, its conventions and its limits are in lindhard::ion::dynamic and the book chapter on dynamic composition.

run.ions is the number of histories of the whole run and fluence_cm2 the fluence they represent, so each ion stands for fluence_cm2 / run.ions ions/cm². The ions are delivered in steps; after each step the grid is updated from that step’s events (a recoil is subtracted where it is created and added where it stops, a stopped beam ion is added, an escaped atom is a loss; the substrate is an immutable reservoir) and relaxed. Primary i of the run always uses the random stream (seed, i), whatever the step sizes and thread count.

With erosion = false the front surface stays at x = 0: swelling moves the interior interfaces and the back face of the slabs, not the front surface (the surface_nm column is that fixed frame, always 0). With erosion = true the lost thickness is removed from the front and the grid is re-anchored so the current surface is again x = 0; surface_nm is then the cumulative recession R in nm and a depth x in the output is x + R in the original frame. The recession of a step is the volume of the removed atoms per area under the chosen relaxation (sum Z removed_Z v_Z, or sum removed / n), and removal equals the sputtered counts per element, so no atom is created or lost. If a step sputters more of an element than the slabs hold, the excess is not removed and the element is counted in the clamped column. With a substrate, atoms sputtered from the substrate remove nothing from the slabs, so the recession falls short of Y F / n once the film is thin. dynamic_summary.json totals gain recession_nm only with erosion on. Depths in the output are measured from the front surface of that step. A dynamic run needs E_d for every element that can occur, including the beam species. The Python bindings run static inputs only.

[run]

[tally]

The depth grid (depth_bin_nm, depth_bins, from the front face) is shared by the range histogram, the dual-Pearson fit and the defect profiles. Particles outside any grid are counted in explicit underflow and overflow entries, never dropped. Every count and per-ion value is for the same incident ions.

Electron runs ([electron])

An input with an [electron] table runs the low-energy electron engine (lindhard::electron::transport) instead of the ion BCA, with the full electron tally (lindhard::tally::FullElectronTally). It exposes what the library does and adds no physics; every choice and every data provenance is written to the output, so a result can be reproduced from its own header. The schema types are lindhard::input::electron. [materials] and [target] are the ion run’s tables (above); [beam], [physics], [stopping], [tally] and [dynamic] are not accepted. Example: ../examples/electron/e_10keV_si.toml.

[electron.beam]
energy_ev = 10000.0

[electron.transport]
cutoff_ev = 1.0
cutoff_reference = "vacuum-level"
secondaries = "kieft-bosch"
boundary = "step-barrier"

[electron.elastic]
potential = "thomas-fermi-yukawa"

[electron.inelastic]
model = "penn-single-pole"

[electron.materials.Si]
optical_elf = "si_elf.toml"
band = { kind = "insulator", valence_band_width_ev = 10.0, band_gap_ev = 2.0, affinity_ev = 3.0, provenance = "..." }

[target]
substrate = "Si"

[run]
histories = 1000
seed = 1

[electron.beam]

Primaries start on the front face at y = z = 0 (just outside it with the step barrier).

[electron.transport] (lindhard::electron::transport::TransportConfig)

[electron.elastic]

The corrections are solved per grid energy with the stand-in’s own Poisson density (AtomicElastic::compute_corrected); the elastic table’s model and provenance strings name them and every polarizability with its source.

[electron.inelastic]

[electron.tables]: one log-spaced energy grid shared by the elastic and inelastic tables of every material.

The transport holds the rates of the first and last rows beyond the grid; a grid that starts above the lowest stopping threshold or ends below the largest possible energy warns.

[electron.materials.<name>]: the electron data of each material the target uses, keyed by the name the target gives it (a [materials] key or an element symbol; an inline target material is an error in an electron run). Every name the target uses needs an entry; an unused entry warns.

No optical or band data of any real material is committed (data-provenance.md); the data files are the user’s, and their terms are the user’s concern. Subshell binding-energy tables (SubshellBindingTable) have no key yet: the transport loop does not use inner-shell channels, so there is nothing to feed them to (see “extending” below).

[electron.tally] (lindhard::tally::ElectronTallyConfig)

[run] of an electron run: histories (required), seed (required) and threads (all cores; not echoed, never changes results).

Output files

lindhard run input.toml --out dir/ writes summary.json and up to five CSV profiles into dir/. The first run shows excerpts of each from a real run. The models behind the damage figures are on Displacement damage.

This page is the output section of docs/cli.md, included here so there is one copy.

summary.json

depth_profile.csv

depth_lo_nm,depth_hi_nm,stopped_primaries,fraction_per_nm: stopped primaries per depth bin, and that count per incident ion per nm. The last row’s upper edge is inf (overflow) and its density is empty.

lateral_profile.csv

quantity,lo_nm,hi_nm,count,per_ion_per_nm: stopped primaries by y, z and radial distance (quantity is y, z or radial), for the grid set by lateral_bin_nm and lateral_bins. Each quantity ends with an underflow row (lo_nm = -inf) and an overflow row (hi_nm = inf) with empty density.

damage_profile.csv

depth_lo_nm,depth_hi_nm,vacancies,interstitials,replacements, then the same three per incident ion per nm: the cascade defect counts by depth, on the depth grid. vacancies is displacements minus replacements. The last row, with depth_hi_nm = inf, holds everything deeper than the grid (empty densities). The totals equal results.damage.cascade.

escape_spectra.csv

species_z,symbol,beam,face,spectrum,lo,hi,count,per_ion_per_unit: the energy (spectrum = energy_ev, lo/hi in eV) and polar-angle (polar_deg, degrees from the outward normal) spectra of every species leaving each face (front or back). Each spectrum ends with -inf and inf rows for entries outside the grid, with empty densities. The density is per incident ion per eV or per degree.

ions.csv

index,fate,x_nm,y_nm,z_nm,energy_ev,dir_x,dir_y,dir_z,layer: the final state of every primary, by history index. fate is stopped, backscattered or transmitted; for escaped primaries the position is on the face and the energy and direction are outside the target. x is depth.

Dynamic runs ([dynamic])

A run with a [dynamic] table writes dynamic_summary.json, dynamic_steps.csv and dynamic_composition.csv instead of the static files.

dynamic_summary.json (format lindhard-dynamic-summary): the echoed input (with dynamic), physics-model list models, the species, totals (ions, steps, rejected attempts, fluence, slab count and thickness before and after, yields per ion), files and the timing run object.

dynamic_steps.csv: one row per accepted step; step 0 is the initial target. step, first_index (global index of the step’s first ion), ions, ions_done, fluence_cm2 (delivered so far), attempts (more than 1 if the adaptive bound rejected the step), max_change, clamped, removed_slabs (slabs that emptied), n_slabs, surface_nm (cumulative surface recession in nm; always 0 with erosion = false), thickness_nm (total of the finite slabs), then cumulative counts since the start: cum_backscattered, cum_transmitted, cum_stopped_in_target, cum_stopped_in_substrate, cum_sputtered, sputter_yield (atoms per ion so far), cum_recoils_transmitted and cum_sputtered_<Sym> per element.

dynamic_composition.csv: the slab profile after every step (step 0 is the initial target), one row per step and slab: step, slab (0 is the front), front_nm, back_nm, thickness_nm, then per element atoms_per_cm2_<Sym> and fraction_<Sym> (atom fraction). Slabs that emptied are gone from later steps.

Electron runs: electron_summary.json and electron_*.csv

An electron run writes these instead of the ion files.

electron_summary.json (format {"name": "lindhard-electron-summary", "version": 1}):

results.table_coverage is a numerical diagnostic: one entry per layer (layer), and for its elastic and inelastic table the grid bounds energy_min_ev and energy_max_ev and the counts below (E < energy_min_ev: the first row’s rate and distribution were used), within (both bounds included: interpolated, or a row read exactly) and above (E > energy_max_ev: the last row’s were used), over primaries and secondaries. They count rate evaluations, not collisions: the transport evaluates both tables of the electron’s layer once before every free flight, including flights cut short at a layer face, the flight after a face reflection and flights with zero total rate, so the totals exceed the number of elastic and inelastic events, and a layer’s elastic and inelastic totals are equal. They are not fractions of path length or of deposited energy either. Nonzero below or above counts say that part of the transport used the constant continuation of a table beyond its grid ([electron.tables]); they do not say how much that changed the result. Summaries written before the key existed lack it.

results.event_caps is another numerical diagnostic, of the collision cap ([electron.transport] max_events): secondary_tracks is the number of secondary electrons the cap cut off (one per capped track) and affected_histories the number of primary histories in which the primary or at least one secondary was cut off (each history counted once, however many of its electrons were capped). Both are counts, not energies; the energy the capped electrons still carried is budget.event_cap_ev. The primaries’ own caps stay in fates.event_capped, so fates.event_capped = 0 alone does not show that the secondary cascades ran to completion: check event_caps.affected_histories. A capped track is a truncated one, not a physical fate. Summaries written before the key existed read back with both counts zero.

electron_escape_spectra.csv: face,spectrum,class,lo,hi,count,per_primary_per_unit. For each face (front, back): the energy spectrum of all escaping electrons (spectrum = energy_ev, class = all, eV) and the polar-angle spectrum of each class (polar_deg, slow or fast, degrees from the outward normal). Each spectrum ends with -inf and inf rows for entries outside the grid, with empty densities. The density is per primary per eV or per degree.

electron_deposition_cylindrical.csv (with tally.cylindrical): ir,ix,r_lo_nm,r_hi_nm,depth_lo_nm,depth_hi_nm,energy_ev,ev_per_primary_per_nm3, one row per cell. electron_deposition_cartesian.csv (with tally.cartesian): ix,iy,iz,x_lo_nm,x_hi_nm,y_lo_nm,y_hi_nm,z_lo_nm,z_hi_nm,energy_ev,ev_per_primary_per_nm3. Energy deposited outside a grid is outside_ev in the summary. Like budget.deposited_ev, the cell energies are measured from the band bottom and, with secondaries, include Fermi-sea energy the beam did not supply.

electron_tables.csv: material,energy_ev,elastic_inverse_mfp_per_nm,inelastic_inverse_mfp_per_nm,inelastic_mean_loss_ev,inelastic_stopping_ev_per_nm, the tables the run used, per material and grid energy (the stopping power is λ⁻¹ ⟨W⟩ of the stored loss distribution).

Cross-section table cache (--table-cache)

Building the elastic and inelastic tables is the slow part of a short electron run (tens of seconds with penn-single-pole, far longer with penn-full; see docs/validation.md, “Electron oracles”). With --table-cache DIR, lindhard run looks each table up in DIR (created if missing) and builds and stores only the ones it does not find, so a series of runs that differ only in seed, history count, tallies or transport settings builds its tables once.

Each entry is two files, named by the SHA-256 of a key document: <kind>-<sha256>.toml, the table in the versioned cache form of lindhard::electron::data::CrossSectionTable, read back with that type’s loader, and <kind>-<sha256>.key.json, the key itself. The key spells out every input the table depends on:

• the key schema version and the table cache format_version;
• the build: crate version, git describe --always --dirty, and the SHA-256 of the running lindhard executable, so any rebuild that changes the code (an uncommitted edit included) misses, and a rebuilt binary never reuses an older binary’s table;
• the table kind and the exact energy grid (electron.tables, after the defaults are filled in);
• the material’s name, composition and density;
• elastic: the potential, the exchange and correlation-polarization corrections with all their inputs, the starting probability grid and the refinement tolerance;
• inelastic: the model (electron.inelastic.model), its Fermi energy, the SHA-256 and provenance of the optical ELF file, and the material’s band parameters.

Every f64 is written in shortest round-trip form, so a change in the last bit of any number is a different key. A lookup must find the stored key equal, byte for byte, to the run’s own (a mismatch under the same hash means the file was edited, and is an error naming the differing field); the table file is then loaded and validated by the library, and its axis and energy grid are checked against the run. A file found under the run’s key that fails any of these checks is an error, never a silent rebuild. A table of another cache format_version can never be found, since the version is part of the key. Writes go to a temporary file renamed into place, the table before its key.

A cached table is the built one bit for bit (the TOML cache form round-trips every f64), so outputs do not depend on whether the tables were built or read, nor on the thread count. Only physics.materials.*_table.source and .cache, and the timings in run, differ. lindhard-cli/tests/examples.rs checks a building run, a storing run and reading runs on 1 and 8 threads against each other. Remove the directory to reclaim space; nothing else prunes it.

Reusing an output directory

--out may name an existing directory; the run overwrites the files it writes and creates the directory if needed. The CLI also owns the reserved optional file names of the run’s mode. After a successful run, an optional file the run did not produce is removed if present: ions.csv (without tally.per_ion), and electron_deposition_cartesian.csv or electron_deposition_cylindrical.csv (without the matching deposition grid). A missing file is not an error; a failed removal is, and names the path. The summary is written last and lists only files that exist. Other files in the directory are never touched, and no cleanup happens between ion, electron and dynamic runs. Do not keep your own data under a reserved name.

Reproducibility

Everything in summary.json except the trailing run object, and every CSV file, is a function of the input and the binary only: byte-identical for the same input and seed at any thread count (for a dynamic run: the same three files, with run the only thread-dependent part of dynamic_summary.json). For an electron run the same holds for electron_summary.json (apart from run) and every electron_*.csv: the tables are built bit-identically on any thread count, and histories run in chunks of a fixed size (16, recorded as physics.transport.chunk_size) merged in chunk order. lindhard-cli/tests/examples.rs checks this on 1 and 4 threads (1, 2 and 8 for the dynamic example; 1 and 4 for the electron example, and with --table-cache 1 and 8). Floats are written in shortest round-trip form.

Compatibility and extension

Consumers must ignore keys they do not know. New results are added as new keys, never by changing existing ones:

• New tallies become new objects under results, with their settings as new keys under [tally] and their profiles as new CSV files listed under files. results.range, results.damage, results.sputtering and results.escapes were added this way, without a version bump.
• New model choices become new values of the existing [physics] keys, or new keys with defaults, so existing inputs keep their meaning.
• format.version is bumped only when an existing key is removed or changes meaning.

The electron schema and its output follow the same rules, with lindhard-electron-summary counting its versions separately:

• A new library model becomes a new value of an existing key (electron.inelastic.model, electron.elastic.potential, electron.transport.secondaries, electron.transport.boundary, a band kind, a phonon preset) and a new physics.models entry; existing values keep their meaning and defaults never change.
• New per-material data (for example a subshell binding-energy table, once inner-shell channels reach the transport loop, or a precomputed cross-section cache) becomes a new optional key of [electron.materials.<name>] that names a file. It must be read with the lindhard::electron::data loader of its type, so data without a provenance is refused, and recorded under physics.materials with its path, SHA-256 and provenance.
• The cross-section table cache is the one exception to the rule above, by decision (#168): it is a command-line flag (--table-cache DIR), not an input key, because it never changes a result (as with --threads, the input and its echo stay the same whether or not tables are reused), and a content-addressed directory keyed on every input of the build cannot name a stale file the way a hand-written path can. Its tables are still read with the lindhard::electron::data loader and echoed under physics.materials with their path, SHA-256 and provenance.
• A new tally becomes a key under [electron.tally], an object under results and, for profiles, a new CSV file listed under files.
• Every table keeps deny_unknown_fields, and every default is echoed, so an input written today still means the same thing.