\section{The computation}
\label{sec:computation}

We need the parts of the two-loop dilatation operator that connect two traces. We obtain them from the Hamiltonian of the theory on $\mathbb{R}\times S^3$, by ordinary quantum-mechanical perturbation theory. This section describes the Hamiltonian, the perturbation theory, the two pieces of $q_4$ and how each was evaluated.

\subsection{The Hamiltonian in Coulomb gauge}

We expand every field in spherical harmonics on $S^3$ \cite{ITT}. The coefficients are matrices in colour, and the theory becomes a quantum mechanics of infinitely many matrices. The free part $H_0$ is a sum of oscillators with the frequencies of table~\ref{tab:letters}. We use Coulomb gauge, $\nabla\cdot A=0$. It has two advantages here. Only physical polarisations propagate, so the letters of table~\ref{tab:letters} are the complete set. There are also no ghosts in the Hamiltonian formulation: $A_0$ is a Lagrange multiplier that imposes Gauss's law, and solving Gauss's law produces a determinant that cancels the Faddeev--Popov determinant exactly.

The price is the instantaneous Coulomb interaction. With $\rho$ the colour charge density of all the fields together,
\begin{equation}
H=H_0+g_{\rm YM}V_3+g_{\rm YM}^2V_4+H_C,\qquad H_C=\frac{g_{\rm YM}^2}2\int\!\!\int\rho^a\,K^{ab}[A]\,\rho^b .
\label{eq:Ham}
\end{equation}
Here $V_3$ contains the cubic gauge vertex, the couplings of the gauge field to scalars and fermions, and the Yukawa couplings. $V_4$ contains the quartic gauge vertex, the seagull and the scalar potential. The kernel $K$ is the Coulomb Green's function in the presence of the transverse gauge field. If $M=-\nabla\cdot D$ is the Faddeev--Popov operator, then $K=M^{-1}(-\nabla^2)M^{-1}$. Writing $M=\Delta+g_{\rm YM}\,a$, where $\Delta=-\nabla^2$ and $a$ is linear in $A$, so that $\Delta^{-1}$ is the free Coulomb propagator,
\begin{equation}
K=\Delta^{-1}-2g_{\rm YM}\,\Delta^{-1}a\,\Delta^{-1}+3g_{\rm YM}^2\,\Delta^{-1}a\,\Delta^{-1}a\,\Delta^{-1}-\ldots
\label{eq:K}
\end{equation}
The Coulomb interaction therefore contributes vertices with four, five and six fields,
\begin{equation}
C_2=\half\rho\,\Delta^{-1}\rho,\qquad C_3=-\rho\,\Delta^{-1}a\,\Delta^{-1}\rho,\qquad C_4=\tfrac32\rho\,\Delta^{-1}a\,\Delta^{-1}a\,\Delta^{-1}\rho ,
\label{eq:C234}
\end{equation}
at orders $g_{\rm YM}^2$, $g_{\rm YM}^3$ and $g_{\rm YM}^4$. All three are needed.

There is a question of operator ordering, since $K[A]$ multiplies charge densities that contain the momentum conjugate to $A$. Christ and Lee showed that the correct quantum Hamiltonian is the Weyl-ordered classical one plus terms of order $\hbar^2$ \cite{ChristLee}. At order $g_{\rm YM}^4$ those extra terms are quadratic in the fields. They correct the energy of a single gluon and cannot contribute to $b$. The rule is therefore to Weyl-order $H_C$ and then to normal-order it. The second step leaves remainders, in which two of the six fields of $C_4$ or of the five of $C_3$ are contracted with each other. They do matter, and they are part of the loop piece below.

\subsection{Fourth-order degenerate perturbation theory}

The spectrum of $H_0$ is highly degenerate. All words with the same total engineering dimension have the same energy. The dilatation operator of section~\ref{sec:oneloop} is the effective Hamiltonian that acts within one degenerate level and reproduces the exact spectrum order by order. Let $P$ project onto a level of energy $E_0$, let $Q=1-P$, and let $R=Q/(E_0-H_0)$ be the resolvent. Collect the interaction by powers of the coupling, $\cV_1=V_3$, $\cV_2=V_4+C_2$, $\cV_3=C_3$, $\cV_4=C_4$. Since $\cV_1$ changes the number of letters by one, $P\cV_1P=0$. The effective Hamiltonian at second order is
\begin{equation}
h_2=P\cV_2P+P\cV_1R\cV_1P .
\label{eq:h2}
\end{equation}
This is the literal one-loop operator of section~\ref{sec:anatomy}, up to the factor $g_{\rm YM}^2$. At fourth order, in the hermitian form of Bloch and des Cloizeaux \cite{Bloch,desCloizeaux},
\begin{equation}
\begin{aligned}
h_4={}&P\cV_4P+P\cV_1R\cV_3P+P\cV_3R\cV_1P+P\cV_2R\cV_2P\\
&+P\cV_1R\cV_1R\cV_2P+P\cV_1R\cV_2R\cV_1P+P\cV_2R\cV_1R\cV_1P\\
&+P\cV_1R\cV_1R\cV_1R\cV_1P-\half\big\{P\cV_1R^2\cV_1P,\;h_2\big\}.
\end{aligned}
\label{eq:h4}
\end{equation}
The last term is the familiar wave-function renormalisation of fourth-order perturbation theory. Up to normalisation, $h_2$ is $D_2$ and $h_4$ is $D_4$.

We do not need $h_4$ as an operator. We need its thermal average with four holonomies. It is efficient to obtain it together with the $\cH^2$ term. One expands $\Tr(x^H\hat U)$ directly to order $g_{\rm YM}^4$ in old-fashioned, time-ordered perturbation theory. The integrals over the Euclidean times of the vertices turn the Boltzmann factors into divided differences of $e^{-\beta E}$ over the energies of the successive intermediate states. When two intermediate states are degenerate the divided difference becomes a derivative, and a power of $\beta$ appears. The term linear in $\beta$ is $-\beta\vev{h_4}$ and the term quadratic in $\beta$ is $\half\beta^2\vev{h_2^2}_{\rm conn}$. We checked this rule, including the fermion signs and the degenerate cases, against exact diagonalisation of small matrix models with gauge groups $U(2)$ and $U(3)$, to ten digits.

The fact that the same diagrams give both terms is a strong check. The quadratic term must agree with $Q_4$ computed from the one-loop operator alone, which involves no energy denominators. It does, order by order in $y$, to nine digits.

\subsection{Which diagrams feed \texorpdfstring{$|u_1|^4$}{|u1|4}}

A term in $h_4$ is a product of letters $W$ and derivatives $\check W$ under traces. Its thermal average closes each $\check W$ onto a $W$ with a twisted propagator \eqref{eq:twisted}. We need four index loops with holonomy. Terms that act on a single letter give two. The planar terms that act on two neighbouring letters give three. Two classes give four.

\emph{(a) Three letters to three letters.} Connected tree diagrams with three incoming and three outgoing letters. With the vertices of \eqref{eq:Ham} and \eqref{eq:C234} the possible contents are: four cubic vertices; two cubic vertices and one quartic or $C_2$; two quartic or $C_2$; one cubic and one $C_3$; one $C_4$. Every colour ordering of the six external legs is needed, since the planar orderings give three index loops and only the others give four.

\emph{(b) Two letters to two letters at one loop.} The same contents with one more contraction: boxes, triangles, bubbles, a loop through $C_3$, the self-contraction of $C_4$, and the ordering remainders. The colour factor of a loop of adjoint fields with four external legs is
\begin{equation}
\begin{aligned}
\tr_{\rm adj}(T^aT^bT^cT^d)={}&N\big[\tr(abcd)+\tr(dcba)\big]\\
&+2\big[\tr(ab)\tr(cd)+\tr(ac)\tr(bd)+\tr(ad)\tr(bc)\big]+\ldots,
\end{aligned}
\label{eq:colour}
\end{equation}
where the dots are terms with a single generator under a trace, which vanish for $SU(N)$. The single-trace part is planar. The double-trace part is what we want. After the thermal average only the forward structure ${:}\tr(W_A\check W^A)\tr(W_B\check W^B){:}$ survives with four holonomies. The double-trace part is totally symmetric in the four legs. It multiplies one kinematic function, which is the sum over the three inequivalent orderings of the legs around the loop.

\subsection{Piece (a): the tree kernel}

Every diagram is a product of vertex coefficients, which are integrals of three or four harmonics over $S^3$, times energy denominators, summed over the magnetic quantum numbers and species of all internal lines. For fixed external letters the internal sums are finite, because angular momentum is conserved at each vertex. The result is a polynomial in $y$ for each set of external letters, and the thermal average is a power series in $y$ with rational coefficients. We find
\begin{equation}
q_4^{(a)}(y)=-450\,y^6-2184\,y^7-6246\,y^8-13632\,y^9-\tfrac{231931}{9}\,y^{10}-\tfrac{132776}{3}\,y^{11}-\ldots
\label{eq:q4a-series}
\end{equation}
The series is known to $y^{16}$ and is given in appendix~\ref{app:series}. The leading term comes from three scalars. All coefficients are negative. The coefficients grow like a power of $n$, as they must for a function whose nearest singularity is at $y=1$, and at $y_H=0.268$ the series converges geometrically. The sum is
\begin{equation}
q_4^{(a)}(x_H)=-0.7340(1).
\label{eq:q4a}
\end{equation}

\subsection{Piece (b): the loop kernel and its divergences}

Here one angular momentum runs freely around the loop, and each diagram diverges logarithmically. We cut the sum off at a maximum angular momentum $J$ on one line of the loop. The divergences must cancel in the sum of diagrams, since there is no counterterm in the theory of the form $(\tr\,\Phi^2)^2$. They do. For every pair of external multiplets the coefficients of $\ln J$ contributed by scalars, gauge bosons and fermions in the loop stand in the ratio
\begin{equation}
6:2:-8 .
\label{eq:628}
\end{equation}
These are the numbers of bosonic and fermionic degrees of freedom of the theory, and their sum is zero.

A finite remainder is more delicate, and the pure Yang--Mills calculation of AMMPV shows why. When logarithmic divergences cancel only between diagrams, a sharp cut-off on modes can leave a finite, regulator-dependent term. AMMPV found that a finite double-trace counterterm was needed to restore gauge invariance, a five per cent effect in their answer. We determined this counterterm for a general matter content. The double-trace part of \eqref{eq:colour} is totally symmetric, so it can be computed in a constant commuting background, where the one-loop energy is $\half(n_B-n_F)\sum\omega$. A sharp cut-off leaves a term from the surface of the cut-off sphere. It is proportional to
\begin{equation}
n_B-n_F ,
\label{eq:nBnF}
\end{equation}
the number of bosonic minus fermionic modes at each momentum. For pure glue $n_B-n_F=2$ and we recover AMMPV's counterterm exactly. For $\cN=4$ it is $8-8=0$. No counterterm is needed, provided the cut-off is applied to the same line in the same way in every diagram, time ordering and species. We return to this in section~\ref{sec:cutoff}.

With this understood the loop kernel is again a power series in $y$ with rational coefficients, each obtained by extrapolating $J\to\infty$,
\begin{equation}
\begin{aligned}
q_4^{(b)}(y)={}&-42\,y^4-384\,y^5-\tfrac{2181}{2}\,y^6-2100\,y^7\\
&-\tfrac{118841}{36}\,y^8-\tfrac{128452}{27}\,y^9-\tfrac{306863}{48}\,y^{10}-\ldots
\end{aligned}
\label{eq:q4b-series}
\end{equation}
It starts two orders earlier than piece (a) because two scalars are enough. The direct sum with an estimate of the tail, and five Pad\'e approximants, give values between $-1.4987$ and $-1.5044$. We quote
\begin{equation}
q_4^{(b)}(x_H)=-1.500(3).
\label{eq:q4b}
\end{equation}

\subsection{Independent evaluations}
\label{sec:independent}

An error in one vertex coefficient among several hundred would be invisible in the final number. The computation was therefore organised so that, with the one exception recorded at the end of this subsection, no number is used until two independent routes agree on it.

There are two constructions of the vertices. One follows the harmonic expansion of \cite{ITT} and expresses the integrals over $S^3$ through $9j$ symbols. The other realises the harmonics as homogeneous harmonic polynomials on $\mathbb{R}^4$ and integrates them directly. They share no code. They agree on every cubic, quartic and Coulomb coefficient, and they give the same one-loop series $c_3$, $d_1$ and $Q_4$ to $y^{14}$.

There are three contraction engines for piece (a) and two for piece (b). One works with reduced matrix elements and recoupling coefficients, one with explicit magnetic quantum numbers, and one at the level of words of matrices with mechanical normal ordering. The three agree on piece (a) coefficient by coefficient from $y^6$ to $y^{12}$, and topology by topology to $y^{10}$.

There is a separate route that does not use the effective Hamiltonian at all. It evaluates three-loop vacuum diagrams on $S^3\times S^1$ in the manner of AMMPV, with the time integrals done by divided differences and the contributions classified by the number of lines that wind. It shares the vertex coefficients with the first construction and nothing else. It reproduces piece (a) exactly to $y^{10}$, separately for each vertex content and each power of $\beta$, and it reproduces the first two coefficients of piece (b).

Finally there is an analytic test. If the thermal lines are restricted to the two scalars of the $su(2)$ sector, with fugacities $x_Z$ and $x_\phi$ and unrestricted virtual states, the known two-loop operator of \cite{BKS} predicts both pieces. Piece (a) must be the polynomial $-(x_Z-x_\phi)^2(x_Z^2+x_\phi^2-3x_Z-3x_\phi)$, and on the critical line $x_Z+x_\phi=1$ piece (b) must be $-2(x_Z^2+x_Zx_\phi+x_\phi^2)$. Three engines reproduce the first exactly. For the coefficient $-2$ in the second we find $-2.0000000$ with one engine and $-2.00(1)$ with the other.

One caveat must be recorded. The coefficients of piece (b) from $y^7$ onwards were produced by one engine only. The second engine confirms $y^4$, $y^5$ and the two largest entries of $y^6$. The unconfirmed terms contribute $-0.34$ of the $-1.50$. The margin for the sign of $b$ is fourteen in these units.