N=4 Gaussian Wilson loop

Can a direct matrix calculation reproduce an exact field-theory observable, including its finite-color correction? This laboratory samples the Gaussian integral for the normalized fundamental half-BPS circle. You will compare its mean with the exact finite polynomial, remove the central mode, and measure the difference from the planar answer.

Prepare with the finite-rank Wilson loop, which owns the analytic derivation. The later planar comparison also uses the planar limit and strong coupling. This code evaluates the already localized integral. Sampling agreement tests that integral’s implementation; the localization argument is separate.

The parent is massless four-dimensional N=4 SYM, with sixteen real Poincaré supercharges in flat space. The physical non-Abelian target here is $SU(N)$, $N\geq2$, of gauge rank $N-1$. We first sample a $U(N)$ Hermitian matrix and then remove its trace. Neither the number of random matrices nor their computational precision counts supercharges or gauge rank.

Fix the distribution and insertion

All sampled variables are dimensionless. Keep $\lambda=g_{\rm YM}^2N\geq0$ and $M=2\pi\ell a$, where $\ell$ is the sphere radius and $a$ the dimensionful localization coordinate. For $\lambda>0$,

$$\langle W\rangle_U= \frac{\int dM\,e^{-2N\operatorname{Tr}M^2/\lambda} \,N^{-1}\operatorname{Tr}e^M} {\int dM\,e^{-2N\operatorname{Tr}M^2/\lambda}}. \tag{LAB07.1}$$

The contour is the full space of Hermitian matrices. The trace is the ordinary fundamental trace, divided by $N$ only in the insertion. At zero coupling the distribution collapses to $M=0$ and the continuous answer is one.

The formula and finite polynomial agree with Drukker–Gross (2.25)–(2.27). That paper’s Gaussian proposal and Pestun’s later localization proof have different logical roles. The prerequisite explains the theory and contour translation before this numerical calculation.

Set $v=\lambda/(4N)$. A diagonal entry has variance $v$. Each real and imaginary part of an off-diagonal entry has variance $v/2$: the trace square contains twice its absolute square. Draw independent real standard-normal matrices $X,Y$, put $A=X+iY$, and form

$$M=\frac{\sqrt v}{2}(A+A^\dagger). \tag{LAB07.2}$$

For example, $M_{ii}=\sqrt v X_{ii}$, while for $i<j$ its real part is $\sqrt v(X_{ij}+X_{ji})/2$. The two terms are independent, giving variance $v/2$. Different independent entries remain independent. This constructs the required measure without a Markov chain, equilibration time or autocorrelation correction.

The eigenvalues are not sampled independently. Their repulsion is already generated by diagonalizing the full random matrix. Multiplying these samples by another squared Vandermonde would count that measure factor twice.

Translate matrices into measurements

If $m_i$ are the eigenvalues of one sample, calculate

$$w_U=\frac1N\sum_i e^{m_i},\qquad c=\frac1N\sum_i m_i,\qquad w_{SU}=\frac1N\sum_i e^{m_i-c}. \tag{LAB07.3}$$

Subtracting $cI$ projects the matrix to its traceless part; it does not require a second diagonalization. The same Gaussian action induces the desired measure on that subspace. Keep the original $\lambda$ while projecting.

Mathematical quantity	Code name
Color number $N$	n
‘t Hooft coupling $\lambda$	coupling
Independent sample count $S$	count or --samples
Hermitian matrices and eigenvalues	matrices, eigenvalues
Central coordinate $c$	centers
One normalized insertion	u, su
Precision of analytic calculations	--dps

The complete sampling core is short:

a = rng.standard_normal((size,n,n)) + 1j*rng.standard_normal((size,n,n))
matrices = (a + a.conj().swapaxes(1,2)) * (math.sqrt(v)/2)
centers = np.trace(matrices, axis1=1, axis2=2).real/n
eigenvalues = np.linalg.eigvalsh(matrices)
u = np.exp(eigenvalues).mean(axis=1)
su = np.exp(eigenvalues-centers[:,None]).mean(axis=1)

The released script processes batches of 2,048 matrices and retains scalar measurements. Its random generator is explicitly PCG64 with seed 20261020 in the default run. The seed makes the calculation reproducible; it does not turn a random estimate into an exact result.

Download and reproduce

Download the complete lab (ZIP) with the program, pinned dependencies, reference output and run instructions.

Save the editable calculation, pinned requirements, and expected output in one folder. The tested environment is Python 3.13.7, NumPy 2.1.2 and mpmath 1.3.0. The calculation needs no network after installation and no external data file.

In PowerShell, run from that folder:

python -m venv .venv
.venv/Scripts/python.exe -m pip install -r requirements.txt
.venv/Scripts/python.exe wilson_loop.py

On macOS or Linux, create the environment with python3 -m venv .venv and replace the interpreter path in the other commands by .venv/bin/python. The default run uses $N=4$, $\lambda=4$, and $S=65{,}536$ independent matrices. Use --help to see the input controls. The full_matrix_targets values are computed at the requested precision and formatted as decimal strings, with 25 significant digits for the loop values. The sampling entries, including their target copies, use ordinary double precision. Minor last-digit differences can occur between numerical libraries and platforms.

The recorded run gives:

Observable	Sample mean	Estimated standard error	Exact target
$U(4)$ loop	1.59426797	0.00188056	1.59422774
$SU(4)$ loop	1.54544235	0.00094836	1.54517851
$c^2$	0.06254535	0.00034244	0.0625
$N^{-1}\operatorname{Tr}M^2$	1.00031549	0.00137948	1

The loop discrepancies are about $0.02$ and $0.28$ estimated standard errors. The exact targets therefore agree well with this run. The extra digits in the exact column do not imply that many digits of sampling accuracy.

Attach an uncertainty to the mean

For independent sample values $w_b$, the script uses

$$\bar w=\frac1S\sum_{b=1}^S w_b,\qquad \widehat{\mathrm{SE}}(\bar w)= \sqrt{\frac{\sum_b(w_b-\bar w)^2}{S(S-1)}}. \tag{LAB07.4}$$

This estimates the random error of the mean, rather than the spread of individual values. Under an adequate central-limit approximation, roughly twice this standard error gives a useful approximate interval. It is not a rigorous enclosure, and no particular run is guaranteed to lie within it. Exponentials amplify rare large eigenvalues; at stronger coupling a small sample can underestimate its own variance. Check additional seeds and sample sizes before trusting a small error estimate there.

The default output also reports its first 16,384 samples. Their $U(4)$ standard error is $0.00377611$, approximately twice the final error, as the $S^{-1/2}$ law predicts. Their mean is $1.59912639$. Individual discrepancies need not decrease monotonically, and nested prefixes are correlated; they are not independent replications.

The central-mode relation is

$$W_U=F_c W_{SU},\qquad F_c=e^{\lambda/(8N^2)}. \tag{LAB07.5}$$

It follows from the independent Gaussian coordinate $c$, of variance $\lambda/(4N^2)$, as in Drukker–Gross (5.1). To test it on correlated samples, form $d_b=w_{U,b}-F_cw_{SU,b}$ before computing its standard error. The default mean is $-0.00023199$ with error $0.00160110$, consistent with zero. Adding two error estimates as if the measurements came from separate samples would discard their covariance. A fit of $F_c$ to these data would also be a different test from using its independently predicted value.

Check the implementation by another method

The exact evaluator uses a positive finite sum. Start $t_0=1$ and update

$$\frac{t_{j+1}}{t_j} =\frac{\lambda(N-j-1)}{4N(j+1)(j+2)}, \qquad W_U=e^{\lambda/(8N)}\sum_{j=0}^{N-1}t_j. \tag{LAB07.6}$$

All terms are nonnegative for the stated coupling. This avoids cancellation between alternating polynomial terms. Multiplication by $e^{-\lambda/(8N^2)}$ gives the exact $SU(N)$ target. The script repeats the sum with twenty extra decimal digits as a precision probe. Raising --dps improves this part of the calculation; it does not change the double-precision random matrices or reduce their sampling noise.

An independent deterministic check uses a traceless two-by-two matrix $M=\boldsymbol r\cdot\boldsymbol\sigma$. Its eigenvalues are $\pm r$ and $\operatorname{Tr}M^2=2r^2$. The full three-dimensional Gaussian therefore reduces to

$$W_{SU(2)}= \frac{\int_0^\infty dr\,r^2e^{-8r^2/\lambda}\cosh r} {\int_0^\infty dr\,r^2e^{-8r^2/\lambda}} =\left(1+\frac\lambda{16}\right)e^{\lambda/32}. \tag{LAB07.7}$$

The script integrates the numerator and denominator separately, and independently integrates the $U(2)$ center. It compares these results with the finite sum at $\lambda=1,4,16$. The Abelian $U(1)$ answer $e^{\lambda/8}$ and zero-coupling limits provide further checks. A traceless one-by-one matrix is identically zero; its unit insertion is a mathematical check, not a claim of interacting $SU(1)$ gauge theory.

All 24 analytic comparisons pass at 60 decimal digits, with maximum scaled residual about $2.84\times10^{-61}$ against a $10^{-48}$ tolerance. This is observed high-precision agreement, not an interval-arithmetic quadrature proof. These checks have different errors from random sampling.

Planar comparison

The planar lesson gives $W_0=2I_1(\sqrt\lambda)/\sqrt\lambda$ after taking $N\to\infty$ at fixed $\lambda$. At $\lambda=4$, $W_0=1.59063685$. The deterministic finite-rank differences are:

$N$	$W_U/W_0-1$	$W_{SU}/W_0-1$
4	0.00225751	-0.02857871
8	0.00056407	-0.00722238

The $SU(4)$ central subtraction is appreciably larger than the $U(4)$ genus correction at this coupling. Their signs need not match. The approximately fourfold improvement on doubling $N$ is consistent with the fixed-coupling $N^{-2}$ expansion.

These numbers use exact formulas and contain no sampling noise. Raising the sample count reduces uncertainty around a fixed finite-rank target. Raising $N$ changes that target and its difference from the planar answer. Increasing arithmetic precision changes neither physical effect. Strong coupling is a further approximation and is not needed for this comparison.

Change one input

1. Double the color number

Run wilson_loop.py --n 8, leaving the coupling and sample count at their defaults. Predict which coupling parameter changes and which error is expected to improve. Is this the same as keeping $g_{\rm YM}$ fixed?

Hint 1

Use $g_{\rm YM}^2=\lambda/N$. The number of samples is a separate input.

Hint 2

Compare the exact finite-rank ratios in the table before assessing a random estimate.

Solution

$g_{\rm YM}^2$ changes from one to one-half. The prescribed experiment keeps $\lambda=4$, as required for the stated planar comparison. Holding $g_{\rm YM}$ fixed instead would double $\lambda$ and compare different entries of the planar function.

The recorded $N=8$ run gives $W_U=1.59204040\pm0.00093099$ and $W_{SU}=1.57982477\pm0.00050243$, where the uncertainties are estimated standard errors. The exact values are $1.59153409$ and $1.57914867$. Their deviations are about $0.54$ and $1.35$ standard errors. The deterministic planar error improves as shown above. The sampling variance also changes with the distribution, so equal sample counts do not require equal standard errors.

2. Remove the off-diagonal entries

Run wilson_loop.py --diagonal-only. This deliberately samples a different measure. Predict its $U(N)$ loop and its normalized trace-square mean. Can the central-factor test alone detect the error?

Hint 1

Each retained diagonal entry is an independent normal variable of variance $v$. Its exponential mean is $e^{v/2}$.

Hint 2

The mean of the diagonal entries still has variance $v/N$. Projecting that Gaussian vector removes one independent direction.

Solution

The wrong full-matrix model gives $e^{\lambda/(8N)}$ and $\langle N^{-1}\operatorname{Tr}M^2\rangle=\lambda/(4N)$. At the default inputs its loop is $1.13314845$, compared with the correct $1.59422774$. The recorded negative-control sample gives $1.13252527\pm0.00117789$, close to its own incorrect-model prediction.

After subtracting the mean diagonal entry, each remaining component has variance $v(1-1/N)$, so its loop is $e^{\lambda(N-1)/(8N^2)}$. The ratio of these two incorrect answers is still $F_c$. Thus the central-factor test passes even though the non-Abelian fluctuations are absent. A useful verification combines the trace-square normalization, an independent finite-rank answer and the center test. One successful identity does not certify the whole implementation.

This laboratory establishes a reproducible computation of a specified matrix observable, with distinct arithmetic, statistical and approximation checks. It does not infer a full field-theory spectrum from one loop or replace the supersymmetry argument that selected the integral.