Physics and Numerical Formulation
Liouville Space Representation
The spectral line-shape intensity \(I(\omega)\) of emitted radiation is related to the Fourier transform of the dipole autocorrelation function \(C(t)\) by [baranger]:
In Liouville space notation, the autocorrelation function is a trace over the quantum emitter states:
where \(\vec{d}\) is the electric transition dipole operator, \(\rho_0\) is the equilibrium density matrix of the initial manifold, and \(U(t)\) is the time-evolution propagator.
In the anti-symmetric subspace (ground-to-excited manifold transitions in the no-quenching approximation), the Liouvillian operator \(L\) is constructed from the excited-state Hamiltonian \(H_u\) and ground-state Hamiltonian \(H_l\):
The diagonal elements of \(L\) give the field-free transition frequencies \(\omega_{ij} = (E_i^u - E_j^l)/\hbar\); off-diagonal elements encode the Stark and Zeeman couplings between dressed states.
From propagator to resolvent.
The time-evolution propagator is generated by \(L\):
Substituting into Eq. (2) gives the explicit correlation function
Inserting this into Eq. (1) and evaluating the one-sided Fourier transform converts the time-domain exponential into a frequency-domain resolvent:
Fast electron collisions are encoded as an additional imaginary damping operator \(\Phi(\omega)\) in the same Liouville space (Section 3.8 ), shifting \(L \to L + i\Phi(\omega)\):
Field-dependent Liouvillian and static profile.
Under the quasi-static ion approximation the radiator is exposed to a static ionic microfield \(\vec{F}\) at angle \(\theta\) to \(\vec{B}\). Both manifold Hamiltonians \(H_u\) and \(H_l\) then depend on the field magnitude \(F\) and orientation \(\mu = \cos\theta\) through the Stark perturbation \(V_E\) (Section 3.4 ), making the Liouvillian parametrically field-dependent:
where \(H_A\) is the field-free magnetic Hamiltonian (Section 3.2 ) and \(V_E(F,\mu)\) is the linear Stark interaction (Section 3.4 ). Equations (7) and (8) together give the field-dependent \(q\)-polarized profile at a given configuration \((F,\mu)\):
The total \(q\)-polarized profile is then the average over the ionic microfield distribution \(W(F)\) and orientation:
Because \(I_q(\omega,F,\mu) = I_q(\omega,F,-\mu)\) (the Hamiltonian is even in \(\mu\) via \(F_z = F\mu\) and \(F_x = F\sqrt{1-\mu^2}\)), the integral over \([-1,1]\) reduces to twice the \([0,1]\) range, discretized by Gauss–Legendre quadrature. The observed intensity at angle \(\alpha\) between the line of sight and \(\vec{B}\) is:
Sections 3.2 –3.8 detail each ingredient of Eq. (9); Section 3.10 extends this framework to include ion dynamics via the FFM.
Radiator Hamiltonian in Magnetic Fields
For each principal quantum number shell \(n\), the atomic Hamiltonian \(H_A\) is constructed in the canonical uncoupled hydrogenic basis \(|n, l, m_l, s{=}\tfrac{1}{2}, m_s\rangle\) of dimension \(2n^2\). Basis states are ordered by an outer loop over \(l \in [0, n-1]\), an intermediate loop over \(m_l \in [-l, l]\), and an inner loop over \(m_s \in \{+\tfrac{1}{2}, -\tfrac{1}{2}\}\). Five physical contributions are assembled into a single matrix [bethesalpeter]:
Unperturbed Hydrogenic Energy (\(H_0\))
The diagonal unperturbed energy is corrected for the finite nuclear mass:
where \(M_{\rm nuc}\) is the nuclear mass of the emitting isotope (fixed by its atomic mass number) and \(\mathrm{Ry}_{\rm red}\) the corresponding reduced-mass Rydberg constant. When use_empirical_data=True this analytic diagonal is replaced by NIST tabulated energies (Section 3.3 ).
Spin-Orbit Coupling (\(H_{\rm SO}\))
The spin-orbit interaction \(H_{\rm SO} = \xi_{nl}\,\vec{L}\cdot\vec{S}\) couples orbital and spin angular momenta. The coupling constant for \(l > 0\) is:
Using the ladder decomposition \(\vec{L}\cdot\vec{S} = L_z S_z + \tfrac{1}{2}(L_+ S_- + L_- S_+)\), the matrix elements within a given \(l\)-subshell are:
Diagonal (\(\Delta m_l = 0,\,\Delta m_s = 0\)):
\[\langle l, m_l, m_s | H_{\rm SO} | l, m_l, m_s\rangle = \xi_{nl}\,m_l\,m_s\]Off-diagonal (\(\Delta m_l = \pm 1,\,\Delta m_s = \mp 1\)): for \(s = \tfrac{1}{2}\) the spin ladder factor equals 1 for all allowed transitions:
\[\langle l, m_l \pm 1, m_s \mp 1 | H_{\rm SO} | l, m_l, m_s\rangle = \tfrac{1}{2}\,\xi_{nl}\sqrt{l(l+1) - m_l(m_l \pm 1)}\]
Fine-Structure Relativistic Corrections (\(H_{\rm FS}\))
To restore the Dirac degeneracy (e.g. \(2s_{1/2}\) and \(2p_{1/2}\)), the mass-velocity and Darwin terms are added [bethesalpeter]. These are diagonal in \(|n, l, m_l, m_s\rangle\) and independent of \(m_l\) and \(m_s\):
where \(A_{\rm FS} = Z^4\alpha^2\,\mathrm{Ry}_{\infty}/n^4\). The algebraic identity \(H_{\rm SO} + H_{\rm FS}\) reproduces the exact Dirac fine-structure eigenvalues for both \(j = l \pm \tfrac{1}{2}\).
Linear Zeeman Effect (\(H_Z^{(1)}\))
where \(\mu_B \approx 5.788\times10^{-5}\) eV/T and \(g_s \approx 2.0023192\).
Quadratic (Diamagnetic) Zeeman Effect (\(H_Z^{(2)}\))
This operator couples states with \(\Delta l = 0\) and \(\Delta l = \pm 2\). The matrix elements factor into radial and angular parts:
Radial integrals:
Diagonal (\(l_1 = l_2 = l\)): \(\displaystyle\langle n, l | r^2 | n, l \rangle = \frac{n^2}{2Z^2}\bigl[5n^2+1-3l(l+1)\bigr]\;[a_0^2]\)
Off-diagonal (\(l_2 = l_1 \pm 2\)): evaluated numerically via \(\int_0^\infty R_{nl_1}\,r^4\,R_{nl_2}\,dr\), cached with
lru_cache.
Angular integrals:
Diagonal: for \(l=0\), \(\langle\sin^2\theta\rangle = 2/3\); for \(l>0\),
\[\langle l, m_l | \sin^2\theta | l, m_l \rangle = 1 - \frac{2l^2 + 2l - 1 - 2m_l^2}{(2l-1)(2l+3)}\]Off-diagonal (\(l_2 = l_1 + 2\), \(l_{\rm low} = \min(l_1,l_2)\)): the constant term in \(1 - \cos^2\theta\) vanishes for orthogonal harmonics, leaving:
\[\langle l_1, m_l | \sin^2\theta | l_2, m_l \rangle = -\sqrt{\frac{\bigl[(l_{\rm low}+1)^2 - m_l^2\bigr] \bigl[(l_{\rm low}+2)^2 - m_l^2\bigr]} {(2l_{\rm low}+1)(2l_{\rm low}+3)^2(2l_{\rm low}+5)}}\]
Empirical Field-Free Energies
The analytic diagonal \(H_0 = -Z^2\,\mathrm{Ry}_{\mathrm{red}}/n^2\), together with the spin-orbit and fine-structure terms, reproduces the field-free level positions only to the accuracy of the hydrogenic Dirac formula. For quantitative line centers, StarkZee can instead inject measured field-free energies through the use_empirical_data=True flag (with the element symbol atom). The values are read from a tabulated database (atomic_data.load_levels) holding NIST level energies [nist] in wavenumbers [cm\(^{-1}\)], resolved by \((n, l, j)\). The database is stored in starkzee/data/atomic_levels.json; currently only hydrogen ("H") is included. Each top-level key is an atom symbol; the value is an object with two sections: "fine_structure_true" (entries with fields n, l, j, energy) and "fine_structure_false" (entries with fields n, energy). All energies are in cm\(^{-1}\) above the ground state (NIST vacuum values). Additional atoms can be supported by adding the corresponding entry to this file. To refresh the bundled data from the NIST ASD levels query, run python scripts/update_nist_levels.py H --spectrum "H I". The updater is a development tool: it writes the local JSON database and checks that each requested fine-structure shell contains all hydrogenic (l, j) states before replacing the data. The empirical Hamiltonian is kept in cm\(^{-1}\) throughout this code path so that callers can work entirely in wavenumber units; all Zeeman contributions are also converted from eV to cm\(^{-1}\) before being added, keeping every term on the same scale.
Because the field-free Hamiltonian is degenerate (the Dirac formula makes \(2s_{1/2}\) and \(2p_{1/2}\) coincide), the empirical energies cannot simply be written on the diagonal of the uncoupled \(|n,l,m_l,m_s\rangle\) basis: numpy.linalg.eigh is free to mix the degenerate \(l\) states arbitrarily. StarkZee therefore:
diagonalizes a degeneracy-broken field-free Hamiltonian (unperturbed energy \(+\) spin-orbit only, omitting the mass-velocity/Darwin term so that \(l\) remains a good label), yielding eigenvectors \(V\);
labels each coupled eigenstate \(k\) by its dominant orbital component \(l\) and its total angular momentum \(j = l \pm \tfrac{1}{2}\), the branch being set by the sign of \(\langle \vec{L}\cdot\vec{S}\rangle = (E_k - E_n)/\xi_{nl}\);
assigns the tabulated energy \(D_k = E^{\mathrm{emp}}(l, j)\) [cm\(^{-1}\)] to each eigenstate and reconstructs the field-free Hamiltonian as
\[H_0^{\mathrm{emp}} = V\,\mathrm{diag}(D)\,V^\dagger \quad [\mathrm{cm}^{-1}].\]
The Zeeman terms (\(\mu_B B (m_l + g_s m_s)\) and the diamagnetic correction) are converted from eV to cm\(^{-1}\) and added on top of \(H_0^{\mathrm{emp}}\). Since the tabulated values are absolute level energies (ground state at \(0\)), transition wavenumbers follow directly as differences \(E_u - E_l\) and reproduce the observed NIST Lyman/Balmer line centers.
Full Electron-Radiator Hamiltonian
Under the quasi-static ion approximation, the radiator is subjected to a constant electric microfield \(\vec{F}\) at an angle \(\theta\) relative to the magnetic field \(\vec{B}\). The microfield is decomposed into longitudinal (\(F_z = F\cos\theta\)) and transverse (\(F_x = F\sin\theta\)) components. The Stark interaction is
where \(\hat{\boldsymbol{d}} = e\vec{r}\) is the electric dipole operator. The full electron-radiator Hamiltonian is
\(H_A\) and \(V_E\) are assembled into a single matrix and diagonalized as a whole by numpy.linalg.eigh at every microfield quadrature point (solve_starkzee). The Stark interaction is not treated perturbatively: there is no expansion in powers of \(V_E\). This exact treatment is essential when the Stark energy \(eF\langle r\rangle_n\) is comparable to the fine-structure or Zeeman splittings, yielding the Stark-dressed eigenstates \(|k\rangle\) and transition frequencies \(\omega_k\).
Stark Perturbation Matrix Elements
The selection rule for the within-shell Stark coupling is \(\Delta n = 0\), \(|l_i - l_j| = 1\), \(\Delta m_s = 0\). In energy units [eV]:
Radial matrix element.
For \(l = \max(l_i, l_j)\), Gordon’s formula gives the within-shell radial dipole element exactly:
Angular matrix elements.
The angular factor is the spherical dipole tensor, written \(z/r = T_0^{(1)}\) and \(x/r = (T_{-1}^{(1)} + T_{+1}^{(1)})/\sqrt{2}\). Its non-zero matrix elements (with \(l = \max(l_1, l_2)\)) are:
These same three matrix elements reappear unchanged in the transition dipole operator (Section 3.5 ), since both the Stark coupling and the radiative dipole are built from the rank-1 tensor \(T_q^{(1)}\). The Stark matrices \(M_z\) and \(M_x\) [eV/(V/m)] are precomputed once per \((n, Z)\) and cached (_stark_templates), so that \(V_E = F_z M_z + F_x M_x\) is formed at each quadrature point by a simple scalar multiplication.
Transition Dipole Operator
The transition dipole matrices \(D_q\) (\(q \in \{0,\pm1\}\)) are built in the uncoupled basis \(|n, l, m_l, s, m_s\rangle\). Selection rules require \(\Delta l = \pm 1\), \(\Delta m_l = q\), and \(\Delta m_s = 0\):
Angular Dipole Matrix Elements
The angular factors \(\langle l_l, m_{l,l} | T_q^{(1)} | l_u, m_{l,u} \rangle\) follow from the Wigner-Eckart theorem [edmonds]. They are exactly the rank-1 spherical-tensor matrix elements already written for the Stark coupling in Eqs. (13)–(15), now evaluated between the lower and upper manifolds: \(q=0\) gives the \(\pi\) component (\(\Delta m_l = 0\)), \(q=+1\) the \(\sigma^-\) component (\(\Delta m_l = -1\)), and \(q=-1\) the \(\sigma^+\) component (\(\Delta m_l = +1\)).
Hydrogenic Radial Dipole Integrals (Gordon’s Formula)
For \(n_u \neq n_l\), the inter-shell radial elements \(\langle n_l, l_l | r | n_u, l_u \rangle\) are evaluated exactly using Gordon’s (1929) analytical formula [gordon]:
where \(L = \max(l_u, l_l)\), and \(F_1\), \(F_2\) are terminating Gauss hypergeometric functions:
Because the first argument in both hypergeometric functions is a non-positive integer, each series terminates as a finite sum:
where \((-m)_k\) is the Pochhammer symbol. This exact evaluation carries no integration cutoff error and is cached for reuse.
Diagonalization and Transition Dipole Basis Rotation
At each quadrature point \((F, \mu)\) the total Hamiltonians for both manifolds are:
Diagonalization via numpy.linalg.eigh gives eigenvectors \(V_u\), \(V_l\) and eigenvalues \(E_u^p\), \(E_l^k\):
The transition energy for each dressed-state pair \((p \to k)\) is \(E_{pk} = E_u^p - E_l^k\). The uncoupled dipole matrices \(D_q\) are rotated into the mixed eigenstate bases:
The weight \(S_q(p,k)\) and energy \(E_{pk}\) together define one Stark-Dressed Transition (SDT). The static profile at field configuration \((F,\mu)\) — Eq. (9) — is then realized in practice as a sum over all SDTs weighted by the quadrature weight \(W_{ij}\):
where \(\gamma_e\) is the electron-impact half-width of Section 3.8 .
Plasma Microfield Distributions
The electric microfield is averaged over a probability distribution \(W(F)\). The presence of \(\vec{B}\) breaks spherical symmetry, so the integration must be performed over both the field magnitude \(F = |\vec{F}|\) and the orientation \(\mu = \cos\theta\) between \(\vec{F}\) and \(\vec{B}\).
The Hooper Screened Microfield Distribution
To account for Debye screening of the ion Coulomb fields by surrounding plasma electrons, StarkZee implements the Hooper microfield distribution [hooper]:
where:
\(\beta = F / F_0\) is the reduced electric field strength.
\(F_0\) is the Holtsmark normal field strength, the Coulomb field of a single perturber at the mean inter-particle distance \(r_e\):
\[F_0 = \frac{e}{4\pi\varepsilon_0 r_e^2}, \qquad r_e = \left(\frac{3}{4\pi N_e}\right)^{1/3}\]\(a = r_e / \lambda_D\) is the Debye screening parameter, where \(\lambda_D\) is the multi-species Debye length:
\[\lambda_D = \sqrt{\frac{\varepsilon_0 k_B T_e}{N_e e^2 \bigl(1 + \sum_i X_i Z_i^2\bigr)}}\]with \(X_i = N_i / N_e\) the fractional ion concentration.
\(T(y, a) = \exp\!\bigl(-y^{3/2} S(y,a)\bigr)\) is the screened characteristic function, with Debye screening factor:
\[S(y, a) = \left(1 + \frac{1.5\,a^2}{y^2}\right)^{-3/4}\]
Setting \(a = 0\) yields \(S = 1\) and \(T = \exp(-y^{3/2})\), recovering the unscreened Holtsmark distribution \(W_H(\beta)\) [holtsmark].
The Potekhin Microfield Distribution
The analytical fits of Potekhin et al. (2002) [potekhin] describe the cumulative distribution \(Q(\beta)\) and probability density \(P(\beta)\) for both screened (\(s > 0\)) and unscreened (\(s = 0\)) potentials, for neutral and charged radiators, incorporating the ion-ion coupling parameter \(\Gamma\):
where \(q_n\), \(A, a, \alpha, B, b, \gamma, c\) are tabulated functions of \(\Gamma\) and \(s\), and \(S_N\) is the normalization constant.
Quadrature Discretization
The double integral \(\iint W(F,\mu)\,dF\,d\mu\) (Eq. [eq:microfield_integral] ) is discretized as follows:
Field magnitude \(F\): A uniform grid of \(N_F\) points in reduced field space \(\beta_i \in [0, \beta_{\max}]\) (default \(\beta_{\max} = 10\)):
\[F_i = \beta_i F_0, \qquad W_{F,i} = W(\beta_i, a)\,\Delta\beta, \qquad \sum_i W_{F,i} = 1\]Field orientation \(\mu\): \(N_\mu\)-point Gauss-Legendre quadrature over \([0, 1]\). Legendre roots \(x_j \in [-1,1]\) and weights \(w_j\) are mapped as:
\[\mu_j = \frac{x_j + 1}{2}, \qquad W_{\mu,j} = \frac{w_j}{2}\]Combined weight: \(W_{ij} = W_{F,i}\,W_{\mu,j}\). Points with \(W_{ij} < 10^{-15}\) are skipped.
so that at each active point the field decomposes into the longitudinal and transverse components seen by the radiator,
which are passed directly to the Stark template (Section 3.4 ). The weight \(W(\beta_i, a)\) is drawn from whichever distribution the user selects — Holtsmark, Hooper, or Potekhin (above) — while the orientation quadrature is identical in every case.
Electron Impact Broadening
Fast-moving electrons are treated in the impact (completed-collision) approximation: each collision is completed in a time short compared to both the inverse linewidth and the mean inter-collision time. In this limit, successive collisions are statistically independent and their net effect on the line shape enters through a damping operator \(\Phi(\omega)\) in the Liouville-space resolvent (Eq. [eq:resolvent_broadened] ).
Liouville-space structure of \(\Phi(\omega)\) and connection to \(\langle r^2\rangle\)
In the no-quenching approximation (electron collisions do not mix the upper and lower radiating manifolds), the operator \(\Phi\) in Liouville space factorizes as a sum of independent contributions from each manifold:
where \(\Phi_u\) and \(\Phi_l\) are self-energy operators acting within the upper (\(n_u\)) and lower (\(n_l\)) manifolds respectively.
The GBK semi-classical model [griembaranger] [griem1959] evaluates \(\Phi_n\) from the leading long-range term of the electron–radiator interaction, which is the dipole–dipole coupling \(V_{ee} \propto r_{\rm atom}^2 / r_e^3\). Averaging over the Maxwell–Boltzmann distribution of electron impact parameters and velocities, this yields a self-energy proportional to the operator \(r_n^2\) within shell \(n\):
where \(r_n^2\) is the mean-square-displacement operator within shell \(n\), \(W_0\) is the density-temperature prefactor (below), and \(G(\omega)\) is the frequency-dependent GBK factor (below). The \(r^2\) dependence is the direct quantum-mechanical trace of the dipole-coupling cross-section: a state with larger electronic extent couples more strongly to a passing electron.
Rotating to the SDT basis via the eigenvectors of \(H_n(F,\mu)\) (Section 3.6 ) makes \(\Phi_n\) approximately diagonal, since the dressed-state splitting \(\omega_{kk'}\) is generally large compared to the collision rate. The diagonal element for upper dressed state \(|k\rangle\) is:
The resolvent \((\omega\hat{I} - L(F,\mu) - i\Phi(\omega))^{-1}\) then has eigenvalues \((\omega - \omega_k - i\gamma_k)^{-1}\) in the SDT basis, and Eq. (9) reduces directly to the Lorentzian sum of Eq. (22).
Shell-average approximation.
Computing the per-SDT matrix element \(\langle k | r_u^2 | k \rangle\) exactly requires rotating the \(r^2\) operator to each dressed-state basis at every quadrature point (electron_operator=True). The default approximation replaces \(\langle k | r_u^2 | k \rangle\) by the shell average \(\langle r^2 \rangle_{n_u}\), giving a single homogeneous half-width:
The lower-manifold contribution \(\gamma_k^{(l)}\) is retained in principle through Eq. (16); in practice, for Balmer lines (\(n_l=2\), \(\langle r^2\rangle_2 = 12\,a_0^2\)) it is a factor \({\sim}3\) smaller than the \(n_u=3\) upper term and is dominated by the latter.
Density-Temperature Prefactor
The \(\hbar/e\) factor converts rad s\(^{-1}\) to eV.
Shell-Averaged Mean-Square Radius
This expression scales as \(n^4/Z^2\) and is not replaced by an approximate constant.
Strong-Collision Constant
Values from Ferri, Peyrusse & Calisti [ferri], Table 1:
Dynamical GBK Factor
where \(\Delta E = E - E_{pk}\) [eV] is the detuning from the transition center, \(E_c = \hbar\omega_c\) [eV] is the cutoff energy, \(\mathrm{Ry}_\infty \approx 13.606\) eV is the Rydberg energy (\(e^2/2a_0\), the ionization energy of hydrogen), and \(E_1(x)\) is the first-order exponential integral function defined for \(x > 0\) by:
Cutoff frequency.
The collision integral diverges logarithmically at large impact parameters unless truncated at a maximum \(\rho_{\max}\), or equivalently at a lower cutoff \(\omega_c = v_{\rm th}/\rho_{\max}\) [griem1959]. StarkZee follows the formulation of Ferri, Peyrusse & Calisti [ferri]:
\(\omega_p = \sqrt{N_e e^2/(\varepsilon_0 m_e)}\) — electron plasma frequency: Debye screening suppresses interactions at \(\rho > \lambda_D = v_{\rm th}/\omega_p\).
\(\omega_e = v_{\rm th}/r_e\) — configuration-change frequency: the impact approximation requires each collision to be completed before the surrounding electron configuration changes appreciably. At high density or low \(T_e\) this can be shorter than the Debye timescale (see Section 3.8.6 ).
\(\omega_L = eB/m_e\) — electron Larmor frequency: in \(\vec{B}\), electrons follow helical trajectories, reducing the effective interaction duration.
\(\omega_{\alpha\alpha'}\) — emitter level-splitting frequency: zero for hydrogen (degenerate \(l\)-subshells).
The largest frequency dominates, imposing the tightest bound on \(\rho_{\max}\).
Frequency-dependent vs. resonance-center width.
When frequency_dependent_width=True (default), \(\gamma_e(\Delta E)\) is evaluated at the actual detuning of each transition component using the full \(E_1(y(\Delta E))\) function. Setting this flag to False fixes the width at the line-center value \(\gamma_e(0)\), which reduces computation time at the cost of accuracy in the far wings.
This single half-width \(\gamma_e(\Delta E)\) is the only quantity the electron model hands to the profile builder; how it is applied to each Stark-dressed transition is described in Section 3.9 .
ZEST electron broadening variant (electron_model=’zest’)
Activating electron_model=’zest’ in LineProfile.compute_profile switches to the electron broadening model of the ZEST code [Gilleron2018]. The formula has the same structure as the default GBK expression but replaces the full shell-averaged \(\langle r^2 \rangle_n\) with the intra-shell mean-square radius:
Physical justification — \(\Delta n = 0\) interaction channels.
The ZEST paper (Section 2.1 of [Gilleron2018]) adopts two approximations that together restrict the electron broadening operator to within-shell matrix elements:
No-quenching approximation: perturbers may not induce transitions between the lower and upper manifolds participating in the radiator emission process.
\(\Delta n = 0\) interaction channels: “only states belonging to the same Layzer complex [same principal quantum number \(n\)] may be mixed by Stark and/or Zeeman effects.”
A Layzer complex is the set of all \(n^2\) states sharing principal quantum number \(n\). These two restrictions reduce the formal sum over all intermediate states in the broadening operator to within-shell dipole matrix elements only, directly yielding \(\langle r^2_\mathrm{intra}\rangle_n\) in place of the full \(\langle r^2\rangle_n\). The restriction is also automatic from the GBK dynamical factor: an inter-shell detuning (e.g. the \(n=3 \to n=2\) gap of \(\approx 1.89\) eV at the conditions of Fig. 2) drives \(G(\Delta\omega_\mathrm{inter}) \to 0\) exponentially, so inter-shell contributions vanish regardless of the explicit approximation.
Intra-shell squared radius.
Starting from the intra-shell radial element \(\langle n,l|r|n,l{\pm}1\rangle = \tfrac{3n}{2Z}\sqrt{n^2-(l\pm1)^2}\) and angular weight factors \(C(l,l+1)=(l+1)/(2l+1)\), \(C(l,l-1)=l/(2l+1)\), the per-\(l\) intra-shell sum is
This formula holds for all \(l \in [0, n-1]\), including the boundary cases \(l=0\) (no \(l-1\) neighbor) and \(l=n-1\) (no \(l+1\) neighbor). Averaging over the \(n^2\) spatial states with \((2l+1)\) weight yields the exact closed form:
Numerical values for H (\(Z=1\)): \(n=2\): 13.5 \(a_0^2\); \(n=3\): 81 \(a_0^2\); \(n=4\): 270 \(a_0^2\). These are smaller than the full shell averages (33, 153, 468 \(a_0^2\)) by factors of 0.41, 0.53, and 0.58, respectively.
Maximum wave-number cutoff \(\kappa_m\).
The PPPB model fixes the minimum impact parameter as \(\rho_\mathrm{min} = n^2 a_0/(2Z)\) (Bohr orbit) and uses a combined cutoff \(\omega_c = \max(\omega_p, \omega_e, \omega_L)\). The ZEST convention instead defines a temperature-dependent wave-number upper limit
The geometric branch \(Z/(n^2 a_0)\) dominates above \(T_e \approx \mathrm{Ry}_\infty \approx 13.6\) eV; below that the thermal de Broglie branch takes over. The only cutoff frequency retained in all ZEST G-functions is the electron plasma frequency \(\omega_p\) (no \(\omega_e\) or \(\omega_L\) terms), so \(B\) does not enter the ZEST broadening directly.
G-function choices (electron_model selector).
Three alternatives are available for the dynamical factor \(G(\Delta\omega)\) in Eq. (18), all sharing \(\langle r^2_\mathrm{intra}\rangle_n\), the same \(G_n = C_n\) strong-collision constants, and \(\kappa_m\) as defined above.
zest / zest-gbk (default ZEST). Closed-form exponential integral with a \(\kappa_m\)-derived argument:
This has the same \(E_1\) structure as the PPPB formula but with \(x = \kappa_m\lambda_D\) playing the role of \(\rho_\mathrm{min}/\lambda_D\), and only \(\omega_p\) in the additive term.
zest-lee. Analytical min-blend of the impact (\(\Delta\omega \to 0\)) and wing (\(\Delta\omega \to \infty\)) limits:
\(G_0\) (the line-center plateau) provides the plasma-cutoff floor; the \(E_1\) wing term contains only \(\Delta\omega^2\) (no \(+\omega_p^2\)) because the regularization is already captured by \(G_0\). At large detuning the two forms agree.
zest-dufty. Full random-phase-approximation (RPA) numerical integration accounting for Landau damping:
where the longitudinal RPA dielectric function is
and \(\mathcal{D}(x) = e^{-x^2}\!\int_0^x e^{t^2}\,dt\) is the Dawson function. The lower limit \(\kappa_\mathrm{min} = |\Delta\omega|/(10\,v_\mathrm{th})\) removes the static (\(\kappa \to 0\)) divergence. This is the most accurate of the three models but evaluates one numerical quadrature per frequency point.
In practice the three G-functions give nearly identical line cores; differences appear in the far wings where Landau damping becomes significant.
Profile Accumulation
Each Stark-dressed transition \((p \to k, q)\) is broadened by a Lorentzian of HWHM \(\gamma_e(\Delta E)\):
The three polarization profiles are accumulated by vectorized summation over all microfield quadrature points and all SDTs:
where \(q{=}0 \to P_\pi\), \(q{=}{-1} \to P_{\sigma^+}\) (blue), \(q{=}{+1} \to P_{\sigma^-}\) (red). This direct summation realizes Eq. (10) in practice, with the frequency-dependent half-width \(\gamma_e(\Delta E)\) of Section 3.8 evaluated at the detuning of each \((E, E_{pk})\) pair. The three accumulated channels are the polarization profiles \(I_\pi, I_{\sigma^+}, I_{\sigma^-}\) entering the observed intensity of Eq. (11); for a line of sight at angle \(\alpha\) to \(\vec{B}\) they combine as \(P_\pi\sin^2\!\alpha + \tfrac{1}{2}(P_{\sigma^+}+P_{\sigma^-})(1+\cos^2\!\alpha)\).
Frequency Fluctuation Model (FFM)
The profile assembled in Section 3.9 treats the ionic microfield as frozen during emission. When the ions move appreciably on the emission timescale, that quasi-static average must be replaced by a dynamic one, and the FFM provides this without rebuilding the atomic calculation. Dynamic ion motion causes the static microfield to fluctuate over time; in the FFM [talin] this is modeled as a stationary Markovian process that mixes the Stark-dressed transition components at rate \(\nu_i\).
The Stark-Zeeman Hamiltonian \(H = H_A + V_E\) is still diagonalized at each field configuration \((F,\mu)\), producing dressed-state frequencies \(\omega_k\) and dipole weights \(|d_k|^2\) — the Stark-Dressed Transitions (SDTs). This step is purely static: \(\omega_k\) and \(|d_k|^2\) depend only on the instantaneous ion field. Fast electron collisions add a homogeneous Lorentzian half-width \(\gamma_k\) to each SDT (the GBK width from Section 3.8 ), acting on a timescale short enough that the ion configuration does not change during a single collision. Ion dynamics enter exclusively through \(\nu_i\), estimated as the inverse time for an ion to cross the mean ion spacing at its thermal speed:
The dynamic line profile is then given by the Sherman-Morrison form:
with the static propagator sum:
where \(p_k = |d_k|^2 / r^2\) are the normalized SDT weights and \(r^2 = \sum_k |d_k|^2\). The limit \(\nu_i \to 0\) recovers the static profile; \(\nu_i \to \infty\) gives a single Lorentzian (motional narrowing).
Thermal Doppler Broadening
Thermal motion of the radiating ions causes each photon frequency to be Doppler-shifted by \(\delta\omega = \omega_0\, v_z/c\). Averaging over a Maxwell–Boltzmann velocity distribution yields a Gaussian line shape with \(1/e\) half-width
where \(T_i\) is the ion temperature, \(m_\mathrm{ion}\) the emitter mass, and \(E_0\) the transition energy in eV. For H Balmer-\(\alpha\) at \(T_i = 5\) eV, this gives \(\Delta E_D \approx 0.062\) meV — comparable to the Zeeman splitting at 1 T and negligible relative to the Stark width at \(N_e \sim 10^{23}\) m\(^{-3}\).
The two calculation paths apply Doppler broadening at different stages.
FFM path (calculate_ffm_profile)
When apply_doppler=True (the default), Doppler broadening is applied inside
calculate_ffm_profile on the energy grid, immediately after the Sherman-Morrison
step. The convolution is a zero-padded FFT filter of length \(2N\):
where \(k\) is the discrete-frequency variable (cycles per eV) of the padded rfft
and zero-padding to \(2N\) avoids circular-wrap artifacts. Set
apply_doppler=False to obtain the purely Stark-Zeeman FFM profile.
Static profile path (calculate_static_profile)
calculate_static_profile returns the purely Stark-Zeeman-broadened profile
(quasi-static ion + electron impact only); Doppler is not included. It must be
added afterward as a post-processing step via convolutions.py:
apply_doppler_broadening(wavelengths_nm, profile, Ti_ev, species='H')— convolves with a Gaussian kernel of width \(\Delta\lambda_D\).apply_instrument_broadening(wavelengths_nm, profile, fwhm_nm)— convolves with a Gaussian instrumental slit function.
Both operate on a uniform wavelength grid (convert energy grids to wavelength before calling) and use FFT convolution for efficiency. The Gaussian kernel is:
The complete broadening pipeline for the static path is therefore:
where \(I_\mathrm{SZ}\) is the Stark-Zeeman profile, \(G_D\) is the Doppler Gaussian, and \(G_\mathrm{inst}\) the instrumental Gaussian.