atomic_hamiltonian

class starkzee.atomic_hamiltonian.AtomicState(index, n, l, ml, s=0.5, ms=0.5)[source]

Bases: object

Hydrogenic basis state |n, l, m_l, s=½, m_s⟩.

Parameters:
index: int
n: int
l: int
ml: int
s: float = 0.5
ms: float = 0.5
starkzee.atomic_hamiltonian.build_basis(n)[source]

Return the ordered list of 2n² hydrogenic basis states for shell n.

Each state is an AtomicState labeled by (n, l, m_l, s=½, m_s). The ordering is: outer loop over l ∈ {0, …, n−1}, then m_l ∈ {−l, …, +l}, then m_s ∈ {+½, −½}. This ordering is fixed throughout the package; all Hamiltonian and dipole matrices use the same index convention.

Parameters:

n (int) – Principal quantum number.

Returns:

2n² basis states in the canonical ordering described above.

Return type:

list of AtomicState

starkzee.atomic_hamiltonian.radial_wavefunction(r, n, l, Z)[source]

Return the hydrogenic radial wavefunction R_{nl}(r) [a₀^{−3/2}].

Evaluates the normalized hydrogenic radial wavefunction using the associated Laguerre polynomial representation:

R_{nl}(r) = N_{nl} × exp(−Z r / n) × (2Z r/n)^l × L_{n-l-1}^{2l+1}(2Z r/n)

where r is in units of a₀ and N_{nl} is chosen so that ∫₀^∞ |R_{nl}|² r² dr = 1.

Parameters:

  • r (float or array-like) – Radial coordinate in Bohr radii (a₀).

  • n (int) – Principal quantum number.

  • l (int) – Orbital angular momentum quantum number (0 ≤ l < n).

  • Z (int) – Nuclear charge.

Returns:

Radial wavefunction value [a₀^{−3/2}]. Returns 0 if l ≥ n.

Return type:

float or ndarray

Notes

Uses scipy.special.assoc_laguerre for the polynomial, which handles the high-n case accurately.

starkzee.atomic_hamiltonian.radial_r2_element(n, l1, l2, Z)[source]

Return the radial matrix element ⟨n, l₁ | r² | n, l₂⟩ [a₀²].

Used for off-diagonal elements (abs(l₁ − l₂) = 2) arising from the angular decomposition of r² sin²θ in the quadratic Zeeman term. The diagonal (l₁ = l₂) elements are computed analytically inside build_hamiltonian() using the closed-form formula:

⟨r²⟩_{n,l} = (n²/2Z²) [5n² + 1 − 3l(l+1)]

Off-diagonal elements are computed here by numerical integration:

⟨n, l₁ | r² | n, l₂⟩ = ∫₀^∞ R_{nl₁}(r) r² R_{nl₂}(r) r² dr

Results are cached with functools.lru_cache() to avoid redundant integration during the microfield loop.

Parameters:

  • n (int) – Principal quantum number (same for both bra and ket).

  • l1 (int) – Orbital quantum numbers; the selection rule abs(l₁ − l₂) = 2 must hold.

  • l2 (int) – Orbital quantum numbers; the selection rule abs(l₁ − l₂) = 2 must hold.

  • Z (int) – Nuclear charge.

Returns:

Off-diagonal radial matrix element [a₀²].

Return type:

float

starkzee.atomic_hamiltonian.radial_dipole(n1, l1, n2, l2, Z)[source]

Return the radial transition matrix element ⟨n₁, l₁ | r | n₂, l₂⟩ [a₀].

Computes the radial part of the electric-dipole matrix element by numerical integration:

⟨n₁, l₁ | r | n₂, l₂⟩ = ∫₀^∞ R_{n₁l₁}(r) r R_{n₂l₂}(r) r² dr

The selection rule abs(l₁ − l₂) = 1 is enforced; 0 is returned immediately for forbidden pairs. Results are cached with functools.lru_cache().

Parameters:

  • n1 (int) – Principal quantum number of the first (bra) state.

  • l1 (int) – Orbital quantum number of the bra state.

  • n2 (int) – Principal quantum number of the second (ket) state.

  • l2 (int) – Orbital quantum number of the ket state.

  • Z (int) – Nuclear charge (same for both states; the code is for hydrogen-like ions).

Returns:

Radial dipole matrix element [a₀]. Always ≥ 0 (the sign convention is handled by the angular part angular_dipole_element()).

Return type:

float

Notes

The upper integration limit is set to 150 a₀, which comfortably contains the hydrogenic wavefunction for n ≤ 7 (the largest n used in practice).

starkzee.atomic_hamiltonian.angular_dipole_element(l1, m1, l2, m2, q)[source]

Return the angular part of the dipole matrix element ⟨l₁, m₁ | T_q^(1) | l₂, m₂⟩.

The spherical tensor components T_q^(1) of the position unit vector r̂ are:

  • q = 0 (π): T₀ = cos θ

  • q = +1 (σ+): T₊₁ = sin θ e^{iφ} / √2

  • q = −1 (σ−): T₋₁ = sin θ e^{−iφ} / √2

The matrix elements are expressed via Clebsch-Gordan coefficients and are non-zero only when abs(l₁ − l₂) = 1 and m₁ − m₂ = q. Explicit formulas:

q = 0 (Δm = 0, both states share the same m):

⟨l, m | cos θ | l+1, m⟩ = √((l+1)² − m²) / √((2l+1)(2l+3))

evaluated with l = max(l₁, l₂) − 1.

q = +1 (m₁ = m₂ + 1):

⟨l−1, m | T₊₁ | l, m⟩ = +√((l−m)(l−m−1) / (2(2l−1)(2l+1))) ⟨l, m+1 | T₊₁ | l−1, m⟩ = −√((l+m)(l+m+1) / (2(2l−1)(2l+1)))

q = −1 (m₁ = m₂ − 1):

⟨l−1, m | T₋₁ | l, m⟩ = −√((l+m)(l+m−1) / (2(2l−1)(2l+1))) ⟨l, m−1 | T₋₁ | l−1, m⟩ = +√((l−m)(l−m+1) / (2(2l−1)(2l+1)))

Parameters:

  • l1 (int) – Orbital and magnetic quantum numbers of the bra state.

  • m1 (int) – Orbital and magnetic quantum numbers of the bra state.

  • l2 (int) – Orbital and magnetic quantum numbers of the ket state.

  • m2 (int) – Orbital and magnetic quantum numbers of the ket state.

  • q (int) – Polarization index: 0 (π), +1 (σ+), or −1 (σ−).

Returns:

Angular matrix element (real). Returns 0 if abs(l₁ − l₂) ≠ 1 or m₁ − m₂ ≠ q.

Return type:

float

Notes

The full dipole matrix element between two basis states is −radial_dipole(n1,l1,n2,l2,Z) × angular_dipole_element(l1,m1,l2,m2,q) (the minus sign comes from d = −e r).

starkzee.atomic_hamiltonian.build_hamiltonian(n, Z, B, quadratic_zeeman=True, fine_structure=True, A=1, use_empirical_data=False, atom='H')[source]

Build the (2n²) × (2n²) atomic Hamiltonian matrix in eV.

Constructs the full magnetic Hamiltonian in the uncoupled |n, l, m_l, m_s⟩ basis returned by build_basis(). Four physical contributions are included:

1. Unperturbed energy (diagonal, degenerate across the shell):

E_n = −Z² Ry / n²

2. Spin-orbit coupling L·S (off-diagonal in m_l, m_s):

H_SO = ξ_{nl} L·S, ξ_{nl} = Z⁴ α² Ry / (n³ l (l+½)(l+1))

Written with ladder operators: L·S = L_z S_z + ½(L₊ S₋ + L₋ S₊).

3. Mass-velocity and Darwin corrections (diagonal, only when fine_structure=True):

ΔE_{l=0} = −A (n − ¾) ΔE_{l>0} = −A (n/(l+½) − ¾), A = Z⁴ α² Ry / n⁴

Together with spin-orbit, these reproduce the Dirac fine-structure formula and restore the 2s_{1/2} = 2p_{1/2} degeneracy.

4. Linear Zeeman (diagonal):

H_LZ = μ_B B (m_l + g_s m_s), g_s = 2.0023192

5. Quadratic (diamagnetic) Zeeman (diagonal and off-diagonal in l, only when quadratic_zeeman=True and B > 0):

H_QZ = (e²B²/8m_e) r² sin²θ

Matrix elements in the |n, l, m_l⟩ basis involve ⟨l, m_l | sin²θ | l', m_l⟩ for abs(l − l’) ∈ {0, 2} and the radial elements ⟨n, l | r² | n, l'⟩.

Parameters:

  • n (int) – Principal quantum number.

  • Z (int) – Nuclear charge.

  • B (float) – Magnetic field [T]. B = 0 is fully supported (skips the QZ term).

  • quadratic_zeeman (bool, optional) – Include the diamagnetic quadratic Zeeman term H_QZ (default True). At B ≤ 10 T the QZ shift is typically < 0.1 meV for n ≤ 3 and can be omitted for speed.

  • fine_structure (bool, optional) – Include mass-velocity and Darwin corrections alongside spin-orbit (default True). When False, only H_SO is added, breaking the 2s_{1/2} = 2p_{1/2} degeneracy — useful for isolated unit tests.

  • use_empirical_data (bool, optional) – If True, injects precise empirical field-free energies from atomic_data.

  • atom (str, optional) – Element symbol used to fetch empirical data.

Returns:

Hermitian Hamiltonian matrix in eV.

Return type:

ndarray, shape (2n², 2n²), dtype complex

Notes

The Stark electric-field perturbation is not included here; it is added separately in build_stark_matrix() and combined in solve_starkzee().

At B = 0 the pi/σ decomposition of the resulting profile is physically meaningless (no preferred axis); only the total profile is physically observable.

starkzee.atomic_hamiltonian.diagonalize_hamiltonian(n, Z, B, quadratic_zeeman=True, fine_structure=True, A=1, use_empirical_data=False, atom='H')[source]

Diagonalize the field-free (F = 0) magnetic Hamiltonian for shell n.

Builds the (2n²) × (2n²) Hamiltonian via build_hamiltonian() with no applied electric field and diagonalizes it using numpy.linalg.eigh (real symmetric solver, numerically stable).

Parameters:

  • n (int) – Principal quantum number.

  • Z (int) – Nuclear charge.

  • B (float) – Magnetic field [T]. B = 0 is valid (returns hydrogenic/fine-structure eigenvalues).

  • quadratic_zeeman (bool, optional) – Include the diamagnetic (quadratic Zeeman) term (default True).

  • fine_structure (bool, optional) – Include mass-velocity and Darwin corrections alongside spin-orbit coupling (default True).

  • use_empirical_data (bool, optional) – If True, injects precise empirical field-free energies.

  • atom (str, optional) – Element symbol used to fetch empirical data.

Returns:

  • eigenvalues (ndarray, shape (2n²,)) – Energy eigenvalues in ascending order [eV].

  • eigenvectors (ndarray, shape (2n², 2n²)) – Columns are the corresponding orthonormal eigenstates expressed in the |n, l, m_l, m_s⟩ basis of build_basis().

starkzee.atomic_hamiltonian.dipole_matrix_elements(n_u, n_l, Z, B, quadratic_zeeman=True, fine_structure=True, A=1, use_empirical_data=False, atom='H')[source]

Return transition dipole matrices between the eigenstate bases of n_u and n_l.

Diagonalizes the Hamiltonian for both shells and rotates the uncoupled dipole matrices D_q into the eigenstate basis:

d_q[i, j] = ⟨ψ_l^j | T_q^(1) | ψ_u^i⟩

where ψ_u^i and ψ_l^j are the i-th upper and j-th lower eigenstates.

Parameters:

  • n_u (int) – Upper and lower principal quantum numbers.

  • n_l (int) – Upper and lower principal quantum numbers.

  • Z (int) – Nuclear charge.

  • B (float) – Magnetic field [T].

  • quadratic_zeeman (bool, optional) – Include diamagnetic Zeeman term (default True).

  • fine_structure (bool, optional) – Include MV + Darwin corrections (default True).

  • use_empirical_data (bool, optional) – If True, injects precise empirical field-free energies.

  • atom (str, optional) – Element symbol used to fetch empirical data.

Returns:

  • eigenvalues_u (ndarray, shape (2n_u²,)) – Upper-shell energy eigenvalues [eV].

  • eigenvalues_l (ndarray, shape (2n_l²,)) – Lower-shell energy eigenvalues [eV].

  • d_elements (dict) – Keys are q ∈ {0, 1, −1}. Values are complex arrays of shape (2n_u², 2n_l²) where d_elements[q][i, j] is the dipole matrix element between upper eigenstate i and lower eigenstate j.

starkzee.atomic_hamiltonian.line_strength(n_u, n_l, Z)[source]

Compute the line strength S_ul = Σ_{q,i,j} |⟨l_j|r_q|u_i⟩|² in a₀².

Sums over all polarizations q ∈ {0,±1}, all uncoupled upper states i, and all lower states j, restricted to Δl = ±1 and Δms = 0 (dipole selection rules). Fine structure mixing is neglected (field-free, B = 0 limit with the states labeled by |n, l, ml, ms⟩).

The profile returned by calculate_static_profile integrates to approximately S_ul at B = 0 and Ne → 0. The Gauss-Legendre weights over the microfield angle [0, 1] sum to 1 and the microfield distribution is normalized to 1, so the total quadrature weight is unity.

Returns:

S_ul (float): line strength in a₀²

starkzee.atomic_hamiltonian.oscillator_strength(n_u, n_l, Z)[source]

Compute the weighted absorption oscillator strength gf = g_l × f_lu.

Uses: gf = (2/3) × (ΔE / E_hartree) × S_ul

where E_hartree = 2 × RYDBERG_EV ≈ 27.211 eV, g_l = 2n_l² is the statistical weight of the lower shell, and S_ul is the total line strength (a₀²) summed over all allowed Δl = ±1 pairs.

The result is independent of Z for hydrogen-like ions (the Z-dependence of ΔE ∝ Z² and S_ul ∝ Z⁻² cancel exactly).

Convention matches NIST ASD for full shell-to-shell transitions (e.g. Hα, Hβ). For Lyman series (n_l = 1) only l_u = 1 upper states contribute, so gf is the weighted oscillator strength for the specific 2p/3p/… level.

Returns:

gf (float): dimensionless weighted oscillator strength

starkzee.atomic_hamiltonian.einstein_a(n_u, n_l, Z)[source]

Compute the Einstein A coefficient for spontaneous emission, in s⁻¹.

Uses the atomic-unit formula:

A_ul = (4α³/3) × (ΔE/E_hartree)³ × S_ul / (2n_u²) / τ_au

where α is the fine-structure constant, τ_au = ħ/E_hartree ≈ 2.419 × 10⁻¹⁷ s is the atomic unit of time, and g_u = 2n_u² is the full upper-shell statistical weight (NIST convention: all 2n_u² substates, including any non-radiating ones).

This gives the shell-averaged spontaneous emission rate: the total Ly-α photon emission rate from a uniformly populated n=2 level equals n_{n=2} × einstein_a(2, 1, Z). For NIST’s level-specific A coefficient (e.g. A(2p→1s) using g_k = 6), multiply by 2n_u² / g_participating, where g_participating = Σ_{l_u that have an allowed lower state} 2(2l_u+1).

Returns:

A_ul (float): spontaneous emission rate in s⁻¹