lindhard/ion/stopping/
lindhard_scharff.rs

1//! Lindhard-Scharff electronic stopping.
2//!
3//! J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961): at ion velocities
4//! below `v0 Z1^(2/3)` the electronic stopping is proportional to velocity.
5//!
6//! # Reduced units (Lindhard, Scharff, Schiott, Mat. Fys. Medd. 33 (14), 1963)
7//!
8//! * screening length `a = 0.8853 a0 (Z1^(2/3) + Z2^(2/3))^(-1/2)`
9//!   (Thomas-Fermi, universal-length form; `a0` is the Bohr radius),
10//! * reduced energy `ε = E a M2 / (Z1 Z2 e² (M1 + M2))`, with `e²` meaning
11//!   `e²/(4π ε0)`,
12//! * reduced path length `ρ = N x 4π a² M1 M2 / (M1 + M2)²`,
13//! * reduced stopping `(dε/dρ)_e = k_L ε^(1/2)`, with
14//!   `k_L = ξ_e 0.0793 Z1^(1/2) Z2^(1/2) (A1+A2)^(3/2)
15//!   / ((Z1^(2/3)+Z2^(2/3))^(3/4) A1^(3/2) A2^(1/2))`.
16//!
17//! `ξ_e = Z1^(1/6)` is the Lindhard-Scharff (1961) velocity-proportionality
18//! factor (`ξ_e ≈ Z1^(1/6)`, order 1 to 2). **Deviation from the issue text:**
19//! the formula quoted there omits `ξ_e`; without it the Z1 dependence does not
20//! agree with the dimensional LS form below (off by `Z1^(1/6)`), so it is
21//! included here. The 0.0793 is the rounded coefficient; a test checks the
22//! whole expression against the dimensional form
23//! `S_e = Z1^(1/6) 8π e² a0 Z1 Z2 / (Z1^(2/3)+Z2^(2/3))^(3/2) · v/v0` to 1 %.
24//!
25//! A per-element multiplicative correction `f(Z2)` may be attached
26//! ([`LindhardScharff::with_correction`]); it defaults to 1. It must be finite
27//! and non-negative (zero switches the element's stopping off); anything else
28//! is rejected at evaluation time with [`StoppingError::InvalidParameter`].
29
30use super::{check_energy, target, ElectronicStopping, Ion, StoppingError, ValidityRange};
31use crate::constants::{BOHR_RADIUS, COULOMB_E2};
32use crate::units::J_PER_EV;
33use std::f64::consts::PI;
34
35/// Lindhard-Scharff screening length `a`, m.
36pub fn screening_length(z1: u8, z2: u8) -> f64 {
37    0.8853 * BOHR_RADIUS / (f64::from(z1).powf(2.0 / 3.0) + f64::from(z2).powf(2.0 / 3.0)).sqrt()
38}
39
40/// Reduced energy `ε` for ion energy `energy_ev` on a target of mass `m2_amu`.
41pub fn reduced_energy(ion: &Ion, z2: u8, m2_amu: f64, energy_ev: f64) -> f64 {
42    let a = screening_length(ion.z(), z2);
43    energy_ev * J_PER_EV * a * m2_amu
44        / (f64::from(ion.z()) * f64::from(z2) * COULOMB_E2 * (ion.mass_amu() + m2_amu))
45}
46
47/// Factor converting reduced stopping `dε/dρ` to a cross section in J m².
48pub fn reduced_to_si_factor(ion: &Ion, z2: u8, m2_amu: f64) -> f64 {
49    let a = screening_length(ion.z(), z2);
50    4.0 * PI * a * f64::from(ion.z()) * f64::from(z2) * COULOMB_E2 * ion.mass_amu()
51        / (ion.mass_amu() + m2_amu)
52}
53
54/// The Lindhard-Scharff coefficient `k_L` including `ξ_e = Z1^(1/6)` (see
55/// module docs).
56pub fn k_l(ion: &Ion, z2: u8, m2_amu: f64) -> f64 {
57    let (z1, z2f) = (f64::from(ion.z()), f64::from(z2));
58    let (a1, a2) = (ion.mass_amu(), m2_amu);
59    0.0793 * z1.powf(1.0 / 6.0) * z1.sqrt() * z2f.sqrt() * (a1 + a2).powf(1.5)
60        / ((z1.powf(2.0 / 3.0) + z2f.powf(2.0 / 3.0)).powf(0.75) * a1.powf(1.5) * a2.sqrt())
61}
62
63/// Lindhard-Scharff model with optional per-element correction factors.
64#[derive(Debug, Clone, Default)]
65pub struct LindhardScharff {
66    corrections: Vec<(u8, f64)>,
67}
68
69impl LindhardScharff {
70    /// Plain LS (all corrections 1).
71    pub fn new() -> Self {
72        Self::default()
73    }
74
75    /// Multiply the stopping in target element `z2` by `factor`.
76    pub fn with_correction(mut self, z2: u8, factor: f64) -> Self {
77        self.corrections.retain(|(z, _)| *z != z2);
78        self.corrections.push((z2, factor));
79        self
80    }
81
82    fn correction(&self, z2: u8) -> f64 {
83        self.corrections
84            .iter()
85            .find(|(z, _)| *z == z2)
86            .map_or(1.0, |(_, f)| *f)
87    }
88
89    /// Reduced stopping `dε/dρ` at reduced energy `eps` (without correction).
90    pub fn reduced_stopping(&self, ion: &Ion, z2: u8, m2_amu: f64, eps: f64) -> f64 {
91        k_l(ion, z2, m2_amu) * eps.sqrt()
92    }
93}
94
95impl ElectronicStopping for LindhardScharff {
96    fn name(&self) -> &'static str {
97        "lindhard-scharff"
98    }
99
100    fn stopping(&self, ion: &Ion, target_z: u8, energy_ev: f64) -> Result<f64, StoppingError> {
101        check_energy(energy_ev)?;
102        let m2 = target(target_z)?.atomic_weight;
103        let factor = self.correction(target_z);
104        if !(factor.is_finite() && factor >= 0.0) {
105            return Err(StoppingError::InvalidParameter {
106                name: "correction_factor",
107                value: factor,
108            });
109        }
110        let eps = reduced_energy(ion, target_z, m2, energy_ev);
111        Ok(factor
112            * self.reduced_stopping(ion, target_z, m2, eps)
113            * reduced_to_si_factor(ion, target_z, m2))
114    }
115
116    fn sqrt_energy_coefficient(&self, ion: &Ion, target_z: u8) -> Option<f64> {
117        let m2 = target(target_z).ok()?.atomic_weight;
118        let factor = self.correction(target_z);
119        if !(factor.is_finite() && factor >= 0.0) {
120            return None;
121        }
122        // S = factor k_L sqrt(eps) f, with eps = E eps_per_ev (reduced_energy
123        // is linear in E), so S = [factor k_L sqrt(eps_per_ev) f] sqrt(E).
124        let eps_per_ev = reduced_energy(ion, target_z, m2, 1.0);
125        Some(
126            factor
127                * self.reduced_stopping(ion, target_z, m2, eps_per_ev)
128                * reduced_to_si_factor(ion, target_z, m2),
129        )
130    }
131
132    fn validity(&self, ion: &Ion) -> ValidityRange {
133        // v < v0 Z1^(2/3): E/A below about 25 keV Z1^(4/3).
134        ValidityRange {
135            min_energy_ev: 0.0,
136            max_energy_ev: ion.mass_amu() * super::velocity_scale_energy_per_amu_ev(ion.z(), 1.0),
137        }
138    }
139}
140
141#[cfg(test)]
142mod tests {
143    use super::*;
144    use crate::elements::element;
145
146    #[test]
147    fn reduced_stopping_is_k_sqrt_eps() {
148        let ion = Ion::new(15).unwrap();
149        let m = LindhardScharff::new();
150        let m2 = element(14).unwrap().atomic_weight;
151        let f = reduced_to_si_factor(&ion, 14, m2);
152        let mut ratio0 = None;
153        for e in [1.0e2, 1.0e3, 1.0e4, 1.0e5] {
154            let s = m.stopping(&ion, 14, e).unwrap();
155            let eps = reduced_energy(&ion, 14, m2, e);
156            let r = s / f / eps.sqrt();
157            let r0 = *ratio0.get_or_insert(r);
158            assert!((r / r0 - 1.0).abs() < 1e-12);
159            assert!((r / k_l(&ion, 14, m2) - 1.0).abs() < 1e-12);
160        }
161    }
162
163    #[test]
164    fn matches_dimensional_form_within_one_percent() {
165        for (z1, z2) in [(1u8, 14u8), (15, 14), (11, 29), (33, 14), (5, 79)] {
166            let ion = Ion::new(z1).unwrap();
167            let e = 2.0e3;
168            let s = LindhardScharff::new().stopping(&ion, z2, e).unwrap();
169            let (z1f, z2f) = (f64::from(z1), f64::from(z2));
170            let v = super::super::speed_nonrel(&ion, e) / super::super::bohr_velocity();
171            let dim = z1f.powf(1.0 / 6.0) * 8.0 * PI * COULOMB_E2 * BOHR_RADIUS * z1f * z2f
172                / (z1f.powf(2.0 / 3.0) + z2f.powf(2.0 / 3.0)).powf(1.5)
173                * v;
174            assert!((s / dim - 1.0).abs() < 0.01, "{z1} {z2}: {}", s / dim);
175        }
176    }
177
178    #[test]
179    fn sqrt_energy_coefficient_reproduces_stopping() {
180        for (z1, z2, corr) in [
181            (5u8, 14u8, 1.0),
182            (33, 14, 1.25),
183            (18, 29, 0.0),
184            (1, 79, 1.0),
185        ] {
186            let ion = Ion::new(z1).unwrap();
187            let m = LindhardScharff::new().with_correction(z2, corr);
188            let c = m.sqrt_energy_coefficient(&ion, z2).unwrap();
189            for e in [1.0, 37.0, 5.0e3, 2.0e5] {
190                let s = m.stopping(&ion, z2, e).unwrap();
191                assert!((c * e.sqrt() - s).abs() <= 1e-13 * s, "{z1} {z2} {e}");
192            }
193        }
194        let ion = Ion::new(5).unwrap();
195        // Failing configurations and unknown targets give no coefficient.
196        let bad = LindhardScharff::new().with_correction(14, f64::NAN);
197        assert!(bad.sqrt_energy_coefficient(&ion, 14).is_none());
198        assert!(LindhardScharff::new()
199            .sqrt_energy_coefficient(&ion, 0)
200            .is_none());
201    }
202
203    #[test]
204    fn correction_factor_scales() {
205        let ion = Ion::new(5).unwrap();
206        let base = LindhardScharff::new().stopping(&ion, 14, 1e4).unwrap();
207        let c = LindhardScharff::new()
208            .with_correction(14, 1.25)
209            .stopping(&ion, 14, 1e4)
210            .unwrap();
211        assert!((c / base - 1.25).abs() < 1e-14);
212    }
213
214    #[test]
215    fn rejects_bad_correction_factor() {
216        let ion = Ion::proton();
217        for bad in [-1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
218            let m = LindhardScharff::new().with_correction(14, bad);
219            assert!(matches!(
220                m.stopping(&ion, 14, 1.0e4),
221                Err(StoppingError::InvalidParameter {
222                    name: "correction_factor",
223                    ..
224                })
225            ));
226            // Other target elements are unaffected.
227            assert!(m.stopping(&ion, 8, 1.0e4).unwrap() > 0.0);
228        }
229    }
230
231    #[test]
232    fn zero_correction_factor_gives_zero_stopping() {
233        let m = LindhardScharff::new().with_correction(14, 0.0);
234        assert_eq!(m.stopping(&Ion::proton(), 14, 1.0e4).unwrap(), 0.0);
235    }
236
237    #[test]
238    fn rejects_bad_energy() {
239        let ion = Ion::new(5).unwrap();
240        assert!(LindhardScharff::new().stopping(&ion, 14, 0.0).is_err());
241        assert!(LindhardScharff::new().stopping(&ion, 0, 1.0).is_err());
242    }
243}
lindhard/ion/stopping/lindhard_scharff.rs

lindhard/ion/stopping/
lindhard_scharff.rs
