Tanabe–Sugano calculator for 3d transition-metal ions

Four Jupyter notebooks and a Python package (tanabe_sugano/) that build Tanabe–Sugano diagrams for an arbitrary 3d ion in cubic symmetry, with edrixs as the atomic-data and many-body operator backend:

Notebook	What it computes
notebooks/SingleIon_TanabeSugano.ipynb	Pure single-ion ED on the d-shell (V = 0, no bath). Direct diagonalisation of the 252-state d⁵ subspace gives the canonical atomic Tanabe–Sugano with explicit ⟨S²⟩, ⟨L²⟩ labels and the high-spin → low-spin crossover.
notebooks/SIAM_TanabeSugano.ipynb	Full Single-Impurity Anderson Model with a charge-transfer bath. Sₓ-blocked exact diagonalisation of the 15 504-state Fock space yields true eigenvectors and total ⟨S²⟩ per state. Includes separate cells for sweeping 10Dq, V_eg, or Δ.
notebooks/TS_Progression.ipynb	Atomic → +Δ → +V_eg three-panel progression that disentangles the contributions of the charge-transfer energy and the hybridisation. Panels (a) and (b) cross-validate each other (V = 0 ⇒ d-character SIAM states sit at the atomic energies exactly); panel (c) shows the CT-renormalised SIAM.
notebooks/SIAM_with_Experiment.ipynb	Zoomed SIAM ladder with experimental dd peaks drawn as horizontal FWHM bands. Reads your dd-fit CSV, runs a SIAM 10Dq sweep over the dd-excitation window, and overlays both, a direct check of which theory states fall inside each experimental linewidth.

The notebooks live in notebooks/ and produce CSVs in csv_data/ plus PDF/PNG figures in figures/ (relative to wherever you launch Jupyter). Every output gets tagged with the input parameters so different runs don't overwrite each other.

It is aimed at a spectroscopist or solid-state theorist scanning parameter space to connect a measured dd-excitation ladder to a quantitative model: either the bare atomic spectrum (for sanity-checking) or the CT-renormalised SIAM (for materials with significant ligand hybridisation, e.g. α-MnTe at the Mn L-edge).

A parallel set of Cr³⁺ worked examples (single-ion, SIAM, three-panel progression, and CAPS/CCPS experiment overlays) lives in notebooks_Cr3/: a different ion (d³) on different ligands, so the machinery isn't Mn-specific. Run them to regenerate the figures and CSVs.

---

Quickstart

Python API (5 lines)

import numpy as np
from tanabe_sugano import load_params, solve_siam, ts_4panel_plot

params = load_params('params/bo_champion.json')
df = solve_siam(params, dq_grid=np.linspace(0, 4, 50))
fig = ts_4panel_plot(df, params, savepath='figures/my_run')

Parameter sets are stored as JSON files in params/. The load_params function accepts either a file path or a plain dict, so you can also build parameters inline:

from tanabe_sugano import load_params, solve_atomic_d5
params = load_params({'ION': 'Mn', 'NELEC_D': 5, 'KAPPA_DD': 0.7,
                      'U_DD': 7.18, 'ZETA_3D': None})
df = solve_atomic_d5(params, dq_grid=np.linspace(0, 3, 40))

Computed term labels and spin/covalency (adiabatic continuation)

The assignment module assigns each SIAM eigenstate an O_h term symbol, computed by adiabatic continuation in the hybridisation V rather than guessed. At V=0 the impurity decouples, so the low-lying states are the exact atomic d⁵ multiplets; these are labelled against the analytic Griffith branches by their O_h degeneracy and energy, then threaded to the physical V≠0 state by maximum subspace overlap (robust to the arbitrary basis a degenerate multiplet's eigenvectors come back in). This runs on the real-basis Cluster engine, which preserves the O_h degeneracy; the complex-basis solve_siam does not at V≠0 (see Conventions).

from tanabe_sugano import assign_states, assign_peaks_to_states

states = assign_states('params/sample_4382_GaAs.json')   # needs edrixs; ~10 min
# columns: label, irrep, parent, E_champ, alpha_d5, is_d, covalency, proj,
#          min_overlap, mult, ...
print(states[['label', 'E_champ', 'alpha_d5', 'is_d', 'proj', 'min_overlap']])

# map measured dd-peak centres (eV) to the nearest computed multiplet
peaks = [('exciton', 1.805), ('dd1', 2.30), ('dd2', 2.55)]
print(assign_peaks_to_states(states, peaks))

Each state carries two quality numbers: proj (how much atomic character survives at the physical V; a charge-transfer exciton may keep only ~50 % of its parent atomic multiplet) and min_overlap (how cleanly it threaded; >0.96 ⇒ a clean label). Only the three lowest typically thread cleanly enough to claim individually.

Character classifier (d⁵ configuration weight α): is_d is set by alpha_d5 = the summed |amplitude|² over Fock states with impurity occupation n_d = 5 (the d⁵ configuration weight; α ≥ 0.5 ⇒ d-dominant, α < 0.5 ⇒ ligand/CT). This is more robust than a mean-occupation threshold: the mean n_d can exceed 5.5 while d⁵ is still the single largest configuration. Covalency is reported separately as covalency = n_d − 5.

The companion figures come from plotting: a V-swept Tanabe–Sugano diagram (opacity = α) and its V=0 atomic level diagram.

from tanabe_sugano import vsweep_ts_plot, v0_level_diagram
fig = vsweep_ts_plot('params/sample_4382_GaAs.json', label_n_lowest=3)

A new d⁵ system (single-ion and SIAM paths)

Everything is parameter-driven: point any function at a JSON file or dict. To model a new d⁵ system, copy params/template_d5.json and fill in the ion and Slater scaling and either the single-ion (B, C, 10Dq ≡ TEN_DQ) or the SIAM fields (DELTA, VEG, VT2G, TEN_DQ, TEN_DQ_BATH, KAPPA_DD, KAPPA_DP):

import numpy as np
from tanabe_sugano import (load_params, solve_atomic_d5, solve_siam,
                           assign_states, vsweep_ts_plot)

p = load_params('params/template_d5.json')   # edit a copy for your system
solve_atomic_d5(p, dq_grid=np.linspace(0, 3, 40))   # single-ion TS diagram
solve_siam(p, dq_grid=np.linspace(0, 4, 50))        # full SIAM TS (10Dq sweep)
states = assign_states(p)                            # computed term labels + α
vsweep_ts_plot(p)                                    # V-swept TS diagram

params/sample_4382_GaAs.json and params/sample_2236_SrF2.json (α-MnTe on GaAs / SrF₂, advisor-grid fits) are the worked examples; reuse them as a guide.

Comparison overlays (RIXS / ED / Tanabe–Sugano)

tanabe_sugano.comparison builds overlay figures. The ED eigenvalue sticks and Tanabe–Sugano threads are regenerated from the package's own Cluster engine (so they are self-consistent with assign_states and carry the α character), and overlaid on a full SIAM RIXS spectrum and the experimental fits:

from tanabe_sugano.comparison import rixs_ed_ts_compare, edrixs_siam_overlay

# 3 rows: exp+SIAM RIXS / ED sticks (d-dominant vs ligand by α) / transposed TS
rixs_ed_ts_compare('params/sample_4382_GaAs.json',
                   rixs_npz='examples/data/rixs_sim_GaAs_U5.npz',
                   exp_peaks='examples/data/peak_params_GaAs_corr3.csv', fit_tendq=0.40)

# 3 panels: XAS (exp+SIAM) / RIXS / ED ladder with the 3 lowest term labels
edrixs_siam_overlay('params/sample_4382_GaAs.json',
                    rixs_npz='examples/data/rixs_sim_GaAs_U5.npz',
                    xas_txt='examples/data/xas_0T_GaAs.txt',
                    exp_peaks='examples/data/peak_params_GaAs_corr3.csv', fit_tendq=0.40)

The RIXS spectrum itself (full SIAM Kramers–Heisenberg, θ-averaged) is computed separately by tools/rixs_container.py; it needs edrixs, Docker, and MPI (the only non-portable step). Cached champion rixs_sim_*.npz ship in examples/data/, so the figures render without running it. Regenerate with the Docker command in tools/README.md. The ED/TS panels need edrixs to recompute the engine sweep (a few minutes), or pass a precomputed sweep_df=.

Notebooks

# Activate an environment that has edrixs, scipy, pandas, matplotlib
conda activate edrixs_run    # or your equivalent

# From the package root:
jupyter lab notebooks/SingleIon_TanabeSugano.ipynb
# or
jupyter lab notebooks/SIAM_TanabeSugano.ipynb

In each notebook the only cell you need to touch is § 1 Inputs. Run all cells (Run → Run All Cells) and the figures and CSVs are written to disk under tags built from the parameter choice.

Requirements are in requirements.txt; installation notes in docs/installation.md.

---

Conventions

Basis: complex-harmonic Fock space (EG_CPLX / T2G_CPLX)

All single-particle operators (the crystal field (CF), spin-orbit coupling (SOC), hybridisation, and the S²/L² diagnostics) are constructed in the complex-spherical-harmonic Fock basis produced by edrixs get_fock_bin_by_N. Orbital indices within the 10-state d-shell are:

index  0  1  2  3  4  5  6  7  8  9
m_l   -2 -2 -1 -1  0  0 +1 +1 +2 +2
spin   ↑  ↓  ↑  ↓  ↑  ↓  ↑  ↓  ↑  ↓

In Oh cubic symmetry the symmetry-adapted subsets are:

• e_g (EG_CPLX): m = ±2 → indices [0, 1, 8, 9]
• t₂g (T2G_CPLX): m = -1, 0, +1 → indices [2, 3, 4, 5, 6, 7]

These are exported as tanabe_sugano.EG_CPLX and tanabe_sugano.T2G_CPLX.

Critical for forkers: edrixs cf_cubic_d and atom_hsoc both output matrices in this complex-harmonic basis, so no cb_op rotation should be applied before assembling the many-body Hamiltonian. A previous basis-mismatch bug (rotating CF/SOC into the real-cubic basis before passing to build_opers) shifted the ⁴T₁(G) energy by ~197 meV and is now documented in CHANGELOG.md.

Two engines: solve_siam (fast 10Dq sweeps) vs Cluster (term labels)

There are two SIAM solvers, and they use opposite but internally-consistent bases:

• siam.solve_siam: complex-harmonic basis (above). Fast for 10Dq sweeps and the ts_*_plot figures. Caveat: at V≠0 it applies the e_g/t2g hybridisation as a diagonal on complex-harmonic indices (which are diagonal only in the real basis), so it splits the O_h ⁴T₁ triplet 3 → 2+1. Fine for energies/occupations, not for symmetry labels.
• cluster.Cluster (used by assignment and the V-sweep plots): builds every operator, including the Coulomb tensor (transform_utensor to the real basis), in one consistent real-cubic basis. It preserves the O_h 3-fold ⁴T₁ degeneracy and reproduces the reference ed_siam_fort ⁴T₁ to ~0.02 eV. Use it (via assign_states / vsweep_ts_plot) for anything symmetry-resolved.

A regression test (tests/test_assignment.py) guards the degeneracy and the advisor energy on the Cluster path.

Requirements: edrixs (conda)

Importing the package imports edrixs, so run under an environment that has it (e.g. conda activate edrixs_run; Python ≥3.10). assign_states and the V-sweep plots additionally do real ED, so budget a few minutes. Pure CSV/figure post-processing has no edrixs dependency. For clean term labelling, SOC is turned off (soc=False, the assignment default): ζ_3d ≈ 40 meV shifts dd energies <1 meV but lifts the exact O_h degeneracy the threading relies on.

---

Tests

cd TanabeSugano
conda activate edrixs_run
pytest tests/ -v

The suite (13 tests) covers:

• Atomic d⁵ spectrum vs exact-diagonalisation reference values
• Linear κ-scaling of term energies
• V = 0 SIAM fast-path equals atomic d⁵ exactly
• Per-state observable checks (S², L², multiplicity)
• Regression tests for the basis-mismatch fix (pre-fix: 2243 meV, post-fix: 2440 meV for ⁴T₁(G) at ten_dq = 0.5 eV, κ = 0.7)

The slow SIAM tests (~75 s) run build_static_caches and ARPACK; the fast tests (atomic-limit only) take ~10 s.

---

Adding a new candidate

1. Create params/my_candidate.json following the schema in any existing file. The _metadata block is optional but recommended.
2. Load and validate: load_params('params/my_candidate.json') raises ValueError if required keys are missing.
3. Run the standard sweep:

   params = load_params('params/my_candidate.json')
   df = solve_siam(params, dq_grid=np.linspace(0, 4, 50))
   fig = ts_4panel_plot(df, params)

Required keys for single-ion (no bath): ION, NELEC_D, KAPPA_DD, U_DD, ZETA_3D. Add NBATH_ORB, DELTA, VEG for SIAM mode. See docs/parameters.md for units and ranges.

---

Example outputs

The figures below come directly from the package API. They use one worked example, a d⁵ impurity with a charge-transfer bath (parameters in params/sample_4382_GaAs.json), but nothing in them is specific to that system: point any function at a different parameter file (see Adding a new candidate) to get the same plots for another d-electron count or ligand environment. The single-ion ladder, the three-panel progression, and the RIXS / ED / Tanabe–Sugano comparison come from the notebooks and the other plotting / comparison functions documented above.

Canonical Tanabe–Sugano (10Dq-swept), by total spin

byspin_tendq_plot(...): the classic Tanabe–Sugano control parameter, the crystal-field splitting 10Dq, on the x-axis, with the dd-multiplet ladder split into one panel per total spin. Same overlap-threaded, α-coloured rendering as the V-sweep below (sextet S = 5/2, quartet S = 3/2, doublet S = 1/2, the only families a d⁵ shell allows); the dashed line marks the parameter file's 10Dq and the hybridisation is held fixed at its champion value.

Spin-resolved V-swept Tanabe–Sugano

byspin_vsweep_plot(...): the dd-multiplet ladder vs the e_g hybridisation V_eg, split into one panel per total spin. For a d⁵ impurity only the half-integer families exist (sextet S = 5/2, quartet S = 3/2, doublet S = 1/2; there are no singlet/triplet panels). Each thread is followed continuously through level crossings by eigenvector overlap (not energy-rank sorting, which would kink at every crossing), and coloured by its d⁵ configuration weight α (opaque = d-dominant, faded = ligand / d⁶L̄). V_eg = 0 is the atomic limit; turning hybridisation on pulls the quartet ladder down and mixes in ligand character. The optional fourth panel overlays a SIAM RIXS spectrum on the shared energy-loss axis, lining the computed multiplets up against the measured dd features.

Charge-transfer (Δ-swept) Tanabe–Sugano

byspin_delta_plot(...): the same spin-resolved rendering swept over the charge-transfer energy Δ, at fixed 10Dq and hybridisation. Small Δ is the strong-CT (covalent) limit where ligand / d⁶L̄ character (faded) mixes deep into the multiplets; large Δ approaches the ionic d⁵ atomic spectrum (opaque). The dashed line marks the parameter file's Δ.

---

Documentation

Document	Read this when…
docs/methodology.md	…you want the physics and math behind the code: the Hamiltonian, the S²-from-eigenvectors trick, the Sz-block decomposition, the threading scheme.
docs/usage_guide.md	…you're running the notebook for the first time, debugging a parameter combination, or porting to a different ion. Walk-through of every cell.
docs/parameters.md	…you need a glossary of all input variables, units, and physically reasonable ranges.
docs/outputs.md	…you want to know what every column of the CSVs means and how the figures are laid out.
docs/installation.md	…Python environment, edrixs install, MPI requirements.
docs/atomic_limit_caveat.md	…what "atomic limit" means in the SIAM context; how to identify d-character vs ligand-character eigenstates; why V_eg = 0 is not the same as the single-ion notebook.
docs/SIAM_model.tex / .pdf	…LaTeX reference document for talks: full Hamiltonian, parameter glossary table with α-MnTe values, ⟨S²⟩-from-eigenvectors derivation, numerical scheme.

---

Citation

These notebooks build on the edrixs library (Y. Wang, G. Fabbris, M. P. M. Dean, Comput. Phys. Commun. 243, 151 (2019), doi:10.1016/j.cpc.2019.04.018). If you use the SIAM machinery please also cite the fedrixs Fortran backend even though here we use the Python-side direct ED.

The methodology (direct ED with explicit ⟨S²⟩ from $|\hat S_-|\psi\rangle|^2 + S_z(S_z-1)$, an Sz-blocked Fock basis, and per-state DM diagnostics) is standard material in any many-body textbook; the implementation here is original to this package.

If these notebooks help your paper, please cite this repository and the edrixs paper above. A CITATION.cff is included.

---

Licence

MIT. See LICENSE. edrixs is BSD-3.