The Gaussian Product Theorem is the Parallel Axis Theorem

The Gaussian product theorem (Boys 1950) is the algebraic identity that makes molecular quantum chemistry computationally tractable. It says that a product of two Gaussians centered at different points in space is itself a single Gaussian centered somewhere on the line segment between them, times a known prefactor:

with $\vec{P}=(\alpha \vec{A}+\beta \vec{B})/(\alpha +\beta )$ and $K=\exp (-\alpha \beta |\vec{A}-\vec{B}{|}^2/(\alpha +\beta ))$. This identity is what collapses 2-electron 4-center integrals down to one-center integrals plus algebra — every term you've ever written as (ab|cd) over a Gaussian basis is closed-form because of it. The whole edifice of molecular Gaussian-basis quantum chemistry rests on this one identity.

The thing nobody mentions when they introduce it: it's the parallel axis theorem from classical mechanics, applied to point masses at the Gaussian centers, then exponentiated. The derivation is mechanical, in both senses of the word.

The identity, without the exponential

Strip the exponentials. What's the corresponding identity for the quadratic exponents alone?

Sum of two weighted squared distances from $\vec{r}$ to two fixed points equals total-weight times squared distance from $\vec{r}$ to a third point, plus a constant that depends only on the geometry. The $\vec{r}$-dependence consolidates from two quadratic-in-$\vec{r}$ terms into one. That is what factors — and exponentiating an additive decomposition into a multiplicative one is what the exponential's there for.

The derivation, in three lines

Decompose $\vec{r}-\vec{A}=(\vec{r}-\vec{P})+(\vec{P}-\vec{A})$ and similarly for $\vec{B}$. Expand the squared norms:

Choose $\vec{P}$ to kill the cross-term:

What's left is $(\alpha +\beta )|\vec{r}-\vec{P}{|}^2$ (the $\vec{r}$-dependent part) plus the constant $\alpha |\vec{P}-\vec{A}{|}^2+\beta |\vec{P}-\vec{B}{|}^2$. Using $\vec{P}-\vec{A}=\beta (\vec{B}-\vec{A})/(\alpha +\beta )$ and $\vec{P}-\vec{B}=-\alpha (\vec{B}-\vec{A})/(\alpha +\beta )$, that constant simplifies to $\alpha \beta |\vec{A}-\vec{B}{|}^2/(\alpha +\beta )$. Three lines, identity proved. Exponentiate to recover the Gaussian product theorem; the constant becomes the prefactor $K$.

The parallel axis theorem is exactly this

Now read the identity again, with mechanical eyes. Consider a system of two point masses — mass $\alpha$ at position $\vec{A}$, mass $\beta$ at position $\vec{B}$. Their center of mass is at

The Gaussian-product center is the mechanical center of mass of two point masses at the orbital centers, with the Gaussian exponents playing the role of masses. The moment of inertia of this two-mass system about an arbitrary point $\vec{r}$ is

The moment of inertia about the center of mass is $I_{\text{CM}}=\alpha |\vec{P}-\vec{A}{|}^2+\beta |\vec{P}-\vec{B}{|}^2=\alpha \beta |\vec{A}-\vec{B}{|}^2/(\alpha +\beta )$. And the parallel axis theorem says

where $M=\alpha +\beta$ is the total mass. Substitute:

Same identity, line for line. The exponent of the product Gaussian is the moment of inertia of two fictional point masses, decomposed into a center-of-mass piece plus a fixed-geometry constant. The cross-term that vanishes in the derivation is the same cross-term that vanishes in classical mechanics whenever you decompose $(\vec{r}-\vec{A})$ into motion of the CM plus position relative to the CM — that's why the parallel axis theorem holds. The exponential in the Gaussian product theorem is just transporting this mechanical fact into a product form by converting addition to multiplication.

Verify

# Numerical verification of the Gaussian product theorem AND of the
# claim that it's the parallel axis theorem applied to point masses
# (α at A, β at B), then exponentiated.
import numpy as np

rng = np.random.default_rng(7)
A = rng.normal(size=3)
B = rng.normal(size=3)
alpha = 1.7
beta  = 0.6

# Center of mass / Gaussian-product center
M = alpha + beta
P = (alpha * A + beta * B) / M

# Moment of inertia about the center of mass
I_P_direct = alpha * np.sum((P - A)**2) + beta * np.sum((P - B)**2)
I_P_closed = (alpha * beta / M) * np.sum((A - B)**2)

# Moment of inertia about an arbitrary point r
r = rng.normal(size=3)
I_r_LHS   = alpha * np.sum((r - A)**2) + beta * np.sum((r - B)**2)
I_r_RHS   = M * np.sum((r - P)**2) + I_P_closed

print(f"I_P direct  : {I_P_direct:.12f}")
print(f"I_P closed  : {I_P_closed:.12f}")
print(f"  identity?  diff = {abs(I_P_direct - I_P_closed):.2e}")

print(f"\nI_r LHS  α|r-A|² + β|r-B|²       : {I_r_LHS:.12f}")
print(f"I_r RHS  M|r-P|² + αβ|A-B|²/M    : {I_r_RHS:.12f}")
print(f"  parallel axis identity? diff = {abs(I_r_LHS - I_r_RHS):.2e}")

# Exponentiate both sides — this IS the Gaussian product theorem.
LHS = np.exp(-alpha * np.sum((r - A)**2)) * np.exp(-beta * np.sum((r - B)**2))
K   = np.exp(-(alpha * beta / M) * np.sum((A - B)**2))
RHS = K * np.exp(-M * np.sum((r - P)**2))
print(f"\nexp(-α|r-A|²) · exp(-β|r-B|²)    : {LHS:.6e}")
print(f"K · exp(-M|r-P|²)                : {RHS:.6e}")
print(f"  Gaussian product theorem? rel diff = {abs(LHS - RHS) / LHS:.2e}")

I_P direct  : 0.832764079152
I_P closed  : 0.832764079152
  identity?  diff = 1.11e-16

I_r LHS  α|r-A|² + β|r-B|²       : 4.555635067194
I_r RHS  M|r-P|² + αβ|A-B|²/M    : 4.555635067194
  parallel axis identity? diff = 0.00e+00

exp(-α|r-A|²) · exp(-β|r-B|²)    : 1.050782e-02
K · exp(-M|r-P|²)                : 1.050782e-02
  Gaussian product theorem? rel diff = 3.30e-16

All three identities — the moment-of-inertia closed form, the parallel axis decomposition, and its exponentiation into the Gaussian product theorem — match to machine precision on a random geometry.

Why it's specific to quadratic exponents

The parallel axis theorem holds because squared norms have a clean decomposition under the substitution $\vec{r}-\vec{A}=(\vec{r}-\vec{P})+(\vec{P}-\vec{A})$. Expand $|\cdot {|}^2$ of a sum and you get the squared-norm of each piece plus a cross term that can be cancelled by choosing $\vec{P}$ as the appropriate weighted centroid. The same move on a linear norm $|\vec{r}-\vec{A}|$ (no square) does not decompose this way — distance is sub-additive, not additive, under shifts. The triangle inequality is the obstruction: $|\vec{r}-\vec{A}|\ne |\vec{r}-\vec{P}|+|\vec{P}-\vec{A}|$ in general.

Concretely: try the same identity with linear exponents. Is there a $\vec{P}$ and a constant $C$ such that $\alpha |\vec{r}-\vec{A}|+\beta |\vec{r}-\vec{B}|=(\alpha +\beta )|\vec{r}-\vec{P}|+C$ for all $\vec{r}$?

# What happens if you try the same trick with linear-in-distance exponents
# (STO-style)? The cleanest candidate is the no-constant version:
#   α|r-A| + β|r-B|  ?=  (α+β)|r-P|         (no parallel-axis constant term)
# Sample along a line through the centers.
import numpy as np

A = np.array([0.0, 0.0, 0.0])
B = np.array([2.0, 0.0, 0.0])
alpha, beta = 1.0, 0.5
M = alpha + beta
P = (alpha * A + beta * B) / M

print(f"P = {P}    (weighted centroid)")
print(f"{'x':>5}   α|r-A| + β|r-B|       (α+β)|r-P|          Δ")
for x in np.linspace(-2.0, 4.0, 13):
    r = np.array([x, 0.0, 0.0])
    LHS = alpha * np.linalg.norm(r - A) + beta * np.linalg.norm(r - B)
    RHS = M * np.linalg.norm(r - P)
    print(f"{x:>5.2f}   {LHS:>15.6f}    {RHS:>15.6f}    {LHS - RHS:>+8.4f}")

P = [0.66666667 0.         0.        ]    (weighted centroid)
    x   α|r-A| + β|r-B|       (α+β)|r-P|          Δ
-2.00          4.000000           4.000000     +0.0000
-1.50          3.250000           3.250000     +0.0000
-1.00          2.500000           2.500000     +0.0000
-0.50          1.750000           1.750000     +0.0000
 0.00          1.000000           1.000000     +0.0000
 0.50          1.250000           0.250000     +1.0000
 1.00          1.500000           0.500000     +1.0000
 1.50          1.750000           1.250000     +0.5000
 2.00          2.000000           2.000000     +0.0000
 2.50          2.750000           2.750000     +0.0000
 3.00          3.500000           3.500000     +0.0000
 3.50          4.250000           4.250000     +0.0000
 4.00          5.000000           5.000000     +0.0000

The identity holds exactly outside the segment $[\vec{A},\vec{B}]$ — both sides match to floating-point zero there. Inside the segment it fails by up to one full unit. The reason is geometric: the LHS has corners (kinks in the derivative) at $\vec{r}=\vec{A}$ and $\vec{r}=\vec{B}$ — wherever an absolute value passes through zero. The RHS has a corner only at $\vec{r}=\vec{P}$. Outside the convex hull of $\{\vec{A},\vec{B}\}$, neither side has any corners and both reduce to linear functions of distance with the same slope, so they coincide. Inside, the LHS has two corners and the RHS has one; no choice of $\vec{P}$ can move the LHS corners to overlap with the RHS corner, because there's only one of the latter and two of the former.

Squared norms have no corners — they're smooth everywhere — so the quadratic identity has no obstruction of this kind. The Pythagorean / parallel-axis decomposition holds globally because both sides are smooth quadratics, and a quadratic in $\vec{r}$ is fully determined by its center and an additive constant. Two free parameters in the RHS ($\vec{P}$ and the constant) match the two structural degrees of freedom in the LHS, with the cross-term cancellation fixing $\vec{P}$ uniquely. The linear case has one corner too many.

Slater-type orbitals have linear-in-distance exponents ($e^{-\zeta |\vec{r}-\vec{A}|}$), so they inherit this obstruction: no parallel-axis identity, no product theorem, no closed-form 4-center integrals. This is the structural reason quantum chemistry uses Gaussians despite their being physically wrong (no Kato cusp, wrong asymptotic decay). The Gaussian's quadratic exponent isn't an arbitrary algebraic convenience; it's the form that admits the only known multi-center reduction identity, and that identity is the parallel axis theorem.

What the product theorem buys

Apply the identity twice — once to the electron-1 pair, once to the electron-2 pair — and a 4-center 2-electron integral $(ab|cd)$ becomes

where $g_p,g_q$ are single-center Gaussians on $\vec{P}$ and $\vec{Q}$ respectively. What's left is a two-center Coulomb integral over single-center Gaussians, which evaluates in closed form via the Boys function $F_0(T)={\int }_0^1{e}^{-T{t}^2}dt$. Six floating-point exponentials and one Boys evaluation per integral — that's what makes molecular-scale Hartree-Fock and DFT runnable on laptops.

Without the parallel axis theorem hiding inside the Gaussian, none of this works. STO basis sets sit in the same Hilbert space and have better physics; what they lack is a mechanical identity buried in their algebraic form.

The same identity in nuclear physics

The parallel-axis identity also runs the standard separation of variables in the two-body harmonic oscillator — the foundation of the nuclear shell model — but used as a coordinate change rather than as algebraic factorization. Substituting ${\vec{r}}_1\to \vec{r}-\vec{A}$, ${\vec{r}}_2\to \vec{r}-\vec{B}$ (equal weights) converts the Moshinsky identity ${\vec{R}}^2+{\vec{r}}^2={\vec{r}}_1^2+{\vec{r}}_2^2$ into the equal-weight $\alpha =\beta =1$ case of the parallel-axis identity on this page. Boys uses fixed centers and a varying integration variable to factor a product; Moshinsky uses two varying particle coordinates to rotate into normal modes that separate the 2-body Hamiltonian. Same identity, two costumes. The full Moshinsky-transform treatment — Jacobi coordinates, the coordinate-rotation derivation, what the brackets buy in shell-model calculations — lives on its own page.

Both Boys and Moshinsky rest on the principal axis theorem of linear algebra: any quadratic form on ${\mathbb{R}}^n$ diagonalizes under an orthogonal transformation. The harmonic oscillator IS such a quadratic form on configuration space, which is why it dominates both nuclear and atomic structure — its Hamiltonian is the only potential shape that admits clean orthogonal-transformation separation. Coulomb doesn't, so QC lives with it by using Gaussians (quadratic exponents in the basis) and pushing the separation into the integrals rather than the wavefunctions. The harmonic shell model in nuclei bakes the separation into the wavefunctions directly.

Related on this site

• Two-electron integrals (from scratch) — the GPT applied twice plus the Boys function, derived in code for the H₂ molecule.
• 1s integrals — single-Gaussian special cases where the GPT collapses to identity ($\vec{A}=\vec{B}$).
• STO-nG fitting — how the field fakes STO physics inside Gaussian-product-theorem machinery by least-squares-fitting STOs with N Gaussians.
• Could AAA-fit rationals replace Slater orbitals? — the negative result motivated by trying to find a non-Gaussian product theorem. There isn't one, and this page is the reason why.
• Nuclear shell model — where the Moshinsky form of this identity earns its keep, separating the 2-body harmonic Hamiltonian into CM and relative oscillators so the interaction matrix elements can be evaluated as expansions in relative-coordinate eigenstates.