Module elastic 

Elastic scattering of electrons by a static central potential: the radial Dirac equation, partial-wave phase shifts and the Mott cross sections.

This is step 1 of the Mott elastic model (issue #17): phase shifts, the differential cross section DCS(theta), the total elastic cross section sigma_el, the first transport cross section sigma_tr1 and the Sherman function, all at one kinetic energy and for a potential given as a function. Optional Furness-McCarthy exchange and correlation-polarization corrections (issue #91), off by default, are in corrections (solve_corrected). Energy-grid tables for a material are in table (issue #90); condensed phases are a follow-up issue.

Units

Hartree atomic units inside this module: lengths in bohr, energies in hartree, the electron mass is 1 and the speed of light is c = 1/alpha (alpha from crate::constants::FINE_STRUCTURE). Cross sections come out in bohr^2 (DCS in bohr^2/sr); BOHR2_TO_M2 converts. The public energy argument is the electron kinetic energy in eV.

Method

With P(r) = r g(r) and Q(r) = r f(r) the large and small radial components, the radial Dirac equations for a central potential energy V(r) and Dirac quantum number kappa are (kinetic energy E, c = 1/alpha)

dP/dr = -kappa/r P + (E - V + 2c^2)/c Q
dQ/dr = -(E - V)/c P + kappa/r Q

(the form used in the ELSEPA paper, Salvat, Jablonski & Powell, Comput. Phys. Commun. 165, 157 (2005); here re-derived from the Dirac equation and checked against the closed-form free solution below, not copied from any code). kappa = -(l+1) is the spin-up channel j = l + 1/2 and kappa = l the spin-down channel j = l - 1/2 (l >= 1).

• Free solutions. For V = 0 the regular and irregular solutions are r (u_l(kr), sigma s u_lbar(kr)), with u a spherical Bessel j or y, k = sqrt(E (E + 2c^2))/c, s = c k/(E + 2c^2), sigma = -1 and lbar = l + 1 for kappa < 0, sigma = +1 and lbar = l - 1 for kappa > 0. (Verified by substitution in the equations above and by the tests.)
• Integration. The regular solution is integrated outward from a small radius r_s. Inside the centrifugal barrier the irregular solution decays outward, so an imprecise starting ratio Q/P is harmless as long as r_s is well below the classical turning point: r_s is placed by a WKB criterion on the bisected turning point, so that the irregular admixture is suppressed by exp(-60) there (see wkb_start), but never below a floor of 1e-10/(1+Z) bohr, where the leading-order point-Coulomb or regular power-series start is accurate. (Until #131 the static potentials used r_s = r_t exp(-60/|kappa|), which assumes power-law growth; measured against the WKB start for Au, Salvat DHFS screening, its phase-shift error was below 2e-12 rad at 10 and 20 keV, 6.3e-10 rad at 30 keV, 1.6e-7 rad at 50 keV and 2.5e-4 rad at 100 keV (|kappa| 779, sigma_tr1 off by 2.4e-5 relative).) The integrator is the Gragg-Bulirsch-Stoer modified-midpoint method with polynomial extrapolation in h^2 (Stoer & Bulirsch, Introduction to Numerical Analysis), substep sequence 2, 4, 6, 8 (order 8), with a step set from the local eigenvalue of the system, h = step_scale / (sqrt|q| + (1 + (2 kappa^2)^(1/3))/r) (the cube-root term resolves the Airy region at a turning point), capped at a quarter of step_scale times the potential’s length scale, q = (E-V)(E-V+2c^2)/c^2 - kappa^2/r^2. The amplitude is rescaled when it grows past 1e100; only the ratio Q/P matters.
• Matching. At the matching radius R, where |r V(r)| has fallen below SolverOptions::matching_threshold (a potential that has a finite range ends there exactly), tan(delta) = (Q J_p - P J_q) / (Q Y_p - P Y_q) with J, Y the regular and irregular free solutions. The potentials here are screened, so no Coulomb tail is matched. The neglected tail changes delta by roughly the threshold divided by k.
• Partial-wave sum. f(theta) and g(theta) are the Mott direct and spin-flip amplitudes, f = (1/2ik) sum_l [(l+1)(e^{2 i d_-} - 1) + l (e^{2 i d_+} - 1)] P_l(cos theta), g = (1/2ik) sum_l [e^{2 i d_-} - e^{2 i d_+}] P_l^1(cos theta) with P_l^1 = sin theta dP_l/d cos theta, and DCS = |f|^2 + |g|^2. These are the standard textbook (Mott 1929; Motz, Olsen & Koch 1964) expressions written from memory of the form, not checked digit for digit against a paper in this work; the tests that pin them are the optical theorem, the exact 1 - beta^2 sin^2(theta/2) first-Born limit (which fixes the relative size of g) and the square-well closed form. The Sherman function is S = i (f g* - f* g) / (|f|^2 + |g|^2); its sign convention is not verified against a published definition.
• Cutoff in l. Phase shifts are computed for l = 0, 1, 2, ... and the series is cut after consecutive (5) successive l for which both |delta_-| and |delta_+| are below phase_tolerance (1e-8 rad); PartialWaves::l_max is the highest l computed. When x = kR is below l by enough that the free Bessel functions leave the range of f64, the wave never reaches the potential and its shift is set to zero. At high energy l_max grows into the thousands; this module integrates every partial wave directly (cost linear in l_max times k R). A WKB / asymptotic phase for large l, as in the literature on which the issue draws (Salvat et al. 2005), is not implemented: that paper was not opened for this work (closed access) and the method is not reproduced from memory. See follow-up work.
• Cross sections. sigma_el = (4 pi/k^2) sum_l [(l+1) sin^2 d_- + l sin^2 d_+] (equal to (4 pi/k) Im f(0), the optical theorem). The DCS times (1 - cos theta) is a polynomial of degree at most 2 l_max + 1 in cos theta, so sigma_tr1 is integrated exactly with an (l_max + 2)-point Gauss-Legendre rule.

Potentials

Yukawa and SquareWell are analytically tractable. SalvatDhfs is the analytic Dirac-Hartree-Fock-Slater screening function of Salvat, Martinez, Mayol & Parellada, Phys. Rev. A 36, 467 (1987), doi:10.1103/PhysRevA.36.467, with its Table I coefficients for Z = 1..92 (transcribed from the open repository copy of the paper, #130), see the type and docs/data-provenance.md.

Determinism

Each phase shift is an independent, deterministic computation. Blocks of partial waves are evaluated in parallel with rayon and collected in order, so the result is bit-identical on any thread count (tested).

Re-exports

pub use corrections::solve_corrected;

pub use corrections::CorrectedPotential;

pub use corrections::CorrectionMetadata;

pub use corrections::Corrections;

pub use corrections::CorrelationPolarization;

pub use corrections::ElectronDensity;

pub use corrections::PolarizationCutoff;

Modules

corrections

Optional corrections to the static potential: local exchange and correlation-polarization (issue #91, step 3 of the Mott elastic model).

table

Energy-grid elastic tables: the radial-Dirac solver of the parent module run over a log energy grid for every element of a Material, stored as an electron::data::CrossSectionTable with SamplingAxis::ElasticPolarAngle (issue #90, step 2 of the Mott elastic model).

Structs

ElasticResult

Result of solve: everything the issue asks for at one energy.

ElasticSolver

Radial Dirac solver at one kinetic energy.

PartialWaves

Mott partial-wave data at one energy: phase shifts and the amplitudes built from them.

SalvatDhfs

The three-term analytic screening function of Salvat, Martinez, Mayol & Parellada, Phys. Rev. A 36, 467 (1987), doi:10.1103/PhysRevA.36.467: V(r) = -(Z/r) sum_i A_i exp(-alpha_i r) with sum_i A_i = 1 (neutral atom), r in bohr.

SolverOptions

Numerical settings of the solver.

SquareWell

Spherical square well V = -depth for r < radius, 0 beyond (hartree, bohr). A negative depth is a barrier.

Yukawa

Exponentially screened Coulomb potential V(r) = -Z exp(-r/a)/r (hartree, bohr). Its first-Born cross sections are analytic.

Enums

ElasticError

Errors from the elastic solver.

Constants

BOHR2_TO_M2

Square metres per bohr squared (a0^2, BOHR_RADIUS squared).

HARTREE_EV

Hartree energy in eV (CODATA 2022 E_h over e).

Traits

ScreenedPotential

A static, spherically symmetric potential energy of the electron.

Functions

solve

Solve at one energy and evaluate the observables on a caller-supplied angle grid (thetas, rad).