lindhard/electron/inelastic/
bethe.rs

1//! Nonrelativistic Bethe stopping, the high-energy reference for the
2//! dielectric models.
3//!
4//! The Bethe equation as given by the Particle Data Group, S. Navas et al.,
5//! Phys. Rev. D 110, 030001 (2024), review 34 "Passage of particles through
6//! matter" (D. E. Groom and S. R. Klein, rev. August 2023), eq. (34.5):
7//!
8//! ```text
9//! -dE/dx = K z² (Z/A) (1/β²) [ (1/2) ln(2 m_e c² β² γ² W_max / I²) - β² - δ(βγ)/2 ].
10//! ```
11//!
12//! Here it is used **nonrelativistically and for a unit charge** (`z = 1`,
13//! `γ -> 1`, the `-β²` and density-effect terms dropped, `m_e c² β² = m_e v²
14//! = 2E` for an electron of kinetic energy `E`), and as a linear stopping
15//! power for an electron density `n` (the `K (Z/A) ρ` of the mass stopping
16//! power is `4π n r_e² m_e c²`):
17//!
18//! ```text
19//! S(E) = (4π n (e²/4πε₀)² / (m_e v²)) · (1/2) ln(2 m_e v² W_max / I²)
20//!      = (2π n (e²/4πε₀)² / E) · (1/2) ln(4 E W_max / I²).
21//! ```
22//!
23//! `W_max` is the largest energy loss the model allows. The PDG gives eq.
24//! (34.5) for heavy particles, where `W_max` is the kinematic maximum
25//! transfer to a free electron; for an incident electron the recommended
26//! stopping formula differs (identical-particle exchange). The single-pole
27//! Penn model has no exchange and allows losses up to the full kinetic energy
28//! (`W_max = E`), so its high-energy limit is the expression above with
29//! `W_max = E`, `S = (2π n (e²/4πε₀)² / E) ln(2E/I)`. That is the comparison
30//! this module exists for; it is not a recommended electron stopping power.
31
32use crate::constants::{COULOMB_E2, ELEMENTARY_CHARGE};
33
34/// Nonrelativistic Bethe stopping power of a unit-charge particle of
35/// electron mass (see the module docs), eV/m.
36///
37/// `electron_density_per_m3` is the density of the electrons the mean
38/// excitation energy `mean_excitation_ev` refers to, `energy_ev` the kinetic
39/// energy and `max_loss_ev` the largest energy loss `W_max`. Returns `NaN`
40/// unless all four are finite and positive.
41pub fn stopping_power_ev_per_m(
42    electron_density_per_m3: f64,
43    mean_excitation_ev: f64,
44    energy_ev: f64,
45    max_loss_ev: f64,
46) -> f64 {
47    let ok = [
48        electron_density_per_m3,
49        mean_excitation_ev,
50        energy_ev,
51        max_loss_ev,
52    ]
53    .iter()
54    .all(|v| v.is_finite() && *v > 0.0);
55    if !ok {
56        return f64::NAN;
57    }
58    let e_j = energy_ev * ELEMENTARY_CHARGE;
59    let log =
60        0.5 * (4.0 * energy_ev * max_loss_ev / (mean_excitation_ev * mean_excitation_ev)).ln();
61    // J/m -> eV/m
62    2.0 * std::f64::consts::PI * electron_density_per_m3 * COULOMB_E2 * COULOMB_E2 / e_j * log
63        / ELEMENTARY_CHARGE
64}
65
66#[cfg(test)]
67mod tests {
68    use super::*;
69
70    #[test]
71    fn matches_the_atomic_unit_form() {
72        // In Hartree atomic units S = (Ω²/(2E)) ln(2E/I) with Ω² = 4π n,
73        // for W_max = E. Check the SI evaluation against it.
74        let a0 = crate::constants::BOHR_RADIUS;
75        let eh = crate::constants::HARTREE_ENERGY / ELEMENTARY_CHARGE;
76        let (n_au, e_au, i_au) = (0.05, 300.0, 0.9);
77        let want_au = 4.0 * std::f64::consts::PI * n_au / (2.0 * e_au) * (2.0 * e_au / i_au).ln();
78        let got = stopping_power_ev_per_m(n_au / a0.powi(3), i_au * eh, e_au * eh, e_au * eh);
79        let got_au = got * a0 / eh;
80        // e²/(4πε₀) = 1 hartree bohr holds to the CODATA rounding, ~1e-10.
81        assert!(
82            (got_au / want_au - 1.0).abs() < 1e-9,
83            "{got_au} vs {want_au}"
84        );
85        assert!(stopping_power_ev_per_m(1e29, 0.0, 1e3, 1e3).is_nan());
86    }
87}
lindhard/electron/inelastic/bethe.rs

lindhard/electron/inelastic/
bethe.rs
