The engine of Obara-Saika is that multiplying a Gaussian by a coordinate is the same as differentiating it: $x,e^{-\alpha x^2} = -\tfrac{1}{2\alpha},\tfrac{d}{dx}e^{-\alpha x^2}$. Climbing the angular-momentum ladder is therefore differentiating the Gaussian, and by Rodrigues’ formula the $n$-th derivative is the $n$-th Hermite polynomial:

$$\frac{d^n}{dx^n}e^{-x^2} = (-1)^n H_n(x), e^{-x^2}.$$

So the functions generated rung by rung are the Hermite polynomials, and they obey the Hermite three-term recurrence $x H_n = \tfrac12 H_{n+1} + n H_{n-1}$. The Obara-Saika overlap step

$$S(i{+}1,j) = (P-A),S(i,j) + \frac{1}{2p}\big[,i,S(i{-}1,j) + j,S(i,j{-}1),\big]$$

is that recurrence with the coordinate shifted to the product-Gaussian center $P$ (because two centers are involved). Same skeleton: a multiply-by-coordinate term plus lower-index corrections.

The common root

A three-term recurrence is exactly what “multiply by $x$” produces, because that operator raises degree by one and, written in an orthonormal-polynomial basis, is tridiagonal (a Jacobi matrix). For the Hermite weight $e^{-x^2}$ the off-diagonals are $\beta_n = \sqrt{n/2}$ — the same numbers Lanczos returns when run on that measure. So one structure wears many costumes:

$$\text{multiply-by-}x ;\to; \text{three-term recurrence} ;\to; \text{Jacobi matrix} ;\to; \text{orthogonal polynomials} ;\to; \text{Lanczos / Gauss quadrature} ;\to; \textbf{Obara-Saika}.$$

McMurchie-Davidson makes the link explicit by expanding Cartesian Gaussians in Hermite Gaussians, so OS and MD are the same Hermite recurrence kept two different ways.

Verification

Rodrigues checked to machine zero for $n = 1\ldots4$; the multiply-by-$x$ off-diagonals $\sqrt{n/2}$ match the Lanczos-on-Hermite result. Scripts: scripts/hermite_link.py, scripts/os_overlap.py.