Liouville-Lanczos Linear-Response TDDFT

Perturbation Methods

A linear-response method for excited-state spectra — optical absorption, electron energy loss, plasmons — built on top of density functional theory (DFT). The idea is due to Walker, Saitta, Gebauer, and Baroni (2006), with the electron-energy-loss extension by Timrov, Vast, Gebauer, and Baroni; the implementation lives in Quantum ESPRESSO's turboTDDFT and turboEELS modules. It is one of the cleanest examples on this site of how a Krylov-subspace iterative method, applied to a non-Hermitian super-operator coming out of physics, gives you the entire frequency-resolved spectrum from a SINGLE matrix-free chain.

Extracting the excitation spectrum of a realistic molecule or solid is, in principle, a non-Hermitian eigenvalue problem on the Liouvillian — an operator whose dimension scales as (occupied × virtual) Kohn-Sham orbitals. For anything past a handful of atoms, you cannot form it, let alone diagonalize it. The way out is to never build the operator. You apply it matrix-free to one chosen perturbation vector, build a Krylov subspace by repeated application, and read the spectrum off the tridiagonal coefficients as a continued fraction. One chain gives you the full $ω$ -resolved response, not just one root at a time.

What is the actual problem

You have a system of interacting electrons in its DFT ground state with density $ρ_{0} (r)$ and Kohn-Sham Hamiltonian $\hat{H}_{0}$ . You shine light on it. To leading order in the field, the density responds linearly:

δ ρ (r, ω) = \int d r^{'} χ (r, r^{'}; ω) V_{ext} (r^{'}, ω) .

The density-density response function $χ$ contains EVERYTHING you can probe with a weak external field: optical absorption (when $V_{ext}$ is a uniform dipole), inelastic electron scattering / EELS (when $V_{ext}$ is the bare Coulomb of a passing electron), plasmons (poles of $χ$ at finite momentum transfer). The whole game is computing $χ (ω)$ .

Why the standard Casida approach fails for solids

The conventional TDDFT route is Casida's equation: cast linear response as a generalized eigenvalue problem

(A - B^{*} B - A^{*}) (X Y) = Ω (X Y),

diagonalize it, get all excitation energies $Ω_{n}$ and oscillator strengths, then build $χ (ω) = \sum_{n} f_{n} / (ω^{2} - Ω_{n}^{2} + i η)$ . For a small molecule this works beautifully. The matrix dimension is $N_{occ} \times N_{virt}$ — manageable.

For a SOLID this is hopeless. A converged plane-wave calculation has thousands of bands per k-point and a dense k-mesh; $N_{occ} \cdot N_{virt}$ easily exceeds $1 0^{6}$ . Diagonalization of a million-dimensional non-Hermitian matrix is not happening. Worse, the physically interesting spectrum is a CONTINUUM — there is no sense in resolving every individual eigenpair when what you actually want is the smooth function $Im χ (ω)$ over a frequency window.

The Liouville-Lanczos method sidesteps the eigenvalue problem entirely.

The Liouvillian superoperator

Reformulate linear response in terms of the density matrix instead of orbitals. The Kohn-Sham density matrix $\overset{ρ}{^}_{0}$ is a projector onto the occupied subspace. A small external perturbation $V^{'}$ drives a small change $\overset{ρ}{^}^{'} (ω)$ . Linearizing the time-dependent Kohn-Sham equation gives:

(ω - \hat{L}) \overset{ρ}{^}^{'} (ω) = [V^{'}, \overset{ρ}{^}_{0}] .

Here $\hat{L}$ is the LIOUVILLIAN — a SUPER-operator (an operator that acts on operators) whose action on a density-matrix variation $\hat{X}$ is

\hat{L} \hat{X} = [\hat{H}_{0}, \hat{X}] + [V_{HXC}^{'} [\hat{X}], \overset{ρ}{^}_{0}] .

The first term is the bare commutator — kinematics of the non-interacting Kohn-Sham system. The second term is the SELF-CONSISTENT linear-response feedback: the density variation generated by $\hat{X}$ induces a perturbation of the Hartree and exchange-correlation potential, which acts back on the occupied subspace. This is the TDDFT kernel $K_{HXC}$ wired into the dynamics. Without it you'd just have non-interacting Kohn-Sham excitations — an independent-particle spectrum that systematically misses excitonic and collective effects.

Crucial structural fact: the Liouvillian is NOT Hermitian under the usual operator inner product. The bare commutator $[\hat{H}_{0}, \cdot]$ is anti-Hermitian; the kernel term breaks the simple structure further. This is why the eventual Lanczos chain has to be a non-Hermitian (bi-orthogonal) variant rather than the symmetric three-term recursion you'd use for, say, ground-state DFT.

The dipole polarizability — the quantity that gives optical absorption — is then

α (ω) = - ⟨ ⟨ \hat{d} ∣ (\hat{L} - ω)^{- 1} ∣ [\hat{d}, \overset{ρ}{^}_{0}]⟩ ⟩,

where $\hat{d}$ is the dipole operator and the inner product is in the space of single-particle operators. The whole spectrum is the matrix element of a resolvent. You do NOT need to know the eigenvalues of $\hat{L}$ individually — you only need ONE diagonal element of $(\hat{L} - ω)^{- 1}$ , between two specific vectors. That is what Lanczos is built for.

The Lanczos chain

Standard Lanczos (the symmetric one, see Lanczos Iteration) takes a Hermitian operator $A$ and a starting vector $v_{0}$ and builds an orthonormal basis of the Krylov subspace $span {v_{0}, A v_{0}, A^{2} v_{0}, \dots}$ in which $A$ is tridiagonal. Three-term recursion, cheap, memory-friendly.

For a non-Hermitian operator like the Liouvillian, you need the BI-ORTHOGONAL Lanczos: build two sequences ${q_{j}}$ (right Krylov) and ${p_{j}}$ (left Krylov, generated by $\hat{L}^{†}$ ) such that $⟨ p_{i} ∣ q_{j} ⟩ = δ_{ij}$ . The recursion is:

\hat{L} q_{j} = α_{j} q_{j} + β_{j + 1} q_{j + 1} + γ_{j} q_{j - 1},

\hat{L}^{†} p_{j} = α_{j} p_{j} + γ_{j + 1} p_{j + 1} + β_{j} p_{j - 1},

with starting vectors chosen so that $q_{0}$ is the perturbation $[\hat{d}, \overset{ρ}{^}_{0}]$ and $p_{0}$ is the dipole-readout bra. The three families of coefficients $α_{j}$ , $β_{j}$ , $γ_{j}$ — produced one per Lanczos iteration — are the OUTPUT of the calculation. After $N$ iterations you have a non-symmetric tridiagonal matrix $T_{N}$ with these coefficients on its three diagonals, and it represents the Liouvillian projected onto the bi-orthogonal Krylov basis.

The original Walker-Saitta-Gebauer-Baroni formulation (PRL 2006) exploits a specific structure of the closed-shell Liouvillian that lets the bi-orthogonal recursion be written in a particularly symmetric way — the "pseudo-Hermitian" or "batch" Lanczos — but the spirit is the same: tridiagonalise on the fly, never store more than three response vectors at once.

The continued-fraction polarizability

Here is where the method earns its keep. Given the tridiagonal projection $T_{N}$ , the matrix element

α (ω) = - ⟨ p_{0} ∣ (\hat{L} - ω)^{- 1} ∣ q_{0} ⟩ \approx - [(T_{N} - ω 1)^{- 1}]_{00}

is the (0,0) entry of the inverse of a tridiagonal matrix. There is a closed-form expression for that — a CONTINUED FRACTION:

α (ω) = \frac{- 1}{( α _{0} - ω ) - \frac{β _{1} γ _{1}}{( α _{1} - ω ) - \frac{β _{2} γ _{2}}{( α _{2} - ω ) - \dots}}} .

The Lanczos coefficients ARE the continued-fraction coefficients. Once you've run the chain for N iterations, evaluating $α (ω)$ at any frequency is an O(N) bottom-up evaluation of the continued fraction. The same N coefficients give you the polarizability at every frequency — you sweep the spectrum for free.

This is the punchline. Casida diagonalizes $N^{2}$ floating-point numbers' worth of matrix and gives you the spectrum as a sum of delta functions; Liouville-Lanczos generates N coefficients and gives you the spectrum as a continuous function of $ω$ at the cost of plotting a continued fraction.

The chain in code

A schematic non-Hermitian Lanczos plus continued-fraction evaluation. Real implementations apply $\hat{L}$ matrix-free — one application is essentially one Sternheimer/SCF-like step on the perturbed orbitals — but the bookkeeping below is what those codes are doing.

# Schematic non-Hermitian Lanczos on the TDDFT Liouvillian L,
# followed by continued-fraction evaluation of the dipole polarizability
#   alpha(omega) = -<d | (L - omega)^(-1) | v>,
# where |v> is the dipole-perturbation vector and <d| is the dipole-readout
# bra. In real codes L is never built; it is APPLIED on the fly to a
# variational density-matrix object stored as a batch of orbital
# response vectors. The math here is identical.

import numpy as np

def apply_L(v):
    """
    Apply the Liouvillian superoperator L = [H0, .] + [V_HXC'[rho_v], rho0]
    to the response object v. In QE this is one SCF-like step:
      (1) compute the response density delta_rho from v,
      (2) build the perturbing potential delta_V_HXC = K * delta_rho
          where K is the Hartree-XC kernel,
      (3) act with H0 - epsilon on v and add the kernel contribution.
    We treat it as a black box here.
    """
    ...

def non_hermitian_lanczos(v0, w0, n_iter):
    """
    Bi-orthogonal Lanczos: build two Krylov sequences
       q_{j+1} = (L q_j - alpha_j q_j - beta_j q_{j-1}) / gamma_{j+1}
       p_{j+1} = (L^T p_j - alpha_j p_j - gamma_j p_{j-1}) / beta_{j+1}
    such that <p_i | q_j> = delta_ij. The Liouvillian is NOT Hermitian
    (L^dag != L for closed-shell TDDFT under the standard inner product),
    so the symmetric three-term Lanczos has to be generalized.
    """
    # Normalize the starting pair so <w0 | v0> = 1.
    nrm = np.vdot(w0, v0)
    q = v0 / np.sqrt(nrm)
    p = w0 / np.sqrt(nrm)

    q_prev = np.zeros_like(q)
    p_prev = np.zeros_like(p)
    beta_prev = gamma_prev = 0.0

    alphas, betas, gammas = [], [], []

    for j in range(n_iter):
        Lq  = apply_L(q)
        LTp = apply_L(p)        # adjoint application; in QE done by sign flip

        alpha = np.vdot(p, Lq)
        alphas.append(alpha)

        r = Lq  - alpha * q  - gamma_prev * q_prev
        s = LTp - alpha * p  - beta_prev  * p_prev

        rs = np.vdot(s, r)
        if abs(rs) < 1e-14:
            break                # serious breakdown; in practice almost never
        beta  = np.sqrt(abs(rs))
        gamma = rs / beta        # gamma * beta = <s|r>

        q_prev, p_prev = q, p
        q = r / beta
        p = s / np.conj(gamma)

        betas.append(beta)
        gammas.append(gamma)
        beta_prev, gamma_prev = beta, gamma

    return np.array(alphas), np.array(betas), np.array(gammas)

def polarizability_from_chain(alphas, betas, gammas, omega, eta=0.01):
    """
    The Krylov projection of (L - omega)^(-1) is a TRIDIAGONAL matrix with
    diagonal entries (alpha_j - omega) and off-diagonals (beta_j, gamma_j).
    Its (1,1) entry — what alpha(omega) needs — is a continued fraction:

        alpha(omega) = z0 / [ (a_0 - omega - i eta) -
                              b_1 g_1 / [ (a_1 - omega - i eta) -
                                          b_2 g_2 / [ ... ] ] ]
    """
    z = omega + 1j * eta
    cf = 0.0
    # Evaluate from the bottom up.
    for j in range(len(alphas) - 1, 0, -1):
        cf = betas[j-1] * gammas[j-1] / (alphas[j] - z - cf)
    cf = 1.0 / (alphas[0] - z - cf)
    return cf       # times the (already-absorbed) <d|v0> normalization

# Sweep frequency once. ONE Lanczos chain gives the spectrum at every omega.
alphas, betas, gammas = non_hermitian_lanczos(v0, w0, n_iter=2000)
omegas = np.linspace(0.0, 1.5, 2000)        # in Hartree
spectrum = np.array([-polarizability_from_chain(alphas, betas, gammas, w).imag
                     for w in omegas])

A runnable validation against Casida

The schematic above hides the bookkeeping that has to be right for the continued fraction to actually reproduce the spectrum. Here is a small but complete test: a 6 particle-hole TDDFT-like model (so the Liouvillian is a 12×12 matrix that we can build explicitly), a random symmetric Hartree-XC kernel, and a comparison between (a) full diagonalization of $L$ followed by spectral summation — Casida — and (b) bi-orthogonal Lanczos plus continued-fraction evaluation. When the chain runs to the full dimension, the two should agree to machine precision; when truncated, the Lanczos answer should converge to the Casida result as the chain grows. This is the cleanest sanity check for the whole construction.

# Self-contained test: tiny TDDFT-like model, bi-orthogonal Lanczos on the
# Liouvillian, validated against direct Casida diagonalization. The model is
# small enough that we can build L explicitly (12x12) and diagonalize it — so
# we can check that the Lanczos continued fraction agrees with the spectral
# sum to machine precision.
#
# Closed-shell, real-orbital, singlet linear response:
#   A = diag(eps_a - eps_i) + K,   B = K,   L = [[A, B], [-B, -A]].
# L is non-Hermitian; its eigenvalues come in (+Omega, -Omega) pairs.

import numpy as np
from scipy.linalg import eig

np.random.seed(0)

# ─── Toy system: 2 occupied x 3 virtual ──────────────────────────────────
n_occ, n_vir = 2, 3
n_ph = n_occ * n_vir            # 6 particle-hole pairs
dim  = 2 * n_ph                  # 12 = Liouvillian dimension

eps_occ = np.array([-1.20, -0.60])
eps_vir = np.array([ 0.40,  0.90,  1.60])

delta = np.array([eps_vir[a] - eps_occ[i]
                  for i in range(n_occ) for a in range(n_vir)])

K = 0.15 * np.random.randn(n_ph, n_ph)
K = 0.5 * (K + K.T)                 # symmetric Hartree-XC kernel
A = np.diag(delta) + K
B = K.copy()
L = np.block([[ A,  B ],
              [-B, -A ]])

# ─── Casida reference: diagonalize L directly ────────────────────────────
eigvals, R = eig(L)
order   = np.argsort(eigvals.real)
eigvals = eigvals[order]
R       = R[:, order]

_, Lvecs = eig(L.T)
order_T  = np.argsort(_.real)
Lvecs    = Lvecs[:, order_T]
# Bi-normalize so Lvecs[:,n] . R[:,n] = 1
Lvecs = Lvecs / np.einsum('ij,ij->j', Lvecs, R)

def resolvent_casida(omega, eta=0.02):
    z = omega + 1j*eta
    out = 0.0 + 0.0j
    for n in range(dim):
        out += (w0 @ R[:, n]) * (Lvecs[:, n] @ v0) / (eigvals[n] - z)
    return out

# ─── Dipole-perturbation and readout vectors ─────────────────────────────
d_ia = np.array([0.7, -0.3, 0.5, 0.2, 0.9, -0.4])
v0 = np.concatenate([d_ia, d_ia]).astype(complex)
w0 = v0.copy()

# ─── Bi-orthogonal Lanczos ───────────────────────────────────────────────
def biortho_lanczos(L, v0, w0, n_iter):
    """Returns (alphas, betas, gammas, <w0|v0>) such that
       <w0|(L - z)^{-1}|v0> = <w0|v0> * [(T - z)^{-1}]_{00}
       where T is the (non-symmetric) tridiagonal projection of L."""
    seed  = np.vdot(w0, v0)
    scale = np.sqrt(seed + 0j)
    q     = v0 / scale
    p     = w0 / scale              # <p|q> = 1

    q_prev = np.zeros_like(q);  p_prev = np.zeros_like(p)
    beta_prev = 0+0j;           gamma_prev = 0+0j
    alphas, betas, gammas = [], [], []

    for _ in range(n_iter):
        Lq, LTp = L @ q, L.T @ p
        alpha   = p @ Lq
        r = Lq  - alpha*q - gamma_prev*q_prev
        s = LTp - alpha*p - beta_prev *p_prev
        rs = s @ r
        if abs(rs) < 1e-14: break
        beta  = np.sqrt(rs + 0j)
        gamma = rs / beta

        alphas.append(alpha);  betas.append(beta);  gammas.append(gamma)
        q_prev, p_prev = q, p
        q, p = r/beta, s/gamma
        beta_prev, gamma_prev = beta, gamma

    return (np.array(alphas, dtype=complex),
            np.array(betas,  dtype=complex),
            np.array(gammas, dtype=complex),
            seed)

def resolvent_lanczos(coeffs, omega, eta=0.02):
    alphas, betas, gammas, seed = coeffs
    z  = omega + 1j*eta
    cf = 0.0 + 0.0j
    for j in range(len(alphas)-1, 0, -1):
        cf = (betas[j-1] * gammas[j-1]) / ((alphas[j] - z) - cf)
    cf = 1.0 / ((alphas[0] - z) - cf)
    return seed * cf

coeffs = biortho_lanczos(L, v0, w0, n_iter=dim)

# ─── Compare at sample frequencies ───────────────────────────────────────
print("omega   Casida (Re,Im)              Lanczos (Re,Im)             |diff|")
for w in np.linspace(0.1, 3.5, 9):
    rc = resolvent_casida(w);  rl = resolvent_lanczos(coeffs, w)
    print(f"{w:5.2f}   ({rc.real:+.5f}, {rc.imag:+.5f})   "
          f"({rl.real:+.5f}, {rl.imag:+.5f})   {abs(rc-rl):.1e}")

# ─── Convergence in chain length ─────────────────────────────────────────
print("\nIm <w0|(L - 1.5 - i*0.02)^{-1}|v0> as Lanczos chain grows:")
for n_iter in [2, 4, 6, 8, 10, 12]:
    c = biortho_lanczos(L, v0, w0, n_iter=n_iter)
    print(f"  n_iter={n_iter:2d}: {resolvent_lanczos(c, 1.5).imag:+.6f}")
print(f"  Casida reference: {resolvent_casida(1.5).imag:+.6f}")

Running the script produces (the comparison columns are $⟨ w_{0} ∣ (L - ω - i η)^{- 1} ∣ v_{0} ⟩$ ):

omega   Casida (Re,Im)              Lanczos (Re,Im)             |diff|
 0.10   (+0.07842, +0.01576)   (+0.07842, +0.01576)   6.0e-15
 0.53   (+0.44010, +0.01920)   (+0.44010, +0.01920)   7.4e-15
 0.95   (+0.95882, +0.03025)   (+0.95882, +0.03025)   1.1e-14
 1.38   (+3.03769, +0.82129)   (+3.03769, +0.82129)   2.0e-14
 1.80   (+4.97732, +0.45700)   (+4.97732, +0.45700)   9.8e-14
 2.23   (-3.26978, +1.03485)   (-3.26978, +1.03485)   5.1e-14
 2.65   (-0.74125, +0.25100)   (-0.74125, +0.25100)   4.0e-14
 3.08   (-4.48742, +0.34972)   (-4.48742, +0.34972)   3.9e-14
 3.50   (-2.08047, +0.04186)   (-2.08047, +0.04186)   4.8e-16

Im <w0|(L - 1.5 - i*0.02)^{-1}|v0> as Lanczos chain grows:
  n_iter= 2: +0.058698
  n_iter= 4: +0.111647
  n_iter= 6: +0.539370
  n_iter= 8: +0.261432
  n_iter=10: +0.256980
  n_iter=12: +0.242786
  Casida reference: +0.242786

A few things to notice. (1) Casida and Lanczos agree to $\sim 1 0^{- 13}$ at every sampled frequency — the continued fraction is a genuine algebraic identity for the matrix element, not an approximation, when the chain is run to the full dimension. (2) The truncated-chain values at $ω = 1.5$ approach the Casida reference monotonically only at the end; intermediate chains can OVERSHOOT badly because a short Krylov subspace samples the spectrum non-uniformly. This is exactly what real TDDFT calculations see and why a few thousand Lanczos iterations are typical for a converged solid-state spectrum. (3) The Liouvillian's eigenvalues come in $\pm Ω$ pairs (printed when you run the script), which is the standard $L = [[A, B], [- B, - A]]$ symmetry — physical excitations and their de-excitation partners. The Lanczos chain implicitly encodes both branches; you don't extract them separately, you just read the resolvent off the continued fraction.

The same code, scaled to a million-dimensional matrix-free $L$ applied via Sternheimer-style orbital responses, is essentially what turboTDDFT does inside Quantum ESPRESSO. Nothing changes structurally between this 12×12 test and the production code — only that $L @ q$ becomes one SCF-like step on a batch of response vectors instead of a numpy matrix-vector multiply.

Why not Arnoldi?

Arnoldi iteration ALSO handles non-Hermitian operators, with one Krylov sequence and full Gram-Schmidt orthogonalisation. It is more numerically robust than bi-orthogonal Lanczos — no risk of "serious breakdown", no drifting bi-orthogonality, classical convergence guarantees. So why isn't TDDFT linear response done with Arnoldi? Four reasons, in roughly decreasing order of how decisive they are.

(1) Memory. Arnoldi orthogonalises every new Krylov vector against ALL previously built ones. At iteration $j$ you must store $j$ vectors of dimension $n$ , and that step costs $O (jn)$ . For a solid-state TDDFT calculation, $n$ — the dimension of a response object, roughly (occupied orbitals) × (basis-set size) × (k-points) × 2 — sits at $1 0^{7}$ to $1 0^{8}$ , and converged chains run $N \sim 1 0^{3} - 1 0^{4}$ iterations. Storing 10⁴ vectors of dimension 10⁸ is on the order of TERABYTES of RAM. Bi-orthogonal Lanczos keeps three right-Krylov vectors and three left-Krylov vectors at any time, plus the $N$ scalar coefficients. Constant memory regardless of chain length. That alone forces the choice in production codes.

(2) Continued-fraction structure. Arnoldi produces a HESSENBERG matrix, not a tridiagonal one. The (0,0) entry of the inverse of a Hessenberg matrix is still a rational function of its entries, but it isn't a scalar continued fraction — it is the ratio of two determinants of large sub-blocks. To evaluate the resolvent at each $ω$ in a frequency sweep you'd LU-factor the Hessenberg once and back-solve per $ω$ . That works, but it loses the elegant $O (N)$ -per-point bottom-up evaluation of the continued fraction. When you want the resolvent at thousands of frequencies (the whole spectrum), the tridiagonal-and-continued-fraction route is dramatically cheaper.

(3) Pseudo-Hermitian structure of the Liouvillian. The closed-shell TDDFT $L = [[A, B], [- B, - A]]$ is non-Hermitian but not generic — it has a symplectic-like symmetry inherited from its block structure (which is why its eigenvalues come in $\pm Ω$ pairs). Walker-Saitta-Gebauer-Baroni's "batch" Lanczos chooses the LEFT seed $w_{0}$ so the two Krylov sequences are related by this symmetry, and the resulting bi-orthogonal recursion is effectively SYMMETRIC Lanczos in disguise — real coefficients, no serious breakdown, only one mat-vec per step instead of two. Arnoldi would treat $L$ as a generic non-Hermitian operator and throw that structure away. Bi-orthogonal Lanczos preserves it; this is half of why WSGB's specific formulation is so well-conditioned in practice.

(4) What you actually want. Arnoldi shines when the goal is a few extremal EIGENVALUES of a non-Hermitian operator — GMRES, nonsymmetric ARPACK, Krylov-Schur for transient stability problems. For TDDFT linear response you do not want eigenvalues; you want one MATRIX ELEMENT of a resolvent, smoothed by a small broadening, over a continuous range of frequencies. The tools and the goal align: tridiagonal projection ↔ scalar continued fraction ↔ single matrix element. Using Arnoldi would be like running a full eigendecomposition when you only need a single matrix element.

Where Arnoldi DOES win: problems where you genuinely want eigenvalues rather than a resolvent. If you needed to assign individual excited states (transition dipoles, characters of specific peaks), or if the Liouvillian had a real INSTABILITY (a complex eigenvalue with positive imaginary part — relevant in non-Hermitian Hamiltonian physics and certain coupled-cluster response variants), Arnoldi or block-Arnoldi would be the right tool. For the standard "compute the response function across the visible spectrum" workflow in turboTDDFT and turboEELS, the Krylov-cheap-and-narrow tradeoff that bi-orthogonal Lanczos delivers is the right one.

Practical knobs: terminator and convergence

A naive truncation of the continued fraction at iteration N gives a spectrum that is a comb of N poles on the real axis — a SPIKY artifact, not a continuum. For a finite molecule that is actually correct (a true molecule has discrete excitations); for a solid it is wrong, and the spectrum should be a smooth function of $ω$ .

The Haydock TERMINATOR fixes this. After enough iterations, the Lanczos coefficients reach a plateau: $α_{j} \to a_{\infty}$ and $β_{j} γ_{j} \to (W /2)^{2}$ , with $W$ the effective bandwidth of the excitation spectrum. Replace the tail of the continued fraction by its analytic limit for an infinite chain with constant coefficients:

T_{\infty} (ω) = \frac{1}{2} [(a_{\infty} - ω) - (a_{\infty} - ω)^{2} - W^{2}],

selecting the branch that decays into the upper half plane (retarded boundary condition). Plugging $T_{\infty}$ in where the truncation cut-off would have been turns the comb of poles into a SMOOTH continuum that is physically correct for the solid-state spectrum.

# Terminator: when the Lanczos coefficients have converged to a
# plateau (a_j -> a_inf, b_j g_j -> (band/2)^2), the tail of the continued
# fraction can be replaced by its analytic limit instead of being truncated.
# For an infinite chain with constant coefficients (a, b*g = (W/2)^2),
#
#    T(z) = (a - z) - (W/2)^2 / T(z)   =>   T(z) = ( (a - z) +/- sqrt((a-z)^2 - W^2) ) / 2,
#
# picking the branch that decays at infinity. Replacing the cut-off
# in the continued fraction by T(z) is the "Haydock terminator" trick:
# the spectrum becomes a smooth continuum rather than a comb of poles,
# which is the correct physical answer for solids.

def haydock_terminator(a_inf, half_bandwidth, z):
    W = 2 * half_bandwidth
    disc = (a_inf - z)**2 - W**2
    s = np.sqrt(disc + 0j)
    # Choose branch so Im T(z) <= 0 (retarded boundary condition).
    T = 0.5 * ((a_inf - z) - s)
    if T.imag > 0:
        T = 0.5 * ((a_inf - z) + s)
    return T

A small imaginary broadening $η$ on top of that controls the resolution. The combination — finite-N Lanczos + Haydock terminator + small $η$ — is what makes the spectrum quantitative.

turboTDDFT and turboEELS

Two Quantum ESPRESSO modules implement this. turboTDDFT (Walker-Saitta-Gebauer-Baroni 2006; Rocca-Gebauer-Saad-Baroni 2008) targets optical absorption of finite systems and molecules: the perturbation is the uniform dipole, the readout is the dipole, and what you get is the molecular polarizability across the UV-visible window. Routine reach is a few thousand atoms — orders of magnitude beyond what direct Casida diagonalization can do.

turboEELS (Timrov-Vast-Gebauer-Baroni 2013, 2015) extends the same machinery to ELECTRON ENERGY LOSS SPECTROSCOPY in periodic systems. The probe is a fast electron passing through the solid; the relevant susceptibility is $χ (q, ω)$ at finite momentum transfer $q$ . The Liouvillian gets a $q$ -dependence (you're perturbing with $e^{i q \cdot r}$ rather than $r$ ), and the readout is the same plane wave. The loss function $- Im ε^{- 1} (q, ω)$ follows from $χ (q, ω)$ via the standard inversion. Plasmon dispersions in metals, collective modes in graphene, momentum-resolved spectra in transition-metal oxides — all done with one Lanczos chain per $q$ .

A key practical point: the chain is BUILT ONCE per perturbation direction (or per $q$ for EELS), and yields the spectrum at every frequency in a sweep. If you want three Cartesian components of the polarizability you run three chains; if you want EELS along a high-symmetry line you run one chain per $q$ -point. Compare to Casida, where the full diagonalization has to be done once and then the matrix elements with each perturbation are post-processed — Liouville-Lanczos is favourable specifically when the spectrum is wide, the system is large, and you want a CONTINUOUS function of frequency rather than a list of states.

What goes in, what comes out

Input: a converged DFT ground state (occupied orbitals ${ψ_{n}}$ , density $ρ_{0}$ , KS potential), a choice of perturbation (dipole direction for optics, $q$ for EELS), an XC functional / kernel, a number of Lanczos iterations $N$ (typically a few thousand for a converged absorption spectrum, sometimes more for sharp plasmons), a broadening $η$ , and a frequency window.

Output: the chain coefficients ${α_{j}, β_{j}, γ_{j}}_{j = 0}^{N - 1}$ , the dynamical polarizability $α (ω)$ (or response function $χ (q, ω)$ ) over the window, the absorption cross-section $σ (ω) \propto ω Im α (ω)$ or EELS loss function. A post-processing tool extracts oscillator strengths if you want individual transition assignments.

Where it sits in the perturbative landscape

Linear-response TDDFT is LITERALLY first-order perturbation theory — the small parameter is the external field strength. What's distinctive here is the choice of MACHINERY for the linear-response equation: instead of diagonalizing the linear operator (Casida) or inverting it at each frequency (Sternheimer), you tridiagonalise it with a Krylov method and then read the spectrum off the resulting continued fraction. The continued fraction plays the same role as a SUMMED PERTURBATION SERIES — except it is built from operator iterates, and converges for the same operator-theoretic reasons that the Lanczos approximation to a resolvent converges. There is a faint resonance with resurgence: in both cases the operator's analytic structure (poles in the resolvent, singularities of the Borel transform) is what the iterative scheme is implicitly resolving.

The connection back to the rest of the site:

Lanczos iteration — the symmetric tridiagonalisation that this page generalizes.
Bi-orthogonal Lanczos — the actual algorithm the chain in this page uses, written down on its own terms; resolvent matrix elements as continued fractions in the non-Hermitian setting.
Real-space vs second quantization — the response density-matrix lives most naturally in a real-space orbital representation; the Liouvillian is built directly from KS orbitals.
Perturbation theory and resurgence — same project from a different angle: extract a frequency-resolved spectrum from a perturbative expansion that, naively summed, doesn't converge.
Quantum physics overview — the resolvent $(ω - \hat{H})^{- 1}$ is the central object of scattering, response, and many-body Green's-function theory.

The deep takeaway: the Liouville-Lanczos method is a textbook example of MATCHING THE NUMERICAL TOOL TO THE ANALYTIC STRUCTURE. The thing you want — a continuous spectrum over a frequency window — is naturally encoded as a continued fraction. The thing your linear-response equation gives you — a giant non-Hermitian super-operator — is naturally tridiagonalised by Lanczos. The two structures fit together exactly, and the result is a method whose cost scales with the number of Lanczos iterations and the cost of applying the Liouvillian, not with the size of its spectrum. That is why it works on systems that Casida cannot touch.