Module full_penn 

The full Penn algorithm (FPA): the optical energy-loss function expanded over the Lindhard loss functions of free-electron gases, and the DIIMFP, inelastic mean free path and stopping power built from it.

Sources

The algorithm is Penn’s, D. R. Penn, Phys. Rev. B 35, 482 (1987), doi:10.1103/PhysRevB.35.482, which is closed access and was not read. As for the single-pole model (super::penn), the equations are those of the account 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, PMC5524379, “S2017”), section “Full Penn algorithm (FPA)”, eqs. (4)-(6):

Im[-1/ε(q, ω)] = ∫_0^∞ g(ω_p) Im[-1/ε_L(q, ω; ω_p)] dω_p,   (S2017 eq. (4))
g(ω_p) = (2/(π ω_p)) Im[-1/ε(ω_p)],                         (S2017 eq. (5))

where ε_L is the Lindhard dielectric function of a free-electron gas of plasma frequency ω_p (super::lindhard_gas; ω_p² = 4π n), and the optical ELF Im[-1/ε(ω)] is the user’s table. S2017 eq. (6) splits the result into a plasmon-pole part and a single-electron part,

Im[-1/ε(q, ω)] = Im[-1/ε(q, ω)]_pl + Im[-1/ε(q, ω)]_se.

S2017 defers “the details of the calculation procedures for each term” to its reference 22 (Shinotsuka et al., Surf. Interface Anal. 47, 871 (2015)), which was not read either. The procedure used here is our own, a direct evaluation of eq. (4): for every ω_p the Lindhard loss function is its electron-hole continuum Im ε/|ε|² plus the undamped plasmon delta function of weight π/(∂Re ε/∂ω) (exactly the pole of the SPA, S2017 eq. (7), at q -> 0). Integrating the delta function over ω_p gives a smooth plasmon term g(ω₀) W/|dω_pl/dω_p| at the ω_p = ω₀ whose plasmon sits at the loss ω. Below the plasmon cutoff the plasmon term is the analogue of the SPA’s S2017 eq. (13); above it, and at all q outside the plasmon, the continuum term is the “single-electron excitations” that the SPA neglects. The two are the pl and se of eq. (6) as far as we can tell from the account in S2017; nothing here was checked against Penn’s own formulae.

The kinematics, units and the definitions of the DIIMFP, IMFP and stopping power are those of the single-pole module (S2017 eqs. (2)-(3), nonrelativistic; Hartree units internally, eV and metres in the API).

Sum rules

Each Lindhard gas satisfies ∫ ω Im[-1/ε_L] dω = (π/2) ω_p² at every q (plasmon included; tested in super::lindhard_gas), so with eq. (5) the expanded ELF has ∫ ω Im[-1/ε(q, ω)] dω = ∫ ω_p Im[-1/ε(ω_p)] dω_p, the optical f-sum, at every q (FullPenn::f_sum_ev2_at and the integration tests check this).

How the integrals are done

IMFP and stopping power (FullPenn::imfp_and_stopping) are evaluated with the order of integration ω_p, q, ω:

λ⁻¹ = (1/(π T')) ∫ dω_p g(ω_p) ∫ (dq/q) ∫_0^{ω_c(q)} Im[-1/ε_L] dω,
S   = (1/(π T')) ∫ dω_p g(ω_p) ∫ (dq/q) ∫_0^{ω_c(q)} ω Im[-1/ε_L] dω,
ω_c(q) = min(T, T' - (q - k)²/2),   k = sqrt(2T'),   0 < q < 2k,

(q- <= q <= q+ is ω <= T' - (q - k)²/2, S2017 eq. (2)). The innermost integral is split at the kinks of the Lindhard continuum and includes the plasmon delta function if its energy is below ω_c(q); the q integral is in ln q, split at the kinematic landmarks and at the momenta where the plasmon enters (ω_pl = ω_c) and leaves (cutoff) the allowed region; the outer integral is split at every knot of the table (the ELF is piecewise linear, and a feature narrower than the spacing of a subset of knots could be missed entirely by the adaptive quadrature, so no knot is dropped). All three use the adaptive Gauss-Legendre quadrature of super::quadrature, each with the tolerance of the model (DEFAULT_FULL_TOLERANCE = 1e-4): the result has a relative error of a few times that. The DIIMFP is the same integrals in the order q, ω_p at fixed ω.

The outer (ω_p) integral of both runs on quadrature::integrate_weighted (#256): g(ω_p) is linear between knots, every knot interval still gets its own Gauss-Legendre rule (so the rule above, no knot dropped, holds), and the error of the smooth Lindhard factor is checked on groups of up to 16 neighbouring intervals instead of by bisecting every interval; an interval that fails is bisected as before. The plasmon term finds the plasma frequency whose plasmon sits at ω by bisecting on the sign of Re ε_L (one evaluation per step) instead of solving for the plasmon energy at every step. Neither changes the integrals beyond their tolerance; measured on the Al ELF of Hagemann et al. (148 knots), they cut the time of one DIIMFP evaluation by about 6 and of one IMFP point by about 4.5 (docs/validation.md, cost of a full-Penn build).

The edge resonance and the bound on the ω_p bisection (#298). At fixed (q, ω) the continuum term starts at the plasma frequency x_a whose gas has its continuum edge q k_F + q²/2 at ω. Where that gas’s plasmon is about to enter the continuum there (Re ε_L ≈ 0 at the edge; the plasmon cutoff), Re ε_L crosses zero just above x_a, where Im ε_L is still small, so Im ε_L/|ε_L|² has a narrow, tall resonance (on the Si ELF at q = 1.48 a.u., ω = 160 eV: Re ε_L = 1e-6 at x_a, the zero 1.5e-7 Hartree above it, a peak of height above 1e6 and width about 3e-8 Hartree). It is integrable and holds a finite share of the first knot interval, but the bisection that halves its absolute tolerance at every level then asks the cells on the peak for a relative accuracy below the floating-point noise of Re ε_L (a cancellation of terms of order 1 down to 1e-7) and fills the peak with cells down to the maximum depth: seconds for one evaluation. The bisection of each knot interval therefore also accepts a cell that agrees with its halves to the tolerance relative to its own integral, which keeps the error contract (see quadrature::integrate_weighted), and has a budget of OMEGA_P_INTERVAL_BUDGET rule applications; anything accepted without meeting its tolerance is counted and reported by FullPenn::loss_function_evaluation (LossFunctionEvaluation). The values elsewhere are unchanged within the tolerance (the inverse IMFP and stopping power by up to about 1e-7; measured differences in docs/validation.md, “Electron oracles”). The outer ω_p integral of FullPenn::imfp_and_stopping has the same bisection and budget and counts its unconverged cells (FullPenn::imfp_and_stopping_evaluation, ImfpEvaluation; also FullPenn::imfp_unconverged); its inner integrals over q and ω are not bounded and report nothing.

The momentum profiles of a DiimfpGrid had the same failure one level up. Above the end of a table whose last ELF value is not zero (Si: 199 eV), the plasmon term jumps where the plasmon of the last tabulated gas reaches ω, and its value carries a relative noise of a few 1e-8 (the finite-difference slope); the bisection of a profile located the jump and then filled its neighbourhood with cells chasing the noise (Si, ω = 208.07 eV: over 10⁵ loss evaluations for one profile). The profiles’ bisection (quadrature::integrate_panels) has the same relative acceptance and budget, and counts what did not converge (DiimfpGrid::unconverged_profile_panels), and the loss evaluations inside a profile that did not converge (DiimfpGrid::unconverged_profile_loss_evaluations). The direct DIIMFP counts the loss evaluations of its momentum integral that did not converge, with the same definition (FullPenn::diimfp_evaluation, DiimfpEvaluation); its momentum integral itself is not bounded and reports nothing.

Tables. An inelastic table needs the DIIMFP at some 10⁴ losses per row, which at about 15 ms per evaluation (and with the quadrature noise of a 1e-4 tolerance driving the table’s 1e-5 density refinement to its node cap) is hours per row. The table builder therefore reads the DIIMFP from a DiimfpGrid (FullPenn::diimfp_grid): the q integrand is independent of the electron energy, so it is tabulated once per loss node and integrated over each row’s momentum window, with an interpolation error held to the model tolerance (see DiimfpGrid). The inverse IMFP of every row is still FullPenn::imfp_and_stopping.

No Fermi-energy, exchange or relativistic correction is included, as for the single-pole model.

Comparison with the single-pole model (synthetic Drude fixture)

For the Drude plasmon of the tests and examples (E_p = 20 eV, γ = 5 eV, synthetic, tabulated at 1000 energies from 10 meV to 100 keV; default tolerances; cargo run --release -p lindhard --example penn_full_vs_single_pole), ratios full/single-pole of the stopping power S and of the IMFP λ:

   E / eV     S ratio    λ ratio
     20.0      3.7971     0.2830
     30.0      1.6561     0.5489
     50.0      0.7988     1.1365
    100.0      0.9605     1.0197
    200.0      0.9860     1.0067
    500.0      0.9953     1.0028
   1000.0      0.9977     1.0018
  10000.0      0.9998     1.0010
  50000.0      1.0000     1.0008

Above 10 keV the two agree to 1e-3 in both quantities (inside the 1 % requirement of the issue); both reach the Bethe limit, which fixes S through the f-sum and the mean excitation energy of the optical ELF, and their IMFP constants differ only by the single-electron continuum, a 1e-3 effect at these energies. At 100 eV to 1 keV the differences are at most 4 % (S) and 2 % (λ). Below about 100 eV they are large and of both signs: at 20 to 30 eV the full model’s stopping power and inverse IMFP are several times those of the single pole (the full model has single-electron excitations that the single pole does not, which is presumably the cause; it was not isolated), while at 50 eV the single pole is the larger. That the full-algorithm IMFP is smaller than the single-pole one at low energies is the direction reported for water in S2017 (section on Fig. 2); the sign change near 50 eV here depends on the shape of the synthetic ELF, and no statement about a real material follows. Nothing was tuned to make the two agree: these are the numbers of the implementation as it stands, and the regression bounds of tests/penn_full.rs follow them.

Structs

DiimfpEvaluation

One evaluation of the full Penn DIIMFP (FullPenn::diimfp_evaluation): the value and how many of the loss-function evaluations of its momentum integral did not converge (#368).

DiimfpGrid

The full Penn DIIMFP of one table build, tabulated so that each row of the table costs interpolations instead of nested integrals (issue #256).

FullPenn

The full Penn model built from one optical ELF. See the module docs.

ImfpEvaluation

One evaluation of the full Penn inverse IMFP and stopping power (FullPenn::imfp_and_stopping_evaluation): the point and the cells of its outer ω_p integral accepted without meeting their tolerance (#298, #368). The inner integrals over q and ω are not bounded and report nothing.

LossFunctionEvaluation

One evaluation of the full Penn loss function Im[-1/ε(q, ω)] (FullPenn::loss_function_evaluation): the value, the work of its adaptive ω_p integral, and whether every cell of that integral met its tolerance.

Constants

DEFAULT_FULL_TOLERANCE

Default relative tolerance of the integrals of the full Penn model.

OMEGA_P_INTERVAL_BUDGET

Budget of Gauss-Legendre rule applications (5 evaluations each) for the bisection of one knot interval of the ω_p integral of the loss function (#298); see LossFunctionEvaluation.