Saddle points and eigenvalue densities

Diagonalization replaced a matrix integral by an eigenvalue integral with a squared Vandermonde. It did not evaluate that integral. When a large parameter concentrates the weight, a saddle approximation can organize the remaining calculation. Its usefulness depends on specifying what is large, what is fixed, and which region of the integration contour contributes.

We first do this for one real variable, where an error can be bounded explicitly. Then we apply the same idea to many interacting eigenvalues. These are dimensionless illustrative integrals applicable to methods used in three- and four-dimensional field theory; they require no supersymmetry. The large-rank lesson supplies optional context about choosing a limiting sequence.

First calculation: correct a Gaussian saddle

For real $n>0$ and $g\geq0$, define

$$I_n(g)=\sqrt{\frac n{2\pi}} \int_{\mathbb R}dx\, \exp\!\left[-n\left(\frac{x^2}{2}+\frac{g x^4}{4}\right)\right]. \tag{T17.1}$$

The parameter $n$ need not be an integer. We take $n\to\infty$ at fixed $g$. The action has a unique real minimum at $x=0$, with second derivative one. Its quadratic term confines the important region to $x$ of order $n^{-1/2}$. Set $y=\sqrt n\,x$:

$$I_n(g)=\frac1{\sqrt{2\pi}}\int_{\mathbb R} e^{-y^2/2}\exp\!\left[-\frac{g y^4}{4n}\right]dy .$$

The leading Gaussian answer is one because of our explicit prefactor. Let $\langle\cdot\rangle_0$ denote expectation with this unit-variance Gaussian. Integration by parts gives $\langle y^{2p}\rangle_0=(2p-1)\langle y^{2p-2}\rangle_0$: differentiate $y^{2p-1}e^{-y^2/2}$ and use its vanishing endpoints. In particular $\langle y^4\rangle_0=3$ and $\langle y^8\rangle_0=105$. Expanding the second exponential therefore gives

$$I_n(g)\sim 1-\frac{3g}{4n}+\frac{105g^2}{32n^2}+\cdots . \tag{T17.2}$$

The negative first correction has a physical interpretation within this integral: a positive quartic term suppresses fluctuations relative to the Gaussian. This is not an exact equality to a truncated polynomial.

Here we can control the truncation globally. For $u\geq0$, Taylor’s integral remainder implies $0\leq e^{-u}-1+u\leq u^2/2$. Apply it with $u=gy^4/(4n)$ before doing the Gaussian integral:

$$\boxed{\quad 0\leq I_n(g)-\left(1-\frac{3g}{4n}\right) \leq\frac{105g^2}{32n^2}. \quad} \tag{T17.3}$$

For example, $g=1$, $n=20$ gives the certified interval $0.9625\leq I_{20}(1)\leq0.970703125$. Direct numerical integration gives $I_{20}(1)\approx0.9685121885$, inside it. This bound is on the normalized integral’s absolute error. A relative-error bound additionally needs a positive lower bound on the exact answer.

The all-order coefficient is $(-g)^p(4p-1)!!/(4^p p!\,n^p)$, with the $p=0$ coefficient one. Its growth makes the series divergent for any fixed nonzero $g/n$. Nevertheless each fixed truncation is an asymptotic approximation as $n\to\infty$. Convergence of an infinite series and accuracy of a finite asymptotic expansion are different questions.

Know what can invalidate the local calculation

The contour in (T17.1) is the whole real line. For $g<0$ the integral diverges there, even though a Gaussian expansion near zero can still be written. If instead $g$ grows proportionally to $n$, the quartic term in the rescaled integral remains of order one; the first correction no longer becomes small.

For a general real integral $\int e^{-nS(x)}dx$, a nondegenerate interior minimum $x_*$ with $S''(x_*)=a>0$ contributes $e^{-nS(x_*)}\sqrt{2\pi/(na)}$ at leading order. This local formula needs a convergent integral and control of the rest of its domain. Another accessible minimum with action larger by $\Delta S>0$ is exponentially suppressed only when $n\Delta S$ is large. Equal-depth minima must both be included. A minimum at an endpoint can contribute only part of a Gaussian; a vanishing second derivative changes the fluctuation scale.

For a complex integral, a steepest-descent contour runs through a stationary point along a direction in which the real part of the exponent decreases. The original contour must be deformable to the proposed one without crossing poles or changing endpoint conditions. Solving the stationary equation alone does not determine which saddles contribute.

This distinction matters in the assigned paper: Marino–Putrov, (3.13)–(3.16) use a saddle to invert a generating function. Their later (5.84)–(5.86) explicitly revisit the finite contour and corrections lost by extending it. We will not replace that physical contour problem by our positive real integral.

Keep the finite matrix model explicit

Let $M$ be an unrestricted Hermitian $N$-by-$N$ matrix, with ordinary trace and metric $ds^2=\operatorname{Tr}(dM^2)$. Retain its central direction: this is the $U(N)$ Gaussian, not the traceless $SU(N)$ model. For fixed dimensionless $s>0$, define

$$Z_N(s)=\int dV(M)\, e^{-N\operatorname{Tr}M^2/(2s)} =\frac{C_N}{N!}\int_{\mathbb R^N} \prod_i d\lambda_i\,\Delta(\lambda)^2 e^{-N\sum_i V(\lambda_i)},\qquad V(x)=\frac{x^2}{2s}. \tag{T17.4}$$

$C_N$ is the angular constant fixed by this metric; the preceding lesson derived it explicitly at $N=2$. It and $1/N!$ affect absolute normalization but not the eigenvalue stationary equations. There are $N^2$ orthonormal real matrix coordinates, so direct Gaussian integration gives the exact check $Z_N(s)=(2\pi s/N)^{N^2/2}$.

For distinct eigenvalues the effective action is

$$S_N(\lambda)=N\sum_iV(\lambda_i) -2\sum_{i<j}\log|\lambda_i-\lambda_j|.$$

The second term expresses the Vandermonde’s suppression of collisions. Differentiating gives the finite-rank saddle equations

$$\frac{\lambda_i}{s} =\frac2N\sum_{j\ne i}\frac1{\lambda_i-\lambda_j}. \tag{T17.5}$$

The factor two comes from the squared Vandermonde. A permutation quotient cannot change it. At fixed $s$, the confining and pair terms are both of order $N^2$ for an eigenvalue distribution of order-one width.

Within an ordered chamber this saddle is a stable minimum. For a real displacement $u_i$, the Hessian obeys

$$u^T(\nabla^2S_N)u =\frac Ns\sum_i u_i^2 +2\sum_{i<j}\frac{(u_i-u_j)^2}{(\lambda_i-\lambda_j)^2}>0 \quad(u\ne0).$$

Confinement and the collision divergence keep the minimum away from the chamber’s boundaries. This finite-variable statement does not make evaluation of the full integral at its minimum exact.

Vary a normalized density

For one configuration the empirical density is the distribution $\rho_N(x)=N^{-1}\sum_i\delta(x-\lambda_i)$. It has integral one and consists of spikes. Its ensemble average is a different object. A smooth limiting density describes the leading distribution as $N$ grows, rather than replacing either finite-rank object by an exact compact curve.

At leading order the discrete action suggests the functional

$$\mathcal E[\rho] =\int V(x)\rho(x)\,dx -\iint \rho(x)\rho(y)\log|x-y|\,dx\,dy, \qquad \rho\geq0,\quad\int\rho\,dx=1 . \tag{T17.6}$$

The double integral is over a smooth density. One must not insert the empirical delta functions literally: their diagonal self-interactions would contain $\log0$, whereas the original sum excludes $i=j$. Treating that short-distance region and fluctuations is part of calculating corrections, not part of this leading replacement.

Vary $\rho$ while imposing normalization with a multiplier. At equilibrium,

$$V(x)-2\int\rho(y)\log|x-y|\,dy=\ell_0 \quad\hbox{on the occupied support}. \tag{T17.7}$$

The same left side must be at least $\ell_0$ outside: moving a little positive density to an empty region must not lower the energy. Differentiating inside a smooth part of the support yields

$$\frac{x}{s} =2\,\operatorname{PV}\!\int\frac{\rho(y)}{x-y}\,dy . \tag{T17.8}$$

The Cauchy principal value removes a symmetric interval $|y-x|<\epsilon$ and then takes $\epsilon\to0$. Without this prescription the two sides of the pole cannot be assigned unrelated divergent integrals. Normalization, positivity and the outside inequality are additional conditions; the singular-integral equation alone is not the complete problem.

Solve and check the semicircle

Try $\rho(x)=C\sqrt{A^2-x^2}$ for $|x|\leq A$ and zero elsewhere, with $A>0$. The half-disc area gives $1=C\pi A^2/2$, hence $C=2/(\pi A^2)$.

To evaluate its principal value, introduce the resolvent $G(z)=\int\rho(x)/(z-x)\,dx$ for $z$ off the support. For real $z>A$, substitute $x=A\cos\theta$. The elementary integral

$$\int_0^\pi\frac{d\theta}{z-A\cos\theta} =\int_0^\infty\frac{2\,dt}{(z-A)+(z+A)t^2} =\frac{\pi}{\sqrt{z^2-A^2}}$$

uses $t=\tan(\theta/2)$. Writing the numerator $\sin^2\theta$ in terms of $z-A\cos\theta$ then gives

$$G(z)=\frac2{A^2}\left(z-\sqrt{z^2-A^2}\right). \tag{T17.9}$$

Continue the square root with $\sqrt{z^2-A^2}\sim z$ at infinity. This branch choice ensures $G(z)\sim1/z$, the normalization condition; using the positive real root at large negative $z$ would fail it.

On the cut, $G(x\pm i0)=2(x\mp i\sqrt{A^2-x^2})/A^2$. The real part is the principal value, while the imaginary jump recovers the density. Equation (T17.8) therefore requires $x/s=4x/A^2$. Thus

$$\boxed{\displaystyle \rho_*(x)=\frac1{2\pi s}\sqrt{4s-x^2}, \qquad |x|\leq2\sqrt s .} \tag{T17.10}$$

The radius follows from the force equation, not normalization alone. This agrees with the Gaussian density in Drukker–Marino–Putrov, (2.32) after identifying their width parameter $\lambda=4s$. That comparison does not identify our illustrative integral with a Wilson-loop observable.

There is also an outside check. The derivative of the left side of (T17.7) is $\operatorname{sgn}(x)\sqrt{x^2-4s}/s$ for $|x|>2\sqrt s$. It increases in either outward direction, as required. The normalized moments are

$$\int x\,\rho_*(x)\,dx=0,\qquad \int x^2\rho_*(x)\,dx=s,\qquad \int x^4\rho_*(x)\,dx=2s^2 .$$

These follow by the same trigonometric substitution. Independently, differentiating the exact matrix normalization gives $\langle\operatorname{Tr}M^2/N\rangle=s$ at every $N$. Agreement is a useful test of the width and factor two.

We can see the finite-rank distinction explicitly at $N=2$. Integrating out one eigenvalue uses $\int (x-y)^2e^{-y^2/s}dy=\sqrt{\pi s}(x^2+s/2)$. After normalization, the exact averaged density is

$$\langle\rho_2(x)\rangle =\frac{e^{-x^2/s}}{2\sqrt{\pi s}} \left(1+\frac{2x^2}{s}\right).$$

It integrates to one and has second moment $s$, but is positive on the whole real line. The semicircle’s compact support is a limiting statement, not an exact finite-rank boundary.

The calculation constructs the leading equilibrium description. It supplies neither a general rigorous convergence theorem nor uniform error estimates near the endpoints. Recovering absolute free energies, finite-rank corrections or another model’s saddle requires retaining that model’s measure, contour and limiting sequence.

Exercises

1. Include the second saddle

For fixed real $h$, find the leading large-$n$ behavior of $J_n(h)=\int_{\mathbb R}\exp[-n(x^2-1)^2/4+hx]\,dx$. What goes wrong if only the positive minimum is retained at $h=0$? Why is the answer not automatically uniform when $h$ grows proportionally to $n$?

Hint 1

Find both minima of $(x^2-1)^2/4$ and the second derivative at each.

Hint 2

For fixed $h$, the factor $e^{hx}$ is a smooth amplitude on each shrinking Gaussian region.

Solution

The minima are $x=\pm1$, both with action zero and second derivative two. Each Gaussian width is of order $n^{-1/2}$. The two leading contributions are $e^h\sqrt{\pi/n}$ and $e^{-h}\sqrt{\pi/n}$, so

$$J_n(h)\sim2\cosh(h)\sqrt{\frac{\pi}{n}} \qquad(n\to\infty,\ h\ \hbox{fixed}).$$

At $h=0$, discarding either minimum loses half the leading answer. The stationary point at zero is a maximum of the action, not a third dominant real minimum.

If $h=n\eta$ with nonzero fixed $\eta$, the action becomes $(x^2-1)^2/4-\eta x$. Its stationary equation is $x(x^2-1)=\eta$, so the locations and depths must be reconsidered. A formula derived with a slowly varying fixed amplitude cannot simply be reused.

2. Shift the center without changing the width

Replace the Gaussian potential by $V_b(x)=(x-b)^2/(2s)$ with real fixed $b$. Determine the density, its support, and its first two moments. Check its force equation.

Hint 1

Translate every eigenvalue by the same amount. The Vandermonde depends only on differences.

Hint 2

Use $u=x-b$ in both the density and its principal-value integral.

Solution

The translated density is

$$\rho_b(x)=\frac{\sqrt{4s-(x-b)^2}}{2\pi s}, \qquad b-2\sqrt s\leq x\leq b+2\sqrt s .$$

Translation preserves normalization. The odd centered moment vanishes and the centered second moment is $s$, so $\int x\rho_b\,dx=b$ and $\int x^2\rho_b\,dx=b^2+s$.

Its force is $2\,\operatorname{PV}\int\rho_b(y)/(x-y)\,dy=(x-b)/s=V_b'(x)$. The width is unchanged; the center moved. This argument uses the unrestricted Hermitian model. Imposing tracelessness would forbid an arbitrary common shift and change the constrained problem.

3. Compare a finite-rank saddle with an average

Multiply (T17.5) by $\lambda_i$ and sum to find $\operatorname{Tr}M^2/N$ at the finite-$N$ eigenvalue saddle. Compare with its exact ensemble expectation. For $N=2$, find the two saddle eigenvalues explicitly.

Hint 1

Group the double sum into unordered pairs. Their numerator is $\lambda_i-\lambda_j$.

Hint 2

Summing the unmultiplied saddle equations first shows that the two eigenvalues have opposite signs.

Solution

Each unordered pair contributes $\lambda_i/(\lambda_i-\lambda_j)+\lambda_j/(\lambda_j-\lambda_i)=1$. Therefore

$$\frac1s\sum_i\lambda_i^2 =\frac2N\frac{N(N-1)}2=N-1,\qquad \left.\frac{\operatorname{Tr}M^2}{N}\right|_{\mathrm{saddle}} =s\left(1-\frac1N\right).$$

For $N=2$, write $\lambda_1=-\lambda_2=a>0$. The force equation is $a/s=1/(2a)$, so $a=\sqrt{s/2}$ and the saddle moment is $s/2$.

The exact average is instead $s$ at every rank, from differentiating $Z_N(s)$. Evaluating at the most probable eigenvalue configuration omits fluctuations and is not an exact expectation. The two results have the same leading limit at fixed $s$, but even their first corrections differ. This also explains why simply discretizing the semicircle is not a method for computing all finite-rank corrections.

The next time a matrix integral produces a proposed saddle, the checks are concrete: identify its contour and large parameter, retain the Vandermonde force, normalize the density, test the support conditions, and distinguish the approximation from the original finite integral. The physical applications add their own interactions and contours to this method.