Gaussian Variational Methods

Gaussian trial states turn many infinite-dimensional variational calculations into finite-dimensional matrix problems. Their normalization, Fourier transforms, derivatives, overlaps, and polynomial moments are analytic. A positive-definite width matrix can represent anisotropy and coordinate correlations, while linear combinations of Gaussians can be enlarged systematically.

The convenience has a cost. A single Gaussian is nodeless, has a Gaussian tail, and has no Coulomb cusp. It can be exact for quadratic Hamiltonians and useful near smooth minima, but it can be structurally wrong for weak binding, tunneling between distant wells, singular interactions, or excited states. The variational theorem preserves an upper bound; it does not erase those ansatz limitations.

Required background. The Variational Principle gives the energy bound, Variational Parameters explains optimization, and Gaussian Distributions and Spectral Decomposition supply the matrix integrals. Gaussian States treats Wigner functions and symplectic spectra; Time-Dependent Variational Principle develops packet evolution.

The complete one-dimensional quartic calculation, including optimization, upper-bound diagnostics, and numerical comparison, is Anharmonic Oscillator by Variational Methods.

When the Gaussian width is used as an adjustable reference frequency inside an order-by-order re-expansion, the method becomes Variational Perturbation Theory. Its first order reproduces a Gaussian Rayleigh quotient; higher orders are reorganized perturbative approximants and need not remain upper bounds.

General Real Gaussian Ansatz

Let

$$\mathbf x \in \mathbb R^d$$

collect $d$ Cartesian coordinates, or linear internal coordinates such as Jacobi coordinates with flat measure $d^dx$, and let $A$ be a real symmetric positive-definite $d\times d$ matrix. With

$$\mathbf y = \mathbf x-\mathbf q,$$

define

$$\psi_{A,\mathbf q}(\mathbf x) ={} \left( \frac{\det A}{\pi^d} \right)^{1/4} \times \exp\left[ -\frac12 \mathbf y^{\mathsf T}A\mathbf y \right].$$

The center is $\mathbf q$. The matrix $A$ is a precision matrix for the wavefunction amplitude; its eigenvalues have dimensions of inverse length squared. Positive definiteness ensures decay in every coordinate direction and makes the state normalizable. Curvilinear coordinates with nonconstant measure or kinetic metric require modified integrals.

Then the probability density is

$$\lvert\psi_{A,\mathbf q}(\mathbf x)\rvert^2 = \frac{\sqrt{\det A}}{\pi^{d/2}} e^{-\mathbf y^{\mathsf T}A\mathbf y}.$$

The multidimensional Gaussian integral

$$\int_{\mathbb R^d} e^{-\mathbf y^{\mathsf T}A\mathbf y} \,d^d y = \frac{\pi^{d/2}}{\sqrt{\det A}}$$

confirms normalization.

Width Matrix and Covariance

The position mean is

$$\langle\mathbf X\rangle = \mathbf q,$$

and the position covariance is

$$\Sigma_X = \left\langle (\mathbf X-\mathbf q) (\mathbf X-\mathbf q)^{\mathsf T} \right\rangle = \frac12A^{-1}.$$

Diagonalize the width matrix as

$$A = R \operatorname{diag}(a_1,\ldots,a_d) R^{\mathsf T},$$

where $R$ is orthogonal and every $a_j\gt0$. The columns of $R$ are the principal directions of the probability ellipsoid. Along principal direction $j$, the position standard deviation is

$$\sigma_j = \frac{1}{\sqrt{2a_j}}.$$

Off-diagonal entries of $A$ therefore have physical content: they rotate the principal axes and encode coordinate correlations. They are not merely inconvenient matrix elements to be discarded.

For $A=R\operatorname{diag}(a_1,a_2)R^{\mathsf T}$, the probability contours are ellipses whose principal standard deviations are $\sigma_j=(2a_j)^{-1/2}$. An off-diagonal entry in laboratory coordinates rotates the ellipse and permits correlations that a diagonal product ansatz cannot represent.

Width conventions

Three related matrices are easy to confuse:

Object	Matrix in the exponent	Position covariance
wavefunction amplitude	$A/2$	$A^{-1}/2$
probability density	$A$	$A^{-1}/2$
conventional normal density	$\Sigma_X^{-1}/2$	$\Sigma_X$

Since $\Sigma_X=A^{-1}/2$, the probability density can also be written

$$\lvert\psi\rvert^2 ={} \frac{1}{ \sqrt{(2\pi)^d\det\Sigma_X} } \times \exp\left[ -\frac12 \mathbf y^{\mathsf T} \Sigma_X^{-1} \mathbf y \right].$$

Stating the convention is more reliable than calling every parameter simply a width.

Momentum Covariance and Kinetic Energy

For the real Gaussian,

$$\nabla\psi = -A\mathbf y\,\psi.$$

The momentum mean vanishes,

$$\langle\mathbf P\rangle = \mathbf 0,$$

and direct differentiation gives

$$\left\langle P_iP_j \right\rangle = \frac{\hbar^2}{2}A_{ij}.$$

Thus the momentum covariance matrix is

$$\Sigma_P = \frac{\hbar^2}{2}A.$$

In the principal basis,

$$\sigma_{x,j} = \frac{1}{\sqrt{2a_j}}, \qquad \sigma_{p,j} = \hbar\sqrt{\frac{a_j}{2}},$$

so every principal pair saturates

$$\sigma_{x,j}\sigma_{p,j} = \frac{\hbar}{2}.$$

Consider a kinetic energy with a symmetric positive-definite mass matrix $M$,

$$T = \frac12 \mathbf P^{\mathsf T} M^{-1} \mathbf P.$$

Its expectation is

$$\langle T\rangle_A = \frac{\hbar^2}{4} \operatorname{Tr} \left( M^{-1}A \right).$$

Large eigenvalues of $A$ describe narrow coordinate directions and increase the kinetic energy. This is the matrix version of the $b^{-2}$ localization cost in the harmonic-oscillator variational estimate.

Potential Expectations

For a coordinate-space potential,

$$\langle V\rangle_{A,\mathbf q} = \int_{\mathbb R^d} V(\mathbf x) \lvert\psi_{A,\mathbf q}(\mathbf x)\rvert^2 \,d^d x.$$

If $V$ is sufficiently regular and its Gaussian expectation exists, this is a Gaussian smoothing of the potential. For analytic $V$, it can be organized formally as

$$\langle V\rangle_{A,\mathbf q} = \left. \exp\left[ \frac14 \nabla^{\mathsf T} A^{-1} \nabla \right] V(\mathbf x) \right|_{\mathbf x=\mathbf q}.$$

The first terms are

$$\langle V\rangle_{A,\mathbf q} ={} V(\mathbf q) + \frac14 \operatorname{Tr} \left[ A^{-1} \nabla\nabla V(\mathbf q) \right] + \text{higher even derivatives}.$$

Odd centered moments vanish. Fourth moments obey the Gaussian pairing identity

$$\langle y_i y_j y_k y_\ell\rangle ={} (\Sigma_X)_{ij}(\Sigma_X)_{k\ell} + (\Sigma_X)_{ik}(\Sigma_X)_{j\ell} + (\Sigma_X)_{i\ell}(\Sigma_X)_{jk}.$$

This is the finite-dimensional Gaussian form of Wick or Isserlis pairing. It makes polynomial anharmonicities straightforward to evaluate, although their width optimization can still be nonlinear.

The center obeys a useful stationarity condition. Translating the normalized Gaussian gives

$$\frac{\partial E}{\partial q_i} = \left\langle \frac{\partial V}{\partial X_i} \right\rangle.$$

At an optimized center,

$$\langle\nabla V\rangle = \mathbf0.$$

This is an averaged force-balance equation, not generally the same as placing $\mathbf q$ at a pointwise minimum of $V$.

Exact Solution for Coupled Quadratic Hamiltonians

The matrix Gaussian becomes exact for a stable quadratic Hamiltonian. Let

$$H ={} \frac12 \mathbf P^{\mathsf T}M^{-1}\mathbf P + V_0 + \frac12 (\mathbf X-\mathbf q_0)^{\mathsf T} K (\mathbf X-\mathbf q_0),$$

where both $M$ and $K$ are real symmetric positive-definite matrices.

The center is minimized at $\mathbf q=\mathbf q_0$. The Gaussian energy is

$$E(A) ={} V_0 + \frac{\hbar^2}{4} \operatorname{Tr}(M^{-1}A) + \frac14 \operatorname{Tr}(KA^{-1}).$$

For a symmetric variation $\delta A$,

$$\delta A^{-1} = -A^{-1}(\delta A)A^{-1}.$$

Therefore

$$\delta E ={} \frac14 \operatorname{Tr} \left[ \hbar^2M^{-1}\delta A \right] - \frac14 \operatorname{Tr} \left[ A^{-1}KA^{-1}\delta A \right].$$

Stationarity for every symmetric $\delta A$ requires

$$\hbar^2M^{-1} = A^{-1}KA^{-1},$$

or equivalently

$$\hbar^2 A M^{-1}A = K.$$

Define the mass-weighted force matrix

$$C = M^{-1/2}KM^{-1/2}.$$

Its positive square root gives the unique positive-definite stationary width,

$$A_\star = \frac{1}{\hbar} M^{1/2} C^{1/2} M^{1/2}.$$

The eigenvalues of $C^{1/2}$ are the normal-mode frequencies $\omega_j$. At the optimum,

$$\langle T\rangle_\star = \langle V-V_0\rangle_\star = \frac{\hbar}{4} \operatorname{Tr}(C^{1/2}),$$

and

$$E_\star = V_0 + \frac{\hbar}{2} \operatorname{Tr}(C^{1/2}) = V_0 + \frac{\hbar}{2} \sum_{j=1}^{d} \omega_j.$$

This is the exact coupled-oscillator ground energy. The full width matrix automatically performs the normal-mode rotation and assigns the correct width to every mode.

If $K$ has a zero or negative direction, the assumptions fail. A zero mode has no normalizable oscillator ground state in that coordinate, while a negative direction signals an unstable quadratic potential. The positive-definite solution should not be continued blindly through either case.

Worked Example: Two Coupled Coordinates

Consider equal masses with

$$H ={} \frac{P_x^2+P_y^2}{2m} + \frac{m}{2} \left[ \omega^2(x^2+y^2) + 2gxy \right],$$

where

$$\lvert g\rvert \lt \omega^2$$

ensures stability. The normal coordinates are

$$q_+ = \frac{x+y}{\sqrt2}, \qquad q_- = \frac{x-y}{\sqrt2},$$

with frequencies

$$\omega_+ = \sqrt{\omega^2+g}, \qquad \omega_- = \sqrt{\omega^2-g}.$$

The exact Gaussian width matrix in the $(x,y)$ coordinates is

$$A_\star = \frac{m}{2\hbar} \begin{pmatrix} \omega_++\omega_- & \omega_+-\omega_- \\ \omega_+-\omega_- & \omega_++\omega_- \end{pmatrix}.$$

For $g\ne0$, the off-diagonal entry is nonzero. The exact ground state is correlated in the original coordinates even though it factorizes in the normal coordinates.

The exact energy is

$$E_0 = \frac{\hbar}{2} (\omega_++\omega_-).$$

Now restrict the trial family to an axis-aligned isotropic product,

$$\psi_a(x,y) = \left( \frac{a}{\pi} \right)^{1/2} e^{-a(x^2+y^2)/2}.$$

Because $\langle xy\rangle_a=0$, its energy is

$$E_{\mathrm{prod}}(a) = \frac{\hbar^2a}{2m} + \frac{m\omega^2}{2a}.$$

Optimization gives

$$a_\star = \frac{m\omega}{\hbar}, \qquad E_{\mathrm{prod},\star} = \hbar\omega.$$

For nonzero coupling,

$$\frac{\hbar}{2} \left( \sqrt{\omega^2+g} + \sqrt{\omega^2-g} \right) \lt \hbar\omega.$$

The restricted product misses correlation energy. Adding one off-diagonal width parameter is enough to recover the exact quadratic result.

Positive-Definite Parameterizations

Directly optimizing independent entries of $A$ can step outside the positive-definite cone, making the trial state nonnormalizable. A parameterization should enforce admissibility by construction.

Cholesky factor

Write

$$A = LL^{\mathsf T},$$

where $L$ is lower triangular and its diagonal entries are positive. Parameterize those diagonals as

$$L_{ii} = e^{\eta_i}.$$

This is efficient and guarantees $A\gt0$. The coordinates depend on the ordering and scaling of the original variables, so preprocessing still matters.

Matrix exponential

For any real symmetric matrix $S$,

$$A = e^S$$

is positive definite. This provides global unconstrained coordinates, although differentiating the matrix exponential is more expensive than differentiating a Cholesky factor.

Eigenvalues and rotations

The decomposition

$$A = R \operatorname{diag}(e^{\eta_1},\ldots,e^{\eta_d}) R^{\mathsf T}$$

is physically interpretable: the $\eta_j$ set principal widths and $R$ sets orientation. It becomes coordinate-singular when eigenvalues coincide because rotations inside a degenerate eigenspace do not change $A$.

Parameterization	Main advantage	Main caution
$A=LL^{\mathsf T}$	efficient and stable	coordinate-order dependence
$A=e^S$	unconstrained symmetric parameters	costly matrix derivatives
eigenvalues plus rotations	direct geometric meaning	redundant rotations at degeneracy

The normalization depends on $A$. The identity

$$\delta\log\det A = \operatorname{Tr} (A^{-1}\delta A)$$

is essential when differentiating normalized Gaussian states. When differentiating a normalized expectation, freezing the determinant prefactor produces the wrong gradient. Alternatively, an unnormalized Gaussian and the full Rayleigh quotient give the same answer because its denominator supplies that normalization derivative.

Complex Widths and Phase-Space Correlations

A more general pure Gaussian wave packet uses a real positive-definite amplitude matrix $A$, a real symmetric chirp matrix $B$, a center $\mathbf q$, and a mean momentum $\mathbf p$:

$$\psi_{A,B,\mathbf q,\mathbf p}(\mathbf x) ={} \left( \frac{\det A}{\pi^d} \right)^{1/4} \times \exp\left[ -\frac12 \mathbf y^{\mathsf T}(A+iB)\mathbf y \right] \times \exp\left[ \frac{i}{\hbar} \mathbf p^{\mathsf T}\mathbf y \right].$$

The position covariance remains

$$\Sigma_X = \frac12A^{-1},$$

while the momentum covariance becomes

$$\Sigma_P = \frac{\hbar^2}{2} \left( A + B A^{-1}B \right).$$

The symmetrized position–momentum covariance block is

$$\Sigma_{XP} = -\frac{\hbar}{2} A^{-1}B.$$

Thus $B$ tilts the Gaussian in phase space and represents a quadratic phase or chirp. It is required for generic exact Gaussian dynamics even though it does not alter the position density at one instant.

For the static scalar Hamiltonian used above,

$$\langle T\rangle ={} \frac12 \mathbf p^{\mathsf T}M^{-1}\mathbf p + \frac{\hbar^2}{4} \operatorname{Tr} \left[ M^{-1} \left( A+B A^{-1}B \right) \right].$$

The $\mathbf p$ and $B$ contributions are nonnegative and do not change a coordinate-only potential expectation. In a time-reversal-invariant ground-state problem with no vector potential, their static optimum is therefore

$$\mathbf p_\star = \mathbf0, \qquad B_\star = 0.$$

Magnetic fields, imposed currents, angular momentum, and real-time propagation change that conclusion. The equations of motion for complex widths are derived in the Time-Dependent Variational Principle, while covariance-matrix state classification is developed in Gaussian States Preview.

Squeezed-State Connection

A centered real Gaussian with a width different from a reference oscillator vacuum is a squeezed pure Gaussian state. In one dimension, if the reference precision is $a_0$ and

$$a = a_0e^{2r},$$

then

$$\frac{(\Delta X)^2}{(\Delta X)_0^2} = e^{-2r}, \qquad \frac{(\Delta P)^2}{(\Delta P)_0^2} = e^{2r}.$$

One quadrature narrows while its conjugate broadens. A multidimensional $A$ combines rotations and mode-dependent squeezing. This connection explains why Gaussian variational widths have a direct quantum-state interpretation, but it does not mean that every Gaussian variational problem should be reformulated in quantum-optics language. The canonical operator and noise discussion is in Squeezed States: First Encounter.

Linear Combinations of Gaussians

A single Gaussian is one nonlinear trial manifold. A Gaussian basis enlarges the state to

$$\Psi(\mathbf x) = \sum_{n=1}^{N} c_n \phi_n(\mathbf x),$$

where each $\phi_n$ may have its own center, width matrix, polynomial prefactor, and symmetry projection.

For fixed nonlinear Gaussian parameters, optimizing the coefficients gives the generalized eigenproblem

$$Hc = ESc.$$

An outer optimization then changes centers and width matrices. Separating the exact inner coefficient solve from the nonlinear outer search is usually more stable than treating all variables identically.

For two centered normalized real Gaussians, the overlap is

$$S(A,B) = \frac{ 2^{d/2} (\det A\det B)^{1/4} }{ \sqrt{\det(A+B)} }.$$

As $A$ and $B$ become nearly equal, their overlap approaches one. Large Gaussian bases can therefore develop near-linear dependence and an ill-conditioned overlap matrix.

Gaussian product theorem

For unnormalized one-dimensional functions,

$$g_{\alpha,A}(x) = e^{-\alpha(x-A)^2}.$$

Similarly,

$$g_{\beta,B}(x) = e^{-\beta(x-B)^2}.$$

Their product is

$$g_{\alpha,A}(x)g_{\beta,B}(x) = C_{AB} e^{-(\alpha+\beta)(x-P)^2},$$

where

$$P = \frac{\alpha A+\beta B}{\alpha+\beta}.$$

The separation-dependent factor is

$$C_{AB} = \exp\left[ -\frac{\alpha\beta}{\alpha+\beta} (A-B)^2 \right].$$

Products centered on different points reduce to one Gaussian centered at a weighted average. This identity is a major reason multicenter molecular integrals remain tractable.

Symmetry and Coordinate Choices

The basic centered Gaussian is even and nodeless. It is naturally suited to bosonic or spatial ground states, but not by itself to odd parity, nonzero angular momentum, or fermionic antisymmetry.

Common extensions include:

• multiplying by polynomials to create nodes and angular structure;
• projecting onto parity or rotational quantum numbers;
• antisymmetrizing products or using Slater determinants;
• combining displaced Gaussians to describe multiple centers;
• using explicitly correlated coordinates such as interparticle separations;
• separating center-of-mass and internal coordinates before optimizing widths.

For a translation-invariant many-particle Hamiltonian, a positive-definite Gaussian in all laboratory coordinates artificially localizes the free center of mass. Use Jacobi or other internal coordinates, or factor the center-of-mass motion explicitly. A singular width matrix in the full coordinate space does not define an ordinary normalized wavefunction there.

Coordinate scaling matters as well. If one coordinate is measured in ångströms and another represents a collective mode with a very different natural length, optimizing raw entries of $A$ can create severe conditioning problems. Mass weighting and nondimensionalization should precede nonlinear optimization.

Limitations and Failure Modes

Wrong asymptotic tail

A Gaussian decays as

$$e^{-\alpha r^2},$$

faster than a typical short-range bound-state tail

$$e^{-\kappa r}.$$

This mismatch can strongly affect weak binding, tunneling amplitudes, polarizabilities, and large-distance observables even when the energy looks good.

Missing Coulomb cusp

At an electron–nucleus or electron–electron coalescence, the exact wavefunction satisfies a cusp relation. A smooth Gaussian centered at the coalescence has zero radial derivative there and cannot satisfy a nonzero cusp. No finite sum of smooth Gaussians reproduces an exact nonzero spherical cusp at the origin. An increasing basis can approximate the wavefunction and energy in appropriate norms, but convergence of contact properties may be slow without explicit correlation or cusp-corrected factors.

Nodes and sign structure

A single real Gaussian is positive everywhere. Excited states, fermionic states, and angular-momentum sectors require polynomial factors, symmetry projections, or signed combinations. Optimizing a nodeless Gaussian cannot discover a nodal surface absent from the ansatz.

Multimodal states

A state localized in two distant wells is poorly represented by one ellipse. A sum of displaced Gaussians can represent both lobes and their relative phase; a single covariance matrix cannot.

Nonlinear optimization and redundancy

Centers and widths create a nonconvex optimization problem. Permuting identical Gaussians leaves the state unchanged, coincident functions create nearly flat directions, and extreme widths can make matrix elements ill-conditioned. Report overlap spectra, gradient norms, multiple-start checks, and refinement stability.

Upper bounds do not validate observables

For an admissible normalized state and exact matrix elements, the energy remains an upper bound to the ground energy. A good bound does not certify tails, contact densities, transition amplitudes, or entanglement. Check observables that probe the physics the Gaussian family may miss.

Exercises

Derive the kinetic trace

For

$$T = \frac12 \mathbf P^{\mathsf T}M^{-1}\mathbf P,$$

show that a centered real Gaussian satisfies

$$\langle T\rangle = \frac{\hbar^2}{4} \operatorname{Tr}(M^{-1}A).$$

Solution

Differentiation gives

$$\partial_i\psi = -(A\mathbf y)_i\psi,$$

and

$$\partial_i\partial_j\psi = \left[ (A\mathbf y)_i(A\mathbf y)_j - A_{ij} \right] \psi.$$

Using $\langle\mathbf y\mathbf y^{\mathsf T}\rangle=A^{-1}/2$,

$$\left\langle (A\mathbf y)_i(A\mathbf y)_j \right\rangle = \frac12A_{ij}.$$

Since $P_iP_j=-\hbar^2\partial_i\partial_j$,

$$\langle P_iP_j\rangle = \frac{\hbar^2}{2}A_{ij}.$$

Contracting with $(M^{-1})_{ij}/2$ gives the trace formula.

Quantify the missed correlation energy

For the two-coordinate example, expand the exact ground energy for small

$$\delta = \frac{g}{\omega^2}.$$

Compare it with the optimized isotropic product energy through order $\delta^2$.

Solution

The exact energy is

$$E_0 = \frac{\hbar\omega}{2} \left( \sqrt{1+\delta} + \sqrt{1-\delta} \right).$$

Using

$$\sqrt{1\pm\delta} = 1 \pm \frac{\delta}{2} - \frac{\delta^2}{8} \pm \frac{\delta^3}{16} + O(\delta^4),$$

the odd terms cancel:

$$E_0 = \hbar\omega \left[ 1 - \frac{\delta^2}{8} + O(\delta^4) \right].$$

The optimized isotropic product gives $E_{\mathrm{prod}}=\hbar\omega$. Its leading excess is therefore

$$E_{\mathrm{prod}}-E_0 = \frac{\hbar\omega}{8} \delta^2 + O(\delta^4).$$

The first missed contribution is quadratic because the cross expectation $\langle xy\rangle$ vanishes in the uncorrelated state.

Show that a static chirp costs kinetic energy

For the complex-width Gaussian, prove that the $B$-dependent kinetic contribution is nonnegative and determine when it vanishes.

Solution

The extra term is

$$\Delta T_B = \frac{\hbar^2}{4} \operatorname{Tr} \left( M^{-1}B A^{-1}B \right).$$

Define

$$D = A^{-1/2}BM^{-1/2}.$$

Then cyclicity of the trace gives

$$\Delta T_B = \frac{\hbar^2}{4} \operatorname{Tr} \left( DD^{\mathsf T} \right).$$

This is a Frobenius norm squared:

$$\Delta T_B = \frac{\hbar^2}{4} \left\lVert A^{-1/2}BM^{-1/2} \right\rVert_F^2 \geq 0.$$

Since $A$ and $M$ are invertible, it vanishes exactly when $B=0$.

References

• S. F. Boys, “Electronic wave functions. I. A general method of calculation for the stationary states of any molecular system”, Proceedings of the Royal Society A 200, 542–554 (1950), a foundational Gaussian-orbital construction.
• T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory (Wiley, 2000), for Gaussian basis functions, integral technology, and conditioning in electronic-structure calculations.
• E. J. Heller, “Time-dependent approach to semiclassical dynamics”, Journal of Chemical Physics 62, 1544–1555 (1975), for evolving multidimensional Gaussian packets and complex width parameters.
• C. Lubich, From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis (European Mathematical Society, 2008), for variational Gaussian wave-packet dynamics and geometric numerical structure.
• R. Shankar, Principles of Quantum Mechanics, 2nd ed. (Springer, 1994), for the variational principle, oscillator, and Gaussian wave packets.
• Y. Suzuki and K. Varga, Stochastic Variational Approach to Quantum-Mechanical Few-Body Problems (Springer, 1998), for correlated Gaussian bases and nonlinear parameter selection.