sc-SQ: Self-Consistent Spectral Quadrature

Python implementation of the sc-SQ method for computing many-body Green functions, as described in

S. Yu. Kruchinin, Self-Consistent Spectral Quadrature Approach to Many-Body Green Functions, submitted to PRB.

The method maps spectral moments onto a Gauss–Christoffel quadrature via the Lanczos recurrence, obtains a pole decomposition $G_{[N]}(z) = \sum_i w_i / (z - \varepsilon_i)$, and drives it to self-consistency through an equation-of-motion loop.

---

Repository layout

sc-sq/
├── sc_sq/
│   ├── core.py          # Hankel matrix, SVD rank, Lanczos, poles/weights
│   ├── anderson.py      # Anderson single-impurity model (AIM) solver
│   ├── ed.py            # Exact-diagonalization (Caffarel-Krauth) impurity solver
│   ├── dmft.py          # DMFT + sc-SQ / sc-GW for Hubbard / Bethe lattice
│   ├── spectral.py      # Spectral reconstruction: LB, BRI, SEC, terminator
│   └── bethe.py         # Bethe-lattice benchmark helpers (branch selection, broadening)
├── benchmarks/
│   ├── 01_anderson/
│   │   ├── oneshot.py               # Fig. 1: one-shot sc-SQ, semielliptic bath
│   │   └── mixed_valence.py         # Fig. 2: mixed-valence self-consistency
│   └── 02_hubbard_bethe/
│       ├── hubbard_bethe.py         # Fig. 3: Mott transition, spectral + Z
│       └── convergence_N.py         # Fig. 4: sc-SQ convergence with rank N
├── tests/
│   ├── test_anderson.py
│   ├── test_core.py
│   └── test_dmft.py
└── README.md

---

Module reference

sc_sq/core.py

The algorithmic backbone of sc-SQ. Given an array of spectral moments $\mu_0,\ldots,\mu_{2N-1}$, it assembles the $N\times N$ symmetric Hankel matrix $M_{ij}=\mu_{i+j}$, selects the physical rank $N^*$ via the SVD criterion $\sigma_n/\sigma_1>\tau$, runs the modified Chebyshev (Gautschi) sigma-recurrence to produce Lanczos coefficients ${a_k,b_k}$ for the Jacobi tridiagonal matrix, diagonalises it to obtain Gauss–Christoffel poles $\varepsilon_i$ and weights $w_i$, and evaluates $G_{[N]}(z)$ by either partial-fraction summation or backward continued-fraction recursion. The top-level entry point spectral_quadrature(moments) runs the entire pipeline and returns an SQResult container; lower-level building blocks are also exported for custom workflows.

Public API: hankel_matrix, svd_rank, lanczos_from_moments, poles_and_weights, green_function, green_function_cf, spectral_quadrature, SQResult.

---

sc_sq/ed.py

Exact-diagonalization (Caffarel–Krauth 1994) impurity solver for the single-orbital Anderson model embedded in a Bethe-lattice DMFT loop. ED_BetheLattice.solve() fits n_bath discrete bath levels ${\varepsilon_k, V_k}$ to the current Weiss field on the Matsubara axis, builds the $(N_{\rm bath}+1)$-site Fock-space Hamiltonian in $(N_\uparrow, N_\downarrow)$ sectors, diagonalizes each block with scipy.linalg.eigh, and evaluates the retarded Green function via the Lehmann representation. The returned EDResult exposes A_omega, lehmann_poles, lehmann_weights, and a moments(2N) method that computes exact spectral moments up to order $2N$ — used by convergence_N.py to seed the sc-SQ loop and break EOM-closure traps at higher $N$.

Public API: AndersonED, EDResult, ED_BetheLattice.

---

sc_sq/anderson.py

sc-SQ impurity solver for the Anderson single-impurity model with a semielliptic bath. anderson_moments() computes moments $\mu_0\ldots\mu_7$ analytically from the equations of motion given $\varepsilon_d$, $U$, $\langle n_\sigma\rangle$, $D_{\rm docc}$, and the hybridization parameters returned by wide_band_hybridization(); the $\mu_3$ formula is derived so that all odd moments vanish exactly at particle-hole symmetry for any $\Delta$. AndersonSolver.solve() drives the moment–quadrature–occupancy loop with linear mixing until convergence, computing occupancy as the Lorentzian-broadened integral $\langle n\rangle=\sum_i w_i[\tfrac{1}{2}+\arctan(-\varepsilon_i/\Delta)/\pi]$ and switching between the Hubbard-I closure ($N\le2$) and the Ward-identity closure ($N\ge3$) for the double occupancy. AndersonSolver.one_shot() evaluates the spectral hierarchy at a fixed input state without iteration. Results are returned as AndersonResult dataclasses that expose poles, weights, occupancy, double_occ, converged, and convenience methods green(z), spectral_lb(omega), spectral_bri(omega).

Public API: AndersonParams, AndersonSolver, AndersonResult, wide_band_hybridization, anderson_moments.

---

sc_sq/spectral.py

Four real-frequency spectral reconstruction procedures, all taking a poles/weights pair from core.spectral_quadrature as input. lorentzian_broadening(omega, poles, weights, eta) convolves each delta-function pole with a Lorentzian of half-width $\eta$ (Eq. 25). barycentric_interpolation(omega, poles, weights) constructs a parameter-free smooth envelope via barycentric rational interpolation with alternating-sign weights $\lambda_i=(-1)^i\sqrt{w_i}$ (Eqs. 26–29). self_energy_continuation(omega, poles, weights, G0_inv, eta) extracts $\Sigma=G_0^{-1}-G_{[N]}^{-1}$ and returns the physical Green function $G_{\rm phys}$ with intrinsic bath-driven linewidths (Eqs. 22–23). terminator_green(z, a, b) evaluates the continued-fraction Green function with a Beer–Pettifor semi-elliptic tail (Eq. 30). The utility spectral_from_G(G) converts any retarded Green function to the spectral function $A=-\mathrm{Im},G/\pi$.

Public API: lorentzian_broadening, barycentric_interpolation, self_energy_continuation, terminator_green.

---

sc_sq/dmft.py

DMFT + sc-SQ two-loop solver for the half-filled Hubbard model on the Bethe lattice ($G_0^{-1}(z)=z+\mu_{\rm chem}-t^2G_{\rm loc}(z)$). hubbard_moments_bethe() computes local moments $\mu_0\ldots\mu_3$ from $U$, $t$, and the kinetic amplitude $t^2\langle c^\dagger c\rangle$, updated self-consistently at each DMFT step. BetheLattice.solve() runs an outer DMFT iteration (updating the Weiss field) around an inner sc-SQ loop (converging the impurity moments) and returns a DMFTResult with converged poles, weights, occupancy, double occupancy, kinetic energy, and full iteration history. BetheLattice.quasiparticle_weight() computes $Z=(1-\partial\mathrm{Re},\Sigma/\partial\omega\rvert_0)^{-1}$ by finite-difference differentiation of the self-energy; BetheLattice.spectral_sec() evaluates the physical spectral function via self-energy continuation at the DMFT fixed point. The non-interacting Bethe DOS $\rho_0(\omega)=\frac{2}{\pi W^2}\sqrt{W^2-\omega^2}$ is exposed as bethe_dos().

Public API: BetheLattice, DMFTResult, bethe_dos.

---

benchmarks/01_anderson/oneshot.py

Reproduces Fig. 1 of the paper. For the particle-hole symmetric AIM with a semielliptic conduction band ($\Gamma/D=0.1$, $U/D=1.0$), reproduces the configuration of Žitko & Pruschke, PRB 79, 085106 (2009), Fig. 8. Panel (a) shows the one-shot spectral function at $N=3$ reconstructed via self-energy continuation (SEC) using the analytic semielliptic $\Delta(\omega+i0^+)$, compared with digitized NRG data; optional reference curves from nrg_spectral_semielliptic_U_over_D_{value}.dat are plotted in black when present. Panel (b) shows the Hankel SVD spectrum $\sigma_n/\sigma_1$ with $N^*=3$ selected at the $\tau=10^{-8}$ threshold. Writes oneshot.pdf and oneshot.h5; accepts --show.

---

benchmarks/01_anderson/mixed_valence.py

Reproduces Fig. 2 of the paper. Scans $\varepsilon_d/U$ from $-0.5$ (particle-hole symmetry) to $0$ (empty-orbital limit); panel (a) compares the one-shot occupancy with the sc-SQ self-consistent fixed point at $N=3$. Panel (b) shows the spectral function at $\varepsilon_d/U=-0.2$ via SEC, with three curves: one-shot, sc-SQ, and the Luttinger-corrected sc-SQ (shift $\mu_L$ enforcing $\mathrm{Re},\Sigma_{\rm corr}(0)=0$). Accepts --N and --show. Writes mixed_valence.pdf and mixed_valence.h5.

---

benchmarks/02_hubbard_bethe/hubbard_bethe.py

Reproduces Fig. 3 of the paper. Runs BetheLattice.solve() for $U/D\in{1.0, 2.94}$ (correlated metal and Mott critical point) and GW_BetheLattice for the sc-GW reference on a denser $U/D$ grid. Panel (a) shows Lorentzian-broadened spectral functions $A(\omega)$, with optional NRG overlays from nrg_spectral_U_over_D_{X.X}.dat. Panel (b) plots the quasiparticle weight $Z$ vs $U/D$ for both sc-SQ and sc-GW, with an optional DMFT+NRG reference from nrg_Z_vs_UoverD.dat. Results are cached in hubbard_bethe.h5; use --recalc sc-sq, --recalc sc-gw, or --recalc all to force recomputation. Accepts --N and --show. Writes hubbard_bethe.pdf and hubbard_bethe.h5.

---

benchmarks/02_hubbard_bethe/convergence_N.py

Reproduces Fig. 4 of the paper. Shows sc-SQ convergence with pole rank $N\in{3,5,7}$ for the half-filled Hubbard model on the Bethe lattice, seeding each run with exact Lehmann moments from a 3-bath ED solution (ED_BetheLattice) to break EOM-closure traps at higher $N$. Panel (a): $U/D=1.0$ (correlated metal). Panel (b): $U/D=1.6$ (strongly correlated, near the Mott transition). Broadening uses $\eta_{\rm qp}=0.10,D$ for quasiparticle poles and $\eta_{\rm sat}=0.22,D$ for Hubbard satellites; optional NRG reference from nrg_spectral_U_over_D_{value}.dat. Writes convergence_N.pdf.

---

Installation

Prerequisites

Package	Version	Notes
Python	≥ 3.11
NumPy	≥ 1.24
SciPy	≥ 1.10	optional but recommended
Matplotlib	≥ 3.7	required for the examples
h5py	≥ 3.8	optional, for HDF5 output
TRIQS	≥ 3.2	optional, for GfReFreq I/O

Install from source (editable)

git clone https://github.com/polariton/sc-sq.git
cd sc-sq
pip install -e .

For example/benchmark scripts, install optional plotting and HDF5 extras:

pip install -e ".[examples]"

Or, without a setup.py, add the repo root to your PYTHONPATH:

export PYTHONPATH="$PWD:$PYTHONPATH"

WSL / Ubuntu quick-start

# Inside Ubuntu (WSL)
cd ~/src/sc-sq
python -m pip install numpy scipy matplotlib h5py
python -c "from sc_sq.core import spectral_quadrature; print('OK')"

---

Running the tests

cd sc-sq
python -m pytest tests/ -v

All 51 tests should pass.

---

Reproducing the paper figures

Each script writes a PDF in the current working directory and, if h5py is installed, an HDF5 archive with all numerical data. Pass --show to open an interactive Matplotlib window.

Figure 1 — One-shot sc-SQ, semielliptic bath

Panels: (a) Spectral function at $U/D=1.0$ via SEC, compared with NRG (Žitko & Pruschke 2009) (b) Hankel SVD spectrum $\sigma_n/\sigma_1$ showing $N^*=3$

cd benchmarks/01_anderson
python oneshot.py

Output: oneshot.pdf, oneshot.h5

---

Figure 2 — Mixed-valence self-consistency

Panels: (a) $\langle n_\sigma \rangle$ vs $\varepsilon_d / U$: one-shot vs sc-SQ (b) Spectral function at $\varepsilon_d / U = -0.2$: one-shot, sc-SQ, and Luttinger-corrected

cd benchmarks/01_anderson
python mixed_valence.py

Accepts --N to change the pole rank.

Output: mixed_valence.pdf, mixed_valence.h5

---

Figure 3 — Mott metal-insulator transition (Bethe lattice)

DMFT + sc-SQ and sc-GW for the half-filled Hubbard model.

Panels: (a) Local spectral function $A(\omega)$ at $U/D = 1.0$ and $2.94$ (b) Quasiparticle weight $Z$ vs $U/D$ for sc-SQ and sc-GW

cd benchmarks/02_hubbard_bethe
python hubbard_bethe.py

Output: hubbard_bethe.pdf, hubbard_bethe.h5

---

Figure 4 — sc-SQ convergence with rank $N$

Panels: (a) $A(\omega)$ for $N = 3, 5, 7$ at $U/D = 1.0$ (b) $A(\omega)$ for $N = 3, 5, 7$ at $U/D = 1.6$

cd benchmarks/02_hubbard_bethe
python convergence_N.py

Output: convergence_N.pdf

---

Minimal library usage

import numpy as np
from sc_sq.core import spectral_quadrature
from sc_sq.anderson import AndersonSolver
from sc_sq.spectral import lorentzian_broadening

# Anderson impurity model: U/Delta = 4, particle-hole symmetric
solver = AndersonSolver(eps_d=-2.0, U=4.0, Delta=1.0, D=10.0)

# Self-consistent solution at N=3
result = solver.solve(N=3)
print(result)

# Lorentzian-broadened spectral function
omega = np.linspace(-12, 12, 2000)
A = lorentzian_broadening(omega, result.poles, result.weights, eta=0.02)

# One-shot (no self-consistency) hierarchy
sq = solver.one_shot(N=3, n_sigma=0.5, D_docc=0.25)
print(sq)

---

Key formulas

Symbol	Meaning
$\mu_n = \langle d \vert H^n \vert d \rangle$	spectral moments
$M_{ij} = \mu_{i+j}$	Hankel moment matrix
$N^* = \max{n : \sigma_n/\sigma_1 > \tau}$	physical rank (SVD criterion)
$G_{[N]}(z) = \sum_i w_i / (z - \varepsilon_i)$	sc-SQ Green function
$A_{LB}(\omega) = \pi^{-1} \sum_i w_i \eta / [(\omega-\varepsilon_i)^2+\eta^2]$	Lorentzian broadening

---

License

MIT