lindhard/electron/inelastic/sum_rules.rs
1//! Sum rules of an optical energy-loss function.
2//!
3//! Two checks of the internal consistency of an optical ELF, from Shinotsuka,
4//! Da, Tanuma, Yoshikawa, Powell, Penn, Surf. Interface Anal. 49, 238 (2017),
5//! doi:10.1002/sia.6123 (open copy: PMC5524379), in Hartree atomic units
6//! there and in eV here:
7//!
8//! - **f-sum** (eq. (18)): the effective number of electrons per target unit
9//! (atom or molecule) taking part in losses up to `W`,
10//! `N_eff(W) = (2 / (π (ħΩ_a)²)) ∫_0^W W' ELF(W') dW'`, with
11//! `(ħΩ_a)² = ħ² n_a e² / (ε₀ m_e)` for the target-unit density `n_a`.
12//! Without a target density, the report gives the effective electron
13//! density `n_eff(W) = (ε₀ m_e / ħ²) (2/π) ∫_0^W W' ELF dW'` instead
14//! (`N_eff = n_eff / n_a`).
15//! - **Perfect-screening (Kramers-Kronig) sum** (eq. (19)):
16//! `P_eff = (2/π) ∫ ELF(W)/W dW`. Eq. (19) adds `n(0)⁻²`, the inverse square
17//! of the static refractive index, which is zero for a metal and is not
18//! known from the ELF alone; the report gives the integral only, and the
19//! caller adds `n(0)⁻²` for an insulator. For a complete ELF `P_eff -> 1`
20//! (text below eq. (19)).
21//!
22//! The report also gives the mean excitation energy `I` from the same ELF,
23//! `ln I = ∫ W ELF ln W dW / ∫ W ELF dW` (the definition the Bethe limit uses;
24//! see [`super::penn`]).
25//!
26//! # Integration
27//!
28//! The integrals run over the tabulated range only (nothing is extrapolated;
29//! `OpticalElf` interpolates linearly and the ELF is taken as zero outside the
30//! table). They are exact for the piecewise-linear interpolant: the f-sum by
31//! Simpson's rule (exact for the quadratic `W ELF(W)` on each segment), `P_eff`
32//! in closed form (`∫ (α + βW)/W dW = α ln(W₁/W₀) + β (W₁ - W₀)`), and the
33//! `ln W`-weighted integral by an 8-point Gauss-Legendre rule per segment
34//! (the integrand is smooth on each segment; relative error below 1e-12).
35
36use std::fmt;
37
38use super::quadrature::GaussLegendre;
39use crate::constants::{ELECTRON_MASS, ELEMENTARY_CHARGE, HBAR, VACUUM_PERMITTIVITY};
40use crate::electron::data::{ElectronDataError, OpticalElf};
41
42/// Electron density (m⁻³) whose free-electron plasma energy squared is
43/// `x_ev2` (eV²): `n = (ħΩ)² ε₀ m_e / (ħ² e²)` with `(ħΩ)²` in J².
44fn density_from_plasma_energy_sq(x_ev2: f64) -> f64 {
45 let x_j2 = x_ev2 * ELEMENTARY_CHARGE * ELEMENTARY_CHARGE;
46 x_j2 * VACUUM_PERMITTIVITY * ELECTRON_MASS
47 / (HBAR * HBAR * ELEMENTARY_CHARGE * ELEMENTARY_CHARGE)
48}
49
50/// The sum rules of one [`OpticalElf`]. See the module docs for the
51/// definitions. Built by [`SumRuleReport::new`]; printed by its `Display`
52/// implementation.
53#[derive(Debug, Clone, PartialEq)]
54pub struct SumRuleReport {
55 /// Material identity of the ELF.
56 pub material: String,
57 /// Provenance of the ELF.
58 pub provenance: String,
59 /// The integration range, the tabulated range of the ELF, eV.
60 pub energy_range_ev: (f64, f64),
61 /// `∫ W ELF(W) dW` over the table, eV².
62 pub f_sum_ev2: f64,
63 /// The plasma energy of all the electrons the ELF accounts for,
64 /// `ħΩ_p = sqrt((2/π) ∫ W ELF dW)`, eV.
65 pub plasma_energy_ev: f64,
66 /// The effective electron density `n_eff` of the whole table, m⁻³.
67 pub electron_density_per_m3: f64,
68 /// `P_eff = (2/π) ∫ ELF(W)/W dW` over the table (no `n(0)⁻²` term).
69 pub p_eff: f64,
70 /// Mean excitation energy `I`, eV: `ln I = ∫ W ELF ln W dW / ∫ W ELF dW`.
71 /// `NaN` if the f-sum is zero.
72 pub mean_excitation_energy_ev: f64,
73 /// The target-unit (atom or molecule) density used for `N_eff`, m⁻³, if
74 /// one was given.
75 pub target_density_per_m3: Option<f64>,
76 /// The tabulated energies `W_k`, eV.
77 pub energy_ev: Vec<f64>,
78 /// `n_eff(W_k)`, the cumulative effective electron density up to each
79 /// tabulated energy, m⁻³.
80 pub cumulative_electron_density_per_m3: Vec<f64>,
81 /// The tabulated ELF values, for partial segments.
82 elf_values: Vec<f64>,
83}
84
85impl SumRuleReport {
86 /// The sum rules of `elf`, without a target density (`N_eff` is then not
87 /// available, `n_eff` is).
88 pub fn new(elf: &OpticalElf) -> Self {
89 let w = elf.energy_ev();
90 let e = elf.elf_values();
91 let gl = GaussLegendre::new(8);
92 let mut f_cum = Vec::with_capacity(w.len());
93 f_cum.push(0.0);
94 let (mut f_sum, mut p_sum, mut ln_sum) = (0.0, 0.0, 0.0);
95 for k in 0..w.len() - 1 {
96 let (w0, w1, e0, e1) = (w[k], w[k + 1], e[k], e[k + 1]);
97 let h = w1 - w0;
98 let wm = 0.5 * (w0 + w1);
99 let em = 0.5 * (e0 + e1);
100 f_sum += h / 6.0 * (w0 * e0 + 4.0 * wm * em + w1 * e1);
101 f_cum.push(f_sum);
102 let beta = (e1 - e0) / h;
103 p_sum += (e0 - beta * w0) * (w1 / w0).ln() + (e1 - e0);
104 ln_sum += gl.integrate(&mut |x: f64| [x * (e0 + beta * (x - w0)) * x.ln()], w0, w1)[0];
105 }
106 let two_over_pi = 2.0 / std::f64::consts::PI;
107 let mean_excitation_energy_ev = if f_sum > 0.0 {
108 (ln_sum / f_sum).exp()
109 } else {
110 f64::NAN
111 };
112 Self {
113 material: elf.material().to_string(),
114 provenance: elf.provenance().to_string(),
115 energy_range_ev: elf.energy_range_ev(),
116 f_sum_ev2: f_sum,
117 plasma_energy_ev: (two_over_pi * f_sum).sqrt(),
118 electron_density_per_m3: density_from_plasma_energy_sq(two_over_pi * f_sum),
119 p_eff: two_over_pi * p_sum,
120 mean_excitation_energy_ev,
121 target_density_per_m3: None,
122 energy_ev: w.to_vec(),
123 cumulative_electron_density_per_m3: f_cum
124 .iter()
125 .map(|f| density_from_plasma_energy_sq(two_over_pi * f))
126 .collect(),
127 elf_values: e.to_vec(),
128 }
129 }
130
131 /// The sum rules of `elf` with the density of target units (atoms or
132 /// molecules, m⁻³), so that `N_eff` per unit is available. The density
133 /// must be finite and positive.
134 pub fn with_target_density(
135 elf: &OpticalElf,
136 target_density_per_m3: f64,
137 ) -> Result<Self, ElectronDataError> {
138 if !(target_density_per_m3.is_finite() && target_density_per_m3 > 0.0) {
139 return Err(ElectronDataError::Invalid {
140 what: "target density",
141 reason: format!("must be finite and positive, got {target_density_per_m3} m^-3"),
142 });
143 }
144 let mut r = Self::new(elf);
145 r.target_density_per_m3 = Some(target_density_per_m3);
146 Ok(r)
147 }
148
149 /// `n_eff(W)`, m⁻³: the effective electron density for losses up to `W`
150 /// (eV), exact for the interpolant. Zero below the table and the whole
151 /// table's value above it.
152 pub fn electron_density_up_to(&self, w_ev: f64) -> f64 {
153 let w = &self.energy_ev;
154 let c = &self.cumulative_electron_density_per_m3;
155 if w_ev.is_nan() || w_ev <= w[0] {
156 return 0.0;
157 }
158 if w_ev >= w[w.len() - 1] {
159 return c[c.len() - 1];
160 }
161 let i = w.partition_point(|&x| x <= w_ev).clamp(1, w.len() - 1);
162 let (w0, w1) = (w[i - 1], w[i]);
163 let (e0, e1) = (self.elf_values[i - 1], self.elf_values[i]);
164 // ∫_{w0}^{x} W ELF(W) dW with the segment's linear ELF, by Simpson's
165 // rule (exact for the quadratic integrand).
166 let x = w_ev;
167 let ex = e0 + (e1 - e0) * (x - w0) / (w1 - w0);
168 let (xm, em) = (0.5 * (w0 + x), 0.5 * (e0 + ex));
169 let part = (x - w0) / 6.0 * (w0 * e0 + 4.0 * xm * em + x * ex);
170 c[i - 1] + density_from_plasma_energy_sq(2.0 / std::f64::consts::PI * part)
171 }
172
173 /// `N_eff(W)`: effective electrons per target unit for losses up to `W`
174 /// (eV), if a target density was given.
175 pub fn effective_electrons_up_to(&self, w_ev: f64) -> Option<f64> {
176 self.target_density_per_m3
177 .map(|n| self.electron_density_up_to(w_ev) / n)
178 }
179
180 /// `N_eff` of the whole table, if a target density was given.
181 pub fn effective_electrons(&self) -> Option<f64> {
182 self.target_density_per_m3
183 .map(|n| self.electron_density_per_m3 / n)
184 }
185}
186
187impl fmt::Display for SumRuleReport {
188 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
189 let (lo, hi) = self.energy_range_ev;
190 writeln!(f, "Sum rules of the optical ELF of {}", self.material)?;
191 writeln!(f, " provenance: {}", self.provenance)?;
192 writeln!(f, " integration range: {lo} to {hi} eV")?;
193 writeln!(f, " f-sum, int W ELF dW: {:.6e} eV^2", self.f_sum_ev2)?;
194 writeln!(
195 f,
196 " plasma energy (f-sum): {:.6} eV",
197 self.plasma_energy_ev
198 )?;
199 writeln!(
200 f,
201 " effective electron density n_eff: {:.6e} m^-3",
202 self.electron_density_per_m3
203 )?;
204 if let (Some(n), Some(neff)) = (self.target_density_per_m3, self.effective_electrons()) {
205 writeln!(f, " target-unit density: {n:.6e} m^-3")?;
206 writeln!(f, " N_eff per target unit: {neff:.6}")?;
207 }
208 writeln!(
209 f,
210 " P_eff = (2/pi) int ELF/W dW: {:.6} (add n(0)^-2 for an insulator)",
211 self.p_eff
212 )?;
213 writeln!(
214 f,
215 " mean excitation energy I: {:.6} eV",
216 self.mean_excitation_energy_ev
217 )
218 }
219}