Module penn 

The single-pole Penn algorithm: the momentum-dependent energy-loss function, the differential inverse inelastic mean free path (DIIMFP), the inelastic mean free path (IMFP) and the stopping power, built from an optical ELF.

Sources

The algorithm is Penn’s, D. R. Penn, Phys. Rev. B 35, 482 (1987), doi:10.1103/PhysRevB.35.482 (the “simple Penn algorithm” or single-pole approximation, SPA). That paper is closed access and was not read for this implementation. The equations below are taken from the account of the SPA by Penn and co-authors in H. Shinotsuka, B. Da, S. Tanuma, H. Yoshikawa, C. J. Powell, D. R. Penn, Surf. Interface Anal. 49, 238 (2017), doi:10.1002/sia.6123 (NIST author manuscript, open copy PMC5524379), and every equation number in this module refers to that paper (“S2017”). The plasmon-pole idea for extending an optical ELF into momentum transfer goes back to Lindhard, Mat. Fys. Medd. Dan. Vid. Selsk. 28, no. 8 (1954) (the dielectric function each pole stands for) and Ritchie and Howie, Philos. Mag. 36, 463 (1977), doi:10.1080/14786437708244948 (neither read here; cited for attribution only).

Hartree atomic units (m_e = e = ħ = 1) are used internally, as in S2017; the public API is in eV and metres.

The model

The optical ELF Im[-1/ε(ω)] is written as a superposition of free-electron (Lindhard) loss functions of plasma frequency ω_p with weight g(ω_p) = (2/(π ω_p)) Im[-1/ε(ω_p)] (S2017 eqs. (4)-(5)), and in the SPA each Lindhard loss function is replaced by a single plasmon pole (S2017 eq. (7)),

Im[-1/ε_L(q, ω; ω_p)] ≈ (π/2) (ω_p² / ω_q(ω_p)) δ(ω - ω_q(ω_p)),

with Penn’s dispersion (S2017 eq. (8))

ω_q²(ω_p) = ω_p² + (1/3) (k_F(ω_p) q)² + q⁴/4,    k_F(ω_p) = (3π/4)^(1/3) ω_p^(2/3).

The ω_p integration then gives (S2017 eqs. (9)-(13))

Im[-1/ε(q, ω)] = Im[-1/ε(ω₀)] / (1 + π q² / (6 k_F(ω₀))),

where ω₀(q, ω) is the plasma frequency whose pole sits at ω, the real positive root x = ω₀^(2/3) of the cubic (S2017 eq. (11)) x³ + a(q) x² + b(q, ω) = 0, a = (π²/48)^(1/3) q², b = q⁴/4 - ω². There is a root only for b < 0, i.e. q < sqrt(2ω); otherwise the loss function is zero. (SinglePolePenn::loss_function.)

Kinematics are nonrelativistic throughout: the relativistic factors of S2017 eq. (2) are set to 1 (c -> ∞). Relativistic kinematics is a later option. The DIIMFP is (S2017 eq. (2), nonrelativistic)

p(T, ω) = (1 / (π T')) ∫_{q-}^{q+} (dq/q) Im[-1/ε(q, ω)],
q± = sqrt(2T') ± sqrt(2(T' - ω)),

with T the kinetic energy measured from the Fermi level, T' = T + E_F, and losses 0 < ω <= T (S2017 eq. (3), ω_max = T' - E_F). The Fermi energy E_F defaults to zero, so T = T' is the kinetic energy. The IMFP is λ = [∫ p dω]⁻¹ (S2017 eq. (3)) and the stopping power S = ∫ ω p dω.

Low energies

The SPA is an approximation to the full Penn algorithm (FPA) that holds at high energy and fails at low energy. S2017 computed IMFPs of liquid water from one ELF with the FPA, the SPA, a simplified SPA and the Mermin method. In its section “Calculated IMFPs from the four algorithms” (repeated in the Summary; open copy PMC5524379) it reports that the FPA and SPA IMFPs “agree well … for energies over 50 eV. For energies less than 30 eV, the SPA IMFPs become larger than the FPA IMFPs with decreasing electron energy. The smaller FPA IMFPs in this energy range must be due to the contributions of single-electron excitations to the ELF, as shown in Eqn (6), that were neglected in the SPA.” S2017 shows its values below 50 eV only for the trend, as less reliable, and computed water, not metals. In eq. (7) above each ω_p contributes only on its pole ω = ω_q(ω_p) >= ω_p, so a loss ω draws on ELF weight at ω_p <= ω only: the single-electron (electron-hole) continuum that the full Lindhard function of a plasmon has at losses below ω_p is dropped. An electron whose allowed losses lie below the main ELF peak then has almost nothing to lose.

Measured in this code (#173; lindhard-cli/examples/inelastic_low_energy.rs, mode imfp), λ in nm for the Hagemann-Gudat-Kunz (1975) ELFs, model Fermi energy 0, for the single pole, the full Penn algorithm (super::full_penn) and Mermin (super::mermin):

E (eV)	Al SPA	Al full	Al Mermin	Au SPA	Au full	Au Mermin	Cu SPA	Cu full	Cu Mermin
12	24.6	1.07	1.11	2.31	1.30	1.46	2.02	1.01	1.77
19	9.83	0.720	0.739	1.66	1.00	1.04	1.10	0.735	1.17
30	0.363	0.512	0.512	1.33	0.810	0.805	0.743	0.586	0.889
100	0.402	0.401	0.433	0.528	0.533	0.469	0.507	0.503	0.628
1000	1.96	1.97	2.13	1.29	1.30	1.30	1.63	1.65	1.93

The SPA agrees with the full algorithm to 1 % at 100 eV and above, and below 30 eV becomes longer with decreasing energy, as S2017 reports for water: up to 23 times longer below 20 eV for Al (whose ELF is dominated by the 15 eV plasmon), 1.5 to 1.8 times for Au at 12 to 50 eV, and up to 2 times for Cu at 12 to 19 eV. (At 30 eV the Al SPA value is the shorter one.) These are the energies of the secondary electrons that make up the yield δ: swapping only the table rows below 30 eV for Mermin rows takes the single-pole Al δ at 400 eV from 6.9 to 1.4, and Au keeps most of its single-pole excess over Mermin in the rows below 50 eV (docs/validation.md, “Secondary-electron yield δ(E)”). This is the published low-energy behaviour of the approximation; no error of this implementation was found (at 100 eV and above it agrees with the full algorithm of super::full_penn to about 1 %).

Optional exchange correction (Born-Ochkur)

ExchangeCorrection makes the primary and the struck electron indistinguishable. The source is the open account in P. de Vera, S. Taioli, P. E. Trevisanutto, M. Dapor, I. Abril, S. Simonucci and R. Garcia-Molina, “Energy Deposition around Swift Carbon-Ion Tracks in Liquid Water”, Int. J. Mol. Sci. 23, 6121 (2022), doi:10.3390/ijms23116121 (open access, PMC9181504), Section 2.4.3, eqs. (32)-(35) (“dV2022” below). It cites Ochkur and Rudge for the factor (their refs. 122-125) and its own earlier work (their ref. 42) for the setup; none of these was read. From dV2022, for a loss E = W + B_α that ionises a shell of binding energy B_α and emits an electron of kinetic energy W:

• the exchange term of the DIIMFP is (1/T) ∫ dk/k F Im[-1/ε(k, E)] (eq. (32), last two terms), with the Born-Ochkur factor F = -x + x², x = (k²/2m) / (T - W) (text below eq. (32); ħ = m = 1 here). The denominator is T - W = T - E + B_α: it contains the emitted energy W, not the loss E;
• its k limits are sqrt(2mT) ± sqrt(2m(T - E)) (eq. (34)), the limits of the plain first Born DIIMFP at loss E;
• the emitted energy is limited by indistinguishability to W <= (T - B_α)/2 (text below eq. (35)), i.e. the loss to E <= (T + B_α)/2.

Here that term is (1/(π T')) ∫ dq/q F Im[-1/ε(q, ω)] with ω = E and x = (q²/2) / (T' - ω + B) (SinglePolePenn::diimfp_with_binding_per_m_ev), and the loss is limited to ω <= (T' + B)/2. B is the binding energy of the channel the loss is attributed to (super::inner_shell). The plain model (SinglePolePenn::diimfp_per_m_ev, SinglePolePenn::imfp_and_stopping) and the valence channel use B = 0, so x = (q²/2)/(T' - ω) and ω <= T'/2, which is also the ω_max of the SSPA* variant of S2017. The SPA loss function vanishes for q² >= 2ω (S2017 eq. (11)), and ω <= (T' + B)/2 gives T' - ω + B >= ω, so 0 <= x <= 1 wherever the integrand is nonzero and -1/4 <= F <= 0 for every channel.

Where this departs from dV2022 eq. (32), deliberately:

• Direct term. dV2022 applies a Coulomb-field correction (their ref. 126, not read) to the direct term only: T -> T + 2B_α in its prefactor 1/(T + 2B_α) and in its k limits (eq. (33)). We do not apply it: the direct term here is the plain first Born (S2017 eq. (2)) DIIMFP with prefactor 1/T' and the eq. (34) limits for every channel, which is what dV2022 eq. (32) reduces to for B_α = 0. So for an inner shell only the exchange term, its denominator and the loss limit follow eq. (32); the direct term does not. Adding the correction would also break the property that the channels add up to the plain-model total.
• Valence binding energy. dV2022 uses an energy-dependent mean binding energy B(T) of the outer shells, fitted for liquid water in their ref. 42 (not read). No such function is available for a general target, so the valence channel here has B = 0: its exchange denominator is T' - ω and its limit ω <= T'/2.
• Fermi energy. dV2022 has none; here T' = T + E_F stands for its T (zero by default).

The correction is optional and applied only for kinetic energies below the energy given to ExchangeCorrection::new; at and above it the model is the plain SPA (no ω_max cut either), so λ(E) has a step of the size of the correction at that energy. The total is then evaluated by direct numerical integration over ω of the DIIMFP (no closed form), with the outer tolerance 100 × rel_tol since each DIIMFP value is itself an integral.

How the integrals are done

DIIMFP (SinglePolePenn::diimfp_per_m_ev): the q integral of S2017 eq. (2) is done in ln q over [q-, q+], restricted to the q where ω₀(q, ω) lies inside the tabulated range of the optical ELF (the ELF is zero outside it: nothing is extrapolated), and split at the q where ω₀ crosses a table knot, so that each piece is smooth. ω₀ is the root of the S2017 eq. (11) cubic, found by Newton’s method from x = ω^(2/3) (the cubic is convex and increasing for x > 0, so the iteration converges monotonically).

IMFP and stopping power (SinglePolePenn::imfp_and_stopping): the double integral over (ω, q) is done in the other order. Because the SPA puts the whole weight of each ω_p on one curve ω = ω_q(ω_p), inserting S2017 eqs. (4), (5), (7) into ∫∫ dω (dq/q) gives exactly

λ⁻¹ = (1/(π T')) ∫ dω_p Im[-1/ε(ω_p)] ω_p ∫_{qa}^{qb} dq / (q ω_q(ω_p)),
S   = (1/(π T')) ∫ dω_p Im[-1/ε(ω_p)] ω_p ln(qb / qa),

where [qa, qb] is the set of q whose pole ω = ω_q(ω_p) is allowed: q- <= q <= q+ is equivalent to ω <= q k - q²/2 (k = sqrt(2T')), and ω <= T to q <= q_T with ω_q(q_T) = T. The function φ(q) = q k - q²/2 - ω_q(q) is concave (the numerator of d²ω_q/dq² is c ω_p² + (3/2) q² ω_p² + (3/4) c q⁴ + q⁶/8 > 0, c = k_F²/3), so the allowed set is one interval; its ends are found by Newton’s method from q = 0 and q = k, each converging monotonically. The inner q integral of λ⁻¹ is elementary (y = q², ∫ dy/(y sqrt(y²/4 + c y + ω_p²)) = -(1/ω_p) ln((2ω_p² + c y + 2 ω_p ω_q)/y)). This reduction (the order of integration and the closed-form inner integral) is our own algebra on the published model; the tests check it against the direct ∫ p(T, ω) dω to 1e-5.

The outer integrals use adaptive Gauss-Legendre quadrature split at the table knots (super::quadrature): the estimated relative error of each integral is at most twice the tolerance (default DEFAULT_RELATIVE_TOLERANCE = 1e-7), plus the error of the Newton roots (relative 1e-15).

High-energy (Bethe) limits

As T -> ∞, qa -> ω_p/k and qb -> k. The stopping power tends to S -> (Ω_p²/(2T)) ln(2T/I), with Ω_p² = (2/π) ∫ ω Im[-1/ε] dω and ln I = ∫ ω Im[-1/ε] ln ω dω / ∫ ω Im[-1/ε] dω: the nonrelativistic Bethe formula with W_max = T (see super::bethe). The inverse IMFP tends to

λ⁻¹ -> (1/(2π T')) ∫ dω_p Im[-1/ε(ω_p)] ln(8 T' / (ω_p + k_F(ω_p)²/3)),

which is the Bethe form T/λ = A ln(B T) of S2017 eq. (27) (with α(T) = 1): its slope A = (1/(2π)) ∫ Im[-1/ε] dω = Ω_p² M_tot² / 2 (Ω_p and M_tot² as in S2017 eq. (29)) is fixed by the optical ELF alone (S2017 eqs. (28)-(29)), and its constant ln B = ln 8 - ⟨ln(ω_p + k_F²/3)⟩ (average weighted by the ELF) by the SPA dispersion. That limit is our own expansion of the model (SinglePolePenn::bethe_inverse_imfp_per_m); the tests check that the model approaches it.

Structs

ExchangeCorrection

Optional exchange correction (module docs, “Optional exchange correction”), applied for kinetic energies strictly below applies_below_ev.

InelasticPoint

One point of an IMFP and stopping-power table.

SinglePolePenn

The single-pole Penn model built from one optical ELF. See the module docs for the equations, the units and the integration tolerance.

Constants

DEFAULT_RELATIVE_TOLERANCE

Default relative tolerance of the IMFP, stopping and DIIMFP integrals.

Functions

born_ochkur_factor

The Born-Ochkur exchange factor F = -x + x² of dV2022 (text below eq. (32); module docs), with x = (ħ²q²/2m) / (T - W) and W = E - B the kinetic energy of the emitted electron, for a recoil energy ħ²q²/2m = recoil_ev, a kinetic energy T = energy_ev (use T' = T + E_F with a Fermi energy), a loss E = loss_ev and a binding energy B = binding_ev, all in eV. This is the factor the DIIMFP integrand of SinglePolePenn::diimfp_with_binding_per_m_ev is multiplied by (1 + F). NaN unless T - W > 0.