Green-Function Numerics

A numerical Green function is often a selected matrix element of a resolvent rather than an entire inverse matrix. Its broadening, pole structure, and spectral sum rules provide strong checks. This page treats the computational resolvent; many-body Green functions explains the operator and ensemble definitions.

Helpful background. Sparse Hamiltonians supplies matrix action; correlations supplies operator-applied source vectors.

Retarded resolvent by linear solves

For Hermitian $H$, let $z=E+i\eta$ with $\eta>0$ and define

$$G(z)=(zI-H)^{-1}.$$

To obtain $\langle a|G(z)|b\rangle$, solve $(zI-H)x=b$ and contract $a^\dagger x$. This avoids an unnecessary dense inverse. Check the relative residual $\|(zI-H)x-b\|/\|b\|$.

The spectral representation is

$$\langle b|G(z)|b\rangle =\sum_n\frac{|\langle n|b\rangle|^2}{E+i\eta-E_n}.$$

Its imaginary part is nonpositive. Consequently

$$A_b(E)=-\frac1\pi\operatorname{Im}\langle b|G(E+i\eta)|b\rangle =\sum_n|\langle n|b\rangle|^2 \frac{\eta/\pi}{(E-E_n)^2+\eta^2}$$

is nonnegative and integrates to $\|b\|^2$ over the full energy axis. A finite integration interval misses Lorentzian tails and must be included in the sum-rule error budget.

There is also a local, frequency-by-frequency identity:

$$-\operatorname{Im}\langle b|G(E+i\eta)|b\rangle =\eta\|G(E+i\eta)b\|^2.$$

It follows from $G-G^\dagger=-2i\eta G^\dagger G$ for Hermitian $H$. Compare the two sides using the solved vector; this checks the retarded sign and conjugation without a frequency integral. For unequal source and test vectors, the imaginary part of $\langle a|G|b\rangle$ need not have a definite sign.

If eigenvalues are available, the expected weight in a finite interval is known exactly:

$$\int_{E_{\min}}^{E_{\max}} A_b(E)\,dE =\sum_n\frac{|\langle n|b\rangle|^2}{\pi} \left[ \arctan\frac{E_{\max}-E_n}{\eta} -\arctan\frac{E_{\min}-E_n}{\eta} \right].$$

This distinguishes frequency-quadrature error from deliberately excluded tails. Comparing a short-window integral only with the full norm incorrectly attributes the latter to the solver.

Lanczos continued fraction

Normalize $q_0=b/\|b\|$ and Lanczos-tridiagonalize $H$ with diagonal $\alpha_j$ and off-diagonal $\beta_{j+1}$. Then the projected resolvent is

$$\langle b|G(z)|b\rangle\approx \frac{\|b\|^2}{ z-\alpha_0-\dfrac{\beta_1^2}{ z-\alpha_1-\dfrac{\beta_2^2}{z-\alpha_2-\cdots}}}.$$

Evaluate the finite continued fraction from the bottom upward. It is the $(0,0)$ element of the tridiagonal resolvent, obtainable by successive Schur complements.

For $b=0$, return the zero matrix element without normalizing. In exact arithmetic, an $m$-step Lanczos projection matches the source moments $\langle b|H^k|b\rangle$ through $k=2m-1$, provided the process has not already terminated on an invariant Krylov subspace. Testing these low moments is more informative than comparing only one resolvent value. Truncating the fraction defines a finite projected model; appending a bath or continuum terminator introduces an additional modeling assumption.

The source norm carries total spectral weight. Loss of Lanczos orthogonality can introduce duplicated or distorted poles; compare with explicit diagonalization on small matrices and monitor moments. Increasing Krylov dimension and reducing $\eta$ are distinct refinements.

Two-site spectral benchmark

For $H=\begin{pmatrix}0&-t\\-t&0\end{pmatrix}$ and $b=(1,0)^{\mathsf T}$,

$$G_{00}(z)=\frac{z}{z^2-t^2} =\frac12\left(\frac1{z-t}+\frac1{z+t}\right).$$

There are equal-weight poles at $\pm t$. Lanczos terminates after two vectors and reproduces this expression exactly up to arithmetic. Test the linear solve, spectral sum, and continued fraction against it.

At $t=0$ the spectrum becomes a single Lorentzian of unit weight. If the sign of $\eta$ is reversed, the spectral density becomes negative; this is a useful retarded-versus-advanced check.

Broadening and conditioning

For Hermitian $H$, $\|G(E+i\eta)\|\leq1/\eta$. As $\eta$ shrinks, poles sharpen and the solve becomes more sensitive near eigenvalues. A residual bound converts to an absolute solution-error bound of at most $\|r\|/\eta$; a small residual at coarse broadening does not certify a much smaller $\eta$.

For a selected matrix element, the corresponding absolute error is bounded by $\|a\|\|r\|/\eta$. This can set a solve tolerance from the required observable accuracy, although it may be conservative far from poles. Neither positivity nor a small residual eliminates the separate broadening and finite-size errors.

In a finite closed system $\eta$ is a numerical resolution parameter unless a physical lifetime has been derived. For a continuum limit, increase system size sufficiently that the level spacing is smaller than the intended resolution before reducing broadening. Reversing the order exposes isolated finite-size poles.

For many-body particle addition and removal, the relevant Hamiltonians act in neighboring particle-number sectors and energy shifts involve the reference state. A single-particle resolvent of a many-body matrix is not automatically the physical fermion Green function.

Exercise

Derive the first two large-$z$ terms and explain their use as checks.

Solution

$(zI-H)^{-1}=I/z+H/z^2+O(z^{-3})$, so $\langle b|G(z)|b\rangle=\|b\|^2/z+\langle b|H|b\rangle/z^2+\cdots$. These coefficients test spectral normalization and the first energy moment independently of detailed peak shapes.

References

• Haydock, R., Heine, V., and Kelly, M. J., “Electronic structure based on the local atomic environment for tight-binding bands,” Journal of Physics C 5, 2845–2858 (1972), doi:10.1088/0022-3719/5/20/004.
• Weiße, A., Wellein, G., Alvermann, A., and Fehske, H., “The kernel polynomial method,” Reviews of Modern Physics 78, 275–306 (2006), arXiv:cond-mat/0504627.