Three representations of 1/r₁₂ (and what they each compute)

The Coulomb kernel $1/r_{12}$ shows up in every 2-electron integral in quantum chemistry, in the Hartree-Fock Fock matrix, in every correlation calculation, in DFT exchange-correlation functionals. It's singular at $r_{12}=0$ and slowly decaying at infinity — neither tractable as a direct integrand nor in any standard function space (not $L^1$, not $L^2$). To compute integrals over it, you replace it with one of three integral representations. The choice of representation determines the algorithm.

Each representation comes from a different identity:

All three are equal as functions of $r_{12}$ away from the singularity at zero. They differ in how they encode that singularity and where they're convenient. Picking which one to substitute into an integral IS picking the algorithm:

• Multipole → fast multipole methods (FMM), used in classical molecular dynamics, periodic Coulomb sums, and hierarchical matrix factorizations of the Coulomb operator.
• Laplace → the Boys function, and through it the Obara-Saika and McMurchie-Davidson recurrences. This is the workhorse representation for molecular Gaussian-basis calculations.
• Fourier → plane-wave DFT, density fitting in $k$-space, Ewald summation, the resolution-of-identity Coulomb-metric methods.

The Laplace representation: where the Boys function comes from

The Laplace identity follows from a Gaussian-moment calculation. Start with the Gaussian integral ${\int }_0^{\infty }{e}^{-s^2{r}^2}ds=\sqrt{\pi }/(2r)$ (for $r>0$), rearrange, and you have $1/r=(2/\sqrt{\pi }){\int }_0^{\infty }{e}^{-s^2{r}^2}ds$. Trivial as identities go. What's remarkable is what it lets you do: if you substitute this into a 2-electron Coulomb integral against two Gaussian densities, the integrand becomes a product of three Gaussians (the two electron densities plus the Laplace-rep Gaussian $e^{-s^2{r}_{12}^2}$). All three are functions of Cartesian coordinates with quadratic exponents, so the 6D spatial integral collapses via the Gaussian product theorem applied twice. What's left is a single integral over the parameter $s$:

where $T=\eta R^2$, $\eta =pq/(p+q)$, and $R=|\vec{P}-\vec{Q}|$ is the distance between the Gaussian-product centers of the two electron-pair densities. The Boys function is a smooth, well-behaved 1D function — everything you compute via the Laplace representation eventually grounds out in $F_n(T)$ evaluations and the angular-momentum recurrences of OS or MD. The "magic" of the Boys function is really the magic of the Laplace representation: it turns a 6D singular integral into a 1D smooth integral.

The Fourier representation needs regularization

The Fourier identity says $1/r_{12}=(1/(2{\pi }^2))\int d^{3}k{e}^{i\vec{k}\cdot {\vec{r}}_{12}}/k^2$, but this is not a classical Fourier transform. $1/r$ is not in $L^1$ (the integral over ${\mathbb{R}}^3$ diverges at infinity, since $r^2\cdot 1/r=r$ doesn't decay) and not in $L^2$ either (same reason but worse). The naive Fourier integral $\int d^{3}r(1/r)e^{-i\vec{k}\cdot \vec{r}}$ doesn't converge. So where does $4\pi /k^2$ come from?

The Yukawa regularization

Multiply $1/r$ by a screening factor $e^{-\epsilon r}$ for some small $\epsilon >0$. The resulting Yukawa potential $V_{\epsilon }(r)=e^{-\epsilon r}/r$ decays exponentially at infinity, so it's in $L^2({\mathbb{R}}^3)$ for any $\epsilon >0$, and its Fourier transform exists classically. Direct calculation in spherical coordinates with the polar axis along $\vec{k}$:

This is everywhere finite. The singularity at $k=0$ is regularized — at $k=0$ we get $4\pi /{\epsilon }^2$, which is large but finite. Now take $\epsilon \to 0$:

The limit is well-defined pointwise for $k\ne 0$ and blows up at $k=0$. The pointwise limit is not a function in any classical sense — but it IS a tempered distribution. So $\mathcal{F}[1/r]=4\pi /k^2$ is meant distributionally, as the $\epsilon \to 0^{+}$ limit of classical Fourier transforms.

# Yukawa regularization in action: the Fourier transform of e^(-εr)/r is
# 4π/(k² + ε²), finite everywhere. As ε → 0 we recover 4π/k², singular at 0.
import math
import numpy as np
from scipy.integrate import quad

# Compute the FT of e^{-εr}/r at a few k values for shrinking ε:
ks = [0.1, 0.5, 1.0, 2.0]
for eps in [1.0, 0.1, 0.01, 0.001, 0.0]:
    fts = []
    for k in ks:
        if eps == 0.0:
            # The ε=0 case: 4π/k² (distributional)
            fts.append(4 * math.pi / k**2)
        else:
            fts.append(4 * math.pi / (k**2 + eps**2))
    print(f"ε = {eps:>6}:  F[e^{{-εr}}/r](k) at k = {ks}")
    print(f"            = {[f'{v:.4f}' for v in fts]}")

ε =    1.0:  F[e^{-εr}/r](k) at k = [0.1, 0.5, 1.0, 2.0]
            = ['12.4444', '10.0531', '6.2832', '2.5133']
ε =    0.1:  F[e^{-εr}/r](k) at k = [0.1, 0.5, 1.0, 2.0]
            = ['628.3185', '48.3322', '12.4444', '3.0807']
ε =   0.01:  F[e^{-εr}/r](k) at k = [0.1, 0.5, 1.0, 2.0]
            = ['1244.4337', '50.2401', '12.5663', '3.1417']
ε =  0.001:  F[e^{-εr}/r](k) at k = [0.1, 0.5, 1.0, 2.0]
            = ['1256.5163', '50.2655', '12.5664', '3.1416']
ε =    0.0:  F[e^{-εr}/r](k) at k = [0.1, 0.5, 1.0, 2.0]
            = ['1256.6371', '50.2655', '12.5664', '3.1416']

Watch the regularized FT approach its limit. At large $k$ every $\epsilon$ gives essentially the same answer ($4\pi /k^2$ dominates over ${\epsilon }^2$). At small $k$, the IR behavior depends sharply on $\epsilon$: as $\epsilon$ shrinks the value at $k=0.1$ climbs from 12 to 1257, blowing up eventually. The singularity at $k=0$ is real but it's the limit of finite quantities, not a problem to be diagnosed and fixed.

Why "infinite singularity, finite answer" works

Now the question that's been hovering: how can integrating against a singular kernel produce finite numbers? Three nested reasons, each operator-theoretic.

The $k^2$ volume element saves you in 3D

Look at the Fourier-rep Coulomb integral against two Gaussian pair densities ${\rho }_{ab}(r_1),{\rho }_{cd}(r_2)$:

The integrand has a $1/k^2$ singularity at the origin. But the volume element in 3D spherical coordinates is $d^{3}k=k^{2}dkd\Omega$ — the $k^2$ Jacobian exactly cancels the $1/k^2$ in the kernel. What's left is a radial integrand that's smooth at $k=0$:

The singularity is removable by 3D-spherical geometry. This is also why the Coulomb potential is logarithmic in 2D (where $d^{2}k=kdkd\theta$, only one factor of $k$ to cancel $1/k^2$, leaving $1/k$ which integrates logarithmically) and linear in 1D (no $k$ in the volume element, $1/k^2$ stays singular). Coulomb in 3D is geometrically blessed: the kernel's IR singularity is exactly the strength needed for the Jacobian to remove it.

The Coulomb kernel as a distribution acting on test functions

Distributions don't have pointwise values; they have action on test functions. The kernel $4\pi /k^2$ is an "improper" function (singular at zero), but it's a perfectly good distribution: paired with a Schwartz function $\varphi (\vec{k})$, the result $\int d^{3}k(4\pi /k^2)\varphi (\vec{k})$ is finite as long as $\varphi$ is smooth enough at the origin. Gaussian-derived $\tilde{\rho }(\vec{k})$s are Schwartz functions — exponentially fast-decaying at infinity, smooth and bounded at the origin — so the pairing is finite by construction.

Operator-theoretically: $1/(4\pi r)$ is the integral kernel of $(-\Delta {)}^{-1}$ on $L^2({\mathbb{R}}^3)$. The Laplacian $-\Delta$ has continuous spectrum $[0,\infty )$, with no eigenvalues. The inverse $(-\Delta {)}^{-1}$ is unbounded (it sends the spectral edge at zero to infinity), but it's well-defined as a map between appropriate function spaces.

The Yukawa regularization shifts the spectrum: $-\Delta +{\epsilon }^2$ has spectrum $[{\epsilon }^2,\infty )$, bounded below by ${\epsilon }^2>0$, so $(-\Delta +{\epsilon }^2{)}^{-1}$ is a bounded operator with norm $1/{\epsilon }^2$. Its integral kernel is exactly $e^{-\epsilon r}/(4\pi r)$, the Yukawa potential. Removing the regularization means letting the spectral edge return to zero, which is when the resolvent picks up its long-range Coulomb tail. The "singularity" of $1/r$ at large $r$ is the spectral edge of the Laplacian sitting at zero — a structural fact about the operator, not an algebraic accident.

The IR singularity is about charge neutrality

The $4\pi /k^2$ singularity at $k=0$ is the Fourier-space encoding of the $1/r$ tail in position space. It only causes trouble when a calculation involves non-neutral charge distributions — because $\tilde{\rho }(0)=\int \rho (r)d^{3}r$ is the total charge. For neutral systems ($\tilde{\rho }(0)=0$), the $1/k^2$ singularity is multiplied by zero at the problematic point, and everything stays finite.

This is exactly why the Madelung sum for a periodic ionic crystal is conditionally convergent term by term but finite overall: the long-range Coulomb tail is suppressed by neutrality of the unit cell. Ewald summation makes this manifest by splitting the Coulomb interaction into a short-range part (handled in real space, convergent locally) and a long-range part (handled in $k$-space, where the $\tilde{\rho }(0)=0$ condition kills the $k=0$ singularity for neutral systems). Molecular QC dodges this whole story because molecules sit in vacuum and have fast-decaying densities, but the operator-theoretic substrate is identical.

Is this a resummation?

Strictly, no — there's no divergent series being tamed. But the Yukawa regularization is structurally the same move as Abel summation: introduce a smoothing parameter ($\epsilon$, or $x$ for Abel: $\sum a_n\to \sum a_n{x}^n$ for $x<1$), compute the well-defined regularized object, then take the limit ($\epsilon \to 0^{+}$, or $x\to 1^{-}$) while extracting the finite part. The "result" is whatever the regularized object converges to as the regulator vanishes, and this is independent of which regulator you use — Yukawa ($e^{-\epsilon r}$), Gaussian ($e^{-\epsilon r^2}$), or analytic continuation in a dimensional parameter all give the same Coulomb FT.

This regulate-then-take-limit pattern is everywhere in physics:

• Casimir energy via $\zeta$-regularization. The vacuum energy ${\sum }_{n}n$ between conducting plates diverges term by term, but $\zeta (-1)=-1/12$ via analytic continuation of the $\zeta$ function gives the right Casimir force. Same pattern: regularize (in this case, by analytic continuation of the series in its index), compute, take limit. The $-1/12$ isn't summation in the literal sense; it's the finite part of an Abel-style regularization.
• Dimensional regularization in QFT. Compute integrals in $d=4-ϵ$ dimensions, where they're finite by power-counting; analytically continue back to $d=4$, subtract the $1/ϵ$ pole; the finite remainder is the renormalized answer. Direct kin of Yukawa regularization, just with $ϵ$ as a dimensional parameter rather than a mass.
• Ewald summation. Split $1/r=\text{erf}(\alpha r)/r+\text{erfc}(\alpha r)/r$. The $\text{erfc}$ piece is short-range (decays Gaussian-fast at large $r$, no IR problem), the $\text{erf}$ piece is smooth at the origin (no UV problem in $k$-space). Sum each over periodic images separately, recombine. This is a genuine resummation — the slowly-convergent real-space Coulomb lattice sum gets reorganized into two fast-convergent pieces. The intermediate parameter $\alpha$ is unphysical; only the sum matters.
• Watson resummation in scattering theory. Replace a slowly-convergent partial-wave sum ${\sum }_{\ell }{f}_{\ell }(\theta )$ with a contour integral that picks up the leading Regge poles. Discrete sum becomes continuous integral; the same Poisson-summation move that converts the multipole expansion of $1/r$ into the Fourier representation.
• Borel summation. A divergent $\sum a_n$ with factorial growth gets replaced by ${\int }_0^{\infty }{e}^{-t}B(t)dt$ where $B(t)=\sum a_n{t}^n/n!$ has a positive radius of convergence. Smooth out by dividing by $n!$, integrate, get a finite answer. Different mathematical machinery from Yukawa, but structurally the same pattern: introduce a regularization, take a limit, extract the finite part.

The unifying principle in all of these: when a "naive" computation diverges, the divergence is often spurious — the limit of a sequence of well-defined regularized objects gives a finite answer that's regulator-independent. The Yukawa derivation of $\mathcal{F}[1/r]=4\pi /k^2$ isn't deeper than this pattern; it's an instance of it. So calling it "a kind of resummation" isn't crazy — it's the right intuition, just in the integral-representation flavor rather than the summed-series flavor.

Numerical verification: Laplace = Fourier

A concrete check: compute the same $(ss|ss)$ integral two ways. Use two Gaussian pair densities centered at $\vec{P}=(0,0,0)$ and $\vec{Q}=(0,0,1.5)$ with combined exponents $p=0.8$ and $q=1.0$. The Laplace route is the closed-form Boys-function expression; the Fourier route is direct radial quadrature of the $k$-space integrand.

# Compute the (ss|ss) Coulomb integral two ways: Laplace and Fourier.
# Both should produce the same number to machine precision.
import math
import numpy as np
from scipy.special import hyp1f1
from scipy.integrate import quad

# Two Gaussian pair densities ρ(r) = exp(-α|r - X|²) centered at P, Q
p, q = 0.8, 1.0
P = np.array([0.0, 0.0, 0.0])
Q = np.array([0.0, 0.0, 1.5])
R = np.linalg.norm(P - Q)
eta = p * q / (p + q)
T = eta * R ** 2

def boys(n, T):
    return (1.0 / (2 * n + 1)) * hyp1f1(n + 0.5, n + 1.5, -T)

# (1) Laplace rep:  closed form via Boys function (this is the OS/MD route)
I_laplace = (2 * math.pi ** 2.5) / (p * q * math.sqrt(p + q)) * boys(0, T)

# (2) Fourier rep:  ∫ d³k (1/k²) ρ̃*(k) ρ̃(k) integrated in spherical coords;
#                   the d³k = k² dk dΩ Jacobian cancels the 1/k² singularity.
sigma2 = 1.0 / (4 * eta)
radial = quad(lambda k: math.exp(-sigma2 * k**2)
                       * (math.sin(k * R) / (k * R) if k > 0 else 1.0),
              0.0, np.inf, limit=400)[0]
I_fourier = (math.pi / p) ** 1.5 * (math.pi / q) ** 1.5 * (2 / math.pi) * radial

print(f"T = η R² = {T}")
print(f"F_0(T) = {boys(0, T):.10f}")
print(f"I (Laplace)  = {I_laplace:.12e}")
print(f"I (Fourier)  = {I_fourier:.12e}")
print(f"diff         = {abs(I_laplace - I_fourier):.2e}")

T = η R² = 1.0
F_0(T) = 0.7468241328
I (Laplace)  = 2.434427140244e+01
I (Fourier)  = 2.434427140244e+01
diff         = 7.11e-15

Same number to 14 digits, computed via completely different machinery — one uses the Boys-function table, the other does a 1D radial $k$-space integral with the singular kernel $1/k^2$ handled implicitly by the volume element. Demonstrates that the choice of representation is operational, not mathematical: the answer is the same; the work is reorganized.

The multipole expansion is the Fourier rep, in disguise

To close the loop: the three representations aren't independent facts. The multipole expansion follows from the Fourier representation by doing the angular part of the $d^{3}k$ integral first using the spherical-wave expansion of the plane wave:

Substituting into the Fourier representation of $1/r_{12}$ and doing the angular $d{\Omega }_k$ integrals (which picks up $Y_{\ell m}(\hat{k})Y_{\ell m}^{\ast }(\hat{k})\to {\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}$) gives a radial $k$-integral over spherical Bessel functions, which evaluates to the $r_{<}^{\ell }/r_{>}^{\ell +1}$ multipole structure term by term. So the Fourier representation is the multipole expansion's discrete angular sum turned into a continuous spherical-harmonic completeness relation. Two costumes of the same decomposition.

Similarly, the Laplace representation is the Fourier representation with the spatial $d^{3}k$ integral done in a different order (radial first). All three are operationally equivalent ways to encode the same singular kernel; the choice is which order to do the integrals in, and that choice is the algorithm.

Related on this site

• Gaussian Product Theorem — what makes the Laplace representation compute cleanly against Gaussian basis functions.
• Obara-Saika and McMurchie-Davidson — both built on the Laplace rep + Boys function.
• Recursion strategies overview — where this page sits in the broader molecular-integral algorithmic landscape.
• Borel-Padé resummation — the divergent-series cousin of the Yukawa-regularization pattern on this page.
• Padé approximants — the original Abel-summation-flavored resummation, structurally akin to Yukawa.

References

• Boys, S. F. (1950). Electronic wavefunctions. I. A general method of calculation for the stationary states of any molecular system. Proc. Roy. Soc. London A, 200(1063), 542-554.
• Ewald, P. P. (1921). Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annalen der Physik, 369(3), 253-287.
• Jackson, J. D. (1999). Classical Electrodynamics, 3rd ed. Ch. 3 (Boundary-Value Problems with Spherical Coordinates) and Ch. 9 (Radiation: angular distributions and the multipole expansion).
• Reed, M. & Simon, B. (1975). Methods of Modern Mathematical Physics II: Fourier Analysis, Self-Adjointness. Ch. IX (Fourier transforms) for the operator-theoretic foundation of $1/r$ as a distributional FT.
• Helgaker, T., Jørgensen, P., Olsen, J. (2000). Molecular Electronic-Structure Theory, ch. 9 §9.6 (Hermite Gaussian functions and the McMurchie-Davidson scheme).