Bootstrap for Finite N Lattice Yang-Mills Theory

Vladimir Kazakov, Zechuan Zheng

► To cite this version:

Vladimir Kazakov, Zechuan Zheng. Bootstrap for Finite N Lattice Yang-Mills Theory. Journal of High Energy Physics, 2025, 2025 (03), pp.099. 10.1007/JHEP03(2025)099. hal-04570608

HAL Id: hal-04570608

https://hal.science/hal-04570608v1

Submitted on 28 Apr 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Bootstrap for finite N lattice Yang-Mills theory

Vladimir Kazakov ID a and Zechuan Zheng ID b,a

E-mail: Vladimir.Kazakov@ens.fr, zechuan.zheng.phy@gmail.com

**ABSTRACT:** We introduce a comprehensive framework for analyzing finite N lattice Yang-Mills theory and finite N matrix models. Utilizing this framework, we examine the bootstrap approach to SU(2) Lattice Yang-Mills Theory in 2,3 and 4 dimensions. The SU(2) Makeenko-Migdal loop equations on the lattice are linear and closed exclusively on single-trace Wilson loops. This inherent linearity significantly improves the efficiency of the bootstrap approach by leveraging the problem's convexity, permitting the inclusion of Wilson loops up to length 24. The exact upper and lower margins for the free energy per plaquette, derived from our bootstrap method, demonstrate good agreement with Monte Carlo data, achieving precision within 0.1% for the physically relevant range of couplings in both three and four dimensions. Additionally, our bootstrap data provides estimates of the string tension, in qualitative agreement with existing Monte Carlo computations.

**KEYWORDS:** Algorithms and Theoretical Developments, Nonperturbative Effects, Wilson, 't Hooft and Polyakov loops, Field Theories in Higher Dimensions

ARXIV ePRINT: 2404.16925

<footnote><sup>a</sup>Laboratoire de Physique de l'École normale supérieure,
ENS, Université PSL, CNRS, Sorbonne Université,
Université de Paris, 24 rue Lhomond, 75005 Paris, France</footnote>

<footnote><sup>b</sup>Perimeter Institute for Theoretical Physics,
Waterloo, ON N2L 2Y5, Canada</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>1.1</td><td>Main results</td><td>2</td></tr><tr><td>2</td><td>Loop equations</td><td>5</td></tr><tr><td>2.1</td><td>Finite N trace identities</td><td>6</td></tr><tr><td>2.2</td><td>Loop equations for the one-matrix model</td><td>7</td></tr><tr><td>2.3</td><td>Multi-trace Makeenko-Migdal loop equations</td><td>8</td></tr><tr><td>2.4</td><td>Makeenko-Migdal loop equations in the continuum</td><td>13</td></tr><tr><td>3</td><td>Positivity</td><td>15</td></tr><tr><td>3.1</td><td>Hermitian conjugation</td><td>16</td></tr><tr><td>3.2</td><td>Reflection Positivity</td><td>17</td></tr><tr><td>3.3</td><td>Symmetry Reduction</td><td>19</td></tr><tr><td>4</td><td>Bootstrap procedure</td><td>21</td></tr><tr><td>4.1</td><td>Initiate the bootstrap</td><td>22</td></tr><tr><td>4.2</td><td>Hierachy</td><td>23</td></tr><tr><td>4.3</td><td>Higher dimensional bound and string tension</td><td>25</td></tr><tr><td>5</td><td>Discussion and Prospects</td><td>28</td></tr><tr><td>A</td><td>Exact results for the unitary one-matrix model</td><td>31</td></tr><tr><td>B</td><td>Loop equations for SU(N)/U(N) one-matrix model</td><td>32</td></tr><tr><td>C</td><td>More bootstrap of the one-matrix model</td><td>33</td></tr></table>

## 1 Introduction

A great deal of our knowledge of the non-perturbative properties of the gauge theories in physical three and four space-time dimensions is based on the Lattice Monte Carlo calculations of the lattice Yang-Mills theories [1, 2]. The current state of the art in this area, thanks to the intensive use of supercomputers and modern algorithms, allows to achieve remarkable precision for the computation of spectra of glueballs and mesons (when the quark fields are included), with the size of lattices now attaining hundreds of lattice units [3]. On the other hand, relying solely on this method is not entirely satisfactory from both theoretical and practical perspectives. The search for other universal non-perturbative approaches did not lose its relevance.

Recent advancements in the field [4–8] have given rise to an innovative bootstrap method for the analysis of matrix models and lattice theories. This approach is unique in that it not only adheres to standard assumptions, like unitarity and global symmetries, but also
incorporates constraints among physical observables that are dictated by the equations of motion. This methodology promises a more rigorous framework for the exploration of matrix model dynamics. The broad utility of this novel bootstrap approach has been rapidly demonstrated across various domains, including lattice field theories [4, 8–10], matrix models [5, 6, 11–15], quantum systems [16–51], and even classical dynamical systems [52–55]. This extensive range of applications highlights the versatility and robustness of the method, underlining its significance in advancing theoretical frameworks across these fields.

In our previous work [8], we employed the bootstrap approach to study the SU(N) lattice Yang-Mills theories at N = ∞ in various dimensions, using the lattice Makeenko-Migdal equation (MME) for the Wilson loop [56]. Initial developments in this area are due to [9]. In this planar limit, the loop equations are inherently non-linear, leading in the context of bootstrap to a non-convex optimization problem. To address this complexity, we implemented a relaxation procedure in [7].

In this paper, we generalize the bootstrap procedure for the matrix models and lattice gauge theories in the planar limit to the finite N situation. First, we notice that, while the multi-trace averages are important to formulate a closed system of MMEs, we can explore a general property of the matrix-valued theories that the so-called trace identities help to reduce the number of traces in the multi-trace products of matrix fields. That is, the product of k > N traces (or k > N − 1, depending on the symmetries) of arbitrary products of matrices can be represented as the sum of products of smaller numbers of traces. Higher trace variables are described through trace identities specific to the corresponding matrix ensemble.

In this paper, we concentrate on systems with variables SU(2) and SU(3), where the loop equations appear to be linear in the space of single-trace and double-trace averages, respectively. For the SU(3) case, we discuss the general framework of linear double-trace loop equations, including the lattice gauge theory case, and solve via bootstrap the one-matrix model for illustrative purposes. The rest of our paper is concentrated on the study and bootstrap solution of the loop equation for SU(2) lattice Yang-Mills theory. We extensively use the fact that, in any dimension, this equation closes exclusively on the single trace Wilson loops and is fully linear.<fnref n="1" /> This linearity obviates the need for the relaxation procedures typically required in the large N limit. With the aid of the hierarchy of Wilson loops discussed in this article, along with the specialized program we developed, we have maximized the capabilities of the bootstrap procedure. The linearity of MMEs helped us to exploit the full power of the utilized hardware and of the advanced optimization software, allowing us to bootstrap the Wilson loops up to a maximum length of 24 lattice units.

We believe that the linear SU(2) loop equations in the single-trace loop space which we
derive in detail in this paper, both on the lattice and in the continuum, can provide important
insights into the nature and properties of the Yang-Mills theory. We will mostly concentrate
here on the practical aspects and advantages of such linearity for the bootstrap procedure.

1.1 Main results

This paper proposes a general framework for bootstrapping finite N lattice Yang-Mills theory and finite N matrix models, highlighting the crucial role of multi-trace variable truncation.

<footnote n="1">¹As noticed already in [57].</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The results of the bootstrap calculations in three and four dimensions for the plaquette average $u_P = \text{tr}U_\square$ as a function of the coupling $\lambda = \frac{2N^2}{\beta}$ are displayed graphically. Blue lines represent interpolations of the upper and lower margins for the maximum length of loops $L_{\text{max}} \le 16$ used in the bootstrap procedure for loop equations. Red lines and dots illustrate the same for $L_{\text{max}} \le 24$. The gray area indicates the allowed values of $u_P$ for $L_{\text{max}} \le 8$. For additional details regarding these plots, refer to section 4.3. The dashed lines in the figure represent the strong and weak coupling expansions, extracted from [58]. Additionally, the results from Monte Carlo simulations are depicted using purple dots [3, 59] and yellow dots [58], each set representing data from different sources. For additional details regarding these plots, refer to section 4.3.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2. Illustration of the different terms in the loop equation (1.4).</figcaption></figure>

Specifically, the number of independent multi-trace variables is bounded, with higher-number trace variables expressed as linear combinations of a subset of multi-trace variables. These truncations, allowing reduction of any multi-trace SU(N) variables to a linear combination of only N − 1, N − 2, . . . , 1-trace variables, are demonstrated here through the following trace identities for SU(2) and SU(3):

$$
\mathrm{tr} U \mathrm{tr} V = \mathrm{tr} UV + \mathrm{tr} UV^\dagger, \quad \forall U, V \in \mathrm{SU}(2) \qquad (1.1)
$$

$$
\begin{equation}
\begin{split}
\mathrm{tr} U \mathrm{tr} V \mathrm{tr} W = -\mathrm{tr} W V U - \mathrm{tr} U V W &+ \mathrm{tr} U \mathrm{tr} V W + \mathrm{tr} U V \mathrm{tr} W + \mathrm{tr} U W \mathrm{tr} V \\
& + \mathrm{tr} U^\dagger V \mathrm{tr} U^\dagger W - \mathrm{tr} U^\dagger V U^\dagger W, \quad \forall U, V, W \in \mathrm{SU}(3)
\end{split}
\tag{1.2}
\end{equation}
$$

We then explored the dynamics of multi-trace Wilson loops, captured by multi-trace loop equations, both on the lattice and in the continuum. For SU(2), the trace identities simplify the representation of every double-trace Wilson loop to a sum of single-trace Wilson loops defined as an averaged trace of the product of group variables $U_{l_ν}$ on the links $l_ν$ along the lattice contour $C$:<fnref n="2" />

$$
W(C) = \langle \mathrm{Tr} \prod_{l_\nu \in C} U_{l_\nu} \rangle. \qquad (1.3)
$$

The specific form of the SU(2) Makeenko-Migdal loop equations for the Wilson lattice
action (2.15) is compactly represented as:

$$
\frac{1}{\lambda} \sum_{\mu \neq \pm \nu} [W(\sigma_{\nu\mu}^l \circ C_{l\nu}) - W(\sigma_{\mu\nu}^l \circ C_{l\nu})] + \frac{1}{2} \sum_{l'_\nu \sim l_\nu \in C} (-1)^{\frac{1+\nu\cdot\nu'}{2}} \left( W(C_{x,\nu}^{(1)} \circ \overline{C_{y,\nu'}^{(2)}}) + \frac{1}{2}W(C) \right) = 0 \quad (1.4)
$$

where $C \circ C'$ denotes the composition of two contours (joint at a given common link), $\sigma_{\nu\mu}^l$ is
a plaquette which includes the link $l$ and has the orientation $\nu, \mu$, where any Greek index

<footnote n="2">²Here and in the whole article, we use the notation that Tr = 1/N tr, so the normalization of the identity equation has Tr[1] = 1.</footnote>
takes the values ±1, ±2 ..., ±d (sign is related to the orientation of the link w.r.t. gauge variable), C means the contour C with the flipped direction (conjugation of the corresponding Wilson path). The notation $\sum_{l'_\nu \sim l_\nu}$ denotes the summation over links in the contour C that coincide with $l_\nu$ on the lattice. Importantly, this summation explicitly includes $l_\nu$ itself. The variations and specific instances of this summing process are visually depicted in figure 2, while a more detailed and explicit mathematical formulation can be found in section 2.

We also derive the finite N Makeenko-Migdal equation in the continuum, focusing here on the SU(2) case, where the continuous loop equation is [57]:<fnref n="3" />

$$ \frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)} W(C_x) = \frac{\lambda}{2} \oint_C dy_{\nu} \delta^{(d)}(x-y) \left[ \frac{1}{2} W(C_x) + W \left( C_{xy}^{(1)} \circ \overline{C_{yx}^{(2)}} \right) \right] \quad (1.5) $$

This situation is simpler than that of the previously studied 't Hooft limit, where double-trace Wilson loops are products of single-trace Wilson loops due to large N factorization. There is no need for relaxation in the bootstrap scheme, which is typically required due to the non-convex nature of the large N loop equations. The inherent linearity of the SU(2) loop equations renders the entire bootstrap procedure convex from the very beginning.

The paper concludes with a detailed examination of SU(2) lattice Yang-Mills theory across two, three, and four dimensions. A major technical achievement of this work is the significant extension of the maximum length of Wilson loops involved in our bootstrap procedure — up to 24 lattice units for the 2D, 3D and 4D cases. The hierarchical system, which categorizes Wilson lines by level and selects operators based on this categorization, has markedly increased the precision of computations for Wilson averages of small loops. Notably, the margins for the plaquette average in three and four dimensions within the physically relevant coupling range achieve a precision of the order of 0.1%. These findings are depicted in figure 1, showcasing results for three and four dimensions, respectively. Additionally, using interpolation of the Creutz ratio, we extracted string tension data. Although these results do not match the precision of Monte Carlo simulations due to the limited loop size (maximum 2 × 3), they are qualitatively sound.

**Note added.** While we were finalizing this article, we became aware of the work [60], which focuses on the Abelian gauge theory and parallels some of our discussions on the one-matrix model.

## 2 Loop equations

Our bootstrap approach is aimed at the expectation values of multi-trace quantities within the given matrix ensemble. Our method differs from the schemes involving the 't Hooft limit, where, due to the large N factorization of multi-trace averages, the loop equations close on the single-trace space<fnref n="4" /> [4, 8].

Although at finite N we cannot factorize the multi-trace averages into products of averages of single traces, we can make use of trace identities to reduce such averages to a

<footnote n="3">³Derivation of such equations on the lattice and in continuum, and in particular of the SU(3) loop equation, is deferred to section 2.</footnote>

<footnote n="4">⁴For more complex physical quantities, such as glueball correlation functions, multi-trace averages are necessary.</footnote>
limited number of traces inside the averages. Several of the trace identities discussed in this section were originally derived in [61, 62]. These identities effectively limit the space of multi-trace averages to those containing no more than N traces, with the specific reduction contingent upon the random matrix ensemble being studied.<fnref n="5" />

The trace identities strongly influence the kinematics of loop variables at finite N and their structure depends on the specific value of N in a non-trivial way. As the first step of our approach, we derive the multi-trace loop equations, which concern the dynamics of the loop variables. We will focus on two models: initially, the SU(3) one-matrix model, which will serve as an illustrative example, followed by the SU(2) lattice Yang-Mills theory employing the Wilson action.

## 2.1 Finite N trace identities

The trace identities we consider here are quite general. We begin with a discussion of the case of general matrices in GL(N). An useful observation follows:

$$ \epsilon^{a_1 a_2 \dots a_{N+1}} \epsilon^{b_1 b_2 \dots b_{N+1}} M_1^{a_1 b_1} \dots M_{N+1}^{a_{N+1} b_{N+1}} = 0, \quad \forall M_1, \dots, M_{N+1} \in \text{GL}(N) \quad (2.1) $$

Additionally, we can express the product of the $\epsilon$ tensors in terms of the Kronecker delta symbols [65]:

$$ \epsilon^{a_1 a_2 \dots a_{N+1}} \epsilon^{b_1 b_2 \dots b_{N+1}} = \delta_{a_1 b_1} \delta_{a_2 b_2} \dots \delta_{a_{N+1} b_{N+1}} + \sum_{P \neq \text{id}} (-1)^P \delta_{a_1 P(b_1)} \delta_{a_2 P(b_2)} \dots \delta_{a_{N+1} P(b_{N+1})} \quad (2.2) $$

The sum runs over all permutations of the b indices, excluding the identity permutation. Combining these formulas, the first term in Equation 2.2 leads to a N + 1 trace variable, while the remainder includes terms with N traces or fewer:

$$ \text{tr} M_1 \text{tr} M_2 \dots \text{tr} M_{N+1} = - \sum_{P \neq \text{id}} (-1)^P \delta_{a_1 P(b_1)} \delta_{a_2 P(b_2)} \dots \delta_{a_{N+1} P(b_{N+1})} M_1^{a_1 b_1} \dots M_{N+1}^{a_{N+1} b_{N+1}} \quad (2.3) $$

As an example, for GL(2), we get:

$$ \text{tr} M_1 \text{tr} M_2 \text{tr} M_3 = - \text{tr} M_3 M_2 M_1 - \text{tr} M_1 M_2 M_3 + \text{tr} M_1 \text{tr} M_2 M_3 + \text{tr} M_1 M_2 \text{tr} M_3 + \text{tr} M_1 M_3 \text{tr} M_2 \quad (2.4) $$

This indicates that in the case of GL(2), the trace variables generally truncate at double traces.

More interestingly, for the subgroups of GL(N), there is evidence suggesting that the truncation of multi-trace variables can be further refined. For SU(N) and SO(N), the unimodularity condition further applies:

$$ 1 = \text{Det} U = \text{Det} U^\dagger = \frac{1}{N!} \epsilon^{a_1 a_2 \dots a_N} \epsilon^{b_1 b_2 \dots b_N} U^{\dagger a_1 b_1} \dots U^{\dagger a_N b_N} \quad (2.5) $$

We multiply both sides of the previous equation by the following term

$$ \frac{1}{N!} \epsilon^{c_1 c_2 \dots c_N} \epsilon^{d_1 d_2 \dots d_N} U^{c_1 d_1} V_1^{c_2 d_2} \dots V_{N-1}^{c_N d_N} \quad (2.6) $$

<footnote n="5">⁵It will be beneficial to further clarify the relationship between the trace identities presented in our study and those associated with finite N discussed in [63, 64].</footnote>
and finally, we obtain:

$$ \epsilon^{c_1 c_2 \dots c_N} \epsilon^{d_1 d_2 \dots d_N} U^{c_1 d_1} V_1^{c_2 d_2} \dots V_{N-1}^{c_N d_N} \\ = \frac{1}{N!} (\epsilon^{a_1 a_2 \dots a_N} \epsilon^{c_1 c_2 \dots c_N}) (\epsilon^{b_1 b_2 \dots b_N} \epsilon^{d_1 d_2 \dots d_N}) U^{c_1 d_1} V_1^{c_2 d_2} \dots V_{N-1}^{c_N d_N} U^{\dagger a_1 b_1} \dots U^{\dagger a_N b_N} \quad (2.7) $$

For similar reasons as outlined in (2.2), the left-hand side of (2.7) contains only one N-trace variable. Upon contracting the $\epsilon$ tensors, $\epsilon^{a_1 a_2 \dots a_N} \epsilon^{c_1 c_2 \dots c_N}$ and $\epsilon^{b_1 b_2 \dots b_N} \epsilon^{d_1 d_2 \dots d_N}$, the right-hand side contains at most $(N-1)$-trace variables. This reduction occurs because $U$ and $U^\dagger$ always contract to yield a trivial trace, tr1. As examples, for SU(2) and SU(3), we obtain the identities (1.1) and (1.2) which will be extensively used in this paper.

Therefore, unlike the general case for GL(N) where the trace variables truncate at N traces, for SU(N) and SO(N), the trace variables truncate at $N-1$ traces.

## 2.2 Loop equations for the one-matrix model

For illustration purposes, let us first consider the “single plaquette” toy model for the SU(3) matrix ensemble:<fnref n="6" />

$$ Z = \int \mathcal{D}U e^{\frac{N}{\lambda}(\text{tr}U^{\dagger} + \text{tr}U)} \quad (2.8) $$

This model represents the simplest non-trivial unitary matrix model, characterized by a real action. Despite its apparent simplicity, the model holds practical physical significance as it describes the dynamics of a single plaquette Wilson loop in a two-dimensional lattice Yang-Mills theory [66, 67].

In light of the trace relation (1.2) derived in the last section, the space of loop variables truncates at double-trace, and the triple-trace operators are a linear sum of the double-trace operators:

$$ \begin{aligned} \mathrm{Tr} U^i \mathrm{Tr} U^j \mathrm{Tr} U^k = & -\frac{2}{9} \mathrm{Tr} U^{i+j+k} + \frac{1}{3} \mathrm{Tr} U^i \mathrm{Tr} U^{j+k} + \frac{1}{3} \mathrm{Tr} U^{i+j} \mathrm{Tr} U^k + \frac{1}{3} \mathrm{Tr} U^{i+k} \mathrm{Tr} U^j \\ & + \frac{1}{3} \mathrm{Tr} U^{j-i} \mathrm{Tr} U^{k-i} - \frac{1}{9} \mathrm{Tr} U^{j+k-2i}, \quad \forall U \in \mathrm{SU}(3) \end{aligned} \quad (2.9) $$

Here, the derivation of the double-trace loop equation stems from the condition that the shift in the probability measure must vanish ($n > 0$):

$$ \int \mathcal{D}U \, \delta_{\epsilon} \left( U_{ab}^{n} \mathrm{Tr} U^{m} e^{\frac{N}{\lambda} (\mathrm{tr} U^{\dagger} + \mathrm{tr} U)} \right) = 0 \quad (2.10) $$

with the definition of the variation:

$$ \delta_{\epsilon}(U_{ab}) = U_{ab} + i\epsilon_{ac}U_{cb}, \quad \delta_{\epsilon}(U_{ab}^{\dagger}) = U_{ab}^{\dagger} - iU_{ac}^{\dagger}\epsilon_{cb} \quad (2.11) $$

Here $\epsilon$ is an arbitrary traceless Hermitian 3 × 3 matrix. Expanding (2.10), we get:

$$ \epsilon_{cd} \left\langle \frac{m}{N} U_{ab}^{n} U_{dc}^{m} + \sum_{i=0}^{n-1} U_{ac}^{i} U_{db}^{n-i} \mathrm{Tr} U^{m} + \frac{N}{\lambda} U_{ab}^{n} \mathrm{Tr} U^{m} (U_{dc} - U_{dc}^{\dagger}) \right\rangle = 0 \quad (2.12) $$

<footnote n="6">⁶We gather a concise summary of the exact solution of this model for SU(N)/U(N) in appendix A, and provide a general treatment of the SU(N)/U(N) multi-trace loop equations in appendix B.</footnote>
We observe that if $\epsilon_{cd}X_{dc} = 0$ for any traceless $\epsilon$, then it follows that $X_{dc} = \frac{1}{N}\delta_{dc}X_{ff}$, from which we can derive:

$$ \begin{aligned} & \left\langle \frac{m}{N} U_{ab}^n U_{dc}^m + \sum_{i=0}^{n-1} U_{ac}^i U_{db}^{n-i} \mathrm{Tr} U^m + \frac{N}{\lambda} U_{ab}^n \mathrm{Tr} U^m (U_{dc} - U_{dc}^\dagger) \right\rangle \\ & = \delta_{dc} \left\langle \frac{m}{N} U_{ab}^n \mathrm{Tr} U^m + \frac{n}{N} U_{ab}^n \mathrm{Tr} U^m + \frac{N}{\lambda} U_{ab}^n \mathrm{Tr} U^m (\mathrm{Tr} U - \mathrm{Tr} U^\dagger) \right\rangle \end{aligned} \quad (2.13) $$

Upon contracting both sides with $\frac{1}{N^2}\delta_{ac}\delta_{bd}$, we arrive at the final form of the loop equation:

$$ \begin{aligned} & \left\langle \frac{m}{N^2} \mathrm{Tr} U^{n+m} + \sum_{i=0}^{n-1} \mathrm{Tr} U^i \mathrm{Tr} U^{n-i} \mathrm{Tr} U^m + \frac{1}{\lambda} \mathrm{Tr} U^m (\mathrm{Tr} U^{n+1} - \mathrm{Tr} U^{n-1}) \right\rangle \\ & = \left\langle \frac{m+n}{N^2} \mathrm{Tr} U^n \mathrm{Tr} U^m + \frac{1}{\lambda} \mathrm{Tr} U^n \mathrm{Tr} U^m (\mathrm{Tr} U - \mathrm{Tr} U^\dagger) \right\rangle \end{aligned} \quad (2.14) $$

The current double-trace loop equation is generally applicable to the SU(N) matrix model. For the U(N) matrix model, however, the left-hand side of the equation vanishes. By incorporating the trace identity (2.9), which is specific to SU(3), we get the loop equation that truncates at the double-trace variable level.

## 2.3 Multi-trace Makeenko-Migdal loop equations

In this section, we derive the general multi-trace Makeenko-Migdal loop equations for SU(N) and U(N) lattice gauge theories employing the Wilson action. The outcomes of this analysis, when integrated with the trace identities previously discussed, demonstrate that the loop equations for SU(N) lattice gauge theory close among $(N-1)$-trace variables. Conversely, for U(N) lattice gauge theory, closure occurs among N-trace variables.

The partition function for SU(N) and U(N) lattice Yang-Mills theories across different spacetime dimensions is given by:

$$ Z = \int \mathcal{D}U \exp(-S), \quad \text{where } S = -\frac{N}{\lambda} \sum_p \mathrm{tr}U_p \quad (2.15) $$

Here, $U_p$ represents the product of four group variables along four links that form a plaquette, with the summation extending over all orientations of all plaquettes. In two spacetime dimensions, our normalization ensures that the dynamics of the single plaquette Wilson loop, as specified by (2.15), are equivalent to those in the one-matrix model (2.8).<fnref n="7" /> Furthermore, our normalization relates to $\beta$, as defined in [3], by:

$$ \beta = \frac{2N^2}{\lambda} \quad (2.16) $$

The multi-trace loop equation corresponds to the vanishing of the following variation:

$$ 0 = \int \mathcal{D}U \delta_{\epsilon} \left( \mathcal{W}[C(x, \mu, \nu)]_{ab} \prod_{i=1}^{m} (\mathrm{Tr} \mathcal{W}[C_i]) \exp(-S) \right). \quad (2.17) $$

<footnote n="7">⁷Which does not apply of course to more complex loops in 2d, see [68–70].</footnote>
Here for the goal of deriving the loop equation, we are considering the following Wilson path:

$$W[C(x, \mu, \nu)]_{ab} = (U_{x,\mu} \dots U_{y_{-},\eta_{-}} U_{y,\eta} U_{y_{+},\eta_{+}} \dots U_{x,\nu}^{\dagger})_{ab} \quad (2.18)$$

That is a specific Wilson path with open indices, starting at $x$ with $\mu$ direction and ending with $-\nu$ direction, at the same point $x$. Here $\delta_{\epsilon}$ is defined by taking the following variation on $U_{x,\mu}$,<fnref n="8" /> whereas all the other link variables are unchanged:

$$U_{x,\mu} \rightarrow (\mathbb{I} + i\epsilon)U_{x,\mu}, \quad U_{x,\mu}^{\dagger} \rightarrow U_{x,\mu}^{\dagger}(\mathbb{I} - i\epsilon), \quad \epsilon^{\dagger} = \epsilon \quad (2.19)$$

If we are talking about the SU(N) group instead of U(N), we need also to impose the tracelessness condition $\text{tr}\epsilon = 0$, which guarantees that the group element remains unimodular. Under this variation, we notice that the action is changed by:<fnref n="9" />

$$\delta_{\epsilon}(S) = -i \frac{N}{\lambda} \epsilon_{cd} \sum_{\nu \neq \pm \mu} (U_p(x, \mu, \nu) - U_p(x, \mu, \nu)^{\dagger})_{dc} \quad (2.21)$$

Applying these variations to (2.17) we can get the loop equation in the following form:

$$\epsilon_{cd} \langle A_{dc,ab} \rangle = 0 \quad (2.22)$$

where $A_{dc,ab}$ reads:

$$\prod_{i=1}^{m} (\text{Tr}W[C_i]) \left( \delta_{ac} W[C(x,\mu,\nu)]_{db} + \frac{N}{\lambda} \sum_{\eta \neq \pm \mu} (U_p(x,\mu,\eta) - U_p(x,\mu,\eta)^{\dagger})_{dc} W[C(x,\mu,\nu)]_{ab} + \sum_{U_{y,\eta}=U_{x,\mu}} W[C_{y1}(x,\mu,-\mu_{-})]_{ac} W[C_{y2}(x,\mu,\nu)]_{db} - \sum_{U_{y,\eta}=U_{x,\mu}^{\dagger}} W[C_{y1}(x,\mu,\mu)]_{ac} W[C_{y2}(x,\mu_{+},\nu)]_{db} \right) + \frac{1}{N} \left( \sum_{U_{y,\eta}=U_{x,\mu}} W[C_n(x,\mu,-\eta_{-})]_{dc} - \sum_{U_{y,\eta}=U_{x,\mu}^{\dagger}} W[C_n(x,\mu_{+},\mu)]_{dc} \right) W[C(x,\mu,\nu)]_{ab} \prod_{i=1, i \neq n}^{m} (\text{Tr}W[C_i]) \quad (2.23)$$

We further notice the fact that due to the arbitrariness of $\epsilon$, we have:

$$\langle A_{dc,ab} \rangle = s_G \frac{\delta_{dc}}{N} \langle A_{ee,ab} \rangle \quad (2.24)$$

Here $s_G$ is given by

$$s_G = \begin{cases} 1, & G = \text{SU(N)} \\ 0, & G = \text{U(N)} \end{cases} \quad (2.25)$$

which takes into account that $\epsilon$ must be traceless for SU(N) gauge group.

$$\begin{align} U_p(x, \mu, \nu)_{dc} &= U_{x,\mu} U_{x+\mu,\nu} U_{x+\mu+\nu, -\mu} U_{x+\nu, -\nu} \tag{2.20} \\ U_p(x, \mu, \nu)_{dc}^{\dagger} &= U_{x,\nu} U_{x+\nu,\mu} U_{x+\mu, -\nu}^{\dagger} U_{x,\mu}^{\dagger} = U_p(x, \nu, \mu)_{dc}. \nonumber \end{align}$$

<footnote n="8">⁸We notice that in this notation, $U_{x,\mu}^{\dagger} = U_{x+\mu, -\mu}$.</footnote>

<footnote n="9">⁹To help set up our notation, we further explain that here:</footnote>
We further contract both sides of (2.24) by $\frac{1}{N^2}\delta_{ac}\delta_{bd}$, and finally get the loop equation in the form:

$$A_{\text{id}} + A_{\text{var}} + A_{\text{split}} + A_{\text{join}} = 0 \qquad (2.26)$$

where we denoted:

$$A_{\text{id}} = \left( 1 - s_G \frac{1 + n_+ - n_- + \sum_{i}^{m} (n_{i+} - n_{i-})}{N^2} \right) \langle \text{Tr} \mathcal{W}[C(x, \mu, \eta)] \prod_{i=1}^{m} \text{Tr} \mathcal{W}[C_i] \rangle, \qquad (2.27)$$

$$A_{\text{var}} = \frac{1}{\lambda} \left\langle \prod_{i=1}^{m} \text{Tr} \mathcal{W}[C_i] \sum_{\eta \neq \pm \mu} \left( \text{Tr} \mathcal{W}[U_p(x, \mu, \eta)C] - \text{Tr} \mathcal{W}[U_p(x, \mu, \eta)^{\dagger}C] \right) \right\rangle \\ - \frac{s_G}{\lambda} \left\langle \prod_{i=1}^{m} \text{Tr} \mathcal{W}[C_i] \sum_{\eta \neq \pm \mu} \left( \text{Tr} \mathcal{W}[U_p(x, \mu, \eta)] \text{Tr} \mathcal{W}[C] - \text{Tr} \mathcal{W}[U_p(x, \mu, \eta)^{\dagger}] \text{Tr} \mathcal{W}[C] \right) \right\rangle \qquad (2.28)$$

$$A_{\text{join}} = \frac{1}{N^2} \left\langle \left( \sum_{U_{y,\eta}=U_{x,\mu}} \text{Tr} \mathcal{W}[C_n(x,\mu, -\eta_-) \circ C] - \sum_{U_{y,\eta}=U_{x,\mu}^{\dagger}} \text{Tr} \mathcal{W}[C_n(x,\eta_+, \mu) \circ C] \right) \prod_{i=1, i \neq n}^{m} \text{Tr} \mathcal{W}[C_i] \right\rangle, \qquad (2.29)$$

$$A_{\text{split}} = \left\langle \prod_{i=1}^{m} \text{Tr} \mathcal{W}[C_i] \left( \sum_{U_{y,\eta}=U_{x,\mu}} \text{Tr} \mathcal{W}[C_{y1}] \text{Tr} \mathcal{W}[C_{y2}] - \sum_{U_{y,\eta}=U_{x,\mu}^{\dagger}} \text{Tr} \mathcal{W}[C_{y1}] \text{Tr} \mathcal{W}[C_{y2}] \right) \right\rangle. \qquad (2.30)$$

The term $A_{\text{id}}$ is proportional to the multi-trace Wilson loop under consideration. Here the terms $n_+$ and $n_-$ denote the number of overlaps with the link $U_{x,\mu}$ in the contour $C$, depending on whether their direction aligns with $\mu$. Similarly, $n_{i,+}$ and $n_{i,-}$ represent the number of overlaps with the link $U_{x,\mu}$ in the contour $C_i$. $A_{\text{var}}$ is the only part that is dynamical, in the sense that it depends on the content of the Wilson action. We notice that in the second line of (2.28) vanishes for SU(2) due to the trace identities (1.1). $A_{\text{split}}$ and $A_{\text{join}}$ describe the splitting and joining of the Wilson loops. These terms are illustrated in the figure 3. We stress that the loop equation we derived so far is independent of the specific value of N. Practically we will set $m = N - 2$ for SU(N) group and $m = N - 1$ for U(N), substituting the finite N trace identities we discussed in section 2.1 to the term $A_{\text{split}}$ and $A_{\text{var}}$, finally getting the loop equation that is closed on the $(m+1)$-trace Wilson loops.

For the specific case of SU(2) lattice gauge theory, which is extensively discussed in this article, we observe that the loop equation derived above can be further simplified. This simplification is achieved by substituting the trace identity for SU(2), as shown in (1.1):

$$\frac{1}{\lambda} \sum_{\mu \neq \pm \nu} \left[ W(\sigma_{\nu\mu}^l \circ C_{l\nu}) - W(\sigma_{\mu\nu}^l \circ C_{l\nu}) \right] + \frac{1}{2} \sum_{l'_{\nu} \sim l_{\nu} \in C} (-1)^{\frac{1+\nu \cdot \nu'}{2}} \left( W(C_{x,\nu}^{(1)} \circ \overline{C_{y,\nu'}^{(2)}}) + \frac{1}{2} W(C) \right) = 0 \qquad (2.31)$$

Here, we specifically emphasize that the sum over $l'_{\nu} \sim l_{\nu} \in C$ includes $l_{\nu}$ itself.
<figure><img src="image_5.png" /><figcaption>Figure 3. Illustration of the different terms in the loop equation (2.26).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. The lowest order single-trace loop equations.</figcaption></figure>

For illustration purposes, we present here the lowest-order single-trace loop equations in a two-dimensional lattice, which correspond to the variation of the figure 4:

$$0 = \begin{array}{c} \Box \\ \Box \end{array} - 1 + \begin{array}{c} \Box \\ \Box \end{array} - \begin{array}{c} \Box \\ \Box \end{array} + (1 - \frac{s_G}{N^2})\lambda \begin{array}{c} \Box \\ \Box \end{array} \\ -s_G \begin{array}{c} \Box \\ \Box \end{array} - \begin{array}{c} \Box \\ \Box \end{array} + \begin{array}{c} \Box \\ \Box \end{array} - \begin{array}{c} \Box \\ \Box \end{array} \tag{2.32}$$

### 2.3.1 Backtracks and redundancies

Backtracks play a rather subtle role in the formulation of the loop equations. They are defined by a string of Wilson lines (or even a backtrack “tree” growing from a vertex of the contour) that are equivalent to group identity due to unitarity. Backtrack loop equations
<figure><img src="image_7.png" /><figcaption>Figure 5. Lowest order back-track loop equation.</figcaption></figure>

originated from a variation on a “cut” Wilson line, i.e. such that it has open indices inside it. Figure 5 shows the lowest order backtrack loop equation in a two-dimensional lattice:

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node (0) at (0,0) {\begin{tikzpicture}[baseline=(current bounding box.center)] \draw (0,0) -- (0,1); \draw (0,1) -- (1,1); \draw (1,1) -- (1,0); \draw (1,0) -- (0,0); \node at (0.5,0.5) {0}; \node at (0.5,1.5) {1}; \node at (1.5,0.5) {1}; \node at (1.5,1.5) {0}; \end{tikzpicture}}; \node at (0.25,0.5) {−}; \node at (0.25,1.5) {−}; \node at (0.75,0.5) {−}; \node at (0.75,1.5) {−}; \node at (1.25,0.5) {+}; \node at (1.25,1.5) {+}; \end{tikzpicture} = 0 \qquad (2.33) $$

More general situations contain a backtrack longer than the length one as shown in figure 5, the backtrack can also overlap with the non-backtrack part of the Wilson loop. In our study of SU(2) lattice Yang-Mills theory, approximately one-third to half of the loop equations are of the backtrack type. They play an essential role in constraining the space of (multi-trace) Wilson loops.

Another subtlety of the backtrack comes from the fact that it is equivalent to identity, so it can be inserted at any point of the Wilson loop. Combined with (1.1), we can expand the following double-trace Wilson loop to a combination of single-trace Wilson loops. Take as an example the double-trace Wilson loop corresponding to two loops:

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node (0) at (0,0) {\begin{tikzpicture}[baseline=(current bounding box.center)] \draw (0,0) -- (0,1); \draw (0,1) -- (1,1); \draw (1,1) -- (1,0); \draw (1,0) -- (0,0); \node at (0.5,0.5) {0}; \node at (0.5,1.5) {1}; \node at (1.5,0.5) {1}; \node at (1.5,1.5) {0}; \end{tikzpicture}}; \node at (0.25,0.5) {−}; \node at (0.25,1.5) {−}; \node at (0.75,0.5) {−}; \node at (0.75,1.5) {−}; \node at (1.25,0.5) {+}; \node at (1.25,1.5) {+}; \end{tikzpicture} \qquad (2.34) $$

There are many different ways to insert here the backtracks and reduce them to a combination of single traces, e.g.:

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node (0) at (0,0) {\begin{tikzpicture}[baseline=(current bounding box.center)] \draw (0,0) -- (0,1); \draw (0,1) -- (1,1); \draw (1,1) -- (1,0); \draw (1,0) -- (0,0); \node at (0.5,0.5) {0}; \node at (0.5,1.5) {1}; \node at (1.5,0.5) {1}; \node at (1.5,1.5) {0}; \end{tikzpicture}}; \node at (0.25,0.5) {2}; \node at (0.25,1.5) {2}; \node at (0.75,0.5) {2}; \node at (0.75,1.5) {2}; \node at (1.25,0.5) {+}; \node at (1.25,1.5) {+}; \end{tikzpicture} = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (0) at (0,0) {\begin{tikzpicture}[baseline=(current bounding box.center)] \draw (0,0) -- (0,1); \draw (0,1) -- (1,1); \draw (1,1) -- (1,0); \draw (1,0) -- (0,0); \node at (0.5,0.5) {0}; \node at (0.5,1.5) {1}; \node at (1.5,0.5) {1}; \node at (1.5,1.5) {0}; \end{tikzpicture}}; \node at (0.25,0.5) {2}; \node at (0.25,1.5) {2}; \node at (0.75,0.5) {2}; \node at (0.75,1.5) {2}; \node at (1.25,0.5) {+}; \node at (1.25,1.5) {+}; \end{tikzpicture} \qquad (2.35) $$

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node (0) at (0,0) {\begin{tikzpicture}[baseline=(current bounding box.center)] \draw (0,0) -- (0,1); \draw (0,1) -- (1,1); \draw (1,1) -- (1,0); \draw (1,0) -- (0,0); \node at (0.5,0.5) {0}; \node at (0.5,1.5) {1}; \node at (1.5,0.5) {1}; \node at (1.5,1.5) {0}; \end{tikzpicture}}; \node at (0.25,0.5) {2}; \node at (0.25,1.5) {2}; \node at (0.75,0.5) {2}; \node at (0.75,1.5) {2}; \node at (1.25,0.5) {+}; \node at (1.25,1.5) {+}; \end{tikzpicture} = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (0) at (0,0) {\begin{tikzpicture}[baseline=(current bounding box.center)] \draw (0,0) -- (0,1); \draw (0,1) -- (1,1); \draw (1,1) -- (1,0); \draw (1,0) -- (0,0); \node at (0.5,0.5) {0}; \node at (0.5,1.5) {1}; \node at (1.5,0.5) {1}; \node at (1.5,1.5) {0}; \end{tikzpicture}}; \node at (0.25,0.5) {2}; \node at (0.25,1.5) {2}; \node at (0.75,0.5) {2}; \node at (0.75,1.5) {2}; \node at (1.25,0.5) {+}; \node at (1.25,1.5) {+}; \end{tikzpicture} \qquad (2.36) $$

These insertions of the backtracks indicate that there are some kinematic redundancies in the space of Wilson loops. Similar kinematic redundancies also appear for double-trace Wilson loops for SU(3) lattice Yang-Mills theory. In addition to the redundancy arising from back-tracking, we observe inherent redundancies in the SU(3) trace identities (1.2). Given that the left-hand side of the equation remains invariant under permutations of U, V, and W, it is straightforward to apply permutations to the right-hand side. This process enables us to delineate the equations that identify redundancies among double-trace variables:

$$ \mathrm{tr} U^\dagger V \mathrm{tr} U^\dagger W - \mathrm{tr} U^\dagger V U^\dagger W = \mathrm{tr} V^\dagger U \mathrm{tr} V^\dagger W - \mathrm{tr} V^\dagger U V^\dagger W, \quad \forall U, V, W \in \mathrm{SU}(3) \qquad (2.37) $$
## 2.4 Makeenko-Migdal loop equations in the continuum

In this section, we study the Makeenko-Migdal loop equations [56] for the U(N) and SU(N) groups, within the context of Yang-Mills theory. The action for this theory is defined as:

$$S = -\frac{1}{4g^2} \int d^d x \operatorname{tr}(F_{\mu\nu} F^{\mu\nu}), \quad \text{where } F_{\mu\nu} = \partial_{\nu} A_{\mu} - \partial_{\mu} A_{\nu} + [A_{\mu}, A_{\nu}]. \qquad (2.38)$$

Our focus will be on multi-trace Wilson loops, which are key observables in this theoretical framework:

$$W(C_1, \dots, C_n) = \langle \operatorname{Tr} \mathcal{P} \exp \left( \oint_{C_1} A_\mu dx^\mu \right) \dots \operatorname{Tr} \mathcal{P} \exp \left( \oint_{C_n} A_\mu dx^\mu \right) \rangle. \qquad (2.39)$$

While our primary emphasis will be on the single-trace and double-trace Wilson loops, the principles discussed can be extended to higher-trace loops.

We begin by examining the scenario for U(N) and subsequently explore the more intricate cases of SU(N), specifically SU(2) and SU(3). Our analysis will demonstrate that for these groups, the loop equations close respectively on single-trace and double-trace Wilson loops, illustrating the nuanced dynamics of these gauge theories.

From the fundamental principle that the functional integral of a full derivative vanishes,<fnref n="10" /> we derive the U(N) Makeenko-Migdal Loop Equation, details of which are extensively discussed in the literature [57, 61, 71]:

$$\frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu \nu}(x)} W(C_x) = \lambda \oint_C dy_{\nu} \delta^{(d)}(x-y) W(C_{xy}^{(1)}, C_{yx}^{(2)}) \qquad (2.40)$$

In this equation, $\delta\sigma_{\mu\nu}$ represents the area derivative, a concept crucial to the formulation of loop equations in gauge theories. The operator $\frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta\sigma_{\mu\nu}(x)}$ and further details on loop kinematics can be found in [61, 72, 73].

A simple generalization of the single-loop MME allows to derive the following U(N) loop equation for the double-loop Wilson average, in parallel to its lattice version (2.31):

$$\frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu \nu}(x)} W(C_x, C') = \\ = \lambda \oint_C dy_{\nu} \delta^{(d)}(x-y) W(C_{xy}^{(1)}, C_{yx}^{(2)}, C') + \frac{\lambda}{N^2} \oint_C dy_{\nu} \delta^{(d)}(x-y) W(C_{xy} \circ C'_{yy}) \qquad (2.41)$$

In this equation, $C_{xy} \circ C'_{yy}$ denotes the recomposition of two contours into a single loop at the point where they intersect, $x = y$. This result illustrates that the U(N) loop equations naturally extend to encompass multi-trace configurations, effectively demonstrating that these equations close on N-trace Wilson loops.

Our exploration of the SU(N) loop equations begins with the observation that the U(N) Wilson average can be decomposed into two independent components when considering the decomposition $U(N) \to SU(N) \times U(1)$:

$$W(C) = \tilde{W}(C) \cdot w(C) \qquad (2.42)$$

<footnote n="10">¹⁰For the vanishing fields at space-time infinity as the boundary conditions.</footnote>
Here, $W(C)$ satisfies the Makeenko-Migdal loop equation as expressed in (2.40), and $w$ adheres to the abelian loop equation:

$$ \frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)} w(C_x) = \frac{g^2}{N} w(C) \oint_C dy_{\nu} \delta^{(d)}(x-y) \quad (2.43) $$

Applying the area derivative, which acts as a linear operator, we leverage the geometric derivative $\frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)}$ to both parts of the decomposed Wilson loop. This application yields the following combined equation for SU(N):

$$ \frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)} W(C_x) = w(C) \frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)} \tilde{W}(C_x) + \frac{g^2}{N} \tilde{W}(C) w(C) \oint_C dy_{\nu} \delta^{(d)}(x-y) \quad (2.44) $$

Upon applying the U(N) loop equations (2.40) and utilizing the factorization (2.42), we simplify the right-hand side (r.h.s.) and, by cancelling $w(C)$ on both sides, derive the SU(N) loop equation for the single-loop Wilson average:

$$ \frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)} \tilde{W}(C_x) = \lambda \oint_C dy_{\nu} \delta^{(d)}(x-y) \left[ \tilde{W}(C_{xy}^{(1)}, C_{yx}^{(2)}) - \frac{1}{N^2} \tilde{W}(C) \right] \quad (2.45) $$

This equation indicates a refined understanding of the non-abelian contributions specific to SU(N) without the influence of the U(1) component.

Similarly, for the connected double-loop Wilson average for U(N) gauge theory we also have the factorization:

$$ W_{U(N)}(C, C') = \tilde{W}_{SU(N)}(C, C') \cdot w_{U(1)}(C \circ C') \quad (2.46) $$

We apply the same procedure to the factorized double loop as we did for the single loop, leading to the following SU(N) double-loop loop equation:

$$ \frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)} \tilde{W}(C_x, C') = \lambda \oint_C dy_{\nu} \delta^{(d)}(x-y) \left[ \tilde{W}(C_{xy}^{(1)}, C_{yx}^{(2)}, C') - \frac{1}{N^2} \tilde{W}(C, C') \right] + \\ + \lambda \oint_{C'} dy_{\nu} \delta^{(d)}(x-y) \left[ \tilde{W}(C_{xx} \circ C'_{yy}) - \frac{1}{N^2} \tilde{W}(C, C') \right] \quad (2.47) $$

Utilizing the trace identities elaborated in section 2.1, we arrive at the general conclusion that the SU(N) loop equations close on $(N-1)$-trace variables. This principle significantly simplifies the mathematical framework for gauge theories by reducing the complexity of multi-loop interactions to more manageable forms.

Specifically for the SU(2) case, the application of the trace identity (1.1)—expressed here as:

$$ W(C_1, C_2) = W(C_1 \circ C_2) + W(C_1 \circ \overline{C_2}), \quad (2.48) $$

transform the loop equation (2.45) into a closed form within the single-loop space. The revised loop equation is:

$$ \frac{\partial}{\partial x_{\mu}} \frac{\delta}{\delta \sigma_{\mu\nu}(x)} W(C_x) = \frac{\lambda}{2} \oint_C dy_{\nu} \delta^{(d)}(x-y) \left[ \frac{1}{2} W(C_x) + W(C_{xy}^{(1)} \circ \overline{C_{yx}^{(2)}}) \right] \quad (2.49) $$
In this equation, the notation ◦ denotes the composition of two loops at open indices into a single-trace loop, and the overbar (¯) indicates the inversion of the ordered product, akin to the Hermitian conjugation of the Wilson line.

For SU(3), the loop equations can similarly be closed on the double-loop space, thanks to the identity (1.2), which for this scenario is expressed as follows:

$$
\begin{aligned}
W(C_1, C_2, C_3) = & \frac{1}{3} W(C_1, C_2 \circ C_3) + \frac{1}{3} W(C_2, C_1 \circ C_3) + \frac{1}{3} W(C_3, C_1 \circ C_2) - \\
& - \frac{1}{9} W(C_1 \circ C_2 \circ C_3) - \frac{1}{9} W(C_1 \circ C_3 \circ C_2) + \frac{1}{3} W(C_1 \circ \bar{C}_3, C_2 \circ \bar{C}_3) \\
& - \frac{1}{9} W(C_1 \circ \bar{C}_3 \circ C_2 \circ \bar{C}_3)
\end{aligned}
$$

Applying this identity to the first term in the second line of (2.47), we derive a refined expression for the double-loop interaction:

$$
\begin{aligned}
\tilde{W}(C_{xy}^{(1)}, C_{yx}^{(2)}, C') = & \frac{1}{3}\tilde{W}(C_{xy}^{(1)}, C_{yx}^{(2)} \circ C') + \frac{1}{3}\tilde{W}(C_{xy}^{(2)}, C_{yx}^{(1)} \circ C') + \frac{1}{3}\tilde{W}(C_{xy}^{(1)} \circ C_{yx}^{(2)}, C') - \\
& - \frac{1}{9}\tilde{W}(C_{xy}^{(1)} \circ C_{yx}^{(2)} \circ C') - \frac{1}{9}\tilde{W}(C_{xy}^{(2)} \circ C_{yx}^{(1)} \circ C') \\
& + \frac{1}{3}\tilde{W}(C_{xy}^{(1)} \circ \bar{C}', C_{yx}^{(2)} \circ \bar{C}') - \frac{1}{9}\tilde{W}(C_{xy}^{(1)} \circ \bar{C}' \circ C_{yx}^{(2)} \circ \bar{C}')
\end{aligned}
$$

Although contour $C'$ does not necessarily intersect the point $x = y$, we can conceptually connect it by introducing a backtrack path $\Gamma_{xz} \otimes \Gamma_{zx}$, leading to the redefined interaction:

$$
\tilde{W}(C_{xy}^{(1)}, C_{yx}^{(2)}, C') \equiv \tilde{W}(C_{xy}^{(1)}, C_{yx}^{(2)}, \Gamma_{xz} \circ C'_{zz} \circ \Gamma_{zx}) \quad (2.50)
$$

Then, applying the identities outlined above, we can effectively recast the three-loop average into a combination of one- and two-loop terms, as shown on the right-hand side of the last identity. This substitution leads to the SU(3) loop equation for the Wilson averages, effectively closing the loop equations on the 1- and 2-loop space. Here, we will not make this substitution for the sake of brevity.

3 Positivity

In the bootstrap method under consideration, the assertion of positivity is contingent upon
the specific definition of the inner product within the operator space or subspace. More
precisely, consider an expansion of operators in a given basis:

$$
\mathcal{O} = \sum_i \alpha_i \mathcal{O}_i \qquad (3.1)
$$

Assuming a well-defined inner product in the operator space, the resulting matrix of the
quadratic form of the scalar product is guaranteed to be positive semidefinite:

$$
\langle \mathcal{O} | \mathcal{O} \rangle = \langle \mathcal{O}^{\dagger} \mathcal{O} \rangle = \alpha^{*\mathrm{T}} \mathcal{M} \alpha \ge 0, \forall \alpha \Leftrightarrow \mathcal{M} \succeq 0. \quad (3.2)
$$

Here, the matrix $\mathcal{M}$ is defined by the entries:

$$
\mathcal{M}_{ij} = \langle \mathcal{O}_i^\dagger \mathcal{O}_j \rangle \qquad (3.3)
$$

Note that in this context, † denotes the adjoint operators under the specified inner product.
## 3.1 Hermitian conjugation

This positivity condition is particularly evident when the measure is positive (i.e., devoid of a sign-problem). It is also the most well-understood category of positivity utilized in the bootstrap approach. In addition to the inherent positivity of the measure itself, this category encapsulates crucial information about the measure space [7, 74]. The resolution of the Hamburger moment problem exemplifies this principle.

Concretely, the positivity condition for this category is expressed as follows:

$$ \langle \mathcal{O} | \mathcal{O} \rangle = \int d\mu \mathcal{O}^* \mathcal{O} \geq 0 \qquad (3.4) $$

where * denotes complex conjugation. Despite its seemingly straightforward appearance — as a result of integrating a positive function over a positive measure — this positivity condition can be remarkably potent, given its validity for any $\mathcal{O}$. The convergence of the bootstrap approach, predicated solely on the positivity of Hermitian conjugation, has been demonstrated in several model examples [7, 74].

In the context of the matrix model, we articulate the following positivity condition:

$$ \left\langle \frac{1}{Nm} \left( \sum_{\{O_i\}} \alpha_{\{O_i\}} O_{1,a_1 b_1} \dots O_{m,a_m b_m} \right)^{\dagger} \left( \sum_{\{O_i\}} \alpha_{\{O_i\}} O_{1,a_1 b_1} \dots O_{m,a_m b_m} \right) \right\rangle \geq 0 \qquad (3.5) $$

where $O_i$ represents specific matrix operators and $a_i, b_i$ are their respective matrix indices, which are summed over. The coefficients list $\alpha_{\{O_i\}}$ is arbitrary here, as in (3.2). The condition (3.5) encapsulates $m$-traces operators. Intuitively, it appears sufficient to consider $m = N - 1$ for the SU(N) matrix model and $m = N$ for both the U(N) and more generally the GL(N) matrix models, as these configurations correspond to the truncation points of the loop equations. Further numerical analysis of the SU(N) toy model (2.8) indicates that the imposition of N-traces positivity conditions does not introduce additional constraints beyond those imposed by $(N - 1)$-traces positivity conditions. A more systematic exploration of this truncation phenomenon remains an open question.

As an example, for the SU(3) version of the toy model (2.8), with a truncation $\text{Tr}U^a \text{Tr}U^b$, $|a| + |b| \leq 2$, we have the following positivity condition:

$$ \left. \begin{array}{c} U_{a_1 b_1}^{\dagger} \mathbb{1}_{a_2 b_2} \quad \mathbb{1}_{a_1 b_1} U_{a_2 b_2}^{\dagger} \quad \mathbb{1}_{a_1 b_1} \mathbb{1}_{a_2 b_2} \quad \mathbb{1}_{a_1 b_1} U_{a_2 b_2} \quad U_{a_1 b_1} \mathbb{1}_{a_2 b_2} \\ U_{b_1 a_1} \mathbb{1}_{b_2 a_2} \quad \left\langle \text{Tr} U^{\dagger} \text{Tr} U \right\rangle \quad \left\langle \text{Tr} U \right\rangle \quad \left\langle \text{Tr} U^2 \right\rangle \quad \left\langle \text{Tr} U^2 \right\rangle \\ \mathbb{1}_{b_1 a_1} U_{b_2 a_2} \quad \left\langle \text{Tr} U^{\dagger} \text{Tr} U \right\rangle \quad 1 \quad \left\langle \text{Tr} U \right\rangle \quad \left\langle \text{Tr} U^2 \right\rangle \quad \left\langle \text{Tr} U^2 \right\rangle \\ \mathbb{1}_{b_1 a_1} \mathbb{1}_{b_2 a_2} \quad \left\langle \text{Tr} U \right\rangle \quad \left\langle \text{Tr} U \right\rangle \quad 1 \quad \left\langle \text{Tr} U \right\rangle \quad \left\langle \text{Tr} U \right\rangle \\ \mathbb{1}_{b_1 a_1} U_{b_2 a_2}^{\dagger} \quad \left\langle \text{Tr} U^2 \right\rangle \quad \left\langle \text{Tr} U^2 \right\rangle \quad \left\langle \text{Tr} U \right\rangle \quad 1 \quad \left\langle \text{Tr} U^{\dagger} \text{Tr} U \right\rangle \\ U_{b_1 a_1}^{\dagger} \mathbb{1}_{b_2 a_2} \quad \left\langle \text{Tr} U^2 \right\rangle \quad \left\langle \text{Tr} U^2 \right\rangle \quad \left\langle \text{Tr} U \right\rangle \quad \left\langle \text{Tr} U^{\dagger} \text{Tr} U \right\rangle \quad 1 \end{array} \right\rangle \geq 0. \qquad (3.6) $$

A more intricate example can be found in the following scenario, which outlines the positivity condition for the 2-dimensional lattice Yang-Mills theory. Here, all open Wilson lines initiate and terminate at the origin, marked by a red dot, and are subject to a length cutoff = 4.
This configuration corresponds to a 9 × 9 matrix that includes 6 independent Wilson loops:

<figure><img src="image_8.png" /></figure>

## 3.2 Reflection Positivity

The study of reflection positivity was initiated in the context of constructive quantum field theory [75, 76], where it comes as part of the condition to guarantee that a Euclidean field theory defines a consistent relativistic quantum field theory through Wick rotation. The conformal bootstrap [77, 78]<fnref n="11" /> and the infrared bound [81] are among the most fruitful applications of reflection positivity.

As illustrated in figure 6, it is widely accepted that there are three known types of reflection positivity in lattice Yang-Mills theory on the square lattice: site-reflection, link-reflection [82], and diagonal-reflection [8]. The inner product defined by the reflection is straightforward in the subspace of operators located on the positive side of the reflection plane:<fnref n="12" />

$$ \langle \mathcal{O}^+ | \mathcal{O}^+ \rangle = \langle \mathcal{O}^+ \Theta \mathcal{O}^+ \rangle \qquad (3.8) $$

In the infrared (IR) regime of the theory, this same inner product undergoes a transformation via Wick rotation to become the inner product of the physical Hilbert space [82].

<footnote n="11">¹¹For recent developments in this field, see [79] and [80], which provide comprehensive updates and insights into current research trends.</footnote>

<footnote n="12">¹²The positive side of the reflection plane can be chosen arbitrarily. It is marked by “+” superscript for the corresponding operators.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Here are examples of site-reflection, link-reflection, and diagonal-reflection. It is important to note that the intersection of the Wilson path with the reflection plane $\pi$ does not necessarily occur at the same point, as illustrated in the figure.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** An example of plaquette who intersect with the reflection plane $\pi$.</figcaption></figure>

As an example, below is a closer investigation of the diagonal reflection positivity. It treats the diagonal of the square lattice as the surface of quantization. The reflection operation defined on a lattice site with lattice coordinates $(x_1, x_2, \dots, x_n)$ is formalized as follows:

$$ \Theta(x_1, x_2, \dots, x_n) = (x_2, x_1, \dots, x_n) \qquad (3.9) $$

This operation swaps the first two coordinates while leaving the remaining coordinates unchanged, reflecting the site across the line or plane that interchanges these dimensions.

For a unitary matrix associated with a link on the lattice, the reflection transformation is defined as follows:

$$ \Theta U_{x \to y} = U_{\Theta x \to \Theta y}^{\dagger} \qquad (3.10) $$

The proof of the reflection positivity is straightforward. We begin by noticing the following fact:

$$ \int \mathcal{D}U \mathcal{O}^{+} \Theta \mathcal{O}^{+} = \int \mathcal{D}U^{+} \mathcal{O}^{+} \int \mathcal{D}U^{\dagger} \Theta \mathcal{O}^{+} = \| \int \mathcal{D}U^{+} \mathcal{O}^{+} \|^{2} \geq 0 \qquad (3.11) $$

Here $\mathcal{O}^{+}$ can be arbitrary operators located in the positive side of the reflection plane, or in our case as in figure 7: $x_2 > x_1$.
To consider the expectation value with the Wilson action, we notice that the main obstacle is the terms $S_Δ$ in the action that are represented by the plaquettes intersected (along their diagonals, as in figure 7) by the reflection plane $π$:

$$ \int \mathcal{D}U \mathcal{O}^+ \Theta \mathcal{O}^+ \exp(-S) = \int \mathcal{D}U \mathcal{O}^+ \Theta \mathcal{O}^+ \exp(-S_+ - \Theta S_- + S_Δ) \\ = \int \mathcal{D}U \mathcal{O}^+ \exp(-S_+) \Theta(\mathcal{O}^+ \exp(-S_+)) \exp(-S_Δ) \qquad (3.12) $$

But this is appears not a problem once we notice that:

$$ \exp(-S_Δ) = \prod_{p \in Δ} \exp(\frac{N}{λ} \mathrm{Tr} U_p) = \prod_{p \in Δ} \exp(\frac{N}{λ} \mathrm{Tr} U_p^+ \Theta U_p^+) \qquad (3.13) $$

Here the product of $p \in Δ$ means the product of all the plaquettes intersecting the reflection plane $π$. Expanding the exponent and the infinite product, we see that this term is a positive sum of terms of the form $\mathcal{O}^+ \Theta \mathcal{O}^+$, which means that the whole expression (3.12) is positive.

## 3.3 Symmetry Reduction

The primary advantage of the bootstrap approach is that it naturally takes into account the symmetries of the theory. This feature has been observed in such subjects as the conformal bootstrap and S-matrix bootstrap [83–86]. In the context of the bootstrap approach discussed in this project, the symmetries of the model under consideration are often manifested in the inner products defined as follows:

$$ \langle (g \circ \mathcal{O}_1) | (g \circ \mathcal{O}_2) \rangle = \langle \mathcal{O}_1 | \mathcal{O}_2 \rangle, \forall g \in G \qquad (3.14) $$

Here, $g \circ \mathcal{O}$ denotes the transformation of the operator $\mathcal{O}$ under the symmetry group $G$, typically representing a reducible representation of this group.

Thanks to Schur's lemma, we know that due to group-theoretical reasons, the operators belonging to different irreducible representations of the symmetry group must be orthogonal. Consequently, the positive semidefinite matrix defined in (3.2) becomes block diagonal, where each block corresponds to an irreducible representation. However, these blocks can be further decomposed using the methods of Invariant Semidefinite Programming.<fnref n="13" /> More specifically, suppose the space of operators forms a reducible representation of the symmetry group $G$ with the following expansion into irreducible representations:

$$ V = \bigoplus_{k=1}^{D} \mathrm{Rep}_k \oplus^{m_k}, \qquad (3.15) $$

where $m_k$ is the multiplicity of a given irreducible representation $\mathrm{Rep}_k$. The block corresponding to this irreducible representation then has dimensions $m_k \times m_k$. It is noteworthy that the final dimension of the positivity matrix is $\sum_k m_k$, rather than $\mathrm{dim}V = \sum_k d_k m_k$.<fnref n="14" />

<footnote n="13">¹³For example, the interested reader could refer to [87].</footnote>

<footnote n="14">¹⁴This can be understood from the classical example of a Hamiltonian with SU(2) symmetry, where states with the same $l$ but different $m$ possess identical energy eigenvalues. In such cases, the lesson is that positivity among the highest weight vectors already captures all necessary information regarding the positivity.</footnote>
<table><caption>**Table 1.** Invariant groups of correlation and reflection matrices.</caption><thead><tr><th>Dimension</th><th>Hermitian<br>Conjugation</th><th>site&amp;link<br>reflection</th><th>diagonal<br>reflection</th></tr></thead><tbody><tr><td>2</td><td>B<sub>2</sub> × Z<sub>2</sub></td><td>Z<sub>2</sub> × Z<sub>2</sub></td><td>Z<sub>2</sub> × Z<sub>2</sub></td></tr><tr><td>3</td><td>B<sub>3</sub> × Z<sub>2</sub></td><td>B<sub>2</sub> × Z<sub>2</sub></td><td>Z<sub>2</sub><sup>3</sup></td></tr><tr><td>4</td><td>B<sub>4</sub> × Z<sub>2</sub></td><td>B<sub>3</sub> × Z<sub>2</sub></td><td>B<sub>2</sub> × Z<sub>2</sub><sup>2</sup></td></tr></tbody></table>

Implementing this symmetry reduction is analogous to projecting the physical state with respect to spin and parity in conformal or S-matrix bootstrap. The systematic steps to obtain the projector are outlined as follows [7]:

1. Identify a specific realization of every irreducible representation (irrep) of the invariant group using GAP software [88].

2. Apply the algorithm proposed in [89] to find an equivalent real representation, if the irrep provided by GAP is complex.

3. Decompose into such irreps using the projector to $Rep_k$ [90]:

$$p_{\alpha\alpha,k} = \frac{\dim(\mathrm{Rep}_k)}{\dim G} \sum_{g \in G} r_{\alpha\alpha}(g^{-1})g \quad (3.16)$$

Here, $r_{\alpha\beta}$ represents the matrix elements of a real representation identified at step 2, with $\alpha, \beta = 1, 2, \dots, \dim(\mathrm{Rep}_k)$. By setting $\alpha = 1$, $P_k = p_{11,k}$ serves as a projector to $\mathrm{Rep}_k$.

Some remarks follow:

* At the second step, the preference for a real representation stems from the advantage that real Semidefinite-Programming typically offers big advantages w.r.t. complex Semidefinite-Programming.<fnref n="15" /> Fortuitously, for the lattice Yang-Mills theory under investigation, all the groups possess real representations. A systematic method to determine whether a finite group has real irreducible representations is documented in [90].

* At the third step, selecting $\alpha = 1$ for the final projector might appear arbitrary. However, the key insight is that choosing different $\alpha$ values yields SDP blocks that are fully equivalent to those obtained with $\alpha = 1$.

For the model under consideration, specifically SU(N) or U(N) lattice Yang-Mills theory, we observe the global symmetries listed in table 1. The notation for the group $B_n$ is adopted from the nomenclature used on the Wikipedia page for the Hyperoctahedral Group (Hyperoctahedral Group). In more standard terms, $B_2$ and $B_3$ correspond to the dihedral and octahedral groups, respectively:

$$B_2 \simeq D_4, \quad B_3 \simeq O_h \quad (3.17)$$

<footnote n="15">¹⁵We note that all multi-trace Wilson loops are real, a property stemming from the unimodularity of the corresponding Haar measure.</footnote>
The $\mathbb{Z}_2$ in table 1 represents the charge conjugation $C$, which acts by inverting the direction of a Wilson loop.

The discussion above parallels the considerations in the planar limit as presented by the author [8]. However, in the case of multi-trace operators, an additional symmetry emerges — specifically: the permutation group among different traces. Consequently, there is an $S_m$ group that complements the global symmetry group listed in table 1.

To illustrate this principle, we conduct the symmetry reduction of (3.6). Besides the permutation group previously discussed, the corresponding model (2.8) also exhibits charge conjugation symmetry, where $U \to U^\dagger$. Therefore, we can categorize the operators according to their representations under the $\mathbb{Z}_2$ group:

$$ 
\begin{aligned}
(+, +) : \quad & 1, \quad U_{a_1 b_1}^\dagger \mathbb{1}_{a_2 b_2} + \mathbb{1}_{a_1 b_1} U_{a_2 b_2}^\dagger + \mathbb{1}_{a_1 b_1} U_{a_2 b_2} + U_{a_1 b_1} \mathbb{1}_{a_2 b_2} \\
(-, +) : \quad & (U_{a_1 b_1}^\dagger \mathbb{1}_{a_2 b_2} + \mathbb{1}_{a_1 b_1} U_{a_2 b_2}^\dagger) - (\mathbb{1}_{a_1 b_1} U_{a_2 b_2} + U_{a_1 b_1} \mathbb{1}_{a_2 b_2}) \\
(+, -) : \quad & (U_{a_1 b_1}^\dagger \mathbb{1}_{a_2 b_2} - \mathbb{1}_{a_1 b_1} U_{a_2 b_2}^\dagger) + (\mathbb{1}_{a_1 b_1} U_{a_2 b_2} - U_{a_1 b_1} \mathbb{1}_{a_2 b_2}) \\
(-, -) : \quad & (U_{a_1 b_1}^\dagger \mathbb{1}_{a_2 b_2} - \mathbb{1}_{a_1 b_1} U_{a_2 b_2}^\dagger) - (\mathbb{1}_{a_1 b_1} U_{a_2 b_2} - U_{a_1 b_1} \mathbb{1}_{a_2 b_2})
\end{aligned} 
\quad (3.18) $$

Consequently, the positivity condition detailed in (3.6) is further refined to the following specifications:

$$ 
\begin{aligned}
& \left( \begin{array}{c} 1 \\ \langle \mathrm{Tr} U \rangle \\ \langle \mathrm{Tr} U \rangle \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U^\dagger \rangle + \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U \rangle + \frac{1}{4} \langle \mathrm{Tr} U^2 \rangle + \frac{1}{4} \end{array} \right) \succeq 0, \\
& \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U^\dagger \rangle - \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U \rangle - \frac{1}{4} \langle \mathrm{Tr} U^2 \rangle + \frac{1}{4} \succeq 0 \\
& - \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U^\dagger \rangle - \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U \rangle + \frac{1}{4} \langle \mathrm{Tr} U^2 \rangle + \frac{1}{4} \succeq 0 \\
& - \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U^\dagger \rangle + \frac{1}{4} \langle \mathrm{Tr} U \mathrm{Tr} U \rangle - \frac{1}{4} \langle \mathrm{Tr} U^2 \rangle + \frac{1}{4} \succeq 0
\end{aligned} 
\quad (3.19) $$

The symmetry reduction of the matrix described in (3.7) has been thoroughly detailed in the supplementary material of [8]. To avoid redundancy, we will not repeat the discussion here.

In concluding this discussion, it is important to note that while our analysis has primarily focused on finite groups, the results are readily generalizable to Lie groups, particularly compact Lie groups [91]. This extension is supported by established correspondences between compact Lie groups and finite groups, suggesting that our findings may have broader applications across various mathematical structures in theoretical physics. This potential for generalization not only underscores the versatility of our theoretical framework but also highlights its relevance for ongoing research in the field of quantum mechanics and field theory, where Lie groups play a crucial role.

## 4 Bootstrap procedure

Our previous discussion of the loop equations in section 2 and the positivity conditions in section 3 facilitate the straightforward bounding of several observables within the corresponding theory. In the subsequent subsection, we will present some immediate results for the bootstrap. However, for the lattice Yang-Mills theory, formulating a large-scale bootstrap
<figure><img src="image_11.png" /><figcaption>**Figure 8.** This figure presents the bounds for the SU(3) one-matrix model. On the left side of the figure, we depict both the upper and lower bounds for $\langle \text{Tr}U \rangle$ at different values of $\Lambda$: 1 (gray), 2 (blue), and 3 (red). The black dashed line represents the exact value of the model. On the right side of the figure, we illustrate the difference between the upper bound and the exact value on a logarithmic scale, which demonstrates a clear exponential convergence.</figcaption></figure>

problem is not as straightforward as it might initially appear. The primary challenge lies in the selection of operators for our bootstrap analysis. In section 4.2, we propose a hierarchy that assigns each operator (open Wilson line) and loop variable (Wilson loop) a specific level, and truncates the infinite space of Wilson lines according to its level. At last, we present the results we have achieved so far, mostly concerned with the free energy per plaquette and string tension, which we compare with the existing Monte Carlo calculations, as well as strong and weak coupling expansions.

## 4.1 Initiate the bootstrap

The bootstrap of the one-matrix model (2.8) is straightforwardly implemented. Here, we present the SU(3) case of equation (2.8) as an example.<fnref n="16" /> The fundamental components are already summarized in section 2 and section 3. The set of loop equations combines those from (2.14) with the trace identities from (2.9). To solve this system, the **Reduce** function in Mathematica is employed, revealing that all double-trace variables are linear functions of two variables: $\langle \text{Tr}U \rangle$ and $\langle \text{Tr}U \text{Tr}U^\dagger \rangle$.

Concerning the positivity condition, as exemplified in (3.6), it arises from the inner product of double trace operators $U_{a_1 b_1}^m U_{a_2 b_2}^n$. We truncate the space of operators such that $|m| + |n| \le \Lambda$.<fnref n="17" /> Increasing $\Lambda$ yields progressively tighter dynamical bounds on the observables.

The results of the bootstrap up to $\Lambda = 3$ are illustrated in figure 8. Both the upper and lower bounds are observed converging towards the exact value.<fnref n="18" /> In the right panel of figure 8, the convergence of the upper bound to the exact value for $\Lambda = 1, 2, 3$ is depicted, demonstrating a clear exponential trend.

<footnote n="16">¹⁶In appendix C, we also present the bootstrap results for SU(2), U(2), and U($\infty$).</footnote>

<footnote n="17">¹⁷The positivity condition (3.6) corresponds to $\Lambda = 1$.</footnote>

<footnote n="18">¹⁸The exact value for the current one-matrix model is documented in appendix A. For illustrative Mathematica code used to generate the plot in figure 8 and for numerical examples referenced in appendix C, please visit the GitHub page https://github.com/Canonical111/Unitary-Matrix-Model-Bootstrap-Example.</footnote>
Before exploring the details of the hierarchy discussed in the next section, it is important to note that from the equations and positivity conditions available for the lattice Yang-Mills theory, we can already derive non-trivial upper bounds on the observables. Specifically, for SU(2) and general U(N) models, we integrate the loop equations (2.32) and (2.33) with the positivity condition (3.7). This setup constitutes a bootstrap problem involving six variables, which encompass all variables below level 2 in the hierarchy discussed in the subsequent section. By minimizing the plaquette average, we obtain non-trivial bounds for SU(2) and uniform bounds for U(N) (including U(∞)). For instance, at λ = 2, the bounds are:

$$ \langle \mathrm{Tr} U \rangle \le \begin{cases} 0.693 & \text{for } \mathrm{U}(N) \\ 0.75755 & \text{for } \mathrm{SU}(2) \end{cases} \qquad (4.1) $$

## 4.2 Hierarchy

Our framework for the bootstrap program applied to the finite N matrix models and lattice Yang-Mills theories can be summarized as follows:

$$ \begin{array}{ll} \min / \max & \llap{\boxed{\phantom{000}}}, \\ \text{subject to} & \mathrm{MM~loop~equations}, \\ & \mathrm{HerM}^{\mathrm{irrep}} \ge 0, \\ & \mathrm{RefM}^{\mathrm{irrep}} \ge 0, \times 3 \end{array} \qquad (4.2) $$

Here, the objective function can, in principle, minimize any convex function of the Wilson loops. We enforce the (multi-trace) Makeenko-Migdal loop equations and apply trace identities to eliminate higher multi-trace Wilson loops. Additionally, four types of positivity conditions are imposed, three of which are reflection positivities, each factorizable into different irreducible representations of the corresponding symmetry group. In the case of large N, we also implemented in [7, 8], in addition to these constraints, a relaxation of the quadratic variables and ensure the positivity of the relaxation matrix $\mathcal{R} \ge 0$.

One critical feature remains to add: to effectively formulate the bootstrap problem (4.2), we must implement a finite truncation of the Wilson loops. This truncation is essential for the efficiency of the bootstrap process. Below, we will illustrate this truncation method (hierarchy) using lattice Yang-Mills theory as an example. We believe that this principle can be universally applied to bootstrap methods within this category of problems.

Our practical selection of positivity generalizes (3.7) by focusing only on the positivity of the inner products of Wilson lines that initiate and terminate at the origin point. This selection is based on the fact that such a subset of Wilson lines forms a representation of the largest possible symmetry group. We then derive all necessary equations among the Wilson loops within the matrix for the inner products of these Wilson lines.

For the selection of Wilson lines as the basis of the inner product, we employ a hierarchical system, where each Wilson line starting and ending at zero is assigned an integer number representing its level. This hierarchy is defined recursively: we begin with the identity operator at level 0 and a single plaquette at level 1. Subsequent operators are derived from the Schwinger-Dyson variations, as illustrated in the first line of figure 3. For instance,
in the context of SU(2) lattice gauge theory in two dimensions, the following defines a level 2 Wilson line:

<figure><img src="image_12.png" /></figure>

The red dot in our diagrams represents the origin point, where the Wilson line has an open index. We also define the hierarchy of Wilson loop expectations, which is determined by the sum of the levels of the Wilson lines whose inner products generate these loops. It is important to note that if multiple levels could be defined recursively for a Wilson line/loop, the smallest one is selected as its level.

As an illustrative example, we list all level two Wilson loops for the two-dimensional SU(2) lattice gauge theory:

<figure><img src="image_13.png" /></figure>

These loops, along with the single plaquette Wilson loop, comprise all the variables that appear in the matrix (3.7).

This hierarchically defined structure ensures that the variables in our positivity conditions are interrelated through the loop equations. When constructing the inner product matrix, if we consider Wilson lines up to level $\Lambda$, then the longest Wilson loop in our problem will be $8 \times \Lambda$. In the current project, we typically consider up to $\Lambda = 3$, corresponding to a maximum loop length of 24.

Utilizing the hierarchical system, it becomes straightforward to derive relatively tight bounds on observables in two-dimensional lattice gauge theory. Up to $\Lambda = 3$, which includes Wilson loops of maximum length 24, we analyze a total of 8335 Wilson loops. These loops are subject to 14591 loop equations.<fnref n="19" /> After processing, it is found that 7291 of these loops can be linearly represented by the remaining 1044 loops. The positivity condition for this setup involves 18 blocks, with sizes given by:

$$54, 52, 46, 45, 98, 45, 46, 52, 53, 98, 30, 27, 21, 24, 16, 11, 8, 12 \qquad (4.5)$$

This configuration presents a complex, yet manageable, bootstrap problem. With the aid of highly optimized software, the entire computation for the $\Lambda = 3$ level, which includes 30 pairs of bounds, is completed within one minute.<fnref n="20" /> The results of the bounds and their convergence are depicted in figure 9.

<footnote n="19">¹⁹This indicates a significant level of linear dependency within the space of loop equations.</footnote>

<footnote n="20">²⁰The solver for SDP we were using for the entire project is MOSEK.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 9.** The results of our bootstrap calculations in two dimensions for the plaquette average $u_P = \langle \text{Tr}U_P \rangle$ as a function of the coupling $\lambda = \frac{2N^2}{\beta}$ are depicted in the following manner: blue lines represent interpolations of the upper and lower margins for the maximum length of loops $L_{\text{max}} \le 16$ i.e. $\Lambda = 2$ used in the bootstrap procedure for loop equations. Red lines illustrate the same for $L_{\text{max}} \le 24$ or $\Lambda = 3$. The gray area denotes the allowed values of $u_P$ for $\Lambda = 1$. The black dotted line indicates the exact value, the explicit form of which can be found in appendix A. The right plot illustrates the difference between the upper bound and the exact value, demonstrating an exponential convergence behavior.</figcaption></figure>

## 4.3 Higher dimensional bound and string tension

The entire framework developed for the two-dimensional scenario can, in principle, be extended to three and four dimensions without substantial conceptual challenges. However, practical implementation in higher dimensions poses significantly greater complexity. The complexity of extending our framework to three and four dimensions arises from two primary factors. Firstly, as the spacetime dimension increases, the number of Wilson loops grows exponentially with our truncation level, leading to a significantly larger computational space. Secondly, the loop equations for the three and four-dimensional lattice Yang-Mills theory become considerably more complex. This complexity manifests as a greater number of independent Wilson loops that must be considered once the loop equations are applied. Due to these practical difficulties, specific adaptations were necessary for conducting the $\Lambda = 3$ calculations, where $L_{\text{max}} = 24$. The actual strategy we adopted involved a further truncation of the problem. Specifically, we considered only half of the $\Lambda = 3$ Wilson lines in our calculations. Moreover, in the four-dimensional analysis, we focused exclusively on the reflection positivity aspect. This selective consideration was essential to manage the computational complexity and ensure the feasibility of the calculations within available resources.

We provide a detailed benchmark of our computational framework for $\Lambda = 3$ bootstrap in different spacetime dimensions in table 2.<fnref n="21" /> This benchmark highlights the performance metrics such as computational time and memory usage, contextualized by the specific demands of each dimensional analysis:<fnref n="22" />

<footnote n="21">²¹The bootstrap calculations at the $\Lambda = 2$ level are completed almost instantaneously, regardless of the spacetime dimension under consideration.</footnote>

<footnote n="22">²²The time consumption figures are based on executions performed using MOSEK on a machine equipped with two AMD EPYC 7282 processors.</footnote>
<table><caption>**Table 2.** The scale of the $\Lambda = 3$ bootstrap problem we addressed in this project. Here #var is the number of all the Wilson loops, #LE is the number of loop equations we considered, #SDPvar is the number of the free variables after imposing the loop equation. #Block and Blocksize indicate the size of the semidefinite programming, four example, in four dimensions, it means 40 positive semidefinite matrices with size around 300 × 300. The last two columns are the time and memory consumption using MOSEK to solve the optimization problem.</caption><thead><tr><th>Dim</th><th>#var</th><th>#LE</th><th>#SDPvar</th><th>#Block</th><th>Blocksize</th><th>Memory</th><th>Time</th></tr></thead><tbody><tr><td>2</td><td>8335</td><td>14591</td><td>1044</td><td>18</td><td>~ 50</td><td>~ 100 MB</td><td>~ 1 s</td></tr><tr><td>3</td><td>174387</td><td>98032</td><td>93561</td><td>38</td><td>~ 200</td><td>220 GB</td><td>~ 1 days</td></tr><tr><td>4</td><td>343851</td><td>152149</td><td>211912</td><td>40</td><td>~ 300</td><td>1 TB</td><td>~ 20 days</td></tr></tbody></table>

* #var: represents the total number of Wilson loops considered.

* #LP: denotes the number of loop equations incorporated into the analysis.

* #SDPvar: indicates the number of free variables remaining after the application of loop equations.

* #Block and Blocksize: these specify the size and number of blocks in the semidefinite programming; for instance, in the four-dimensional case, this involves 40 positive semidefinite matrices, each approximately 300 × 300 in size.

* The final two columns detail the computational resources required for a single optimization, specifically the time and memory consumption, using the MOSEK optimization software to solve the problem.

In the three-dimensional analysis presented in figure 10, we observe that even in the worst-case scenario, the allowed region for our results remains within around 0.1% of the exact value. Furthermore, the Monte Carlo data from [59] aligns perfectly with our bootstrap bounds, underscoring the robustness and accuracy of our computational approach.

It is commonly accepted that for the Wilson averages $W[T, L]$ of rectangular Wilson loops of considerable size $T \times L$, particularly in the regime where $T \gg L \gg 1$, the asymptotic behavior is expressed as:

$$ \log W[T, L] = -\sigma TL + \mathcal{O}(1/L). \qquad (4.6) $$

The first term in this equation is indicative of the confinement “area law”, characterized by the string tension $\sigma(\lambda)$. The string tension $\sigma(\lambda)$ is routinely extracted from Monte Carlo (MC) data employing the Creutz ratio [2]:

$$ a^2\sigma = -\log \frac{W(T, L)W(T-1, L-1)}{W(T-1, L)W(T, L-1)} \qquad (4.7) $$

Although the loops measured in our bootstrap procedure are relatively small compared to those used in various MC simulations for extracting $\sigma$, it is believed (as discussed in [92]) that the regime (4.6) is attained from very small physical distances, of the order of 0.5 fm.
<figure><img src="image_15.png" /><figcaption>**Figure 10.** The results of our bootstrap calculations in three dimensions for the plaquette average $u_P = \langle \text{Tr}U_P \rangle$ as a function of the coupling $\lambda = \frac{2N^2}{\beta}$ are depicted in the following manner: blue lines represent interpolations of the upper and lower margins for the maximum length of loops $L_{\text{max}} \le 16$ i.e. $\Lambda = 2$ used in the bootstrap procedure for loop equations. Red lines illustrate the same for $L_{\text{max}} \le 24$ or $\Lambda = 3$. The gray area denotes the allowed values of $u_P$ for $\Lambda = 1$. The purple points are the results of Monte Carlo simulation from [59]. For a more detailed view of the current plot, please refer to figure 1(a).</figcaption></figure>

Consequently, we attempt here to extract, at least qualitatively, $\sigma(\lambda)$ from our bootstrap data for loops of size $(T, L) = (3, 2)a$. In our analysis, we select the optimal solution corresponding to the lower bound of the one-plaquette Wilson loop and use it to extract the expectation values for the loop configurations (2, 3), (2, 2), (3, 1), (2, 1). While the choice to focus on the lower bound might appear arbitrary, it is important to note that the expectation values exhibit minimal variation when the optimal solution for the upper bound is used instead. This stability between the lower and upper bounds is corroborated by the similar allowed margins for larger Wilson loops, relative to those for the one-plaquette Wilson loop, as demonstrated in figure 10. This consistency across different bounds and loop sizes not only validates our methodological approach but also reinforces the reliability of our bootstrap predictions. This approach ensures that our results are robust, minimizing the impact of the specific choice of lower or upper bounds on the final outcomes.

The comparison of our results with those from Monte Carlo simulations [59] is depicted in figure 11. A notable observation is that for larger values of $\beta$, where the physical scale of lattice spacing $a$ is relatively small, our agreement with the Monte Carlo data is not perfect. Conversely, for smaller $\beta$, which corresponds to a larger lattice spacing $a$, there is a decent match with the Monte Carlo results. This outcome aligns with expectations, as precise estimations of the string tension from our bootstrap method are anticipated only when the product $(3, 2) \times a$ falls within the asymptotic region described by (4.6). This comparison not only illustrates how bootstrap data correlates with established Monte Carlo simulations
<figure><img src="image_16.png" /><figcaption>**Figure 11.** The results of the string tension calculations from the (2,3) rectangular Wilson loop are presented here, with a reminder that $\lambda = \frac{2N^2}{\beta}$. The brown dots represent Monte Carlo simulation results from [59], whereas the blue dots depict string tension extracted from bootstrap expectations for the loop configurations (2,3), (2,2), (3,1), (2,1).</figcaption></figure>

but also provides qualitative validation of our bootstrap method for extracting the string tension, even with relatively small loop sizes.

The bootstrap results for the four-dimensional lattice Yang-Mills theory are depicted in figure 12. We compare these results with perturbation theory [58] and Monte Carlo data from various sources [3, 58]. The comparison shows that both the Monte Carlo and perturbation theory results fall well within the region allowed by our bootstrap method, affirming the validity of our approach.

Unlike the scenarios in two and three dimensions, for $\Lambda = 3$ or $L_{max} = 24$, we derived only three pairs of bounds in the four-dimensional bootstrap. This limitation reflects the maximum computational capacity achievable with our current software and hardware setups. Moreover, we observed a slight lack of numerical precision in the lower bounds, suggesting that merely increasing computational resources might not significantly enhance the results. This observation underscores the need for possibly refining our computational strategies or methodologies to overcome these challenges in higher-dimensional settings.

# 5 Discussion and Prospects

Although our current exploration is limited to SU(2) lattice gauge theory, we believe that our approach applies to a more general framework encompassing SU(N) and U(N) gauge theories, both theoretically and practically. The theoretical applicability is evident and further substantiated by our explorations of SU(3) and U(2) one-matrix models. Practically, the feasibility of broader application is supported by a simple yet significant observation: within a fixed truncation level in our hierarchy, the single trace Wilson loop consistently dominates. The most advanced result achieved in our project so far is at truncation level $\Lambda = 6$, leading
<figure><img src="image_17.png" /><figcaption>**Figure 12.** The results of our bootstrap calculations in four dimensions for the plaquette average $u_P = \langle \text{Tr}U_P \rangle$ as a function of the coupling $\lambda = \frac{2N^2}{\beta}$, ($N = 2$), in the SU(2) lattice gauge theory are depicted in the following manner: blue lines represent interpolations of the upper and lower margins for the maximum length of loops $L_{max} \le 16$ i.e. $\Lambda = 2$ used in the bootstrap procedure for loop equations. Three pairs of red points illustrate the upper and lower bounds of $u_P$ for $L_{max} \le 24$ or $\Lambda = 3$ at three different values of the coupling $\lambda = 3.19, 3.48, 3.64$. The gray area denotes the allowed values of $u_P$ for $\Lambda = 1$. The blacked-dashed lines in the figure represent the strong and weak coupling expansions, extracted from [58]. Additionally, the results from Monte Carlo simulations are depicted using purple dots [3] and yellow dots [58], each set representing data from different sources. For a more detailed view of the current plot, please refer to figure 1(b).</figcaption></figure>

us to anticipate the possibility of establishing a non-trivial bound for SU(6) or U(5). Here, the lowest non-trivial five-trace loop equation corresponds to level 6 in our hierarchy. For higher values of $N$, while valid bounds can still be obtained, their dependence on $N$ tends to be trivial, as even the lowest order trace identity does not manifest below level $\Lambda = 6$.

Our current project has not overcome the main bottleneck facing the state of the art in lattice bootstrap, as it was already manifest in previous studies [8–10]: our analysis remains confined to a local region of the infinite lattice, despite notable technical achievements. To significantly advance beyond this limitation, we foresee the potential for implementing a coarse-graining procedure analogous to the Density Matrix Renormalization Group (DMRG) used in quantum spin chain studies [93]. This approach has recently shown substantial progress in related fields [36, 94]. Additionally, a more detailed understanding of the close-to-null operators in lattice gauge theory could lead to the development of specific functionals based on these operators. Recent advances in conformal bootstrap suggest that such strategies could considerably enhance analytic efficiency [95].

It would be exciting to extend our framework from the statistical ensemble to quantum systems involving finite $N$ matrix models. Particularly in scenarios like the finite SU(2) matrix quantum mechanics, where even the ground state remains inadequately understood [96], our approach could provide significant insights. Given the rigor of the bounds established by
our method, it could offer valuable complementary information to the numerical studies of matrix quantum mechanics, especially those related to gravity duals of such models [97]. Such an extension could enhance our understanding of quantum behavior in matrix models and potentially contribute to a broader field of quantum gravity research.

In addition, we aim to elucidate the effectiveness of reflection positivity, which has remained somewhat enigmatic. Previous studies, such as [7] on the matrix model and [74] on the lattice Ising model, have clarified the role of positivity from Hermitian conjugation; it is closely associated with the absence of a sign problem in the model under consideration. It is generally believed that merely imposing Hermitian conjugation positivity ensures the convergence of the bootstrap procedure. However, reflection positivity extends beyond this by guaranteeing the possibility of defining a physical Hilbert space on the spatial reflection plane, which suggests the potential for a Hamiltonian formulation of the theory in the infrared regime. This observation also hints that the bootstrap method on the lattice could potentially transcend the limitations imposed by the sign-problem, which is prevalent in the Monte Carlo simulation methods. This suggests a promising avenue for advancing non-perturbative studies in lattice theories without the computational burdens typically associated with sign issues.

The convergence of the bootstrap based on the reflection positivity is also non-conclusive. In a previous study of the large $N$ Lattice Yang-Mills theory [8], it was observed that imposing reflection positivity as the sole type of positivity condition did not yield any non-trivial upper bounds for $\langle \text{Tr}U \rangle$. However, in our current project and in the previous study of the lattice Ising model across two and three spacetime dimensions [10], we have observed that reflection positivity exhibits a markedly greater influence than Hermitian conjugation positivity. This observation suggests that reflection positivity may play a more pivotal role in shaping the dynamics and theoretical outcomes of such models.

The fact that SU(2) loop equations are linear and closed on the single-loop space (and similarly, the linear SU(3) loop equations are closed on the two-loop space) can have, in principle an important impact on the understanding of the structure of non-abelian gauge theories in general. It would be important to understand what are the redundancies for these equations in the loop space and what kind of boundary conditions (or maybe asymptotics on large loops) we should impose to single out the relevant physical solutions. These insights could also significantly improve the current algorithms and may enable an efficient application of neural network methods, to try to extrapolate our results to longer Wilson loops.

An important problem is to include the dynamical quarks into the bootstrap. In principle, the multi-loop equations allow the inclusion of internal quark loops and then sum up their shapes to render the physical results. Moreover, for SU(2) case all such multiloop configurations are in principle generated by single loop Wilson averages. A similar, though more involved statement is true for the two-loop space of SU(3) gauge theory. However, the practical bootstrap computations in this direction may be very involved.

It is also important to note that our formalism allows us to approach the baryonic loop averages, as defined and formulated in [98].

## Acknowledgments

We are grateful for the discussions with Yiming Chen, Minjae Cho, Barak Gabai, Nikolay Gromov, Martin Kruczenski, Henry Lin, Raghu Mahajan, Alexandre Migdal, Colin Oscar
Nancarrow, Wenliang Li, Zhijin Li, Ying-Hsuan Lin, Victor A. Rodriguez, Joshua Sandor, Stephen Shenker, Evgeny Sobko, Douglas Stanford, Yuan Xin and Xi Yin. V.K. thanks the Perimeter Institute for Theoretical Physics and the Simons Center for Geometry and Physics, where a part of this work has been done. Z.Z. would like to express his sincere gratitude to Pedro Vieira, who provided invaluable support and guidance during his stay at ICTP-SAIFR, São Paulo, where the final stage of the project is finished. Z.Z is supported by Simons Foundation grant #994308 for the Simons Collaboration on Confinement and QCD Strings.

## A Exact results for the unitary one-matrix model

In this section, we present some exact results of the one-matrix model discussed in the main text. The partition function $Z$ is defined as:

$$Z = \int dU e^{\frac{N}{\lambda}(\text{tr}U^{\dagger} + \text{tr}U)} \qquad (A.1)$$

where the Haar measure is given by:

$$dU_{U(N)} = \frac{1}{N!} \prod_{1 \le j < k \le N} ((e^{i\theta_j} - e^{i\theta_k})(e^{-i\theta_j} - e^{-i\theta_k})) \prod_{j=1}^{N} \frac{d\theta_j}{2\pi}, \qquad (A.2)$$

$$dU_{SU(N)} = \frac{1}{N!} \sum_q \delta\left(\sum_{j=1}^{N} \theta_j - 2q\pi\right) \prod_{1 \le j < k \le N} ((e^{i\theta_j} - e^{i\theta_k})(e^{-i\theta_j} - e^{-i\theta_k})) \prod_{j=1}^{N} \frac{d\theta_j}{2\pi}.$$

The integration is performed from 0 to $2\pi$ for each $\theta_i$.

The exact values for the partition function are known [99]:

$$Z_{U(N)}(N, \lambda) = \text{Det} \left[ I_{j-k} \left( \frac{2N}{\lambda} \right) \right]_{j,k=1,...,N}, \qquad (A.3)$$

$$Z_{SU(N)}(N, \lambda) = \sum_{p=-\infty}^{\infty} \text{Det} \left[ I_{j-k+p} \left( \frac{2N}{\lambda} \right) \right]_{j,k=1,...,N}. \qquad (A.4)$$

For the special case of SU(2), the infinite sum can be explicitly calculated:

$$Z_{SU(2)}(2, \lambda) = \frac{1}{4} \lambda I_1 \left( \frac{8}{\lambda} \right), \qquad (A.5)$$

where $I_n$ represents the modified Bessel function of the first kind, as defined, for example, in the BesselI function in Mathematica.

The single-trace moments can be derived by taking derivatives of the partition function. It can be directly verified that both $\langle \text{Tr}U \rangle_{SU(N)}$ and $\langle \text{Tr}U \rangle_{U(N)}$ converge to the known analytic result in the 't Hooft limit [66, 67]:

$$\langle \text{Tr}U \rangle_{SU(\infty)} = \begin{cases} 1 - \frac{\lambda}{4}, & \text{for } 0 \le \lambda \le 2 \\ \frac{1}{\lambda}, & \text{for } \lambda \ge 2 \end{cases}. \qquad (A.5)$$

For various values of $N$, $\langle \text{Tr}U \rangle_{SU(N)}$ monotonically decreases, whereas $\langle \text{Tr}U \rangle_{U(N)}$ monotonically increases.
## B Loop equations for SU(N)/U(N) one-matrix model

In this appendix, we provide the derivation of the multi-trace loop equations for the following model:

$$Z = \int dU e^{\frac{N}{\lambda}(\mathrm{tr}U^{\dagger} + \mathrm{tr}U)} \qquad (B.1)$$

This framework is extended to accommodate more general cases of SU(N)/U(N) matrix model.

The derivation of loop equations for this model is crucial for understanding the dynamics within our bootstrap approach and can be represented as follows:

$$\int \mathcal{D}U \, \delta_{\epsilon} \left( U_{ab}^{n} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} e^{\frac{N}{\lambda}(\mathrm{tr}U^{\dagger} + \mathrm{tr}U)} \right) = 0 \qquad (B.2)$$

The loop equations are derived using the concept of small left variations on the group, denoted as $\delta_{\epsilon}$. This variation is defined as follows:

$$\delta_{\epsilon}(U_{ab}) = U_{ab} + i\epsilon_{ac}U_{cb}, \qquad \delta_{\epsilon}(U_{ab}^{\dagger}) = U_{ab}^{\dagger} - iU_{ac}^{\dagger}\epsilon_{cb} \qquad (B.3)$$

Here, $\epsilon$ represents an arbitrary traceless Hermitian matrix in the context of SU(N) gauge theory, and an arbitrary Hermitian matrix for U(N) scenarios. This distinction ensures that the variations are appropriately aligned with the structural properties of each group:

$$\epsilon_{cd} \left\langle \sum_{i=1}^{N-2} \frac{m_i}{N} U_{ab}^n U_{dc}^{m_i} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{i-1}} \mathrm{Tr} U^{m_{i+1}} \dots \mathrm{Tr} U^{m_{N-2}} \right. \\ \left. + \sum_{i=0}^{n-1} U_{ac}^i U_{db}^{n-i} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} + \frac{N}{\lambda} U_{ab}^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} (U_{dc} - U_{dc}^{\dagger}) \right\rangle = 0 \qquad (B.4)$$

We observe that if $\epsilon_{cd}X_{dc} = 0$ for any $\epsilon$, then it implies that $X_{dc} = \frac{s_G}{N}\delta_{dc}X_{ff}$. Here, $s_G$ is a scalar factor specific to the group in question and is defined in (2.25) in the main text. This relationship can be expressed as follows:

$$\begin{aligned} & \left\langle \sum_{i=1}^{N-2} \frac{m_i}{N} U_{ab}^n U_{dc}^{m_i} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{i-1}} \mathrm{Tr} U^{m_{i+1}} \dots \mathrm{Tr} U^{m_{N-2}} \right. \\ & \left. + \sum_{i=0}^{n-1} U_{ac}^i U_{db}^{n-i} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} + \frac{N}{\lambda} U_{ab}^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} (U_{dc} - U_{dc}^{\dagger}) \right\rangle \\ &= s_G \delta_{dc} \left\langle \frac{1}{N} n U_{ab}^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} + \frac{N}{\lambda} U_{ab}^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} (\mathrm{Tr} U - \mathrm{Tr} U^{\dagger}) \right. \\ & \left. + \left( \sum_{i=1}^{N-2} m_i \right) U_{ab}^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} \right\rangle \end{aligned} \qquad (B.5)$$
<figure><img src="image_18.png" /><figcaption>Figure 13. SU(2) bounds and convergences.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 14. U(2) bounds and convergences.</figcaption></figure>

Contracting both side by $\frac{1}{N^2}\delta_{ac}\delta_{bd}$, we get (for $n > 0$):

$$ 
\begin{aligned}
& \left\langle \sum_{i=1}^{N-2} \frac{m_i}{N^2} \mathrm{Tr} U^{n+m_i} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{i-1}} \mathrm{Tr} U^{m_{i+1}} \dots \mathrm{Tr} U^{m_{N-2}} \right. \\
& \qquad \left. + \sum_{i=0}^{n-1} \mathrm{Tr} U^i \mathrm{Tr} U^{n-i} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} + \frac{1}{\lambda} \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} (\mathrm{Tr} U^{n+1} - \mathrm{Tr} U^{n-1}) \right\rangle \\
& = \frac{s_G}{N^2} \left\langle n \mathrm{Tr} U^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} + \frac{N^2}{\lambda} \mathrm{Tr} U^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} (\mathrm{Tr} U - \mathrm{Tr} U^\dagger) \right. \\
& \qquad \left. + \left( \sum_{i=1}^{N-2} m_i \right) \mathrm{Tr} U^n \mathrm{Tr} U^{m_1} \dots \mathrm{Tr} U^{m_{N-2}} \right\rangle
\end{aligned}
\quad (B.6) $$

## C More bootstrap of the one-matrix model

In this appendix, we present additional results from the bootstrap analysis of the one-matrix model (2.8). The corresponding figures, figure 13, figure 14, and figure 15, represent the SU(2), U(2), and U($\infty$) counterparts of figure 8 discussed in the main text. We offer the following remarks on these results:
<figure><img src="image_20.png" /><figcaption>Figure 15. SU($\infty$) bounds and convergences.</figcaption></figure>

• The structure of the loop equations reveals that all moments are functions of $\langle \mathrm{Tr}U \rangle$ for both SU(2) and U($\infty$).<sup>23</sup> In contrast, the U(2) model requires two variables to determine the rest: $\langle \mathrm{Tr}U \rangle$ and $\langle \mathrm{Tr}U \mathrm{Tr}U^\dagger \rangle$, which closely resembles the SU(3) scenario.

• The convergence behavior of the U(2) bootstrap closely parallels that of the SU(3) model discussed in the main text, exhibiting similar rates. However, the SU(2) bootstrap converges significantly faster. Notably, the U($\infty$) bootstrap shows intriguing characteristics: the upper and lower bounds converge at distinctly different rates, particularly in the weak coupling phase where $\lambda \le 2$. Additionally, there is a pronounced spike at the phase transition point $\lambda = 2$.

**Data Availability Statement.** This article has no associated data or the data will not be deposited.

**Code Availability Statement.** This article has no associated code or the code will not be deposited.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

**References**

[1] M. Creutz, *Monte Carlo Study of Quantized SU(2) Gauge Theory*, Phys. Rev. D **21** (1980) 2308 [INSPIRE].

[2] M. Creutz, *Asymptotic Freedom Scales*, Phys. Rev. Lett. **45** (1980) 313 [INSPIRE].

[3] A. Athenodorou and M. Teper, *SU(N) gauge theories in 3+1 dimensions: glueball spectrum, string tensions and topology*, JHEP **12** (2021) 082 [arXiv:2106.00364] [INSPIRE].

[4] P. Anderson and M. Kruczenski, *Loop equation in Lattice gauge theories and bootstrap methods*, EPJ Web Conf. **175** (2018) 11011 [INSPIRE].

<footnote n="23"><sup>23</sup>However, it is important to note that these are polynomials of $\langle \mathrm{Tr}U \rangle$ in the case of U($\infty$).</footnote>
[5] H.W. Lin, *Bootstraps to strings: solving random matrix models with positivity*, JHEP 06 (2020) 090 [arXiv:2002.08387] [INSPIRE].

[6] X. Han, S.A. Hartnoll and J. Kruthoff, *Bootstrapping Matrix Quantum Mechanics*, Phys. Rev. Lett. 125 (2020) 041601 [arXiv:2004.10212] [INSPIRE].

[7] V. Kazakov and Z. Zheng, *Analytic and numerical bootstrap for one-matrix model and “unsolvable” two-matrix model*, JHEP 06 (2022) 030 [arXiv:2108.04830] [INSPIRE].

[8] V. Kazakov and Z. Zheng, *Bootstrap for lattice Yang-Mills theory*, Phys. Rev. D 107 (2023) L051501 [arXiv:2203.11360] [INSPIRE].

[9] P.D. Anderson and M. Kruczenski, *Loop Equations and bootstrap methods in the lattice*, Nucl. Phys. B 921 (2017) 702 [arXiv:1612.08140] [INSPIRE].

[10] M. Cho et al., *Bootstrapping the Ising Model on the Lattice*, arXiv:2206.12538 [INSPIRE].

[11] A. Jevicki, O. Karim, J.P. Rodrigues and H. Levine, *Loop Space Hamiltonians and Numerical Methods for Large N Gauge Theories*, Nucl. Phys. B 213 (1983) 169 [INSPIRE].

[12] A. Jevicki, O. Karim, J.P. Rodrigues and H. Levine, *Loop Space Hamiltonians and Numerical Methods for Large N Gauge Theories. 2*, Nucl. Phys. B 230 (1984) 299 [INSPIRE].

[13] R.M. Koch et al., *Large N optimization for multi-matrix systems*, JHEP 01 (2022) 168 [arXiv:2108.08803] [INSPIRE].

[14] H.W. Lin, *Bootstrap bounds on D0-brane quantum mechanics*, JHEP 06 (2023) 038 [arXiv:2302.04416] [INSPIRE].

[15] K. Mathaba, M. Mulokwe and J.P. Rodrigues, *Large N master field optimization: the quantum mechanics of two Yang-Mills coupled matrices*, JHEP 02 (2024) 054 [arXiv:2306.00935] [INSPIRE].

[16] Y. Aikawa, T. Morita and K. Yoshimura, *Application of bootstrap to a θ term*, Phys. Rev. D 105 (2022) 085017 [arXiv:2109.02701] [INSPIRE].

[17] Y. Aikawa, T. Morita and K. Yoshimura, *Bootstrap method in harmonic oscillator*, Phys. Lett. B 833 (2022) 137305 [arXiv:2109.08033] [INSPIRE].

[18] D. Bai, *Bootstrapping the deuteron*, arXiv:2201.00551 [INSPIRE].

[19] D. Berenstein and G. Hulsey, *Bootstrapping Simple QM Systems*, arXiv:2108.08757 [INSPIRE].

[20] D. Berenstein and G. Hulsey, *Bootstrapping more QM systems*, J. Phys. A 55 (2022) 275304 [arXiv:2109.06251] [INSPIRE].

[21] D. Berenstein and G. Hulsey, *Semidefinite programming algorithm for the quantum mechanical bootstrap*, Phys. Rev. E 107 (2023) L053301 [arXiv:2209.14332] [INSPIRE].

[22] D. Berenstein and G. Hulsey, *Anomalous bootstrap on the half-line*, Phys. Rev. D 106 (2022) 045029 [arXiv:2206.01765] [INSPIRE].

[23] D. Berenstein and G. Hulsey, *One-dimensional reflection in the quantum mechanical bootstrap*, Phys. Rev. D 109 (2024) 025013 [arXiv:2307.11724] [INSPIRE].

[24] J. Bhattacharya et al., *Numerical bootstrap in quantum mechanics*, Phys. Lett. B 823 (2021) 136785 [arXiv:2108.11416] [INSPIRE].

[25] M.J. Blacker, A. Bhattacharyya and A. Banerjee, *Bootstrapping the Kronig-Penney model*, Phys. Rev. D 106 (2022) 116008 [arXiv:2209.09919] [INSPIRE].
[26] W. Ding, Variational Equations-of-States for Interacting Quantum Hamiltonians, arXiv:2307.00812 [inSPIRE].

[27] B.-N. Du, M.-X. Huang and P.-X. Zeng, Bootstrapping Calabi-Yau quantum mechanics, Commun. Theor. Phys. 74 (2022) 095801 [arXiv:2111.08442] [inSPIRE].

[28] J. Eisert, A note on lower bounds to variational problems with guarantees, arXiv:2301.06142 [inSPIRE].

[29] H. Fawzi, O. Fawzi and S.O. Scalet, Entropy constraints for ground energy optimization, J. Math. Phys. 65 (2024) 032201 [arXiv:2305.06855] [inSPIRE].

[30] X. Han, Quantum Many-body Bootstrap, arXiv:2006.06002.

[31] M.B. Hastings and R. O'Donnell, Optimizing strongly interacting fermionic Hamiltonians, in the proceedings of the 54th Annual ACM Symposium on Theory of Computing, Rome, Italy, June 20-24 (2022) [DOI:10.1145/3519935.3519960] [arXiv:2110.10701] [inSPIRE].

[32] M.B. Hastings, Perturbation Theory and the Sum of Squares, arXiv:2205.12325 [inSPIRE].

[33] H. Hessam, M. Khalkhali and N. Pagliaroli, Bootstrapping Dirac ensembles, J. Phys. A 55 (2022) 335204 [arXiv:2107.10333] [inSPIRE].

[34] X. Hu, Different Bootstrap Matrices in Many QM Systems, arXiv:2206.00767 [inSPIRE].

[35] S. Khan, Y. Agarwal, D. Tripathy and S. Jain, Bootstrapping PT symmetric quantum mechanics, Phys. Lett. B 834 (2022) 137445 [arXiv:2202.05351] [inSPIRE].

[36] I. Kull, N. Schuch, B. Dive and M. Navascués, Lower Bounds on Ground-State Energies of Local Hamiltonians through the Renormalization Group, Phys. Rev. X 14 (2024) 021008 [arXiv:2212.03014] [inSPIRE].

[37] W. Li, Null bootstrap for non-Hermitian Hamiltonians, Phys. Rev. D 106 (2022) 125021 [arXiv:2202.04334] [inSPIRE].

[38] W. Li, Taming Dyson-Schwinger Equations with Null States, Phys. Rev. Lett. 131 (2023) 031603 [arXiv:2303.10978] [inSPIRE].

[39] Y. Guo and W. Li, Solving anharmonic oscillator with null states: Hamiltonian bootstrap and Dyson-Schwinger equations, Phys. Rev. D 108 (2023) 125002 [arXiv:2305.15992] [inSPIRE].

[40] T. Morita, Universal bounds on quantum mechanics through energy conservation and the bootstrap method, PTEP 2023 (2023) 023A01 [arXiv:2208.09370] [inSPIRE].

[41] Y. Nakayama, Bootstrapping microcanonical ensemble in classical system, Mod. Phys. Lett. A 37 (2022) 2250054 [arXiv:2201.04316] [inSPIRE].

[42] C.O. Nancarrow and Y. Xin, Bootstrapping the gap in quantum spin systems, JHEP 08 (2023) 052 [arXiv:2211.03819] [inSPIRE].

[43] A. Tavakoli, A. Pozas-Kerstjens, P. Brown and M. Araújo, Semidefinite programming relaxations for quantum correlations, Rev. Mod. Phys. 96 (2024) 045006 [arXiv:2307.02551] [inSPIRE].

[44] S. Tchoumakov and S. Florens, Bootstrapping Bloch bands, J. Phys. A 55 (2022) 015203 [arXiv:2109.06600] [inSPIRE].

[45] W. Fan and H. Zhang, Non-perturbative instanton effects in the quartic and the sextic double-well potential by the numerical bootstrap approach, arXiv:2308.11516 [inSPIRE].

[46] W. Fan, H. Zhang and Z. Li, Unify the Effect of Anharmonicity in Double-Wells and Anharmonic Oscillators, Int. J. Theor. Phys. 63 (2024) 266 [arXiv:2309.09269] [inSPIRE].
[47] R.R. John and K.P. R, *Anharmonic oscillators and the null bootstrap*, arXiv:2309.06381 [INSPIRE].

[48] W. Li, *Principle of minimal singularity for Green's functions*, Phys. Rev. D 109 (2024) 045012 [arXiv:2309.02201] [INSPIRE].

[49] M. Zeng, *Feynman integrals from positivity constraints*, JHEP 09 (2023) 042 [arXiv:2303.15624] [INSPIRE].

[50] H. Fawzi, O. Fawzi and S.O. Scalet, *Certified algorithms for equilibrium states of local quantum Hamiltonians*, Nature Commun. 15 (2024) 7394 [arXiv:2311.18706] [INSPIRE].

[51] W. Li, *φn trajectory bootstrap*, Phys. Rev. D 111 (2025) 045013 [arXiv:2402.05778] [INSPIRE].

[52] D. Goluskin, *Bounding Averages Rigorously Using Semidefinite Programming: Mean Moments of the Lorenz System*, Journal of Nonlinear Science 28 (2017) 621 [arXiv:1610.05335].

[53] D. Goluskin, *Bounding extrema over global attractors using polynomial optimisation*, Nonlinearity 33 (2020) 4878 [arXiv:1807.09814].

[54] I. Tobasco, D. Goluskin and C.R. Doering, *Optimal bounds and extremal trajectories for time averages in nonlinear dynamical systems*, Phys. Lett. A 382 (2018) 382 [arXiv:1705.07096].

[55] M. Cho, *Bootstrapping the Stochastic Resonance*, arXiv:2309.06464 [INSPIRE].

[56] Y.M. Makeenko and A.A. Migdal, *Exact Equation for the Loop Average in Multicolor QCD*, Phys. Lett. B 88 (1979) 135 [Erratum ibid. 89 (1980) 437] [INSPIRE].

[57] R. Gambini and J. Griego, *A geometric approach to the Makeenko-Migdal equations*, Phys. Lett. B 256 (1991) 437 [INSPIRE].

[58] A. Denbleyker et al., *Series expansions of the density of states in SU(2) lattice gauge theory*, Phys. Rev. D 78 (2008) 054503 [arXiv:0807.0185] [INSPIRE].

[59] A. Athenodorou and M. Teper, *SU(N) gauge theories in 2+1 dimensions: glueball spectra and k-string tensions*, JHEP 02 (2017) 015 [arXiv:1609.03873] [INSPIRE].

[60] Z. Li and S. Zhou, *Bootstrapping the Abelian lattice gauge theories*, JHEP 08 (2024) 154 [arXiv:2404.17071] [INSPIRE].

[61] A.A. Migdal, *Loop Equations and 1/N Expansion*, Phys. Rept. 102 (1983) 199 [INSPIRE].

[62] R. Gambini and A. Trias, *Gauge Dynamics in the C Representation*, Nucl. Phys. B 278 (1986) 436 [INSPIRE].

[63] R. Dempsey, I.R. Klebanov, L.L. Lin and S.S. Pufu, *Adjoint Majorana QCD$_2$ at finite N*, JHEP 04 (2023) 107 [arXiv:2210.10895] [INSPIRE].

[64] R. Dempsey, I.R. Klebanov, S.S. Pufu and B.T. Søgaard, *Lattice Hamiltonian for adjoint QCD$_2$, JHEP 08 (2024) 009 [arXiv:2311.09334] [INSPIRE].

[65] Wikipedia contributors, *Levi-civita symbol*, Wikipedia, the free encyclopedia, https://en.wikipedia.org/w/index.php?title=Levi-Civita_symbol&oldid=1243909782.

[66] D.J. Gross and E. Witten, *Possible Third Order Phase Transition in the Large N Lattice Gauge Theory*, Phys. Rev. D 21 (1980) 446 [INSPIRE].

[67] S.R. Wadia, *A Study of U(N) Lattice Gauge Theory in 2-dimensions*, arXiv:1212.2906 [INSPIRE].

[68] V.A. Kazakov, *Wilson loop average for an arbitrary contour in two-dimensional U(N) gauge theory*, Nucl. Phys. B 179 (1981) 283 [INSPIRE].
[69] V.A. Kazakov and I.K. Kostov, *Nonlinear strings in two-dimensional U(infinity) gauge theory*, Nucl. Phys. B 176 (1980) 199 [INSPIRE].

[70] V.A. Kazakov and I.K. Kostov, *Computation of the Wilson Loop Functional in Two-dimensional U(infinite) Lattice Gauge Theory*, Phys. Lett. B 105 (1981) 453 [INSPIRE].

[71] Y. Makeenko, *Methods of Contemporary Gauge Theory*, Cambridge University Press (2002) [DOI:10.1017/9781009402095] [INSPIRE].

[72] A. Migdal, *Statistical equilibrium of circulating fluids*, Phys. Rept. 1011 (2023) 1 [arXiv:2209.12312] [INSPIRE].

[73] A. Migdal, *To the theory of decaying turbulence*, Fractal Fract. 7 (2023) 754 [arXiv:2304.13719] [INSPIRE].

[74] M. Cho and X. Sun, *Bootstrap, Markov Chain Monte Carlo, and LP/SDP hierarchy for the lattice Ising model*, JHEP 11 (2023) 047 [arXiv:2309.01016] [INSPIRE].

[75] K. Osterwalder and R. Schrader, *Axioms for Euclidean Green's functions*, Commun. Math. Phys. 31 (1973) 83 [INSPIRE].

[76] K. Osterwalder and R. Schrader, *Axioms for Euclidean Green's functions*. 2, Commun. Math. Phys. 42 (1975) 281 [INSPIRE].

[77] A.A. Belavin, A.M. Polyakov and A.B. Zamolodchikov, *Infinite Conformal Symmetry in Two-Dimensional Quantum Field Theory*, Nucl. Phys. B 241 (1984) 333 [INSPIRE].

[78] R. Rattazzi, V.S. Rychkov, E. Tonni and A. Vichi, *Bounding scalar operator dimensions in 4D CFT*, JHEP 12 (2008) 031 [arXiv:0807.0004] [INSPIRE].

[79] D. Poland, S. Rychkov and A. Vichi, *The Conformal Bootstrap: Theory, Numerical Techniques, and Applications*, Rev. Mod. Phys. 91 (2019) 015002 [arXiv:1805.04405] [INSPIRE].

[80] S. Rychkov and N. Su, *New developments in the numerical conformal bootstrap*, Rev. Mod. Phys. 96 (2024) 045004 [arXiv:2311.15844] [INSPIRE].

[81] J. Frohlich, B. Simon and T. Spencer, *Infrared Bounds, Phase Transitions and Continuous Symmetry Breaking*, Commun. Math. Phys. 50 (1976) 79 [INSPIRE].

[82] K. Osterwalder and E. Seiler, *Gauge Field Theories on the Lattice*, Annals Phys. 110 (1978) 440 [INSPIRE].

[83] R. Rattazzi, S. Rychkov and A. Vichi, *Bounds in 4D Conformal Field Theories with Global Symmetry*, J. Phys. A 44 (2011) 035402 [arXiv:1009.5985] [INSPIRE].

[84] Y. He, A. Irrgang and M. Kruczenski, *A note on the S-matrix bootstrap for the 2d O(N) bosonic model*, JHEP 11 (2018) 093 [arXiv:1805.02812] [INSPIRE].

[85] L. Córdova and P. Vieira, *Adding flavour to the S-matrix bootstrap*, JHEP 12 (2018) 063 [arXiv:1805.11143] [INSPIRE].

[86] M.F. Paulos and Z. Zheng, *Bounding scattering of charged particles in 1 + 1 dimensions*, JHEP 05 (2020) 145 [arXiv:1805.11429] [INSPIRE].

[87] C. Bachoc, D.C. Gijswijt, A. Schrijver and F. Vallentin, *Invariant semidefinite programs*, arXiv:1007.2905 [DOI:10.1007/978-1-4614-0769-0_9].

[88] The GAP Group, *GAP — Groups, Algorithms, and Programming, Version 4.14.0*, (2024), https://www.gap-system.org.

[89] N. Xu and P. Doerschuk, *Computation of real-valued basis functions which transform as irreducible representations of the polyhedral groups*, SIAM J. Sci. Comput. 43 (2021) A3657.
[90] J.-P. Serre, *Linear representations of finite groups*, vol. 42, Springer (1977).

[91] H.W. Lin and Z. Zheng, *Bootstrapping ground state correlators in matrix theory. Part I*, JHEP 01 (2025) 190 [arXiv:2410.14647] [INSPIRE].

[92] C.A. Doran, L.A. Pando Zayas, V.G.J. Rodgers and K. Stiffler, *Tensions and Luscher Terms for (2+1)-dimensional k-strings from Holographic Models*, JHEP 11 (2009) 064 [arXiv:0907.1331] [INSPIRE].

[93] S.R. White, *Density matrix formulation for quantum renormalization groups*, Phys. Rev. Lett. 69 (1992) 2863 [INSPIRE].

[94] M. Cho et al., *Coarse-grained Bootstrap of Quantum Many-body Systems*, arXiv:2412.07837 [INSPIRE].

[95] K. Ghosh and Z. Zheng, *Numerical conformal bootstrap with analytic functionals and outer approximation*, JHEP 09 (2024) 143 [arXiv:2307.11144] [INSPIRE].

[96] J. Hoppe, *Classical dynamics of SU(2) matrix models*, arXiv:2303.12169 [INSPIRE].

[97] S. Komatsu et al., *Gravity from quantum mechanics of finite matrices*, arXiv:2401.16471 [INSPIRE].

[98] V.A. Kazakov and A.A. Migdal, *Baryons as Solitons in Loop Space*, Phys. Lett. B 103 (1981) 214 [INSPIRE].

[99] R.G. Jha, *Finite N unitary matrix model*, arXiv:2003.00341 [INSPIRE].