# Analytic and numerical bootstrap for one-matrix model and “unsolvable” two-matrix model

Vladimir Kazakov, Zechuan Zheng

► To cite this version:

Vladimir Kazakov, Zechuan Zheng. Analytic and numerical bootstrap for one-matrix model and “unsolvable” two-matrix model. Journal of High Energy Physics, 2022, 2022 (6), pp.30. 10.1007/JHEP06(2022)030. hal-03324265

HAL Id: hal-03324265

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

Submitted on 12 Jul 2024

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
RECEIVED: March 30, 2022

ACCEPTED: May 10, 2022

PUBLISHED: June 7, 2022

# Analytic and numerical bootstrap for one-matrix model and “unsolvable” two-matrix model

Vladimir Kazakov and Zechuan Zheng

Laboratoire de Physique de l'École Normale Supérieure, CNRS,
Université PSL, Sorbonne Universités,
24 rue Lhomond, 75005 Paris, France

E-mail: kazakov@lpt.ens.fr, zechuan.zheng@phys.ens.fr

**ABSTRACT:** We propose the *relaxation* bootstrap method for the numerical solution of multi-matrix models in the large *N* limit, developing and improving the recent proposal of H. Lin. It gives rigorous inequalities on the single trace moments of the matrices up to a given “cutoff” order (length) of the moments. The method combines usual loop equations on the moments and the positivity constraint on the correlation matrix of the moments. We have a rigorous proof of applicability of this method in the case of the one-matrix model where the condition of positivity of the saddle point solution appears to be equivalent to the presence of supports of the eigenvalue distribution only on the real axis and only with positive weight. We demonstrate the numerical efficiency of our method by solving the analytically “unsolvable” two-matrix model with tr[A, B]² interaction and quartic potentials, even for solutions with spontaneously broken discrete symmetry. The region of values for computed moments allowed by inequalities quickly shrinks with the increase of the cutoff, allowing the precision of about 6 digits for generic values of couplings in the case of Z₂ symmetric solutions. Our numerical data are checked against the known analytic results for particular values of parameters.

**KEYWORDS:** 1/*N* Expansion, Matrix Models

ARXIV ePRINT: 2108.04830
# Contents

1 Introduction 1
1.1 Main results 4

2 Hermitian one-matrix model bootstrap 5
2.1 Hermitian one-matrix model in the planar limit and loop equations 5
2.2 Bootstrap method for the large $N$ one matrix model 9

3 Hamburger moment problem and positivity of resolvent 11
3.1 Hamburger problem versus the positivity condition on resolvent 11
3.2 Classification of physical solution of quartic one-matrix model 13
3.3 Comments 15

4 Relaxation bootstrap method 17
4.1 Physical constraints 18
4.2 Relaxation matrix 20

5 Bootstrap for “unsolvable” two-matrix model with Tr[$A, B$]$^2$ interaction 23
5.1 Bootstrap solution for a generic choice of $g, h$ 24
5.2 Demonstration for analytically solvable cases 26
5.3 Phase diagram and convergence rate 29
5.3.1 Rate of convergence 31
5.4 Bootstrapping the symmetry breaking solution 32
5.4.1 Schemes for symmetry breaking bootstrap 33
5.4.2 Numerical results for symmetry breaking solution 34

6 Conclusion and discussion 37

A Analytic solvability of two-matrix model with cubic interactions and arbitrary potentials 39

B Solving the positivity condition of the resolvent 41
B.1 Cuts and zeros 41
B.2 A working example 42

C Dual formulation and relaxation 44
C.1 Dual problem of general optimization problem 45
C.2 Relaxation problem and dual problem 46

D Implementation in SDPA 47
D.1 Choice of SDP solvers 48
D.2 Generation of the loop equations 48
D.3 Solving the loop equations 49
(1/N-expansion appears to be a topological expansion: the “fat” graphs of a given genus g are weighted with the factor $N^{2-2g}$). This allows the counting of such planar graphs [14–16] and enables the introduction and exact solution of statistical mechanical models on random planar dynamical lattices — Ising model on random triangulations [17, 18] and various generalizations [6, 19–21].

The direct analytic computation of a majority of such multi-matrix integrals is virtually impossible, apart from some trivial, albeit important, cases, such as the quadratic potential $V(x, y, z, \dots)$ leading to the gaussian integral.<fnref n="1" />

For a sub-class of such integrals with specific potentials the problem can be reduced to integrations over a smaller number of variables than ~ $N^2$. For example, sometimes the problem can be reduced to the integrations or summation only over ~ $N$ variables, such as eigenvalues of the matrices. Then, in the large $N$ limit, the problem can be reduced to the saddle point calculation, significantly simplifying the problem of computation of that functional integral.<fnref n="2" /> The basic example of such a simplification is the one matrix model

$$Z = \int d^{N^2} A e^{-\text{tr}\mathcal{V}(A)} \qquad (1.3)$$

solvable for any potential $\mathcal{V}(x)$. Once we have two or more matrix integration variables in (1.1) the problem usually gets much more complicated. Generically, such models are unsolvable, i.e. the number of degrees of freedom cannot be efficiently reduced, and the saddle point approximation is inappropriate since the characteristic “energy” and entropy of the integration variables are both of the order ~ $N^2$. Here comes the question whether we can study these integrals at least numerically.

Virtually the only universal general method of numerical computation of functional integrals is the Monte-Carlo method. It has been applied to some matrix integrals with more or less of success. Its main drawbacks are well known: i) the result comes with a statistical error; ii) it is sometimes difficult to reach the numerical equilibrium state in a reasonable time; iii) MC is bad for the systems with sign-changing Boltzmann weights or non-local interactions; iv) the size of the system (the number of integrals) is limited by computational facilities v) The precision is usually rather modest, maximum about 3–4 digits.

Do we have any alternative?

In the 1980’s, in a series of papers [23–25], the authors formulated the problem of large N matrix integral and large N quantum mechanics in the loop space (space of moments). The authors attempted the numerical study for the loop variables by minimizing an effective action. They were the first to stress the importance of positive semi-definiteness conditions for certain matrices of loop variables in getting physically meaningful results.

Recently, an important progress has been made in the computations of multi-point correlators in conformal field theories in various dimensions, due to the conformal bootstrap

<footnote n="1">¹Expansions w.r.t. parameters around the gaussian point lead in the ‘t Hooft limit to the perturbation theory formulated in terms of planar Feynman graphs. It can help to study the model in a specific, narrow domain of the parameter space.</footnote>

<footnote n="2">²An (incomplete) review of such *solvable* matrix models can be found in [22].</footnote>
method [26]. The method uses various properties of correlators, such as crossing and positivity, to “bootstrap” numerically their values and the values of the critical exponents. It appeared to be far more efficient and precise then other numerical approaches, giving the critical exponents of 3d Ising model with the record 6-digits precision [27]. An appealing property of this method is the absence of any statistical error in the results, which are given within rigorously established margins.

Inspired by this success a few authors applied the philosophy of the numerical bootstrap to the computations of various matrix integrals [28, 29] and even of the lattice multi-color QCD and $N = 4$ SYM theory [30]. Instead of the direct study of the matrix integrals they proposed to study the large $N$ Schwinger-Dyson equations which are often also called loop equations, in analogy with their applications to QCD [31]. They are easily obtained by the obvious Ward identities resulting from insertion of the full matrix derivative under the matrix integral:

$$0 = \int d^{N^2} A d^{N^2} B d^{N^2} C \dots \frac{\mathrm{tr}}{N} \left( \frac{\partial}{\partial A} A^m B^n C^k \dots \right) e^{-\mathrm{tr}(A,B,C,\dots)} \quad (1.4)$$

where the matrix derivative inside the trace $\frac{\partial}{\partial A}$ acts on all $A$-matrices, including the potential. All other loop equations correspond to all possible “words” of matrices under the trace and to all insertions of various matrix derivatives at any place in the “words”.<fnref n="3" /> Then the positivity conditions are imposed stating that the inner product<fnref n="4" /> of any operator with itself is positive. Rigorous bounds on the dynamical quantities of the theory can be derived from these positivity conditions and loop equations.

This new approach, compared to the previous work of loop variables [23–25], imposes the large $N$ Schwinger-Dyson equations (loop equations) explicitly, rather than getting loop equations as a result of effective action minimization. In parallel with the philosophy of conformal bootstrap, this approach focuses more on the geometry of the space of loop variables under sensible physical constraints, which guarantees the rigorousness of the bounds on physical quantities.

In the inspiring work of Lin [28] the method was rather successfully applied to the one-matrix model mentioned above, to the exactly solvable two-matrix model with tr($AB$) interaction [15, 16] describing the Ising model on planar graphs [17] as well as to the model with tr($AB^2 + A^2B$) interaction, presented there as a case of “unsolvable” matrix model.<fnref n="5" /> Lin uses the non-linear equations (1.4) to bootstrap the loop averages up to the positive semi-definite matrix of size 45.

This new approach has, in our opinion, a great potential for the precision computations of physically important matrix integrals in the ‘t Hooft limit. But at the same time it is very much perfectible at this stage.

Firstly, the numerical matrix bootstrap approach of [28–30], based on the loop equations and positivity constraint, is not well understood analytically. Its efficiency, and the

<footnote n="3">³In the ‘t Hooft limit, the single trace “words” are enough due to the factorization property which we will describe in the next section.</footnote>

<footnote n="4">⁴Here inner product of an operator $O$ means $\langle \mathrm{tr} O^\dagger O \rangle$.</footnote>

<footnote n="5">⁵We will demonstrate in the appendix A that, in fact, all two-matrix models with cubic interactions, including this one, are solvable in the above-mentioned sense.</footnote>
power of positivity, still looks quite mysterious. It is not even fully understood why we need the positivity condition. Secondly, the matrix bootstrap has a very distinguished feature comparing to most of the other bootstrap problems we dealt with so far: it is in general non-convex. The non-convexity comes from the quadratic terms in the loop equation, which is a result of large $N$ factorization. In optimization theory, this is called Nonlinear SDP (semi-definite programming) and all the solvers for it are not mature enough compared with the highly developed SDP solvers dealing with linear problems. In [28–30], the authors tried to bootstrap the matrix models by the Nonlinear SDP directly, and this non-linearity limited the bootstrap capabilities to very simple models, or to more complex models but only up to very small lengths of operators.

## 1.1 Main results

In this work, we advance the matrix bootstrap approach trying, on the one hand, to understand analytically the role of positivity conditions, and on the other hand, to overcome, at least partially, the above-mentioned limitations of the method.

First, we derive a necessary and sufficient condition for the positivity of bootstrap for large $N$ one-matrix model, to clarify how this method is working. Namely, we show that the positivity is equivalent to the condition for the resolvent to have the cuts only on the real axis, with the positive imaginary part corresponding to the positive density of eigenvalue distribution. This condition actually enables us, in principle, to analytically solve the bootstrap problem for any one-matrix model. For the illustrative purposes, we will apply the new positivity condition to the one-matrix model with quartic potential:

$$V'(x) = \frac{\mu}{2}x^2 + \frac{g}{4}x^4, \quad (\mu = \pm 1) \qquad (1.5)$$

where we normalized the quadratic term to $\pm 1$. We will use the analytic bootstrap to completely classify the admissible set of solutions of the loop equations and positivity conditions, and to locate the critical value of $Z_2$ symmetry breaking.

So far, we could solve exactly a very limited set of bootstrap problems, and most of them correspond to very simple theories, such as Sine-Gordon theory in S-matrix bootstrap [32] and 1d mean field theory in conformal bootstrap [33]. Since this method appears to be applicable to any one-matrix model and generalized to some solvable multi-matrix models, it provides us with a big new family of exactly solvable bootstrap problems. Hopefully these solvable bootstrap models will give us more of intuition about the way the bootstrap method works.

The other new result of this work is a new bootstrap scheme for the study of non-linear SDP for multi-matrix integrals, which appears to be numerically much more efficient than those proposed in the past. The main ingredient of the method is the introduction of relaxation matrix in the place of non-linearity of the loop equation. Namely, we treat the quadratic terms as independent variables and impose the positivity condition on these variables. Surprisingly, it seems enough to bootstrap the region of admissible values of the computed quantity that is quickly shrinking with the increase of the “cutoff” — the maximal length of “words” in the involved operators.
As a particular example of analytically unsolvable matrix integral we will study by this method the following two-matrix model

$$Z = \lim_{N \to \infty} \int d^{N^2} A d^{N^2} B e^{-N \text{tr}(-h[A,B]^2/2 + A^2/2 + gA^4/4 + B^2/2 + gB^4/4)} \quad (1.6)$$

Various versions of this model have been studied in the past in connection to certain $N = 1$ supersymmetric Yang-Mills theories [34]. In the particular case $g = 0$ the model is solvable and it will serve us as an important check of applicability and efficiency of our relaxation bootstrap method. Our results show a very good precision: up to 6 digits with the maximal cutoff equal to 22 for the words under averages. We were also able to establish with a reasonable accuracy the phase structure of the model in the $g, h$ coupling space, i.e. the positions of critical lines corresponding to the convergence radius of planar perturbative expansion, as well as to the spontaneous $Z_2$ symmetry breaking.

The two-matrix model (1.6) considered in this paper serves mostly for the illustration of the power of our method, though it could have in principal some physical applications, such as the statistical mechanics on dynamical planar graphs, in the spirit of [35–37].

This article is organized as follows. The next section 2 serves as a retrospect of the Hermitian matrix integral and the numerical bootstrap technique developed for it so far. Then in section 3 we propose our equivalent condition for the positivity condition described in section 2. This condition will justify the numerical bootstrap method and enable us to analytically solve the corresponding bootstrap problem. in section 4, we will describe the way our relaxation method works for analytically unsolvable large $N$ multi-matrix integrals. We test this relaxation method in section 5 on the concrete unsolvable model (1.6). We will see that our relaxation method is able to largely meet our expectations, with remarkable precision. In the last section, after short conclusions, we will briefly discuss possible applications of our method to some more physical problems, such as the multicolor lattice Yang-Mills theory.

**Note:** the main results of this work are compared with the later Monte Carlo (MC) results [38]. This comparison convinces us that the bootstrap method is more efficient than MC regarding the large $N$ two-matrix model calculation.

## 2 Hermitian one-matrix model bootstrap

In this section we will revisit several basic facts about large $N$ limit Hermitian one-matrix model and the related numerical bootstrap proposed in [28]. We will be mainly focused here on the aspects of this model which are crucial for the theoretical development in the next section and provide us with important intuition. The reader can refer to numerous works and reviews, some already cited above (see e.g. [39] for a good state-of-art description of results on Hermitian one-matrix model).

### 2.1 Hermitian one-matrix model in the planar limit and loop equations

The Hermitian one-matrix model is defined by matrix integral:

$$Z_N = \int d^{N^2} M e^{-N \text{tr} V(M)} \quad (2.1)$$
where the invariant Hermitian measure is $d^{N^2}M = \prod_{i,j=1}^{N} dM_{ij}$. The potential is usually taken polynomial:<fnref n="6" />

$$V(x) = \sum_{k=2}^{d+1} \frac{g_k}{k} M^k. \qquad (2.2)$$

The main “physical observable” is the k-th moment:

$$W_k = \langle \mathrm{Tr} M^k \rangle = \int \frac{d^{N^2} M}{Z_N} \frac{1}{N} \mathrm{tr} M^k e^{-N \mathrm{tr} V(M)}. \qquad (2.3)$$

This model is solvable in the planar limit for arbitrary polynomial potentials [14]. There exist several methods for that: direct recursion relations for planar graphs, orthogonal polynomials, saddle point approximation for the eigenvalue distribution, loop equations (see [2, 40] for a review). The loop equations will play the crucial role in our bootstrap method.

To derive them we simply use the Schwinger-Dyson method by writing

$$0 = \int d^{N^2} M \frac{\mathrm{tr}}{N} \left( \frac{\partial}{\partial M} M^k \right) e^{-N \mathrm{tr} V(M)} \qquad (2.4)$$

since the expression under the integral is a total derivative. The boundary terms are absent assuming that the highest power $d+1$ of the potential is even and its coefficient is positive $t_{d+1} > 0$.<fnref n="7" />

Applying explicitly the matrix derivative in (2.4) we write the loop equation in terms of the moments:<fnref n="8" />

$$\langle \mathrm{Tr} V'(M) M^k \rangle = \sum_{l=0}^{k-1} \langle \mathrm{Tr} M^l \mathrm{Tr} M^{k-l-1} \rangle. \qquad (2.5)$$

In the $N \to \infty$ limit we can use the factorization property:

$$\langle \mathrm{Tr} M^l \mathrm{Tr} M^m \rangle = \langle \mathrm{Tr} M^l \rangle \langle \mathrm{Tr} M^m \rangle + \mathcal{O}(1/N^2). \qquad (2.6)$$

Then the loop equation reduces to

$$\sum_{j=1}^{d} g_j W_{k+j} = \sum_{l=0}^{k-1} W_l W_{k-l+1}. \qquad (2.7)$$

The simplest way to solve (2.7) is to introduce the generating function of moments — the resolvent — as a formal power series in terms of $z^{-1}$:

$$G(z) = \sum_{k=0}^{\infty} z^{-k-1} W_k. \qquad (2.8)$$

<footnote n="6">⁶We believe that our final conclusion can be generalized to non-polynomial potentials, but there may be some subtleties.</footnote>

<footnote n="7">⁷For the “unstable” potentials, which do not satisfy one of these conditions, the matrix integral might still exist with appropriate deformation of the integration contour. The large $N$ solutions can exist even independently of the contour deformation since they correspond to local minima of the effective potential for the eigenvalues.</footnote>

<footnote n="8">⁸Here for the conciseness, we introduce the normalized trace $\mathrm{Tr} = \frac{1}{N} \mathrm{tr}$, so that $\mathrm{Tr} I = 1$.</footnote>
We have not yet assumed anything about the convergence of the series. Multiplying (2.7) by $z^{-k}$ and summing from $k=1$ to $\infty$ we represent the loop equation in a compact form as a quadratic equation for the resolvent.

$$G(z)^2 + P(z) = V'(z)G(z). \qquad (2.9)$$

The function $P(z)$ comes from carefully collecting in the summation the terms with small $k$'s. It can be written compactly as:

$$P(z) = \langle \text{tr} \frac{V'(z) - V'(M)}{z - M} \rangle. \qquad (2.10)$$

This is a polynomial of $z$ and a linear function of $W_k$, ($k=1, \dots, d-1$). For example, if $V'(z) = z + gz^3$, then we have $P(z) = 1 + gz^2 + gW_2 + gzW_1$. We can solve (2.9), picking the relevant branch of the root which reproduces the leading $z^{-1}$ behavior of the resolvent $G(z) = \frac{1}{z} + \mathcal{O}(\frac{1}{z^2})$ at infinity:

$$G(z) = \frac{1}{2} \left( V'(z) - \sqrt{V'(z)^2 - 4P(z)} \right). \qquad (2.11)$$

This result will play an important role in our work, so we make several comments on it:

1. By (2.11), the resolvent is understood as a genuine analytic function, at least in the neighborhood of infinity point. So the formal series defined in (2.8) has a finite radius of convergence. As a consequence, there must be an exponential bound for the moments:<fnref n="9" />

$$W_n \le C_0 R(g_k, W_1, \dots, W_{d-1})^n. \qquad (2.12)$$

2. It is clear from this formula, that all the moments are determined by several low-order moments in $P(z)$ and the couplings:

$$W_n = W_n(g_k, W_1, \dots, W_{d-1}). \qquad (2.13)$$

This is what the loop equation tells us. But the loop equation doesn't tells us how to fix the low-order moments involved in $P(z)$. Those can be fixed only by additional assumptions on the solution, such as for example the single support solution for the eigenvalues (single cut on the physical sheet of $G(z)$). We will see how to classify the solutions which are picked up by the bootstrap method.

To have more intuitive ideas of possible large $N$ solutions it is useful to reduce the matrix integration (2.1) to the integration over the eigenvalues of the Hermitian matrix. Namely, if we represent it as $M = \Omega^\dagger X \Omega$ where $X = \text{diag}(x_1, x_2, \dots, x_N)$ is the diagonal matrix of eigenvalues and $\Omega$ is the diagonalizing unitary matrix, the matrix integral reduces to only $N$ integrations over the eigenvalues [14]:

$$Z_N = \int \prod_{j=1}^{N} (dx_j e^{-NV(x_j)}) \Delta^2(x_1, \dots, z_N) \qquad (2.14)$$

<footnote n="9">⁹Strictly speaking the radius of convergence is the inverse of the module of largest root of the polynomial under the square root in (2.11), unless two of such roots merge. Here it is enough for us that it is bounded exponentially.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Coulomb gas interpretation for the eigenvalue configurations and a typical cut configuration. For a general solution of the SDP (2.15), we can have a complex cut at the maximum of the potential.</figcaption></figure>

where the square of the Vandermonde determinant $\Delta(x_1, \dots, z_N) = \prod_{i>j} (x_i - x_j)$ represents the Jacobian of the change of integration variables (Dyson measure). Here the integrand is of the order $e^{N^2(\dots)}$ whereas the number of variables is reduced to $N$. This allows for the application of the saddle point approximation, giving the BIPZ saddle point equations (SPE) [14]

$$V'(x_j) = \sum_{k(\neq j)} \frac{2}{x_j - x_k}, \quad j = 1, 2, \dots, N. \qquad (2.15)$$

It looks as the condition of electrostatic equilibrium of two-dimensional point-like electric charges (of the same sign) with coordinates $x_j$ on a line, locked in the potential $V(x)$ (see the figure 1).

The possible physical solutions correspond to filling $n$ minima of such a potential with fractions $N_1, N_2, \dots, N_n$ of these charges, such that $\sum_{l=1}^n N_l = N$. The eigenvalues then form a continuous distribution with $n$ finite supports along the real axis.

The SPE describes all extrema of the effective potential

$$V_{\text{eff}} = N \sum_k V(x_k) + \log \Delta^2(x_1, \dots, z_N), \qquad (2.16)$$

not only the minima but also the maxima. For the solutions with the filling of some maxima of the potential the linear supports of distributions around the maxima should inevitably turn into the complex plane, with the complex conjugate endpoints, as shown in figure 1. We will call such solutions “unphysical”. The values of fractions $\nu_j = N_j/N$, $j = 1, 2, \dots, d$ are in one-to-one correspondence with the values of first $d-1$ moments $W_1, \dots, W_{d-1}$ and they completely fix the algebraic curve of the solution [8, 9]. The solutions where we fill only the minima of the effective potential will be called “physical”. The supports for such
solutions will be located only on the real axis, with positive weight for the distribution of the eigenvalues.

In the large N limit, the distribution of eigenvalues converges to a continuous function $\rho(x)$ and the corresponding SPE actually becomes the quadratic equation for resolvent (2.9). In this limit, the resolvent function, the eigenvalue distribution and the series of moments are closely related.

The moments can be computed via the resolvent (2.8) by a simple contour integration formula:

$$W_n = \frac{1}{2\pi i} \oint_{\Gamma} z^n G(z) dz = -\frac{1}{4\pi i} \oint_{\Gamma} z^n C(z) dz \quad (2.17)$$

where the contour $\Gamma$ must encircle all branch points of $G(z)$ (see [8, 9] for the details). We introduced here the “cut-function” — the square root of the discriminant — by the formula

$$C(x) = \sqrt{V'(x)^2 - 4P(x)} = \sqrt{D(x)} = g_{d+1} \sqrt{\prod_{k=1}^{d} (x - a_k)(x - b_k)}. \quad (2.18)$$

The roots $\{a_k, b_k\}$ of the discriminant become the branch points of the cut-function. Their number is always even. For real couplings in the potential the branch points lay only on the real axis or come in complex conjugate pairs. Note that it is easy to relate the eigenvalue fractions to these branch points:

$$\nu_j = \frac{1}{2\pi i} \oint_{\Gamma_j} G(z) dz = -\frac{1}{4\pi i} \oint_{\Gamma_j} C(z) dz \quad (2.19)$$

where the contour $\Gamma_j$ encircles anticlockwise only the branch points $\{a_k, b_k\}$, and $\nu_j$ can be in principal of either sign.

The eigenvalue density is expressed as the discontinuity of the cut function $C(x)$:

$$\rho(x) = \frac{1}{4\pi i} (C(x + i0) - C(x - i0)) = \frac{1}{2\pi} \Im C(x + i0) \quad (2.20)$$

and the moments can be expressed by the eigenvalue density:

$$W_n = \int_{-\infty}^{\infty} x^n \rho(x) dx. \quad (2.21)$$

## 2.2 Bootstrap method for the large N one matrix model

As we stated above, the loop equation (2.7), (2.9) has in general a continuum of solutions of the form (2.11) labeled by a finite number of parameters — the lowest moments $W_1, \dots, W_{d-1}$ which can take a priori arbitrary values. But not all of these solutions are “physical”, i.e. rendering all moments $W_k$, $k \in \mathbb{Z}_+$ real and compatible with the finite N Hermitian matrix ensemble. For example, the physical even moments should be positive, but this condition is not the only one.

A more general physical condition on a solution is the positivity of inner product for the matrix integral. This condition states that, for any operator of the form $O =$
$\sum_{i=0}^{n-1} \alpha_i M^i$ s.t. $\alpha \in \mathbb{R}^n$, and for any $n \in \mathbb{Z}_+$, we have the positive semi-definite quadratic form:<fnref n="10" />

$$ \langle \mathrm{Tr} \mathcal{O}^\dagger \mathcal{O} \rangle = \alpha^\mathrm{T} \mathcal{W} \alpha \ge 0 \quad \forall \alpha \in \mathbb{R}^n. \qquad (2.22) $$

Here we introduced the matrix $\mathcal{W}_{ij} = \mathcal{W}_{i+j-2}$ which will be called below for convenience the correlation matrix. The above condition is equivalent to the positive definiteness of correlation matrix:<fnref n="11" />

$$ \mathcal{W} \succeq 0. \qquad (2.23) $$

The condition (2.22) is obvious for a finite $N$ matrix model with converging integral (2.1), i.e. when $d+1$ is even and $g_{d+1} > 0$ in the potential (2.2), since the moments are just given by the integration of positive definite functions with a positive measure. In addition, at finite $N$ there exists only one solution for the moments.<fnref n="12" /> But it is far from trivial in large $N$ limit, where we have to understand what solutions from the continuum are really physical.

To do numerical bootstrap, we set a finite cutoff $\Lambda$, i.e. the highest moment in the correlation matrix $\mathcal{W}$ and loop equations is $\mathcal{W}_{2\Lambda}$ and the size of the correlation matrix is $(\Lambda+1) \times (\Lambda+1)$. We can use the loop equations to express the higher moments through a certain number of the lower moments and substitute them into the correlation matrix as functions of lower moments. Then the positivity of the correlation matrix provides us with algebraic inequality on these lower moments. In general, we expect to get the inequalities for each of the lower moments both from above and from below, for example:

$$ \mathcal{W}_{2\Lambda}^{\min} \le \mathcal{W}_2 \le \mathcal{W}_{2\Lambda}^{\max} \qquad (2.24) $$

at a cutoff $\Lambda$. In practice, the allowed region $[\mathcal{W}_{2\Lambda}^{\min}, \mathcal{W}_{2\Lambda}^{\max}]$ shrinks fast as we increase $\Lambda$, giving us tight bounds on $\mathcal{W}_2$.

We exemplify this approach on the case of quartic potential $V(z) = \frac{1}{2}z^2 + \frac{g}{4}z^4$. We plot in figure 2 the region for the allowed values of function $\mathcal{W}_2(g)$ under the assumption of $M \Leftrightarrow -M$ symmetry of the solution, i.e. $\mathcal{W}_{2l+1} = 0$. Under this assumption, all the higher moments are polynomials in terms of $\mathcal{W}_2$ and $g$. The positivity of the correlation matrix reduces to a list of algebraic inequalities on $\mathcal{W}_2$.<fnref n="13" />

It is a bit surprising that the bootstrap scheme described above for Hermitian one-matrix model is generally analytically solvable, considering that it is usually non-trivial

<footnote n="10">¹⁰Here we assume that all the expectation values we study are real. This is actually a non-trivial result from the $M \to M^\mathrm{T}$ symmetry of the potential. Since this symmetry is always present for all the models considered in this article, we will implicitly assume this to be always true.</footnote>

<footnote n="11">¹¹Here we slightly abused the notations: $\mathcal{W}$ sometimes means the matrix with finite size, involving only the moment up to a certain order $\Lambda$ (“cutoff”), and sometimes it means the infinite dimensional matrix. But the positive semi-definiteness is always well defined as a positivity of the corresponding quadratic form.</footnote>

<footnote n="12">¹²To define the matrix integral for unstable potentials, when $d+1$ odd or/and $g_{d+1} < 0$ one usually deforms appropriately the integration contours. Then the questions of positivity become less obvious for finite $N$. But we will see on the example of quartic potential that at infinite $N$ we can still have positivity for certain solutions, even for such, globally unstable, potentials.</footnote>

<footnote n="13">¹³To depict the allowed region, maybe the simplest way is the Mathematica RegionPlot function. Although not really numerically efficient, it is already good enough for the simplest quartic one-matrix model.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The allowed region for quartic model $V(z) = \frac{1}{2}z^2 + \frac{g}{4}z^4$ for different cutoff $\Lambda$, compared with the analytic solution. Here we assume $W_1 = 0$, we also note that for $g > 0$, we didn't plot $\Lambda \ge 5$ since they are almost indistinguishable from the analytic solution on the figure.</figcaption></figure>

to solve an infinite series of algebraic inequalities. In the next section we will propose a necessary and sufficient condition for the positivity constraint (2.23) by virtue of a result of solution of Hamburger moment problem. By this condition we can not only generally solve the bootstrap problem analytically but also justify why the numerical bootstrap process excludes the unphysical solutions of SPE described in section 2.1.

# 3 Hamburger moment problem and positivity of resolvent

As it should be clear from the previous section, the bootstrap method for one-matrix model has two main ingredients: 1. The Schwinger-Dyson loop equations for the moments of the random matrix variables; 2. The positive definiteness of the correlation matrix of these moments. We will rigorously prove that the second ingredient, in virtue of the Hamburger moment problem [41], picks up in the planar limit the solutions of loop equations only with real-positive supports of the matrix eigenvalue distribution. We will employ this condition to analytically solve the bootstrap condition.

## 3.1 Hamburger problem versus the positivity condition on resolvent

The loop equation (2.7) renders all possible large $N$, saddle point solutions of the Hermitian one matrix model. Some of them “look” physical, i.e. corresponding to the stable equilibrium in the Coulomb gas picture for the eigenvalues locked in the effective potential (2.16). For the other, the stable solution corresponds to “unphysical” picture when the supports of eigenvalue distribution become complex. What are the solutions captured by our bootstrap procedure?

As we just reviewed, numerical bootstrap of one matrix model consists of two ingredients: loop equation and positivity of correlation matrix. It sets a cutoff on both constraints
and gets a rigorous bound on the physical quantities we are interested in. Analytically all the information contained in the loop equations is encoded in the quadratic equation of resolvent (2.9), which has a simple solution (2.11). It describes the hyper-elliptic algebraic curve parameterized by complex variable $x$. It is natural to ask the question: can the positivity of correlation matrix also be expressed as a simple condition on the resolvent? The answer is, luckily and a bit surprisingly, yes.

For convenience, we make the following definition: a resolvent satisfies the positivity condition if the corresponding eigenvalue density (2.20) is supported on the real axis and is positive on its support. Our main conclusion of this section will be:

$$ \textit{Positivity of correlation matrix} \Leftrightarrow \textit{Positivity of Resolvent} \qquad (3.1) $$

We prove the necessity first. The proof is based on a well-known mathematical conclusion that will play an important role in our demonstration — the result of the solution of the Hamburger moment problem [41]:

For a given series of real numbers {$m_n$}$_{n=0}^{\infty}$, there exists a positive Borel measure $\mu$ such that:

$$ m_n = \int x^n d\mu \qquad (3.2) $$

if and only if the matrix $H_{ij} = m_{i+j-2}$ is positive semi-definite. Moreover, if there exist the constants $C$ and $D$, such that $|m_n| \le CD^n n!$, the measure is unique.

Applying the result of Hamburger momentum problem (3.2), we have for each moment $W_n = \int_{\mathbb{R}} x^n d\mu(x)$. We notice from the exponential bound condition (2.12) that $\mu$ must be supported in a finite region $[-R, R]$, since otherwise, if we have $\epsilon > 0$, such that $\mu((-\infty, -(R+\epsilon)) \cup (R+\epsilon, \infty)) = \mu_0 > 0$, then

$$ W_{2n} = \int_{\mathbb{R}} x^{2n} d\mu(x) > (R + \epsilon)^{2n} \mu_0, \qquad (3.3) $$

which contradicts the exponential bound (2.12).

Consequently for $|z| > R$:

$$ G(z) = \sum_{k=0}^{\infty} z^{-k-1} W_k = \sum_{k=0}^{\infty} z^{-k-1} \int_{[-R, R]} x^k d\mu(x) = \int_{[-R, R]} \frac{d\mu(x)}{z-x}. \qquad (3.4) $$

The exchange of infinite sum and integration is justified by Fubini's theorem. Due to this equation, $G(z)$ is analytic in the region outside of the disk $|z| > R$. The last equality in (3.4) enables us to analytically continue $G(z)$ to the whole region $\mathbb{C}\setminus[-R, R]$. So the function $G(z)$ must be analytical away from the real line, which eliminates the possibility of cuts between complex branch points. Comparing with (2.11) and (2.20), we come to the conclusion that all supports of $\rho(x)$ function must be located on the real line. For the positivity of eigenvalue density, we note that we can extract by contour deformation the coefficient of the series in (2.8):

$$ W_n = \frac{1}{2\pi i} \oint z^n G(z) dz = \int_{[-R, R]} x^n \rho(x) dx. \qquad (3.5) $$
By the uniqueness of solution of the Hamburger momentum problem, we must have $\rho(x)dx = d\mu$ i.e. they are equal in terms of positive measure.<fnref n="14" /> So we have $\rho(x)$ real supported and positive. This concludes our proof of necessity for (3.1).

The proof of sufficiency is straightforward. It is already true because the sufficiency is a part of the result of the solution of Hamburger momentum problem. For a more direct argument, suppose we have a resolvent that satisfies the positivity condition, i.e. (2.21) with $\rho(x) > 0$. We notice that the matrix $(X)_{ij} = x^{i+j-2}$ is trivially positive semi-definite for real $x$, so that if we integrate the matrix $X$ w.r.t. the positive measure $\rho(x)dx$, it stays positive semi-definite as well. The result of the integration is actually our correlation matrix $W$. This concludes the proof of sufficiency and hence of the equivalence (3.1).

It is easy to demonstrate by the direct computation that in the presence of complex branch points in $G(x)$ the corresponding correlation matrix is not positive definite. A simple example is the resolvent for the matrix model with the unstable potential $V(M) = -\frac{1}{2}\text{Tr } M^2$, which is $G(x) = \sqrt{x^2+2} - x = \frac{1}{x} - \frac{1}{2x^3} + O((\frac{1}{x})^4)$. We see that $\langle \text{Tr } M^2 \rangle = -\frac{1}{2}$ so that the correlation matrix is not positive definite.

This suggests the validity of the numerical bootstrap approach at least in the case of the one-matrix model: by imposing the positive semi-definiteness condition on the solutions of loop equations, at least for a finite cutoff $\Lambda$, we exclude the “unphysical” large $N$ solutions with the eigenvalue distributions having complex or negative supports, i.e. violating the hermiticity of the matrix measure.

The result of the present section actually enables us to analytically solve the bootstrap constraints, since the infinite series of inequalities from the positivity of the correlation matrix have been proven to be equivalent to the positivity property of resolvent. As an example, in the next subsection we will present the analytic result of solving the bootstrap problem of quartic one-matrix model.

The application of numerical bootstrap to the multi-matrix models, such as the one studied in section 5, has not as strong theoretical basis as the one presented in this section for the one-matrix model. However the arguments presented here give a good intuition why the numerical bootstrap can work even in the multi-matrix model case. In the following sections we will demonstrate its viability empirically, by showing its numerical efficiency for a specific, “unsolvable” matrix model.

## 3.2 Classification of physical solution of quartic one-matrix model

In figure 2, we saw that as we increase $\Lambda$ the allowed region converges to the analytic solution. One may ask whether the allowed region will ultimately exclude all other solutions as $\Lambda$ increases, or it will stabilize to a very tiny island which will not shrink further. The results of the current section will support the first of these options.

<footnote n="14">¹⁴We note that the uniqueness is not strictly necessary here. To see this, the reader can combine Stone-Weierstrass theorem and the fact that compactly supported continuous function is dense in $L_p$, $1 \le p < \infty$. Then if $\rho(x)dx$ is not positive almost everywhere then the positivity of correlation matrix is violated.</footnote>
In this section, we will apply the positivity of resolvent to the one-matrix model with quartic potential

$$V(x) = \frac{1}{2}\mu x^2 + \frac{1}{4}gx^4 \qquad (3.6)$$

in order to fully classify all physical solutions. This is equivalent to solving the positivity condition of one-matrix bootstrap analytically. Since in the previous subsection we have already formulated this problem as a precise mathematical theorem, we will not present the formal mathematical derivation here. For the details the reader can refer to the appendix B.

For the bootstrap problem we are trying to solve, we will not assume the $\mathbb{Z}_2$ symmetry of the solutions. This symmetry would mean $W_{2k+1} = 0$. We will see that there exist solutions that break this symmetry. In fact, for solutions we find numerically the breakdown or preservation of $\mathbb{Z}_2$ symmetry will be established dynamically and not necessarily imposed as an input. Alternatively, if we assume $\mathbb{Z}_2$ symmetry from the beginning, the numerical efficiency for such solutions considerably increases.

For the specific potential (3.6) the positivity condition for the resolvent

$$G(x) = \frac{1}{2} \left( V'(x) - \sqrt{V'(x)^2 - 4P(x)} \right) \\ = \frac{1}{2} \left( -\sqrt{-4g(W_2 + x(W_1 + x)) + (gx^3 + \mu x)^2 - 4 + gx^3 + \mu x} \right) \qquad (3.7)$$

translates into the condition that it has a only real positive eigenvalue distribution. This condition can be solved rigorously, namely:

1. $\mu = 1$ and $g \ge -\frac{1}{12}$: $W_1 = 0, W_2 = \frac{(12g+1)^{3/2}-18g-1}{54g^2}$.

2. $\mu = 1$ and $g < -\frac{1}{12}$, there is no possible solution.

3. $\mu = -1$ and $g \le 0$, there is no possible solution.

4. $\mu = -1$ and $g \ge \frac{1}{4}$: $W_1 = 0, W_2 = \frac{(12g+1)^{3/2}+18g+1}{54g^2}$.

5. $\mu = -1$ and $0 < g < \frac{1}{4}$: this situation is a bit involved. The bootstrap solution is a curve segment parametrized by $W_1$. Explicitly, the solution is a branch of the algebraic equation:

$$\begin{aligned} 0 = & 11664g^6 W_2^5 + (-27216g^5 - 864g^4) W_2^4 \\ & + W_2^3 (-16200g^5 W_1^2 - 13824g^5 + 19872g^4 + 1440g^3 + 16g^2) \\ & + W_2^2 ((43200g^5 + 33480g^4 + 888g^3) W_1^2 + 23040g^4 - 3232g^3 - 544g^2 - 16g) \\ & + W_2 ((4125g^4 - 22500g^5) W_1^4 + (-65280g^4 - 24568g^3 - 1480g^2 - 16g) W_1^2 \\ & \qquad + 4096g^4 - 8704g^3 - 1072g^2 - 32g) \\ & + 3125g^5 W_1^6 + (18500g^4 - 3925g^3 - 16g^2) W_1^4 \\ & + (-1024g^4 + 22848g^3 + 7096g^2 + 608g + 16) W_1^2 - 4096g^3 - 512g^2 - 16g. \end{aligned} \qquad (3.8)$$
The physical branch of solution is selected by the one passing through $W_1 = 0$, $W_2 = \frac{1}{g}$,<fnref n="15" /> with $-W_{1c} \le W_1 \le W_{1c}$.

For $0 < g \le \frac{1}{15}$,

$$W_{1c} = \frac{2\sqrt{4500g^2 + 75g - 2(1 - 15g)^{3/2}(60g + 1) + 2}}{75\sqrt{5}g^{3/2}}, \quad (3.9)$$

and for $\frac{1}{15} < g < \frac{1}{4}$,

$$W_{1c} = \frac{2\sqrt{12000g^2 + 1200g - \sqrt{3}(20g + 7)^{3/2}(60g + 1) + 102}}{75\sqrt{5}g^{3/2}}. \quad (3.10)$$

This reproduces the exact solution of quartic one-matrix model. In figure 2 we have already compared the exact solution and the numerical bootstrap result for $\mu = 1$. A typical comparison for $\mu = -1$ case is figure 3.

In figure 3 we take a representative from each phase and compare it with the above analytic solutions. We see that the numerical bootstrap results converge quickly to the analytic result. A distinguishable feature of these figures is that the allowed region is not guaranteed to be convex. This is very different from the convex optimization problems which we encountered in CFT bootstrap and S-matrix bootstrap. Generally, the large-scale non-convex problem is hard and usually unsolvable. We will discuss in the next section how to overcome this difficulty.

## 3.3 Comments

Here we present several comments on the results of this section:

1. There may exist certain doubts on particular choices of the positivity condition of the correlation matrix in numerical bootstrap. In the work [28, 30], the authors showed that in some cases one only needs the positivity of even moments $W_{2k} \ge 0$ to make bootstrap converging to the analytically known solution. But in general one should be careful about the choices of the positivity condition. For example, consider the model with $V(x) = -\frac{1}{2}x^2 + \frac{1}{4}x^4$. Under the assumption of the $\mathbb{Z}_2$ symmetry, the loop equations of this model read:

$$W_{2k} = W_{2k-2} + \sum_{l=0}^{k-2} W_{2l}W_{2k-4-2l} \quad k = 2, 3, 4, \dots \quad (3.11)$$

We see that the positivity condition on even moments only provides us with the constraint $W_2 \ge 0$, evident by induction in loop equations. In this situation we can bootstrap the physical solution only with the positivity condition on the full correlation matrix. This fact explains to some extent why the convergence in figure 3 is not as fast as for the model with positive quadratic coefficient.

<footnote n="15">¹⁵Actually for $\mu = -1$ and $0 < g \le \frac{1}{4}$, the $\mathbb{Z}_2$ symmetry preserving solution is just $W_1 = 0$, $W_2 = \frac{1}{g}$. So the first discontinuity of $W_2(g)$ at $g = 1/4$ happens for second derivative.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The comparison of numerical bootstrap with our analytic bootstrap results (the black point or the black curve on the figures) in section 3.2 with $\mu = -1$ in potential (3.6). They are representatives from different phases of the model, with $g = 1, \frac{1}{10}, \frac{1}{30}$ (for the figures from above to below, respectively). We notice even visually that for $g < \frac{1}{4}$ i.e. when the symmetry breaks, the exact solution is a non-convex set.</figcaption></figure>
2. For the one-matrix integral with integration over the unitary matrix instead of the Hermitian matrix, we can establish and justify a similar bootstrap method. This enables us with the analytic solution of such bootstrap problems. The main difference in this case comparing to the Hermitian integral is that the correlation matrix is of the form $W_{ij} = W_{i-j}$. It is called the Toeplitz matrix in linear algebra.<fnref n="16" /> For this correlation matrix, we have the following result of solution, this time for trigonometric moment problem:

For a given series of real numbers {$m_n$}$_{n=-\infty}^{\infty}$ such that $m_{-k} = m_k^*$, there exists a positive Borel measure $\mu$ on $[0, 2\pi]$ such that:

$$m_n = \frac{1}{2\pi} \int_0^{2\pi} \exp(-int)d\mu(t) \qquad (3.12)$$

if and only if the matrix $T_{ij} = m_{i-j}$ is positive semi-definite.

Applying this result to our unitary matrix integral, we come to the conclusion that the positivity of correlation matrix for large $N$ unitary matrix integral is equivalent to the positivity of the eigenvalue density which is supported on the unit circle in the complex plain.

## 4 Relaxation bootstrap method

Now we turn to the discussion of the bootstrap method for multi-matrix models. We will see that a naive generalization of the previous one-matrix model bootstrap will lead to a Non-linear SDP.<fnref n="17" /> But it is widely known that a general large-scale Non-linear SDP cannot be solved efficiently. In this section we will propose a systematic numerical bootstrap procedure to solve the large $N$ multi-matrix models via SDP.

SDP, unlike the Nonlinear-SDP which is directly applicable in the case of large $N$ matrix model bootstrap [28, 30], has a long history in academic research as well as in applied sciences. The standard primal form of SDP is:<fnref n="18" />

$$ \begin{array}{ll} \text{minimize} & \displaystyle\sum_{i=1}^{m} c_i x_i \quad \text{w.r.t.} \ \{x_1, x_2, \dots, x_m\} \in \mathbb{R}, \\ \text{subject to} & \displaystyle\sum_{i=1}^{m} F_i x_i - F_0 \succeq 0, \quad F_i \in \mathcal{S}^n. \end{array} \qquad (4.1) $$

Here $\mathcal{S}^n$ denotes the space of $n \times n$ real symmetric matrices. As long as we can transform our bootstrap problem to the form (4.1), we can get rigorous bounds on the physical quantities of interest — linear functions of {$x_1, x_2, \dots, x_m$} — by efficiently solving the SDP problem (4.1).

<footnote n="16">¹⁶For the Hermitian integral the correlation matrix is of the form of the Hankel matrix.</footnote>

<footnote n="17">¹⁷SDP means semi-definite programming.</footnote>

<footnote n="18">¹⁸There exists also the dual form of these problems, which will be discussed in appendix C. We also note that in some literature different conventions for dual and primal for SDP are used.</footnote>
So the problem reduces to the question how to efficiently transform our matrix integral problem into the constraints of the form (4.1). Then the original physical problem is transformed into a purely numerical SDP problem.

In this section we will describe our relaxation bootstrap method on the example of single trace moments in a large N two-matrix model with the partition function:<fnref n="19" />

$$Z = \lim_{N \to \infty} \int d^{N^2} A \, d^{N^2} B e^{-N \mathrm{tr} V(A, B)} \qquad (4.2)$$

where $V(A, B)$ is assumed to be a so far general polynomial in $A$ and $B$, to make the loop equations more tractable. In the next section we will apply it to a model with a concrete potential, generally unsolvable by the known analytic methods. We will see that our method has four types of constraints: loop equations, global symmetries, positivity of correlation matrix and positivity of relaxation matrix (which will be explained later).

## 4.1 Physical constraints

To make this section as self-contained as possible, we briefly review here the terminology already introduced in the previous sections and show how the constraints of the type (4.1) are specified in the two-matrix model.

The positivity of correlation matrix is still at the heart of our method. Since we are doing numerical analysis, we set the cutoff $2\Lambda$ to the length of operators that we are considering, i.e. to the length of “words” built from two “letters” — the matrices $A$ and $B$: $O = ABBAAAB...$. For any word $O$ of the length $\le \Lambda$, we assume:

$$\langle \mathrm{Tr} O^\dagger O \rangle \ge 0. \qquad (4.3)$$

The set of words with length $\le \Lambda$ is a vector space spanned by all the words constructed from two letters with the length cutoff $\Lambda$. This is a set of $L = 2^{\Lambda+1} - 1$ elements which we denote as $O_i$, where $i$ runs from 1 to $L = 2^{\Lambda+1} - 1$. For example, when $\Lambda = 2$ the basis of this vector space reads:

$$I, A, B, A^2, AB, BA, B^2. \qquad (4.4)$$

We can expand the equation (4.3) w.r.t. this base:

$$\langle \mathrm{Tr} \left( \sum_{i=1}^{L} \alpha_i O_i \right)^{\dagger} \left( \sum_{i=1}^{L} \alpha_i O_i \right) \rangle = \alpha^{\mathrm{T}} \mathcal{M}_L \alpha \ge 0. \qquad (4.5)$$

Let us introduce the correlation matrix $\mathcal{M}_{Lij} = \langle \mathrm{Tr} O_i^\dagger O_j \rangle$ which consists of expectation values of operators with the lengths up to $2\Lambda$. Since (4.3) is true for all operators, the condition (4.5) holds for all $\alpha \in \mathbb{R}^L$, i.e. the semi-definite positivity of correlation matrix is ensured:

$$\mathcal{M}_L \succeq 0. \qquad (4.6)$$

This correlation matrix condition can be directly applied to the two-matrix model. We see that the main difference with the one-matrix model is that the dimension of correlation matrix grows exponentially with $\Lambda$.

<footnote n="19">¹⁹The generalization to multi-matrix models with more matrices is straightforward.</footnote>
Another important ingredient for our bootstrap method is the loop equations. For the two-matrix model it can be schematically represented as:

$$ \int d^{N^2} A d^{N^2} B \operatorname{tr} \left( \partial_M \left( \operatorname{Word} \times e^{-N \operatorname{tr} V(A,B)} \right) \right) = 0, \quad M = \{A, B\} \quad (4.7) $$

where “Word” means the matrix word built by arbitrary finite product of matrices A and B.<fnref n="20" /> The differentiation $\partial_M$ can be either w.r.t. the matrix A or w.r.t. the matrix B.

The loop equations for large N multi-matrix model in general close on all words.<fnref n="21" /> Schematically, they have the following quadratic form:

$$ \langle \operatorname{Tr} (\operatorname{Word}_l \times \partial_M V(A, B)) \rangle = \sum_{l_1=1}^{l} \langle \operatorname{Tr} \operatorname{Word}_{l_1-1}^{(M)} \rangle \cdot \langle \operatorname{Tr} \operatorname{Word}_{l-l_1}^{(M)} \rangle \quad (4.8) $$

which is a direct generalization of (2.7) of the one-matrix model. Here Word$_{l_1-1}$ and Word$_{l-l_1}$ are the words obtained by cutting the word Word$_l$ in two words whenever one has the matrix M on the l₁-th place in Word$_l$. The matrix factor $\partial_M V(A, B)$ in the l.h.s. comes from the derivative of the exponential factor in (4.7), which generically renders a sum over single trace operators with lengths from l₁ to l₁ + d (the degree of polynomial V(A, B) is assumed to be d + 1). So we expect that a loop equation of length l₁ involves quadratic relations of operators with lengths up to l₁ + d. In the next section we will precise all these steps on a particular example of the two-matrix model.

The set of all loop equations can be efficiently generated by applying the derivatives in $M = \{A, B\}$ to any word of the length less than a certain cutoff.<fnref n="22" /> However, the loop equations obtained in this way are not all independent, which means that there may exist linear dependence and/or algebraic dependence among them. It turns out that these redundancies are numerically crucial when applying the SDP solver to the constraints of our system, but they are not important at this stage of explanation. We will discuss these technicalities in appendix D.

If the model has some discrete symmetries, such as $M \to -M$ or $A \leftrightarrow B$, it is not necessary to assume them from the beginning in our bootstrap scheme, but factoring it out will significantly increase our numerical efficiency if we are only interested in the symmetry preserving solution. Generally, the symmetry assumptions not only simplify the loop equation by reducing the number of operators<fnref n="23" /> but in certain cases they make the correlation matrix block diagonal, thus greatly simplifying our problem. We will encounter this situation in the next section for a concrete model.

At last, we identify all the operators related by cyclicity of trace and the reversion of the word. These transformations also reduce considerably the number of unknowns in our scheme.

<footnote n="20">²⁰Note that “word” is not yet traced, so that generically it is not cyclically symmetric: a cyclic transformation gives in general a new word.</footnote>

<footnote n="21">²¹Here we mean that there is generally no infinite subset of loop equations and operators closed among themselves. This fact will be explored further in section F.</footnote>

<footnote n="22">²²We will discuss the detail of the choice of the cutoff in the appendix D.</footnote>

<footnote n="23">²³For example, if the potential has $Z_2$ symmetry $A \leftrightarrow B$, we could identify all the operators identical by $A \leftrightarrow B$ transformation.</footnote>
In summary, for the two-matrix integral (4.2), assuming the global symmetries or not, we set up all the physical constraints. A natural question is what is the solution of these constraints. But this is not a good question since generally, apart from some solvable models where the loop equations close on a very limited subclass of operators (like in a two-matrix model [42, 43] or some n-matrix models [6, 19, 20]), the number of operators grows faster than the number of constraints, which means that the solution is a region in an extremely high dimensional space. A constructive question at this stage is: given a cutoff to the length $\le 2\Lambda$, what is the minimal or maximal possible value of a physical quantity? This amounts to asking what is the allowed interval when the region allowed by the constraints is projected on the linear subspace corresponding to the specific physical quantity.

Rephrasing it in the language of optimization theory, we deal with the problem of the form:

$$ 
\begin{array}{ll}
\text{minimize} & c^T x \\
\text{subject to} & x^T A_i x + b_i^T x + a_i = 0 \quad (i\text{'th loop equation}), \\
\text{and} & M_0 + \sum_{j=1}^{L} M_j x_j \ge 0
\end{array}
\qquad (4.9) $$

where $c$ is a vector defining the dynamical quantity we want to optimize and $x$ is the column vector of all our operator expectations $x_i = \langle \text{Tr} O_i \rangle$, up to the length $2\Lambda$. The quadratic loop equation (in the middle) is written in the vector form, where $A_i$ is the quadratic form encountered in the $i$th equation; linear and constant terms are represented accordingly. The matrix inequality is the expansion of the correlation matrix in terms of the operator expectations. This is certainly not equivalent to the standard SDP which we introduced by (4.1) since the quadratic equations represent non-convex conditions. One of the conventional methods to deal with it is relaxation.

## 4.2 Relaxation matrix

The constraints discussed in the last section define a problem which is called Non-Linear SDP in optimization theory. There are indeed some solvers specialized for it but, from our limited trials, they are not mature enough to solve large-scale problems such as the ones we encountered in matrix bootstrap. To improve the situation, we propose to modify the problem (4.9) by relaxing the non-convex conditions involving the non-linear loop equations, into convex ones. Our intuition here is that we don't really need all of the loop equation constraints for our bootstrap method to converge as $\Lambda$ increases.

To see how our method works, let us begin with a simple example which will provide us with a heuristic argument. Suppose we have only three quadratic “loop equations”:

$$ 
\left\{ 
\begin{array}{l}
x^2 = T_1 \\
y^2 = T_2 \\
xy = T_3
\end{array} 
\right. 
\qquad (4.10) $$

Here $T_i = \sum_j q_i^j w_j$, $(i = 1, 2, 3)$ denote linear combinations of some other variables $w_1, w_2, \dots$. These equations are of course non-convex. But we can relax them to make them
convex by replacing $x^2 = T_1$ with $x^2 \le T_1$ or, in the positive semi-definite matrix form,

$$ \begin{pmatrix} 1 & x \\ x & T_1 \end{pmatrix} \succeq 0. \qquad (4.11) $$

We can do the same thing with the second equation $y^2 = T_2$, to relax it to a convex condition. But the same operation cannot be reproduced for equation $xy = T_3$, since neither $xy \le T_3$ nor $xy \ge T_3$ is convex.<fnref n="24" /> It is tempting to consider the positive semi-definite combinations:

$$ (x + \alpha y)^2 \le T_1 + \alpha^2 T_2 + 2\alpha T_3, \forall \alpha \in \mathbb{R}. \qquad (4.12) $$

It is not very elegant to implement (4.12) by introducing extra parameters like $\alpha$, although numerically this is viable. Can we write instead of (4.12) a condition that does not contain explicitly $\alpha$? In fact yes. Since $T_1 \ge x^2 \lor T_2 \ge y^2$, we only need the discriminant of (4.12) w.r.t. $\alpha$ to be non-positive, to exclude the existence of real solution for $\alpha$ when (4.12) becomes an equality. That means

$$ (T_1 T_2 - T_1 y^2 - T_2 x^2 - T_3^2 + 2 T_3 x y) \ge 0 \qquad (4.13) $$

is equivalent to (4.12) for all $\alpha \in \mathbb{R}$. In its turn, it is equivalent to:

$$ \mathrm{Det} \begin{pmatrix} 1 & x & y \\ x & T_1 & T_3 \\ y & T_3 & T_2 \end{pmatrix} \ge 0. \qquad (4.14) $$

Combining (4.11) and (4.14) we come to the conclusion that:

$$ \begin{pmatrix} 1 & x & y \\ x & T_1 & T_3 \\ y & T_3 & T_2 \end{pmatrix} \succeq 0. \qquad (4.15) $$

This is mathematically more elegant and numerically more efficient.

To apply this relaxation method to the case of our loop equations is a simple generalization of what we just proposed. We make in the loop equation the substitution $\langle \mathrm{Tr} \mathcal{O}_i \rangle \langle \mathrm{Tr} \mathcal{O}_j \rangle = X_{ij}$, or in matrix notations:

$$ X = x x^T \qquad (4.16) $$

where again $x$ is the column vector whose components are $\langle \mathrm{Tr} \mathcal{O}_i \rangle$. Formally, this changes the loop equations in (4.9) to a linear form:

$$ \mathrm{Tr} X \mathcal{A}_i + b_i^T x + a_i = 0. \qquad (4.17) $$

To apply the relaxation method sketched above, we relax (4.16) by imposing the inequality:

$$ (\alpha^T x)^2 \le \alpha^T X \alpha, \quad \forall \alpha \in \mathbb{R}^L \qquad (4.18) $$

<footnote n="24">²⁴Because the bilinear form $xy$ is not positive semi-definite.</footnote>
which is equivalent to:

$$X \succeq xx^T. \qquad (4.19)$$

By Schur's complement, this can be re-arranged into a more compact form:

$$\mathcal{R} = \begin{pmatrix} 1 & x^T \\ x & X \end{pmatrix} \succeq 0. \qquad (4.20)$$

Here we introduced the relaxation matrix by $\mathcal{R}_{ij} = X_{ij}$ and $\mathcal{R}_{0i} = \mathcal{R}_{i0} = \langle \mathrm{Tr} \mathcal{O}_i \rangle = x_i$. This step concludes our translation of the nonlinear bootstrap problem into an SDP. This SDP takes now a numerically much more tractable, convex form:

$$
\begin{array}{ll}
\text{minimize} & c^T x \\
\text{such that} & \mathrm{Tr} X \mathcal{A}_i + b_i^T x + a_i = 0, \\
\text{and} & M_0 + \sum_{j=1}^L M_j x_j \succeq 0, \qquad (4.21) \\
\text{and} & \begin{pmatrix} 1 & x^T \\ x & X \end{pmatrix} \succeq 0.
\end{array}
$$

It has now two types of variables to bootstrap: a column vector variable $x$ and a symmetric matrix variable $X$.

Several comments are in order:

* One of the primary questions to the method is: does the relaxed SDP generate the same bounds as the previous Non-linear SDP problem? Generally, the answer is “no”. It is obvious that when the optimal solution of the relaxed problem satisfies the constraint of the original problem the relaxed problem will generate the same bound as the original one. From our experience, this is not the case for any finite $\Lambda$.<fnref n="25" /> But as we increase the cutoff $\Lambda$, the mismatches for the quadratic conditions are tending to zero. So we are tempted to believe that for infinite $\Lambda$, the relaxed problem and the original problem give the same result for most of the questions we are interested in. This indicates that the non-linear constraints in the loop equations are somehow contained in the positivity conditions for correlation matrix and relaxation matrix.

* One can regard our relaxation scheme as a numerical compromise: doing relaxation we replace the nonlinear equalities by linear inequalities but we can thus explore the correlation matrices of a much higher order since we can significantly increase the length cutoff $\Lambda$. This enable us to embrace more information from correlation matrix. Our numerical results in the next section will show that this is a worthy trade-off.

* There is another point of view on our relaxation formulation (4.21). The problem (4.21) is actually the dual of the dual of the problem of (4.9). Although this

<footnote n="25">²⁵More precisely, if the relaxation is saturated for the optimal solution, we expect that the relaxation matrix will only have one non-zero eigenvalue. But practically, we always observe multiple non-zero eigenvalues for the relaxation matrix.</footnote>
fact is in principle simple to show its proof is quite lengthy, so we put it into the appendix C. In that appendix, we also briefly review the definition and basic facts about the dual formulation. As it is known, the dual problem of any general optimization problem is always convex [44], so the double dual is guaranteed to be convex. In some sense, this point of view is more general and universal.

* We believe that the key condition for the relaxation method to work well is that, under our bootstrap assumption, there is a unique exact solution.<fnref n="26" /> Then since a single point (corresponding to the $\Lambda = \infty$ solution of bootstrap) is convex, our relaxation procedure leading to convex constraints will not make the results too different even for a finite but sufficiently large $\Lambda$. However, we observed in section 3.2 that the set of exact solutions may become non-convex in the presence of a symmetry breaking. In such situation, we need further assumptions to make the exact solution unique. We will further discuss these aspects in the next section when bootstrapping the symmetry breaking solutions.

# 5 Bootstrap for “unsolvable” two-matrix model with $\mathrm{Tr}[A, B]^2$ interaction

In this section, we implement the relaxation bootstrap method described in the previous section to the case of generically unsolvable large $N$ two-matrix model:

$$Z = \lim_{N \to \infty} \int d^{N^2} A d^{N^2} B e^{-N \mathrm{tr}(-h[A, B]^2/2 + A^2/2 + gA^4/4 + B^2/2 + gB^4/4)} \quad (5.1)$$

where the integration goes over Hermitian matrices $A$ and $B$. This model is unsolvable analytically for generic parameters $h$ and $g$, at least with the known methods, such as reduction to eigenvalues or the character expansion. It is still analytically solvable for some particular values: for $g = 0$ it can be reduced to a specific one-matrix model and solved via saddle point method or via the reduction to a KP equation [34, 45]; for $h = 0$ it reduces to two decoupled one-matrix models; for $h = \infty$ we have $[A, B] = 0$ and it reduces again explicitly to another eigenvalue problem. These particular solvable cases are useful to test the power of our numerical method.

The present section is organized as follows: the bootstrap results for the model (5.1) with particular choice of parameters $g = 1$, $h = 1$ (which represent a generic analytically “unsolvable” example) are shown in section 5.1. Then in section 5.2 we compare the bootstrap result for the analytically solvable cases $h = 0$ or $g = 0$ with the corresponding analytic solution, to test our method. In section 5.3, we explore the phase diagram of the this model and make several comments about the convergence rate in different regions. At last, in section 5.4, we investigate the symmetry breaking in the model by our relaxation bootstrap method.

<footnote n="26">²⁶Here exact solution means bootstrap solution with infinite cutoff.</footnote>
## 5.1 Bootstrap solution for a generic choice of g, h

In this subsection, we present the results of the bootstrap for the model (5.1) where we specify, for definiteness, the parameters: $g = 1, h = 1$. We stress that this choice has nothing specific for the properties of the model and it is made mostly for the demonstrative reasons, as an example of generic values of parameters. The method appears to be very efficient almost everywhere in the physical domain of parameters $g, h$, except when we approach the critical lines where it is less efficient. We will discuss in the next subsection the phase structure of the model in the $g-h$ plane.

The symmetry of this model can be described by the Dihedral Group $D_4$,<fnref n="27" /> with generator:

$$ \left\{ \begin{array}{l} A \to -A \\ B \to -B \\ A \leftrightarrow B \end{array} \right. \qquad (5.2) $$

We saw already on the example of the one matrix model that in the large $N$ limit there could be a multitude of saddle point solutions, many of them breaking this kind of symmetries. We begin with the study of $D_4$ symmetric large $N$ solutions. Later we will discuss the solutions with broken symmetries as well.

In the fully $D_4$ symmetric solutions, only the operators with even number of $A$ and even number of $B$ can be non-vanishing, and we should identify the operators under the exchange $A \leftrightarrow B$. Obviously this assumption of $D_4$ symmetry of solution is in principle not necessary for our bootstrap method to work. However, assuming this symmetry we gain a lot in the efficiency since we are left with approximately 1/8 of operators comparing to a general non-symmetric setup. It also happens that the symmetry assumption simplifies the correlation matrix by much. Namely, when constructing the correlation matrix, only the words with the same $Z_2$ parity in both $A$ and $B$ can appear in the inner product for a non-vanishing correlator. So our correlation matrix break into 4 block-diagonal matrices, corresponding to $Z_2$ parities in $A$ and $B$: even-even, even-odd, odd-even, odd-odd. By $A \leftrightarrow B$ symmetry, the even-odd and odd-even blocks are actually the same. So the original correlation matrix can be reduced to three block diagonal matrices: even-even, even-odd, and odd-even.

Here we bootstrap the allowed region for the first two non-vanishing operators $t_2 = \langle \text{Tr}A^2 \rangle$, $t_4 = \langle \text{Tr}A^4 \rangle$.<fnref n="28" /> According to (4.21), this corresponds to setting the objective function of the optimization problem as:

$$ \text{Minimize: } t_2 \cos \theta + t_4 \sin \theta \qquad (5.3) $$

Scanning it in the interval $0 \le \theta < 2\pi$ we can fix the allowed region for these two operators. Using the general method described in the last section we can use the SDP solvers to solve

<footnote n="27">²⁷Actually we implicitly assume the $A \to A^T, B \to B^T$ which basically means that all the moments are real. We will assume throughout this paper that this symmetry cannot be broken. At least intuitively, this is unlikely to happen in our model (5.1).</footnote>

<footnote n="28">²⁸Here we give up the $W_k$ notation for the moments we used in one-matrix model since for two-matrix model the moments cannot be characterized by a single positive number.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4. The allowed region of $t_2 - t_4$ of model (5.1) with parameter $g = 1, h = 1$ for the cutoff $\Lambda = 7, 8, 9, 10, 11$. We recall the definition of $\Lambda$: the longest operators in the correlation matrix and in the loop equations have the length $2\Lambda$.</figcaption></figure>

these problem. The readers interested in the details of the implementation can refer to the appendix D, where we gather all the technical detail of numerical implementations. We also demonstrated in appendix E our numerical procedure explicitly, step by step, on the example of the system with $\Lambda = 4$ cutoff.

Let us demonstrate our results for various values of the length cutoff $\Lambda$. We summarized the allowed regions for the first two correlators $t_2 = \langle \text{Tr}A^2 \rangle$ and $t_4 = \langle \text{Tr}A^4 \rangle$ in figure 4. The regions for $\Lambda = 10$ and $\Lambda = 11$ are too small to be plotted on the figure, so we give here the upper and lower bound of $t_2$ and $t_4$. For $\Lambda = 10$:

$$ \begin{cases} 0.421780275 \le t_2 \le 0.421785491 \\ 0.333339083 \le t_4 \le 0.333343006 \end{cases} \qquad (5.4) $$

and for $\Lambda = 11$:

$$ \begin{cases} 0.421783612 \le t_2 \le 0.421784687 \\ 0.333341358 \le t_4 \le 0.333342131 \end{cases} \qquad (5.5) $$

We see here that for $\Lambda = 11$ we already have a six digits precision at $g = h = 1$. The $\Lambda = 11$ calculation is the largest problem in this work, it is done with SDPA-dd, a solver in SDPA family with the double-double float type. The input to SDPA has 95 variables, with the correlation matrix size: even-even 683, odd-odd 682, even-odd 1365, and with relaxation matrix size 8. We note that this is still within the capability of a single laptop, it only takes 150000s CPU time for a single maximization cycle. We also stress that these inequalities, unlike the Monte Carlo methods, are exact: increasing the cutoff $\Lambda$ we can only improve the margins.
## 5.2 Demonstration for analytically solvable cases

It is instructive to apply our numerical method to the analytically solvable cases $g = 0$ or $h = 0$, which is a good check for our approach, convincing us that it works well indeed even for the generic parameters, where we have no analytic data to compare with. In this part we will firstly review the analytic solution for both cases and then compare it with the numerical results of our relaxation bootstrap method.

As we mentioned, for $h = 0$ this model reduces to two decoupled one-matrix models — the case which we already discussed and studied analytically in section 3.2. Integrating out one of the decoupled matrices, we expect the operator containing only one matrix to have exactly the same expectation value as for the result in section 3.2:

$$g_c = -\frac{1}{12}, \quad t_2 = \frac{(12g + 1)^{3/2} - 18g - 1}{54g^2}. \qquad (5.6)$$

For $g = 0$, this model is already solved analytically in [34, 45]. Here we simply present the analytic solution derived there in our notations and normalization. To have a compact form, we introduce the short-hand notations $E = E(m)$, $K = K(m)$, $\vartheta = E/K$, where K and E are the complete elliptic integrals of first and second kind:

$$K(m) = \int_{0}^{\pi/2} \frac{d\theta}{\sqrt{1 - m^2 \sin^2(\theta)}}, \quad E(m) = \int_{0}^{\pi/2} \sqrt{1 - m^2 \sin^2(\theta)} d\theta. \qquad (5.7)$$

We introduce the new parameter $m$ related with $h$ by:

$$h(m) = \frac{K((m-1) - 2(m-2)\vartheta - 3\vartheta^2)}{6\pi^4}, \qquad (5.8)$$

and we can express $t_2$ as:

$$h(m)t_2(m) = \frac{1}{12} - \frac{K^2(-(m-2)(m-1) + 10(m-2)\vartheta^2 + 2((m-6)m + 6)\vartheta + 10\vartheta^3)}{5\pi^2(-(m-1) + 2(m-2)\vartheta + 3\vartheta^2)}. \qquad (5.9)$$

This formula is valid when $h > 0$. For $h < 0$, we need to analytically continue the solution to the other sheet of Riemann surface of the variable $m$. For that we introduce the analytic continuation of the elliptical integral $K(m)$ and $E(m)$:

$$K_a = K_a(m) = \frac{K\left(\frac{1}{m}\right) + iK\left(1 - \frac{1}{m}\right)}{\sqrt{m}}, \\ E_a = E_a(m) = \frac{-(m-1)K\left(\frac{1}{m}\right) + iK\left(\frac{m-1}{m}\right) + mE\left(\frac{1}{m}\right) - imE\left(\frac{m-1}{m}\right)}{\sqrt{m}}, \qquad (5.10)$$

$$\vartheta_a = E_a/K_a.$$

To make (5.8) and (5.9) valid for $h < 0$, we simply replace all the $K$, $E$, $\vartheta$ by $K_a$, $E_a$, $\vartheta_a$.

The critical point of the smallest possible $h_c$ for $h < 0$ can be defined as the solution of the equation:<fnref n="29" />

$$\frac{dh(m)}{dm} = 0 \qquad (5.11)$$

<footnote n="29">²⁹We thank Nikolay Gromov for sharing with us his computation of $h_c$.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5. Comparison with the exact analytic solution of model (5.1) with $h = 0$, i.e. two decoupled quartic one-matrix model. The lower plot is for $\Lambda = 8$.</figcaption></figure>

which can be numerically solved as:

$$h_c \approx -0.04965775; \quad t_{2c} \approx 1.18960475. \qquad (5.12)$$

The comparison of our numerical results with analytic result is presented on figure 5 and figure 6. Indeed, we see that our numerical results nicely agree with the analytic formula (5.6) and (5.9). An apparent feature of these plots is that when $g < 0$ or $h < 0$, the allowed region is much larger than the one for the positive coupling case, thus giving less of precision. In general we have the worst convergence in the neighborhood of critical value. We will discuss this feature in more details in the next subsection.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of the numerical bootstrap results with the exact analytic solution of the model (5.1) with $g = 0$. The lower plot is for $\Lambda = 8$.</figcaption></figure>

There is another fact which is not obvious from the figure 5. If we compare this figure with figure 2 in section 2 we will find that for same values of $\Lambda$, the non-relaxed one-matrix bootstrap bound for $t_2$ (denoted by $W_2$ in that section) and our relaxation bootstrap bound for $h=0$ case of the 2-matrix model actually coincides within the error bar. This is a very striking feature of our relaxation method since we relaxed all the quadratic equalities to inequalities but we compensated this with many more mixed operators of two decoupled matrices. So the correlation matrix is much larger in the relaxed case and the final results are basically the same. We will see this feature of relaxation again when we discuss later the bootstrapping of the symmetry breaking solutions. We don't have a very clear explanation for these phenomena in general.
## 5.3 Phase diagram and convergence rate

In this part, we will discuss the phase diagram of the matrix model (5.1) and the corresponding convergence rate in different regions of the diagram.

In general, for finite N matrix integral the potential $V(A, B)$ must be bounded from below to define a sensible integral over Hermitian matrices. But this is not necessary for a large N theory, where we only need deep enough local minima to have a stable saddle point solution. Even for the unstable potentials, the tunnelling effects between the local minima, or to the infinity are suppressed exponentially. We saw this in section 5.2, where the bootstrap procedure allowed the existence of solutions with negative values of g and h. This provides us with a possibility to study the boundaries of possible g and h values (we will call the region of possible g and h values the feasible region in the following) even when the corresponding potential is not bounded from below.

Before going deeper into the technicalities of bootstrapping the boundaries of the feasible region, we can get a rough estimate of them by deriving the parameter region of g and h which leads to the matrix potential bounded from below. It is obvious that the domain $(h \ge 0, g \ge 0)$ is one part of the region we are looking for. Another, less obvious part is $(h < 0, g \ge -4h)$, as in this case we should have:

$$ 
\begin{aligned}
\mathrm{tr}V(A, B) &= \mathrm{tr}\left(-\frac{h[A, B]^2}{2} + \frac{A^2}{2} - \frac{h(A^4 + B^4)}{4} + \frac{B^2}{2} + \frac{(g+4h)(A^4/4 + B^4/4)}{2}\right) \\
&= \mathrm{tr}\left(-\frac{h((AB+BA)^2)}{2} + \frac{(A^2-B^2)^2}{2} + \frac{A^2}{2} + \frac{B^2}{2} + \frac{(g+4h)(A^4/4 + B^4/4)}{2}\right) \\
&\ge 0. \qquad (5.13)
\end{aligned}
$$

The union of these two domains represents the maximal region where the matrix potential is bounded from below, since for $(h \ge 0, g < 0)$ and $(h < 0, g < -4h)$ we can always find A, B configurations where the potential is not bounded from below. For $(h \ge 0, g < 0)$, one of these configurations is taking $B = 0$ and $A \to \infty$. For $(h < 0, g < -4h)$, we simply put A and B to be some constants $\alpha$ times generalized Pauli matrices of dimension $N \alpha\sigma_1$ and $\alpha\sigma_2$, where $\alpha$ is a large real number. Then we have:

$$ \mathrm{tr}V(\alpha\sigma_1, \alpha\sigma_2) = N \left(\alpha^2 + (g+4h)\alpha^4\right). \qquad (5.14) $$

This must be unbounded from below when $(h < 0, g < -4h)$.

In conclusion, the region of potential bounded from below is $(g \ge 0 \cap g \ge -4h)$. In addition, the domain $(g \ge 0 \cap g \ge -4h)$ is guaranteed to lie within the feasible region. But due to the large N effects, we expect the feasible region to be a little bigger than that. Specifically, for analytically solvable cases, when $h = 0$ we have $g \ge -\frac{1}{12}$ and when $g = 0$ we have $h \ge -0.04965775$. These facts give us an additional information about the location of the boundary of the feasible region.

To numerically bootstrap the boundary of the feasible region, we can obtain the critical boundary between the allowed and forbidden parameter regions by bisection. Namely, for a given $\Lambda$ we fix $h$ and take two values of $g$, as $g_1$ and $g_2$. Here $g_1$ is a point that is guaranteed to be forbidden for a given $h$, and $g_2$ is a point that is guaranteed to be allowed. Then we test the geometric average value $g_m = \frac{g_1g_2}{g_1+g_2}$. If $g_m$ is allowed, then we make the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The numerical phase diagram of the model (5.1). The gray region is strictly excluded by our relaxation bootstrap method at $\Lambda = 8$. In the blue region the matrix potential is bounded from below and its boundary is located above the straight lines $g = 0$ and $h = -\frac{1}{4}g$. The red and green dots are the critical values for $h = 0$ and $g = 0$.</figcaption></figure>

substitution $g_2 = g_m$, otherwise we take $g_1 = g_m$. In this way we can recursively approach the maximal forbidden value of $g$ at fixed $h$. Then we scan over the values of $h$ and get the plot shown in figure 7.

Some explanations for the plot figure 7 are in order. The gray region is rigorously forbidden as the result of bootstrap at $\Lambda = 8$. On the contrary, the white region is not guaranteed to be allowed for any physical large $N$ solution. As we increase $\Lambda$, the gray region will expand a little. But we have several hints about the position of the exact boundary line:

* We notice the red and green dots on the plot, which are the critical points of the analytic solutions. They are located on the exact boundary of feasible region, i.e. no matter how large is $\Lambda$, the gray curve cannot go beyond these two points. From this fact we convince ourselves that our numerical curve in figure 7 is already very accurate, since the red dot and the green dot are very close to the gray curve.

* The blue region is the region where the potential is strictly bounded from below. It is enclosed by the lines $g = -4h$ and $g = 0$. Its boundary can be considered as the exact solution in the “classical” limit $\hbar \to 0$ for this matrix integral, where $1/\hbar$ is the coefficient put in front of the potential. In this case we scale the couplings as $h \to 1/\hbar$, $g \to 1/\hbar$. The boundary of the gray will coincide for $\hbar \to 0$ with the boundary of the blue area on figure 7. Then inside the blue area we have a well-defined theory even for finite $N$. At large $N$ and finite $h, g$ there is a gap between blue region and gray region, as is visible on the figure 7.
### 5.3.1 Rate of convergence

As the reader may have noticed already in section 5.2, when $g < 0$ or $h < 0$ the convergence is very bad compared to the case $g > 0$ and $h > 0$. From our experience, this is a generic situation when we are outside of the blue region in figure 7, which is defined by the region of parameters yielding a potential bounded from below. For example, figure 8 depicts the allowed region for $t_2$ when we fix $g = 1$ and scan over $h$ in the neighborhood of $h = -1/4$. It is clear from this figure that for $h < -1/4$ there is drop in the rate of convergence. Actually, from numerical data, the difference of the upper bound and the lower bound varies between the orders of magnitude from $10^{-4}$ to around $10^{-2}$ when $h$ varies from $h = -0.25$ to $h = -0.26$.

Nonetheless we can get a rather accurate estimate of physical quantities in the region discussed in the last paragraph. We note that in figure 5 and figure 6, the analytic solution is very close to the lower bound, comparing to the upper bound.<fnref n="30" /> Actually, as we increase $\Lambda$, the lower bound stabilizes already at rather small $\Lambda$. Empirically this is a typical behavior in the unbounded region. Under the assumption that there is a unique solution satisfying the constraint for arbitrarily large $\Lambda$, we expect that the optimization results for the maximum and the minimum of $t_2$ will ultimately converge with increasing $\Lambda$ to the same value. This has been proven for some parameters of the one-matrix model in section 3.2, and we have strong numerical evidence to believe it will hold for our model (5.1) as well. So we can simply bootstrap the physical quantities by the minimization of $t_2$ in this region (in the following, we will call this procedure the minimization scheme as opposed to the maximization scheme). Comparing it to the analytically solvable particular cases we learned that this method can yield especially accurate estimate of physical quantities. However, we lost the rigorous margin in the region with good convergence (blue region in figure 7).<fnref n="31" />

There exists a region in the phase diagram figure 7 where the bootstrap is valid only for very high cutoff $\Lambda$: it is $g < 0$, $|g| \ll h$. The figure 9 shows the allowed region when we fix $h = 1$ and vary $g$. For the lower bound of pink region $\Lambda = 9$, there are some numerical instabilities for $-0.09 < g < 0$. From careful inspection of our data at various values of $\Lambda$ it seems that the lower bound at $\Lambda = 9$ should stabilize in this region at the value $t_2 \approx 0.5$ if no numerical instabilities happened in our SDP solver. We notice a few very distinguishable features of this plot:

1. For a fixed $\Lambda$, there is a region where $t_2$ is slightly larger than 0.5 and not bounded from above. In other words, the dual SPD problem for the upper bound is infeasible. In this bad region of parameter space, the bootstrap with such $\Lambda$ essentially tells us nothing about the right physical values. Luckily, the “bad” region is shrinking when we increase $\Lambda$, and hopefully it will disappear when we have a high enough cutoff.

<footnote n="30">³⁰We believe that the upper bound and the lower bound converge to the same value, but it seems they have rather different convergence behaviors.</footnote>

<footnote n="31">³¹This situation is similar to that of the early days of conformal bootstrap when people used the kink of a plot to estimate the dimension of operators in the 3d Ising model, cf. [46].</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 8. The allowed region for $t_2$ when we fix $g = 1$ and vary $h$. In the region where the potential is bounded from below i.e. $h \ge -\frac{1}{4}$, we have a decent convergence whereas for $h < -\frac{1}{4}$ the convergence gets much slower.</figcaption></figure>

2. As already stated in the last paragraph, when we are not in the “bad region”, the minimization scheme converges much faster than the maximization scheme. So for a reasonable estimate of the operator expectation we should privilege the minimization scheme.

3. We also notice that for the region $g > 0, h > 0$, the convergence is excellent as expected, but there is a huge drop in the rate of convergence in the neighborhood of $h = 0$.

## 5.4 Bootstrapping the symmetry breaking solution

In the previous parts of this section, we always assumed the global symmetry, or in other words, we bootstrapped the symmetry preserving solutions. Here in the following, we will make the first attempt to study the symmetry breaking solutions with our relaxation bootstrap method. Consequently, in this subsection we will not make assumptions on a specific global symmetry of operator expectations. For example, we will assume that it is possible to have:

$$t_1 = \langle \mathrm{Tr} A \rangle = \langle \mathrm{Tr} B \rangle \neq 0 \qquad (5.15)$$

and any other nonzero expectations containing odd number of letter A or B, unlike the solutions with such $A \to -A, B \to -B$ symmetry.

To understand the general features of symmetry breaking solutions, figure 3 in section 3.2 is a good source for our intuition. We see on that figure that the exact solution
<figure><img src="image_10.png" /><figcaption>Figure 9. The allowed region for the fixed $h = 1$ and varying $g$. For the lower bound of the pink region $\Lambda = 9$, there are some numerical instabilities at $-0.09 < g < 0$. From careful inspection of our data at various values of $\Lambda$ it seems that the lower bound at $\Lambda = 9$ should stabilize in this region at the value $t_2 \approx 0.5$ if no numerical instabilities happened in our SDP solver.</figcaption></figure>

is not unique anymore but there is rather a continuous family of solutions parametrized by $t_1$. This is a non-convex set of exact solutions, so we don't expect that our relaxation bootstrap method, as applied in the case of figure 4, will converge to such a non-convex set as $\Lambda$ increases. Namely, if we impose the relaxation bootstrap constraint without the assumption of $A \to -A$, $B \to -B$ symmetry, then minimize the value of

$$t_1 \cos \theta + t_2 \sin \theta \qquad (5.16)$$

and then scan over $\theta$ in $[0, 2\pi)$, we expect to get a convex set instead of the non-convex one, due to the convex nature of the relaxation method.

So to bootstrap the symmetry breaking solutions, new techniques are needed to tackle the non-convexity. We will describe the general method for bootstrapping solutions and then we apply it to the study of our model (5.1).

### 5.4.1 Schemes for symmetry breaking bootstrap

The main problems in the study of symmetry breaking solutions in the multi-matrix model of the type considered here are:

1. How to identify the range of parameters for which the model has a possible symmetry breaking solution?

2. How to numerically bootstrap the symmetry breaking solution?
The answer to the first problem is quiet straightforward. We can establish the relaxed constraint without the symmetry assumption, and bootstrap a dynamical quantity which signals the symmetry breaking. For example, for the $A \to -A$, $B \to -B$ symmetry breaking solution we take the objective function ($c^T x$ in (4.21)) as:

$$t_1 = \langle \mathrm{Tr} A \rangle \qquad (5.17)$$

and for the $A \leftrightarrow B$ symmetry breaking we take the objective function as:<fnref n="32" />

$$\langle \mathrm{Tr} A^2 \rangle - \langle \mathrm{Tr} B^2 \rangle. \qquad (5.18)$$

If the bound of the symmetry breaking expectation is significantly larger than the error bar at the current $\Lambda$ for a given value of parameters, we believe that this is a strong signal of existence of a symmetry breaking solution.

For the second problem, we propose to transform the non-convex set of exact solutions to a convex one, which means that for the case of figure 3 we fix $t_1$ by $t_1 = t_1^{(0)}$ in our bootstrap procedure. For this particular value of $t_1$ we should have at infinite cutoff $\Lambda$ a unique exact solution for $t_2$ and for other higher moments, which is definitely a convex set. Therefore our relaxation bootstrap method with a finite cutoff $\Lambda$ will yield a rigorous upper bound and lower bound for $t_2$. Next we scan over $t_1^{(0)}$ until such values that the problem becomes infeasible. In this way we get the allowed region in $t_1, t_2$ plane.

The above method is easily generalizable to the problem of bootstrapping solutions with the other symmetry breaking patterns. Namely, we establish the bootstrap scheme by fixing the dynamical quantity signaling the symmetry breaking, and then we bootstrap the quantities we are interested in. At this step, we expect that after fixing such dynamical quantity, the exact solution of the bootstrap problem is unique. At the next step, we scan over all possible values of the quantity which was fixed in the previous step. In this way we can bootstrap a non-convex set of solutions.

There is another possible solution for the first problem, i.e. to locate the symmetry breaking region. We can assign to the dynamical quantity signaling the symmetry breaking a specific value and then use the method similar to that of section 5.3, i.e. using a bisection to approach the maximal possible value of expectation signaling the symmetry breaking. In principle this bisection method could have given us a tighter bound than our initially proposed method. However, from our test, the two methods yield basically the same numerical result, so we will not bother to use the bisection method in what follows.

### 5.4.2 Numerical results for symmetry breaking solution

Here we apply the method proposed above to the model (5.1). Our results in this part concern the breaking of the following symmetries:

$$A \to -A, B \to -B \qquad (5.19)$$

<footnote n="32">³²Here in these two situations the dynamical quantity signaling the symmetry breaking is respectively $\langle \mathrm{Tr} A \rangle$ and $\langle \mathrm{Tr} A^2 \rangle - \langle \mathrm{Tr} B^2 \rangle$.</footnote>
and

$$A \leftrightarrow B. \qquad (5.20)$$

In the bootstrap setup, we don't impose the global symmetry assumptions for the corresponding symmetries, i.e. that the non-singlet operator expectations of the $Z_2$ symmetry vanish. Then we pick up the dynamical quantities signaling the symmetry breaking as:

$$\langle \mathrm{Tr} A \rangle \qquad (5.21)$$

and

$$\langle \mathrm{Tr} A^2 \rangle - \langle \mathrm{Tr} B^2 \rangle, \qquad (5.22)$$

respectively and set them as the objective functions in the corresponding bootstrap problem.

As the result, in the feasible region of figure 7 we didn't find any evidence of the existence of a symmetry breaking solution for the model (5.1). We tried several points in different regions of figure 7. The results show that the maximized values are always lying within the error bar (typically $10^{-3}$ and $10^{-4}$, depending on the cutoff $\Lambda$ and the parameters $g$ and $h$). In particular, for $\Lambda = 8$ and some generic values of $g$ and $h$, we have:

$$-10^{-4} \lesssim \langle \mathrm{Tr} A \rangle, \langle \mathrm{Tr} A^2 \rangle - \langle \mathrm{Tr} B^2 \rangle \lesssim 10^{-4}. \qquad (5.23)$$

We believe this to be a strong evidence that the two symmetries we investigated are not spontaneously broken for all the regions in figure 7.

Some other interesting facts:

1. For $g = 0$, i.e. when the quartic coefficient vanishes, the preservation of symmetry is automatic from the loop equation. This fact provides us with the intuition that the commutator square interaction is to some extent not a symmetry-breaking interaction. Regarding that at $h = 0$ the model is not in symmetry breaking phase, since it reduces to two decoupled one-matrix models, intuitively it points on the absence of symmetry breaking phase for model (5.1) (with positive coefficients in front of quadratic terms).

2. For the region $h > 0$ and $g$ slightly smaller than zero, we have a very large upper bound for the exposed quantities, sometimes of order 10, which might signal the symmetry breaking. But we note that the bootstrap convergence is really bad in this region where some bootstrap results for symmetry preserving solution are presented on figure 9, and the error bar here is almost infinitely big. So we believe this cannot be a reliable evidence that there a symmetry breaking takes place in this region.

As we don't find evidence for the existence of symmetry breaking solutions for the model (5.1), we consider the same model but with negative coefficients in front of quadratic terms:

$$Z = \lim_{N \to \infty} \int d^{N^2} A d^{N^2} B e^{-N \mathrm{tr}(-h[A,B]^2/2 - A^2/2 + gA^4/4 - B^2/2 + gB^4/4)}. \qquad (5.24)$$
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The allowed $t_1 - t_2$ region for $\Lambda = 6$ and $\Lambda = 9$. The corresponding parameters in the model (5.24) are $g = \frac{1}{100}$, $h = 0$. The shaded region is the result of one-matrix bootstrap. The black line is the exact analytic solution described in section 3.2. The green and red dots are the upper bounds and lower bounds of our relaxation method from scanning over $t_1$.</figcaption></figure>

We know from the section 3.2 that for $h = 0$ where we have just two decoupled one-matrix models, we have a symmetry breaking phase for $g < 1/4$. At such values of $g$ we can test our method for bootstrapping the symmetry breaking solutions.

In figure 10 we compare the results of our relaxation method described above with the exact results and the one-matrix bootstrap plot at the same cutoff and the parameters $g = \frac{1}{100}$, $h = 0$. We see that our method is indeed able to bootstrap the symmetry breaking solution, even though it is non-convex. It is especially striking that not only our relaxation method converges to the highly non-convex exact solution, but it even coincides with the one-matrix bootstrap at each cutoff $\Lambda$ within the error bars. It seems that, in spite of some loss of information when applying the relaxation method, we recover this information by considering the positivity condition of the mixed operators containing both matrices, such as $\langle \text{Tr}ABAB \rangle$. We don't have yet a good explanation why these two approaches give equal or very close results. We also note that in this case the maximization scheme converges faster than the minimization scheme. Namely, the upper bound (green dots) in the plot is much closer to the exact solution than the lower bound (red dots). This suggests that if we are looking for a good approximation for the exact solution, we should use the upper bound solution as the best approximation.

For generic values of $h$ and $g$ for the model (5.24), the convergence is slower than in the analytically solvable particular case. It would be good to understand whether such a situation for solvable versus unsolvable models is typical. In figure 11 we plot the bootstrap result for $g = \frac{1}{30}$, $h = \frac{1}{15}$. Obviously, it is still a symmetry breaking solution. We expect that taking the upper bound we can get a very accurate estimation of the physical
<figure><img src="image_12.png" /><figcaption>**Figure 11.** The allowed region for $g = \frac{1}{30}$ and $h = \frac{1}{15}$. These numerical results are obtained for $\Lambda = 8$.</figcaption></figure>

quantities. We didn't try to further increase the value of $\Lambda$, being already satisfied to see
that the proposed method works for rather generic values of parameters.

# 6 Conclusion and discussion

In this work, we develop further the matrix bootstrap method pioneered in the papers [28, 30] and propose a crucial improvement — the relaxation procedure — applicable to a large class of multi-matrix problems and allowing to bootstrap them with a much higher precision. The relaxation transforms a Non-linear SDP, with the non-linearity due to the structure of loop equations, to the usual, linear SDP. We demonstrate the efficiency of our approach on the analytically unsolvable two-matrix model and establish its phase structure with rather high precision. The method appears to work well even for the discrete symmetry breaking large *N* solutions.

Actually, the efficiency of the matrix bootstrap, based on the positivity of the correlation matrix, remains an enigma. Nevertheless we managed to theoretically study this question in the one-matrix model and to establish precisely the class of physical solution singled out by such bootstrap.

As it was noticed in [28, 30, 47] the numerical bootstrap for large *N* multi-matrix
models presents a higher degree of difficulty than the bootstrap of the large *N* one-matrix
model. The bootstrap study of multi-matrix models was restricted to low orders in the
length cut-off for the moments (or “words” forming them) and consequently it provided us
with a limited precision. The main reason for this inefficiency is the non-convex formulation
of the problem. Our main task in this work was to overcome these drawbacks of the matrix
bootstrap.

Compared to the cited above papers we achieved a better understanding and efficiency
of the matrix bootstrap in the following aspects:
* In the case of the large *N* Hermitian one-matrix model, we managed to prove that the bootstrap constraints pick up precisely the exact “physical” solutions, with positive measure for the distribution of the eigenvalues (the corresponding resolvent has only the cuts on the real axis). In other words, we established the exact analytic solutions of the bootstrap conditions for the one-matrix model, thus justifying the numerical bootstrap techniques. We don’t have yet the generalization of such analytic argument for the multi-matrix models, which is an interesting question to address.

* Then for the multi-matrix models, we developed the relaxation bootstrap method to overcome the crucial obstacle of non-convexity of the original problem. We demonstrated that this relaxation method was a systematic approach, capable to provide the numerically viable procedure for the large *N* multi-matrix models. We tested this method on a model that is genuinely analytically unsolvable (unlike Lin’s 2-matrix model with cubic interactions in [28]). For particular parameters, when the analytic solution is known, our numerics reproduces extremely well the analytic results. For generic values of these parameters, we bootstrap the physical values with a remarkable precision (6 digits).

* This method is also able to detect critical behaviors, though the precision gets less impressive in the vicinity of critical lines.

* Remarkably, our bootstrap method is also applicable for bootstrapping the symmetry-breaking solutions and transform the non-convex problem to a convex one.

Here we make several comments on the bootstrap method proposed here and sketch out some further directions:

* All the numerical results in this study can be, in principle, reproduced on a single decent laptop in a decent time laps. So it looks very promising to implement it on a big cluster with parallelization. The main technical difficulty of the method is that we used the precision bigger than the machine precision (double-double or quad-double) in our current work. We believe that this is mainly due to the fact that our problem is badly-scaled as an SDP: the involved variables can have very different orders of magnitude. It would be good to find a systematic approach to scale appropriately the variables for very large-scale problems.

* The positivity of the correlation matrix in the matrix bootstrap method must be satisfied for any multi-matrix integrals with a reasonably converging positive measure. It follows from the fact that the integral of a positive function against a positive measure is positive. Contrary to the conformal bootstrap and *S*-matrix bootstrap where unitarity is one of the most important conditions, we don’t know whether the unitarity or reflection positivity can be imposed in the bootstrap method for the matrix models. We also hope that our method can be generalized for bootstrapping non-unitary quantum field theories.
* We expect that the correlation matrix contains a lot of redundancies, i.e. very few of its minors may contain 99% of the information of the whole correlation matrix. This is reminiscent of a similar feature of the conformal bootstrap: we don't impose the positivity condition on all spin channels, rather a very limited number of spin channels are good enough to make the algorithm to converge [26]. At the moment we don't have any scheme to isolate the minors of the correlation matrix that are more important than the other, which would be very beneficial when considering large-scale problems.

* It would be interesting to apply our methods to the Matrix Quantum Mechanics, in the spirit of the work [29], including for the non-singlet states there. Another interesting two-matrix model to study by bootstrap would be the generalization of (1.6) by taking the q-deformed version of interaction: tr[$A, B$]$_q^2$ = tr($qAB - q^{-1}BA$)$^2$. This model interpolates between the solvable cases with tr($A^2B^2$) or tr($BABA$) interactions [35].

* An obvious, and one of the most ambitious possible applications of our relaxation bootstrap methods is the lattice Yang-Mills theory. We have thus good chances to significantly improve on this way the very preliminary results of [30]. A method alternative to the wide-spread Monte-Carlo simulations, even at large $N$ [48], would be extremely welcome for the study and a deeper understanding of QCD. Obvious advantages of the bootstrap method based on Migdal-Makeenko loop equations [31] w.r.t. Monte-Carlo are: i) Exact inequalities on loop averages, no statistical error; ii) absence of finite boundary conditions (the lattice is infinite); iii) One gets some information on all loops at once up to a given length, although with better precision for short loops. That gives access to more of the physical quantities. The obvious drawback is the limited length of Wilson loops. We hope to establish by the future numerical work whether this drawback is crucial indeed.

## A Analytic solvability of two-matrix model with cubic interactions and arbitrary potentials

The Hermitian 2-matrix model with the general cubic interactions between two matrices and general potentials in the action

$$S = \text{tr} \left( h (AB^2 + BA^2) + W(A) + \tilde{W}(B) \right), \qquad (A.1)$$

has been studied in [28] by numerical bootstrap method, as an example of bootstrap approach to an analytically “unsolvable” matrix integral. Here we show that this matrix model is in fact analytically solvable for generic potentials $W$ and $\tilde{W}$, in the sense that the matrix integral can be explicitly reduced to $\sim N$ amount of variables, instead of the original $\sim N^2$ matrix variables, which in principle allows the application of the saddle point method at large $N$. Our derivation will be schematic and we will repeatedly neglect the non-dynamical factors before the integral of partition function. It is unclear whether this
integrability influences the efficiency of Lin's bootstrap method but this is our motivation to choose a different, truly "unsolvable" 2-matrix integral, with the interaction tr[A, B]², as the main example of application of bootstrap in this paper.

Since tr(A + B)³ = tr(A³ + 3A²B + 3AB² + B³) we can always rewrite it as

$$S = \text{tr} \left( (A+B)^3/3 - V(A) - \tilde{V}(B) \right) \qquad (A.2)$$

where $V(A) = A^3/3 - W(A)$ and $\tilde{V}(B) = B^3/3 - \tilde{W}(B)$ and we have set without loss of generality $h = 1$.

We can always reduce it to Itzykson-Zuber-Charish-Chandra (IZC) integral by an extra matrix integration, following the trick similar to proposed in [6] in the context of solution of the Potts model on dynamical planar Feynman graphs. Namely, represent the first term in (A.2) in terms of an extra hermitian matrix integral over $X^{33}$

$$e^{\text{tr} C^3/3} = \int d^{N^2} X e^{\text{tr}(iXC + F(X))}, \quad C = A + B \qquad (A.3)$$

where the function $F(X)$ is defined as the inverse matrix Fourier transform:

$$e^{\text{tr}F(X)} = \int d^{N^2}C e^{\text{tr}(-iXC+C^3/3)} = \int \prod_j dc_j e^{c_j^3/3} \Delta^2(c) \frac{\det_{j,k} e^{-ix_j c_k}}{\Delta(c)\Delta(x)}, \qquad (A.4)$$

i.e. it represents the matrix Airy function. In the last equality we applied the IZC integral. Then we rotate the contour of $X$ integration, change variable $C \to -C$ and write the partition function of cubic 2MM model in the form

$$Z = \int d^{N^2}X \int d^{N^2}C e^{\text{tr}(XC-C^3/3)} \int d^{N^2}A e^{\text{tr}[XA-V(A)]} \int d^{N^2}B e^{\text{tr}[XB-\tilde{V}(B)]} \qquad (A.5)$$

We can compute the function $F(X)$ in terms of the eigenvalues of $X = \Omega^\dagger x \Omega$ where $x = \text{diag}\{x_0, x_1, \dots, x_{N-1}\}$. We compute the angular integral in (A.4) via IZC integrals:

$$e^{\text{tr}F(X)} = \Delta^{-1}(x) \int \prod_j dc_j e^{x_j c_j - c_j^3/3} \Delta(c) = \frac{W[-\text{Ai}(x_0), \dots - \text{Ai}(x_{N-1})]}{\Delta(x)} \qquad (A.6)$$

where in the denominator we have the $N \times N$ "Wronskian" of Airy-type functions:

$$\text{Ai}(x) = \int_C dc \, e^{-xc + \frac{1}{3}c^3}, \qquad (A.7)$$

and:

$$W[f_0(x_0), \dots, f_{N-1}(x_{N-1})] = \sum_{\sigma \in S_n} \text{sgn}(\sigma) \prod_{i=0}^{N-1} \partial_{x_i}^{\sigma(i)} f_i(x_i). \qquad (A.8)$$

Here the complex contour $C$ is usually chosen so that it goes from infinity with the slope $-\pi/3$ and ends up at infinity with the slope $\pi/3$. However, when we study the limit $N \to \infty$ the saddle point configuration of the eigenvalues will adjust itself to the relevant

<footnote n="33">³³We drop here and further all inessential overall factors.</footnote>
distribution on the real axis given by a solution of the integral saddle point equation. Similarly, the bootstrap numerical procedure should single out such solutions.

Then we treat similarly the other two integrals in (A.5) and represent them also in terms of Wronskians of

$$f(x) = \int da e^{xa-V(a)}, \quad \tilde{f}(x) = \int db e^{xb-\tilde{V}(b)}. \qquad (A.9)$$

In this way, we managed to rewrite the cubic 2MM entirely in terms of eigenvalue integral:

$$Z = \int \prod_j dx_j \frac{W[-\mathrm{Ai}(x_0), \dots] W[f(x_0), \dots] W[\tilde{f}(x_0), \dots]}{\Delta(x)}. \qquad (A.10)$$

Hence we reduced the cubic two-matrix integral (A.2) to an explicit integral over N eigenvalues of an auxiliary matrix X. We treat such a matrix model as “solvable” though the further details of the explicit solution can be rather involved. Instead of studying the saddle point in terms of wronskians it is better to apply the method (inspired by [49]) which was proposed by V.Kazakov and I.Kostov for solution of Potts model on random planar graphs [50]; it is well presented in [20]. We will not pursue here this route and we leave it for the future work.

## B Solving the positivity condition of the resolvent

We saw in section 3.1 that the positivity of correlation matrix is equivalent to the positivity of the resolvent. As it was noticed there, this equivalence enables us to analytically solve the bootstrap condition. Here we propose a general method to solve the lower moments from the positivity condition of resolvent. This finishes our analytic solution of the bootstrap problem corresponding to the Hermitian one-matrix model. As a specific example, we also apply this method to one-matrix model with quartic potential, for which the results were summarized in section 3.2.

### B.1 Cuts and zeros

Finding the lower moments from the positivity of the resolvent is a well-posed problem in complex analysis, and abundant mathematical tools can be employed to solve it. Here we study the configuration of cuts and zeros of the cut function defined in the main text (2.18). Due to the polynomiality of the discriminant $D(x) = C(x)^2$ we can give the full classification of all possible configurations of cuts and zeros of the cut function $C(x)$ on the real line. It turns out that the properties of these configurations provide not only necessary but also sufficient condition for the positivity of the resolvent. Generally, if we know a configuration of the cuts and zeros, we are able to fix a few lowest moments which we want to find the solution. This is how we solve the positivity condition of the resolvent. In the following we will be considering a Hermitian one-matrix model with general polynomial potential $V(x)$.

First let us pick a single cut of $C(x)$ on the real axis, namely $[a_i, b_i]$, as shown in figure 12. The positivity condition of the resolvent implies that for $x \in [a_i, b_i]$ we have
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. The left plot shows the sign of $C(x)$ in the neighborhood of a positive cut. The right plot shows there must be at least one zero, or generally odd number of zeros, between two cuts to fix the sign.</figcaption></figure>

$\mathrm{Im}C(x + i0) \propto \rho(x) \ge 0$. From the definition of the cut function, we must have $C(x) > 0$ in the right neighborhood of $b_i$, and $C(x) < 0$ in the left neighborhood of $a_i$. A direct consequence is that we must have at least one zero, or generally odd number of zeros, between two positive cuts to fix the sign. This is also illustrated in figure 12.

There is yet another constraint on the zeros. We notice that $V'(x)$ and $C(x)$ get unbounded at infinity, but $G(z) \xrightarrow{z \to \infty} \frac{1}{z}$, i.e. it is analytic there. As a result, asymptotic behavior of $C(x)$ must match the asymptotics of $V'(x)$. For example, if $\lim_{x \to \infty} V'(x) < 0$, we must add another zero to the right of all the cuts to fix the sign of $C(x)$, preserving both positivity of the resolvent and the asymptotic behavior.

For the zeros of $D(x)$ which are not located on the real axis, there exist roots with even multiplicities since there shouldn't exist complex cuts for $C(x)$. We also note that since $D(x)$ is a polynomial with real coefficients, all its complex roots must come in pairs.

The above analysis gives the way to count the possible number of cuts. Here we list the maximum number of cuts $m$ when the degree of the potential is $d+1$, under the asymptotic behavior $(\pm, \pm):^{34}$

$$ \bullet \quad (-, +): 2m + 2(m - 1) \le 2d \Rightarrow m \le (d + 1)/2, $$

$$ \bullet \quad (+, -): 4 + 2m + 2(m - 1) \le 2d \Rightarrow m \le (d - 1)/2, $$

$$ \bullet \quad (+, +) \text{ and } (-, -): 2 + 2m + 2(m - 1) \le 2d \Rightarrow m \le d/2. $$

## B.2 A working example

In this part, we will use these cuts and zeros considerations solve the bootstrap condition of the model with potential:

$$ V(x) = \frac{1}{2}x^2 + \frac{1}{4}gx^4, \quad g < 0. \qquad (B.1) $$

This potential has a $(+, -)$ asymptotic, so in the minimal case, we must have a positive cut in the middle, and two zeros, one placed on the right and another one on the left of the cut, to fix the asymptotics. Since the polynomial $D(x) = V'(x)^2 - 4P(x)$ is of 6th degree, we have already reached the maximum number of zeros. The conclusion is that the cut configuration in figure 13 is the only possibility for such asymptotic behavior.

To see when this cut configuration is possible, we note that:

$$ D(x) = V'(x)^2 - 4P(x) = (gx^3 - x)^2 - 4gx^2 + 4g - 4g(W_1x + W_2) = D_1(x) - D_2(x). \quad (B.2) $$

<footnote n="34">³⁴For example, $(-, +)$ means $\lim_{x \to -\infty} V'(x) < 0$ and $\lim_{x \to \infty} V'(x) > 0$.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 13. The only possible cut configuration for asymptotic (+, -).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14. The only possible cut configuration for $\mu = 1$ and $g < 0$. The critical case is $g = -\frac{1}{12}$, below that, there is no more bootstrap solution.</figcaption></figure>

Here we split the discriminant into the part depending on $W_1$ and $W_2$:

$$D_1(x) = (gx^3 - x)^2 - 4gx^2 + 4g \qquad (B.3)$$

and the part depending only on $g$:

$$D_2(x) = 4g(W_1x + W_2). \qquad (B.4)$$

$D_2(x)$ is a straight line with negative intercept. As depicted in figure 14, it is possible that a straight line crosses $D_1(x)$ with sufficient number of intersections, but only for $g \ge g_c = -\frac{1}{12}$. For $g < g_c$ the qualitative shape of the graph of $D_1(x)$ disqualifies the only possible cut configuration figure 13, excluding the existence of any bootstrap solution.

From the figure 14, it is obvious that $W_1 = 0$. Otherwise we won't have a resolvent satisfying the positivity condition. We can find $W_2$ from the vanishing of the discriminant of the polynomial $D(x)$, namely:

$$W_2 = \frac{(12g + 1)^{3/2} - 18g - 1}{54g^2} \qquad (B.5)$$

which corresponds to the standard one-cut IBPZ solution [14].

To make our intuitive arguments above more systematic, and applicable to higher degree potentials, let us list the conditions by which the positivity of the resolvent translates into the properties of the polynomial $D(x)$:

1. $D(x)$ has 6 real roots counting multiplicity, two of them are double roots. Two of them are single roots.

2. The simple roots of $D(x)$ lie in between the double roots on the real line.
<figure><img src="image_17.png" /><figcaption>Figure 15. 4 solutions under the condition 1 when $g = -1/15$. The blue line is the function $D_1(x)$ which is independent of moments $W_1, W_2$, whereas other color straight lines are different solution of $D_2(x)$. We see only the green line has the correct configuration of zeros depicted in figure 14.</figcaption></figure>

The property 1 boils down to the condition that the degree of the greatest common divisor of $(D(x), D'(x))$ is 2. We apply the Euclidean algorithm to $D(x)$ and $D'(x)$ to rewrite them in the form:

$$ \left\{ \begin{array}{l} D = q_1 D' + r_1 \\ D' = q_2 r_1 + r_2 \\ r_1 = q_3 r_2 + r_3 \\ r_2 = q_4 r_3 + r_4 \end{array} \right. \qquad (B.6) $$

The remainder $r_4$ is a polynomial of degree 1 in $x$ and we set all its coefficients to zero. In this way we get 2 algebraic equations which, in principle, fix $W_1$ and $W_2$ in terms of $g$. For the reality of the roots, we note that this means the discriminant of $(D(x), D'(x)) = r_3$ and the discriminant of $D(x)/r_3^2$ is non-negative. These are already a lot of conditions. Luckily the mathematica Reduce function can treat these conditions efficiently, leaving us with four solutions. The four solutions correspond to different orders of single zeros and double zeros. We show the corresponding solutions in figure 15. We can take advantage of the condition 2 to select the only physical solution for this case:

$$ \left\{ \begin{array}{l} -\frac{1}{12} \le g < 0 : W_1 = 0, W_2 = \frac{(12g+1)^{3/2}-18g-1}{54g^2}. \\ g < -\frac{1}{12} : \text{No bootstrap solution.} \end{array} \right. \qquad (B.7) $$

## C Dual formulation and relaxation

In this appendix, we review some basic facts about dual formulation in optimization theory and clarify the relationship between our relaxation method introduced in section 4 and the dual formulation. For the readers interested in more details about the optimization theory, the book [44] is a good starting point.
## C.1 Dual problem of general optimization problem

Consider a general optimization problem of the form:

$$ 
\begin{array}{ll}
\min & c(x) \\
\text{subject to} & f_i(x) \le 0, \quad i = 1, \dots, m \\
& h_j(x) = 0, \quad j = 1, \dots, p \\
& x \in \mathbb{R}^n.
\end{array}
\qquad (C.1) $$

For convenience, we denote the optimal value for (C.1) by $p^*$. Whenever $c(x)$, $f_i(x)$ are convex functions and $h_j(x)$ are linear functions the problem (C.1) is defined to be convex, otherwise it is non-convex.<fnref n="35" />

For a general problem of the form (C.1), convex or not, we can construct its dual problem starting with the Lagrangian defined by:

$$ \mathcal{L}(x, u, v) = c(x) + \sum_{i=0}^{m} u_i f_i(x) + \sum_{j=0}^{p} v_j h_j(x), \quad u_i \ge 0 \qquad (C.2) $$

and minimize over $x$:<fnref n="36" />

$$ g(u, v) = \min_{x \in \mathbb{R}^n} \mathcal{L}(x, u, v). \qquad (C.3) $$

It is a simple exercise to show that $g(u, v)$ is concave, since it is a minimization over a family of linear functions in $u, v$. For all $u, v$ under the constraints $u_i \ge 0$, we have:

$$ g(u, v) \le \max g(u, v) \le p^*. \qquad (C.4) $$

This indicates us the formulation of the dual problem:

$$ 
\begin{array}{ll}
\max & g(u, v) \\
\text{s.t.} & u_i \ge 0 \quad i = 1, \dots, m.
\end{array}
\qquad (C.5) $$

We denote the optimal solution of this problem to be $d^*$. This optimization problem is guaranteed to be convex, since it is maximizing a concave function $g(u, v)$. In this sense the dual problem (C.5) is simpler when the primal problem is non-convex. Of course we always have:

$$ d^* \le p^* \qquad (C.6) $$

due to (C.4). The inequality (C.6) is conventionally called weak duality. It would be actually nice to have the equality, i.e. when strong duality holds. For that case we have the well-known Slater's condition:

The strong duality holds if the primal problem is convex and it has a strictly feasible solution.

By definition, a solution $x^*$ is strictly feasible when:

$$ 
\begin{array}{ll}
f_i(x^*) < 0 & i = 1, \dots, m \\
h_j(x^*) = 0, & j = 1, \dots, p.
\end{array}
\qquad (C.7) $$

In general, we don't have strong duality for non-convex problem.

<footnote n="35">³⁵It turns out that usually we can only solve the convex optimization problem efficiently, and non-convex problems are generally NP-hard.</footnote>

<footnote n="36">³⁶We stress here that the minimization is the unconstrained minimization, i.e. $x \in \mathbb{R}^n$.</footnote>
## C.2 Relaxation problem and dual problem

In this section, for completeness we present a proof that the relaxation problem introduced in (4.21) is the dual of the dual problem of the original problem (4.9). The proof is actually trivial but lengthy, so the reader could treat it as an implementation example of the dual formulation introduced in the last section.

First we transform the original problem (4.9) to the following form:

$$
\begin{array}{ll}
\min & c^\mathrm{T}x \\
\text{s.t.} & \mathrm{tr}XA_i + b_i^\mathrm{T}x + a_i = 0, \\
& M_0 + \sum_{j=1}^L M_j x_j \ge 0, \\
& X = xx^\mathrm{T}.
\end{array}
\qquad (C.8)
$$

To take the dual problem of it, we write down its Lagrangian:

$$
\mathcal{L}(x, X, \lambda, Y, Z) = c^{\mathrm{T}}x + \sum_{i} \lambda_{i} (\mathrm{tr}XA_{i} + b_{i}^{\mathrm{T}}x + a_{i}) - \mathrm{tr}\left(Y\left(M_{0} + \sum_{j=1}^{L} M_{j}x_{j}\right)\right) \\
\qquad + \mathrm{tr}\left(Z(xx^{\mathrm{T}} - X)\right). \qquad (C.9)
$$

The matrices $Y$ and $Z$ introduced here are real symmetric matrix variables, satisfying $Y \succeq 0$ due to the inequality condition in (C.8). To minimize the Lagrangian, we collect all the terms involving primal variables $x$ and $X$:

$$
\mathcal{L}(x, X, \lambda, Y, Z) = \operatorname{tr}\left(\left(\sum_{i} \lambda_{i} A_{i}-Z\right) X\right)+x^{\mathrm{T}} Z x+\left(c+d+\sum_{i} \lambda_{i} b_{i}\right)^{\mathrm{T}} x \\
+ \sum_{i} \lambda_{i} a_{i}-\operatorname{tr}\left(Y M_{0}\right)
\qquad (C.10)
$$

where we introduced the vector variable $d_k = -\mathrm{Tr}(YM_k)$ to make the formula more compact. We get $g(\lambda, Y, Z)$ by taking the minimization over $x, X$:

$$
\begin{array}{l}
g(\lambda, Y, Z) = \sum_{i} \lambda_{i} a_{i} - \mathrm{tr}(YM_{0}) - \frac{1}{4} \left( c + d + \sum_{i} \lambda_{i} b_{i} \right)^{\mathrm{T}} Z^{\dagger} \left( c + d + \sum_{i} \lambda_{i} b_{i} \right), \\
\text{if } \sum_{i} \lambda_{i} A_{i} - Z = 0, \quad Z \succeq 0, \quad (I - ZZ^{\dagger}) \left( c + d + \sum_{i} \lambda_{i} b_{i} \right) = 0. \qquad (C.11)
\end{array}
$$

The reader can verify that if the conditions in the second line of (C.11) are not saturated, the minimal value $g(\lambda, Y, Z)$ is $-\infty$, which is irrelevant since we are only interested in its the maximal. Note that $Z^\dagger$ here denotes the pseudo-inverse of the matrix $Z$:

$$
Z^{\dagger} = \lim_{\epsilon \to 0^{+}} (Z^{\mathrm{T}}Z + \epsilon I)^{-1} Z^{\mathrm{T}}. \qquad (C.12)
$$
Introducing an auxiliary variable $\gamma$ by Schur complement, we can formulate the dual problem in a more compact form:

$$ \begin{array}{ll} \max & \gamma \\ \text{s.t.} & \displaystyle\sum_i \lambda_i A_i - Z = 0, \quad Y \ge 0, \\ & \left( \begin{array}{cc} \sum_i \lambda_i a_i - \text{tr}(YM_0) - \gamma & (c+d+\sum_i \lambda_i b_i)^\text{T}/2 \\ (c+d+\sum_i \lambda_i b_i)/2 & Z \end{array} \right) \ge 0. \end{array} \qquad (C.13) $$

To do the second dualization, we introduce again the Lagrangian:

$$ \mathcal{L}_2(x, X, \delta, W, S, \gamma, \lambda, Y, Z) = -\gamma - \text{tr}(SY) + \text{tr}\left(W\left(\sum_i \lambda_i A_i - Z\right)\right) \\ -\text{tr}\left(\begin{pmatrix} \delta & x^\text{T} \\ x & X \end{pmatrix} \begin{pmatrix} \sum_i \lambda_i a_i - \text{tr}(YM_0) - \gamma & (c+d+\sum_i \lambda_i b_i)^\text{T}/2 \\ (c+d+\sum_i \lambda_i b_i)/2 & Z \end{pmatrix}\right) \qquad (C.14) $$

Here we slightly abuse the notation by introducing a dual variable contains $x$ and $X$, which satisfies:

$$ \begin{pmatrix} \delta & x^\text{T} \\ x & X \end{pmatrix} \ge 0. \qquad (C.15) $$

We also note that $S \succeq 0$. The dual problem becomes apparent when we collect the coefficients of $\gamma, Y, Z, \gamma$:

$$ \mathcal{L}_2(x, X, \delta, W, S, \gamma, \lambda, Y, Z) = (\delta - 1)\gamma + \text{tr}\left(\left(M_0\delta - S + \sum_k M_k x_k\right)Y\right) - \text{tr}\left((W+X)Z\right) \\ + \sum_i \lambda_i \left(\text{tr}(WA_i) - a_i\delta - b_i^\text{T}x\right) - c^\text{T}x. \qquad (C.16) $$

Again, all the coefficients of the linear terms must vanish, otherwise the minimal value of $\mathcal{L}_2$ will be $-\infty$ when we take the minimization. So this time we can write down the dual problem as:

$$ \begin{array}{ll} \min & c^\text{T}x, \\ \text{subject to} & M_0 + \sum_k M_k x_k \ge 0 \\ & \text{tr}(XA_i) + b_i^\text{T}x + a_i = 0, \\ & \begin{pmatrix} 1 & x^\text{T} \\ x & X \end{pmatrix} \ge 0 \end{array} \qquad (C.17) $$

which is exactly our original problem with relaxation.

## D Implementation in SDPA

Here we gather the details on the implementation of our methods which are not covered in section 5. We stress that although we are quite satisfied by the current efficiency
of the method, we cannot guarantee that the choices we made are the best possible. In appendix E, we will demonstrate explicitly the importance of all these details on a simple example.

## D.1 Choice of SDP solvers

The SDP solver for all of the numerical results in section 5 is SDPA [51]. In the search for a better performance, we tried a lot of other solvers, including PENLAB, SDPB, SDPA, MOSEK, SDPT3 and SeDuMi. In the early stage of this work, we were mainly working with SDPB [52]. Our experience shows that SDPB gives the most stable convergence in all situations, even when we set the SPDB precision to 64. All other solvers sometimes fail to converge. The disadvantage of SDPB is that it is much slower and uses much more memory compared to all other solvers. This is due to the fact that SDPB is an arbitrary precision solver relying on GMP, which potentially slows the program by more than 1000 times. Due to this reason, when we try to deal with large-scale problems, we switch to SDPA which, according to our limited tests, converges better than all other SDP solvers on large-scale problems.

## D.2 Generation of the loop equations

For loop equations of a given length $l$, the set of all loop equations could be generated by exhaustive method.<fnref n="37" /> This is conceptually simple and numerically easy to implement. In fact, a slight optimization of Lin's code in [28] is efficient enough to generate loop equations up to $\Lambda = 11$ in our work. The tricky part here is to determine what is the sufficient length cutoff of loop equations for our numerics.

Let's take the model (5.1) as an example again. Suppose we are considering a cutoff $\Lambda$, which means that the longest operator in the correlation matrix is $2\Lambda$. We recall that in loop equation of length $l$ in this model, we have the operator expectations of lengths $l-1, l+1, l+3$.<fnref n="38" /> We can consider only the loop equations up to the length $l = 2\Lambda-3$, since they contain only the operators we are interested in. But higher length loop equations could also generate non-trivial relations on the operator expectations of lengths within the cutoff. Namely, as loop equations of higher length contain lower operators that we are interested in, we could eliminate higher operators and get new, generally linear independent loop equations containing only operators up to the length $l = 2\Lambda$. However, empirically for the model (5.1) this doesn't happen: the procedure described above does not produce any new independent loop equations, in addition to the standard ones, constraining the operators of lengths $l \le 2\Lambda - 3$. One exception is the case $h = 0$ and $g = 0$, i.e. when the model is solvable. Then we need to consider the loop equations up to the lengths $l = 2\Lambda + 1$ and $l = 2\Lambda - 1$ to generate all the equations between operators we are interested in. We don't have any proof why this is true. But this is not hard to understand intuitively. For example, when $g = 0$ and we consider the loop equations only up to $l = 2\Lambda - 3$, so we have a smaller number of length $2\Lambda$ operators comparing to the case $g \ne 0$ in the loop equations.

<footnote n="37">³⁷Here exhaustive method means that we generate all the loop equations by taking matrix derivatives on all possible positions of all possible words in Schwinger-Dyson equations.</footnote>

<footnote n="38">³⁸Here we define the length of quadratic terms in loop equations as the sum of the lengths of their factors.</footnote>
But we also note that even for $h = 0$ or $g = 0$, when we consider only the loop equation up to length $l = 2\Lambda - 3$ in our bootstrap assumptions, we still get a bound very close to the one using the higher length loop equations as well. This example supports our general intuition that higher loop equations are expected to have a relatively limited influence on small length operators.

In conclusion, for the model (5.1) with cutoff $2\Lambda$, it is enough to take into account the loop equations up to the length $l = 2\Lambda - 3$.

### D.3 Solving the loop equations

Unlike other solvers, the input form of SDPA is rather demanding. It only admits an input problem strictly of the (4.1) form. As a consequence, we need to solve the loop equation and substitute the solution into the correlation matrix and relaxation matrix. More precisely, by solving the loop equations we mean that there is a reduced set of operator expectations, such that all other operator expectations and all the quadratic terms $R_{ij}$ can be solved by loop equations as linear combinations of them. For later convenience, we denote the vector of this reduced set by $x_{\text{red}}$.

A bad news is that the loop equations generated by exhaustive method are not linearly independent, so solving it with floating point numbers will cause some numerical instabilities. One way out of this difficulty is to remove the linear redundancies in loop equations. Another way is to use exact numbers in Mathematica. The first is hard since we found that for some specific values of $h$ and $g$, there will be coincident degeneracies. So we doubt we could have a systematic way to select the linear independent loop equations. Actually we are using the second way, i.e. only using the exact number computation when generating and solving loop equations. The exact number arithmetics consumes a lot of memory and makes the solving procedure much slower, but up to the most involved case $\Lambda = 11$ in this article, we still have enough of memory(around 200 Gigabytes) and the solving functions are efficient enough.

### D.4 Relaxation matrix

There are some technical details concerning generating the relaxation matrix constraints. First we note that if we are considering loop equations up to $l = 2\Lambda - 3$, they contain quadratic terms only up to length $2\Lambda - 4$. The operators in the bottom right corner of the relaxation matrix $R$ are neither constrained at all, nor contained in the objective function (cf (4.20)). So they are actually irrelevant. Formally these can be shown by Schur complement:<fnref n="39" />

$$ \mathcal{R} = \begin{pmatrix} X_1 & B_1 \\ B_1^T & X_2 \end{pmatrix} \succ 0 \Leftrightarrow X_1 \succ 0 \land X_2 - B_1^T X_1^{-1} B_1 \succ 0. \quad (\text{D.1}) $$

Here $X_1$ is the sub-matrix contain the product of $\langle \text{Tr} O_i \rangle \langle \text{Tr} O_j \rangle$, where the length of $O_i$ plus $O_j$ is less or equal than $\Lambda - 2$. Since $X_2$ is not constrained at all, the second inequality

<footnote n="39">³⁹Numerically there is no difference between strict inequalities and non-strict ones.</footnote>
can always be satisfied. In this way, the positivity condition of $\mathcal{R}$ is equivalent to the positivity of the $X_1$.

The matrix $X_1$ needs to be further reduced. Since we know that the operator expectations are linearly dependent on the loop equations, it is not surprising that the quadratic terms are linearly dependent. We must reduce this linear dependency, otherwise it will induce numerical instabilities for SDPA. To solve this potential problem, we consider the minor $X_{1\text{red}}$ in the form:

$$X_{1\text{red}} = x_{1\text{red}}x_{1\text{red}}^{\text{T}}. \qquad (\text{D.2})$$

The vector $x_{1\text{re}}$ consists of operator expectations $\langle \text{Tr}\mathcal{O}_i \rangle$, such that not only the lengths of $\mathcal{O}_i$ are equal or less than $\Lambda - 2$ but also they are in the set $x_{\text{red}}$ derived by solving the loop equations. Then we impose the positivity condition on the reduced relaxation matrix:

$$\mathcal{R}_{1\text{red}} = \begin{pmatrix} 1 & x_{1\text{red}}^{\text{T}} \\ x_{1\text{red}} & X_{1\text{red}} \end{pmatrix} \succeq 0. \qquad (\text{D.3})$$

## D.5 Feasibility

In section 5.3, in the bisection process we need to test whether the SDP is feasible. However, SDPA doesn't have a built-in option to test the feasibility of a SDP. To deal with this difficulty, we introduce a slack variable $\mu$ in addition to the original variables in the SDP, to transform the original problem (4.21) into the form:

$$\begin{array}{ll}
\min & \mu \\
\text{subject to} & M_0 + \sum_k M_k x_k + \mu I_1 \succeq 0, \\
& \text{and } \text{tr}(XA_i) + b_i^{\text{T}}x + a_i = 0, \\
& \text{and } \begin{pmatrix} 1 & x^{\text{T}} \\ x & X \end{pmatrix} + \mu I_2 \succeq 0.
\end{array}
\qquad (\text{D.4})$$

Here $I_1$ and $I_2$ are the identity matrices of the appropriate size. If the optimal value is negative, then the original problem (4.21) is feasible, otherwise it is infeasible.

## D.6 Normalization

When we are dealing with a large-scale SDP, a common problem is that our input data is badly-scaled. We usually need to adjust our normalization to make all the numbers in the optimization problem to be of similar order of magnitude. In our problem, we make the substitution:

$$\mathcal{O}'_l = g_s^{l/2} \mathcal{O}_l \qquad (\text{D.5})$$

and then we bootstrap the primed operator instead of the original one. Here $l$ is the length of the operator. When $g$ and $h$ are not vanishing we put $g_s = \min\{|g|, 2|h|\}$. When it is vanishing, we take $g_s$ to be the absolute value of the non-vanishing coupling.

This naive choice to scale is not guaranteed to be the best for the problem, though it appears to lead to decent results for our range of cutoffs. It may be more advantageous to scale this problem according to the asymptotic behavior of large length operator expectations, which would be interesting to study analytically.
## D.7 Precision of the solver

Finally, we want to make some comments on the numerical precision. Most of the SDP solvers including SDPA are machine precision solvers. Numerical instabilities may happen when we are too ambitious about the precision of the results or the input data is badly scaled. In this situation, one could try arbitrary precision solvers like SDPB and SDPA-GMP, or other solvers in SDPA family called SDPA-DD and SDPA-QD, which is actually our recommendation. SDPA-DD and SDPA-QD is based on double-double and quad-double data type in QD library. From our tests, 32 digits or 64 digits are always enough for our purpose, and they are still efficient enough to solve the problem.

### E $\Lambda = 4$ example

In this appendix, we demonstrate the numerical implementation of SDPA introduced in the appendix D on the simplest non-linear case $\Lambda = 4$ of the model (5.1), under the assumption of $\mathbb{Z}_2^3$ symmetry.

### E.1 Operators and loop equations

Before generating the loop equations, one of the preparatory work is to generate all the nonequivalent operators and possible quadratic terms up to the cutoff $\Lambda = 4$. As already discussed in section 5, we have only very few operators that are non-equivalent up to identifications due to the symmetries of the problem and of the class of solutions considered here.40 There are 20 operators with the length smaller or equal to 8:

$$ \begin{aligned} & \mathrm{Tr} A^2, \mathrm{Tr} A^4, \mathrm{Tr} A^2 B^2, \mathrm{Tr} ABAB, \mathrm{Tr} A^6, \mathrm{Tr} A^4 B^2, \mathrm{Tr} A^3 BAB, \mathrm{Tr} A^2 BA^2 B, \mathrm{Tr} A^8, \\ & \mathrm{Tr} A^6 B^2, \mathrm{Tr} A^5 BAB, \mathrm{Tr} A^4 BA^2 B, \mathrm{Tr} A^4 B^4, \mathrm{Tr} A^3 BA^3 B, \mathrm{Tr} A^3 BAB^3, \mathrm{Tr} A^3 B^2 AB^2, \\ & \mathrm{Tr} A^2 BABAB^2, \mathrm{Tr} A^2 BAB^2 AB, \mathrm{Tr} A^2 B^2 A^2 B^2, \mathrm{Tr} ABABABABAB. \end{aligned} \quad (E.1) $$

All other operators are identical to those in this list or they vanish under the $\mathbb{Z}_2^3$ symmetry assumption. There is only one quadratic term left under this assumption. For conciseness, we denote it by $\beta$ in this appendix:

$$ \beta = (\mathrm{Tr} A^2)^2 = (\mathrm{Tr} B^2)^2 = \mathrm{Tr} A^2 \mathrm{Tr} B^2. \quad (E.2) $$

<footnote n="40">40 We slightly abuse the notations in this appendix: all single trace operators in this section actually mean their expectation values. For example $\mathrm{Tr} A^2$ is actually $\langle \mathrm{Tr} A^2 \rangle$.</footnote>
To generate the loop equations, we simply apply the exhaustive method and then delete the duplicates. For general couplings $h$ and $g$, we have 14 loop equations left:

$$1 = \mathrm{Tr}A^2 + g\mathrm{Tr}A^4 - h(-2\mathrm{Tr}A^2B^2 + 2\mathrm{Tr}ABAB)$$

$$0 = -2\mathrm{Tr}A^2 + \mathrm{Tr}A^4 - h(2\mathrm{Tr}A^3BAB - 2\mathrm{Tr}A^4B^2) + g\mathrm{Tr}A^6$$

$$0 = -\mathrm{Tr}A^2 + \mathrm{Tr}A^2B^2 - h(-\mathrm{Tr}A^2BA^2B + 2\mathrm{Tr}A^3BAB - \mathrm{Tr}A^4B^2) + g\mathrm{Tr}A^4B^2$$

$$0 = -h(2\mathrm{Tr}A^2BA^2B - 2\mathrm{Tr}A^3BAB) + g\mathrm{Tr}A^3BAB + \mathrm{Tr}ABAB$$

$$\beta = -2\mathrm{Tr}A^4 + \mathrm{Tr}A^6 - h(2\mathrm{Tr}A^5BAB - 2\mathrm{Tr}A^6B^2) + g\mathrm{Tr}A^8$$

$$\beta = -\mathrm{Tr}A^2B^2 + \mathrm{Tr}A^4B^2 - h(-\mathrm{Tr}A^3B^2AB^2 + 2\mathrm{Tr}A^3BAB^3 - \mathrm{Tr}A^4B^4) + g\mathrm{Tr}A^6B^2$$

$$0 = -2\mathrm{Tr}A^2B^2 - h(-\mathrm{Tr}A^2B^2A^2B^2 + 2\mathrm{Tr}A^2BABAB^2 - \mathrm{Tr}A^3B^2AB^2) + \mathrm{Tr}A^4B^2 + g\mathrm{Tr}A^6B^2$$

$$0 = -\mathrm{Tr}A^4 + \mathrm{Tr}A^4B^2 + g\mathrm{Tr}A^4B^4 - h(-\mathrm{Tr}A^4BA^2B + 2\mathrm{Tr}A^5BAB - \mathrm{Tr}A^6B^2)$$

$$0 = \mathrm{Tr}A^3BAB - h(2\mathrm{Tr}A^2BAB^2AB - \mathrm{Tr}A^2BABAB^2 - \mathrm{Tr}A^3BAB^3) + g\mathrm{Tr}A^5BAB - \mathrm{Tr}ABAB$$

$$0 = \mathrm{Tr}A^3BAB + g\mathrm{Tr}A^5BAB - 2\mathrm{Tr}ABAB - h(-2\mathrm{Tr}A^2BABAB^2 + 2\mathrm{Tr}ABABABAB)$$

$$0 = \mathrm{Tr}A^3BAB + g\mathrm{Tr}A^3BAB^3 - h(-\mathrm{Tr}A^3BA^3B + 2\mathrm{Tr}A^4BA^2B - \mathrm{Tr}A^5BAB)$$

$$0 = g\mathrm{Tr}A^3BA^3B + \mathrm{Tr}A^3BAB - h(2\mathrm{Tr}A^3B^2AB^2 - 2\mathrm{Tr}A^3BAB^3)$$

$$0 = -\mathrm{Tr}A^2B^2 + \mathrm{Tr}A^2BA^2B - h(-\mathrm{Tr}A^2BAB^2AB + 2\mathrm{Tr}A^2BABAB^2 - \mathrm{Tr}A^3B^2AB^2) + g\mathrm{Tr}A^4BA^2B$$

$$\beta = \mathrm{Tr}A^2BA^2B + g\mathrm{Tr}A^3B^2AB^2 - h(2\mathrm{Tr}A^3BA^3B - 2\mathrm{Tr}A^4BA^2B). \qquad (E.3)$$

This is a system of 14 linear equations for 21 variables. For generic values of $h$ and $g$, like the one we often chose in this paper $h = g = 1$, they are all linearly independent. So we can express 14 variables including $\beta$ through seven variables of shortest lengths. These seven variables form the subset $x_{\mathrm{red}}$ introduced in appendix D.4. For $h = g = 1$, we can take it as:

$$x_{\mathrm{red}} = (\mathrm{Tr}A^2, \mathrm{Tr}A^4, \mathrm{Tr}A^2B^2, \mathrm{Tr}A^6, \mathrm{Tr}A^8, \mathrm{Tr}A^6B^2, \mathrm{Tr}A^5BAB)^{\mathrm{T}}. \qquad (E.4)$$

The other operators, including $\beta$, can be expressed as linear combinations of these variables:

$$\mathrm{Tr}ABAB = \frac{1}{2}\mathrm{Tr}A^2 + \frac{1}{2}\mathrm{Tr}A^4 + \mathrm{Tr}A^2B^2 - \frac{1}{2}$$

$$\mathrm{Tr}A^4B^2 = \frac{1}{6}\mathrm{Tr}A^2 - \mathrm{Tr}A^2B^2 + \frac{1}{6}\mathrm{Tr}A^6 + \frac{1}{6}$$

$$\mathrm{Tr}A^3BAB = -\frac{5}{6}\mathrm{Tr}A^2 + \frac{1}{2}\mathrm{Tr}A^4 - \mathrm{Tr}A^2B^2 + \frac{2}{3}\mathrm{Tr}A^6 + \frac{1}{6}$$

$$\mathrm{Tr}A^2BA^2B = -\mathrm{Tr}A^2 + \mathrm{Tr}A^4 - \mathrm{Tr}A^2B^2 + \mathrm{Tr}A^6$$
$$ \begin{aligned} \mathrm{Tr} A^4 B A^2 B = & -\frac{8}{3}\mathrm{Tr} A^2 + 9\mathrm{Tr} A^4 - \frac{14}{3}\mathrm{Tr} A^2 B^2 - \frac{1}{3}\mathrm{Tr} A^6 - \frac{8}{3}\mathrm{Tr} A^8 - \frac{16}{3}\mathrm{Tr} A^6 B^2 \\ & + \frac{28}{3}\mathrm{Tr} A^5 B A B + \frac{1}{3} \end{aligned} $$

$$ \begin{aligned} \mathrm{Tr} A^4 B^4 = & \frac{5}{2}\mathrm{Tr} A^2 - 8\mathrm{Tr} A^4 + \frac{17}{3}\mathrm{Tr} A^2 B^2 + \frac{1}{6}\mathrm{Tr} A^6 + \frac{8}{3}\mathrm{Tr} A^8 + \frac{13}{3}\mathrm{Tr} A^6 B^2 \\ & - \frac{22}{3}\mathrm{Tr} A^5 B A B - \frac{1}{2} \end{aligned} $$

$$ \begin{aligned} \mathrm{Tr} A^3 B A^3 B = & -\frac{9}{2}\mathrm{Tr} A^2 + \frac{31}{2}\mathrm{Tr} A^4 - \frac{23}{3}\mathrm{Tr} A^2 B^2 - \frac{2}{3}\mathrm{Tr} A^6 - \frac{14}{3}\mathrm{Tr} A^8 - \frac{28}{3}\mathrm{Tr} A^6 B^2 \\ & + \frac{46}{3}\mathrm{Tr} A^5 B A B + \frac{1}{2} \end{aligned} $$

$$ \mathrm{Tr} A^3 B A B^3 = 2\mathrm{Tr} A^4 - \frac{2}{3}\mathrm{Tr} A^2 B^2 - \frac{2}{3}\mathrm{Tr} A^6 - \frac{2}{3}\mathrm{Tr} A^8 - \frac{4}{3}\mathrm{Tr} A^6 B^2 + \frac{7}{3}\mathrm{Tr} A^5 B A B $$

$$ \begin{aligned} \mathrm{Tr} A^3 B^2 A B^2 = & -\frac{8}{3}\mathrm{Tr} A^2 + 10\mathrm{Tr} A^4 - 5\mathrm{Tr} A^2 B^2 - \frac{2}{3}\mathrm{Tr} A^6 - 3\mathrm{Tr} A^8 - 6\mathrm{Tr} A^6 B^2 \\ & + 10\mathrm{Tr} A^5 B A B + \frac{1}{3} \end{aligned} $$

$$ \begin{aligned} \mathrm{Tr} A^2 B A B A B^2 = & -\frac{14}{3}\mathrm{Tr} A^2 + 14\mathrm{Tr} A^4 - \frac{26}{3}\mathrm{Tr} A^2 B^2 - 4\mathrm{Tr} A^8 - 8\mathrm{Tr} A^6 B^2 \\ & + 14\mathrm{Tr} A^5 B A B + \frac{2}{3} \end{aligned} $$

$$ \begin{aligned} \mathrm{Tr} A^2 B A B^2 A B = & -3\mathrm{Tr} A^2 + 8\mathrm{Tr} A^4 - \frac{17}{3}\mathrm{Tr} A^2 B^2 - \frac{7}{3}\mathrm{Tr} A^8 - \frac{14}{3}\mathrm{Tr} A^6 B^2 \\ & + \frac{26}{3}\mathrm{Tr} A^5 B A B + \frac{2}{3} \end{aligned} $$

$$ \begin{aligned} \mathrm{Tr} A^2 B^2 A^2 B^2 = & -\frac{41}{6}\mathrm{Tr} A^2 + 18\mathrm{Tr} A^4 - \frac{28}{3}\mathrm{Tr} A^2 B^2 + \frac{1}{2}\mathrm{Tr} A^6 - 5\mathrm{Tr} A^8 - 11\mathrm{Tr} A^6 B^2 \\ & + 18\mathrm{Tr} A^5 B A B + \frac{5}{6} \end{aligned} $$

$$ \begin{aligned} \mathrm{Tr} A B A B A B A B = & -\frac{67}{12}\mathrm{Tr} A^2 + \frac{55}{4}\mathrm{Tr} A^4 - \frac{61}{6}\mathrm{Tr} A^2 B^2 + \frac{1}{3}\mathrm{Tr} A^6 - 4\mathrm{Tr} A^8 - 8\mathrm{Tr} A^6 B^2 \\ & + \frac{29}{2}\mathrm{Tr} A^5 B A B + \frac{5}{4} \end{aligned} $$

$$ \beta = -2\mathrm{Tr} A^4 + \mathrm{Tr} A^6 + \mathrm{Tr} A^8 + 2\mathrm{Tr} A^6 B^2 - 2\mathrm{Tr} A^5 B A B. \qquad (E.5) $$

## E.2 Correlation matrix and relaxation matrix

As we discussed in section 5, under the $Z_2^3$ symmetry our correlation matrix decouples into a block-diagonal matrix with three blocks. They are, respectively, the inner product$^{41}$ matrix of even-even words:

$$ I, AA, BB, AAAA, AABB, ABAB, ABBA, BAAB, BABA, BBAA, BBBB \qquad (E.6) $$

odd-odd words:

$$ AB, BA, AAAB, AABA, ABAA, ABBB, BAAA, BABB, BBAB, BBBA \qquad (E.7) $$

<footnote>$^{41}$Here inner product of $O_1$ and $O_2$ is defined to be $\langle \mathrm{Tr} O_1^\dagger O_2 \rangle$.</footnote>
and even-odd words:

$$ B, AAB, ABA, BAA, BBB. \qquad (E.8) $$

For example, the block for the even-odd words reads:

$$ \left( \begin{array}{cccccc} \mathrm{Tr} A^2 & \mathrm{Tr} A^4 & \mathrm{Tr} A^2 B^2 & \mathrm{Tr} A B A B & \mathrm{Tr} A^2 B^2 & \\ \mathrm{Tr} A^4 & \mathrm{Tr} A^6 & \mathrm{Tr} A^4 B^2 & \mathrm{Tr} A^3 B A B & \mathrm{Tr} A^4 B^2 & \\ \mathrm{Tr} A^2 B^2 & \mathrm{Tr} A^4 B^2 & \mathrm{Tr} A^4 B^2 & \mathrm{Tr} A^3 B A B & \mathrm{Tr} A^2 B A^2 B & \\ \mathrm{Tr} A B A B & \mathrm{Tr} A^3 B A B & \mathrm{Tr} A^3 B A B & \mathrm{Tr} A^2 B A^2 B & \mathrm{Tr} A^3 B A B & \\ \mathrm{Tr} A^2 B^2 & \mathrm{Tr} A^4 B^2 & \mathrm{Tr} A^2 B A^2 B & \mathrm{Tr} A^3 B A B & \mathrm{Tr} A^4 B^2 & \end{array} \right) \qquad (E.9) $$

It is easy to construct the relaxation matrix for this example using the explanations of appendix D.4. It is:

$$ \left( \begin{array}{cc} 1 & \mathrm{Tr} A^2 \\ \mathrm{Tr} A^2 & \beta \end{array} \right) \succeq 0. \qquad (E.10) $$

We note that in our current setting the vector $x_{1\mathrm{red}}$ is a single component vector with the component $\mathrm{Tr} A^2$.

To turn the original problem into the form (4.1), we substitute the “solution” of the loop equations (E.5) into the correlation matrix and the relaxation matrix. As for the objective function, we choose it to minimize $\mathrm{Tr} A^2$ and $-\mathrm{Tr} A^2$ to find the minimal and the maximal value of $\mathrm{Tr} A^2$, respectively. In this way, we can get the allowed region of $\mathrm{Tr} A^2$. At the next step we generate the input file for SDPA and solve it. With the appropriate setup, the entire time consumed for generating the input and solving it should last less than 0.1s CPU time. The result for $\Lambda = 4$ bootstrap is then:

$$ 0.393566 \leq \mathrm{Tr} A^2 \leq 0.431148. \qquad (E.11) $$

# F Structure of loop equations and solvable 2-matrix model

Generalizing the results derived in section 3.1 to multi-matrix model is far from straightforward. In one-matrix model, all the information contained in the loop equations and the positivity conditions can be encoded compactly into the resolvent function. On the contrary, for the multi-matrix model, as the loop equations and the correlation matrix both get much more involved, we don't expect to have such an analytic function enclosing the information of all the moments. Due to this complication, the bootstrap problem for a general multi-matrix model is generally not exactly solvable. In this appendix, we will discuss the nature of these complications in the structure of the loop equations and remind, from this point of view, an old result for the simplest 2-matrix model with $\mathrm{tr}(AB)$ interaction when the loop equations greatly simplify [43, 53].

## F.1 Base moments

As demonstrated in section 2 by the loop equations, all higher moments of one-matrix model are fully determined by a fixed number of lower moments (which we will call the base moments). But this is not generally true in multi-matrix model: namely, the number
of such base moments generally grows with the increase of the cutoff $\Lambda$. Let us take the model d (5.1) as an example. The following results can be observed from our numerical investigation:<fnref n="42" />

1. In the simplest case $h = 0$ when the model effectively factorizes into two decoupled one-matrix models, all the moments can be expressed by a polynomial of $t_2 = \langle \mathrm{Tr} A^2 \rangle$ and $g$. The number of base moments is 1 here.<fnref n="43" />

2. For $g = 0$ which is also solvable, all the moments are fixed by $t_{2k} = \langle \mathrm{Tr} A^{2k} \rangle$, $k = 1, 2, \dots$. This is very different from the one-matrix model since here we have to specify the value of infinite number of moments to determine the value of the remaining moments. If we set a finite cutoff $2\Lambda$ to the length of the moments, we have a set of truncated set of base moments of the size $\Lambda$.

3. For the general parameters ($g \neq 0, h \neq 0$), there are much more base moments than the case $g = 0$.

The intuition here is that, for a given multi-matrix model, the number of the base moments is negatively related with the solvability of the multi-matrix model. For a given cutoff, the number of the base moments are relatively easy to calculate. So one of the possible application of this intuitive observation is that one can use the number of the base moments for a finite but high cutoff to predict the solvability of the model.

## F.2 Closed subset of loop equations

Sometimes, the loop equations can be closed on a proper, much reduced subset of all moments. This usually leads to a great simplification of the system of loop equations and potentially makes the bootstrap problem for the multi-matrix model exactly solvable. We demonstrate this on the simplest solvable 2-matrix model, with a long history of stydy and applications [15–17] and used in [28] to demonstrate the matrix bootstrap:

$$Z = \lim_{N \to \infty} \int d^{N^2} A d^{N^2} B e^{-N \mathrm{tr}(-AB+V(A)+V(B))}, \quad V(x) = g_2 x^2 / 2 + g_3 x^3 / 3. \quad (F.1)$$

We notice that the following subset of loop equations:

$$\begin{align}
-g_3 t_{n+1} &= g_2 t_n - t_{n-1,1} - \sum_{j=0}^{n-2} t_j t_{n-2-j}, \nonumber \\
-g_3 t_{n+1,1} &= g_2 t_{n,1} - t_{n-1,2} - \sum_{j=0}^{n-2} t_j t_{n-2-j,1}, \tag{F.2} \\
-g_3 t_{n-1,2} &= g_2 t_{n-1,1} - t_n, \nonumber
\end{align}$$

are closed among the operators $t_n, t_{n,1}, t_{n,2}$, here:

$$t_n = \langle \mathrm{Tr} A^n \rangle, \quad t_{n,m} = \langle \mathrm{Tr} A^n B^m \rangle. \quad (F.3)$$

<footnote n="42">⁴²At some finite but high cutoffs, but we strongly believe that it holds also for arbitrarily high cutoff.</footnote>

<footnote n="43">⁴³We assumed the global symmetry.</footnote>
Summing over the equations by the way in section 2.1, we can get the Master loop equation for model (F.1) [53]:

$$ (Y(z) - V'(z))(z - V'(Y(z))) + P(Y(z), z) = 0, \qquad (F.4) $$

where

$$ \begin{aligned} P(x, y) &= -\langle \text{Tr} \frac{V'(x) - V'(A)}{x - A} \frac{V'(y) - V'(B)}{y - B} \rangle + 1, \\ G(z) &= \langle \text{Tr} \frac{1}{z - A} \rangle = \sum_{i=0}^{\infty} z^{-i-1} t_i, \\ Y(z) &= V'(z) - G(z). \end{aligned} \qquad (F.5) $$

Since this is a closed subset of loop equations, we assume that the sub-correlation matrix defined by $T_{i,j} = t_{i+j-2}$ is positive semi-definite, which is the positivity condition for the minor of the whole correlation matrix consisting of the elements $t_n$. This brings us back to the one-matrix type bootstrap problem considered in section 3.2. Here the positivity condition is equivalent to that the eigenvalue distribution corresponding to $G(z)$ is real and positive, or $Y(z)$ has a negative cut. In principle, this problem is analytically solvable, the complication compared with one-matrix model is that we have a cubic equation instead of a quadratic one (or an equation of nth degree for the potentials of order n). The solution of these loop equations has been found in [54] in terms of an algebraic curve depending on the base moments [55].

## Acknowledgments

We thank M.Paulos for a useful discussion. This work benefits a lot from the help from Walter Landry on SDPB in the early stage.

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

## References

[1] G. 't Hooft, *A Planar Diagram Theory for Strong Interactions*, Nucl. Phys. B 72 (1974) 461 [INSPIRE].

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

[3] F. David, *Planar Diagrams, Two-Dimensional Lattice Gravity and Surface Models*, Nucl. Phys. B 257 (1985) 45 [INSPIRE].

[4] V.A. Kazakov, A.A. Migdal and I.K. Kostov, *Critical Properties of Randomly Triangulated Planar Random Surfaces*, Phys. Lett. B 157 (1985) 295 [INSPIRE].

[5] V.A. Kazakov, *Bilocal Regularization of Models of Random Surfaces*, Phys. Lett. B 150 (1985) 282 [INSPIRE].
[6] V.A. Kazakov, *Exactly Solvable Potts Models, Bond And Tree Like Percolation On Dynamical (Random) Planar*, in *International Symposium on Field Theory of the Lattice*, Moscow, Russia, December 1, (1987), *Nucl. Phys. B Proc. Suppl. 4* (1988) 93 [INSPIRE].

[7] O. Bohigas, M.J. Giannoni and C. Schmit, *Characterization of chaotic quantum spectra and universality of level fluctuation laws*, *Phys. Rev. Lett. 52* (1984) 1 [INSPIRE].

[8] R. Dijkgraaf and C. Vafa, *Matrix models, topological strings, and supersymmetric gauge theories*, *Nucl. Phys. B 644* (2002) 3 [hep-th/0206255] [INSPIRE].

[9] R. Dijkgraaf, S. Gukov, V.A. Kazakov and C. Vafa, *Perturbative analysis of gauged matrix models*, *Phys. Rev. D 68* (2003) 045007 [hep-th/0210238] [INSPIRE].

[10] B. Eynard and N. Orantin, *Invariants of algebraic curves and topological expansion, Commun. Num. Theor. Phys. 1* (2007) 347 [math-ph/0702045] [INSPIRE].

[11] M. Kontsevich, *Intersection theory on the moduli space of curves and the matrix Airy function*, *Commun. Math. Phys. 147* (1992) 1 [INSPIRE].

[12] H.L. Montgomery, *The pair correlation of zeros of the zeta function*, in *Proc. Symp. Pure Math 24* (1973) 181.

[13] T. Eguchi and H. Kawai, *Reduction of Dynamical Degrees of Freedom in the Large N Gauge Theory*, *Phys. Rev. Lett. 48* (1982) 1063 [INSPIRE].

[14] E. Brézin, C. Itzykson, G. Parisi and J.B. Zuber, *Planar Diagrams*, *Commun. Math. Phys. 59* (1978) 35 [INSPIRE].

[15] C. Itzykson and J.B. Zuber, *The Planar Approximation. II*, *J. Math. Phys. 21* (1980) 411 [INSPIRE].

[16] M.L. Mehta, *A Method of Integration Over Matrix Variables*, *Commun. Math. Phys. 79* (1981) 327 [INSPIRE].

[17] V.A. Kazakov, *Exact Solution of the Ising Model on a Random Two-dimensional Lattice*, *JETP Lett. 44* (1986) 133 [INSPIRE].

[18] D.V. Boulatov and V.A. Kazakov, *The Ising Model on Random Planar Lattice: the Structure of Phase Transition and the Exact Critical Exponents*, *Phys. Lett. B 186* (1987) 379 [INSPIRE].

[19] I.K. Kostov, *O(n) Vector Model on a Planar Random Lattice: spectrum of Anomalous Dimensions*, *Mod. Phys. Lett. A 4* (1989) 217 [INSPIRE].

[20] J.-M. Daul, *Q states Potts model on a random planar lattice*, [hep-th/9502014] [INSPIRE].

[21] V.A. Kazakov and A.A. Migdal, *Recent Progress in the Theory of Noncritical Strings*, *Nucl. Phys. B 311* (1988) 171 [INSPIRE].

[22] V.A. Kazakov, *Solvable matrix models*, [hep-th/0003064] [INSPIRE].

[23] 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].

[24] 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].

[25] J.P. Rodrigues, *Numerical Solution of Lattice Schwinger-dyson Equations in the Large N Limit*, *Nucl. Phys. B 260* (1985) 350 [INSPIRE].
[26] 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].

[27] D. Simmons-Duffin, The Lightcone Bootstrap and the Spectrum of the 3d Ising CFT, JHEP 03 (2017) 086 [arXiv:1612.08471] [INSPIRE].

[28] H.W. Lin, Bootstraps to strings: solving random matrix models with positivity, JHEP 06 (2020) 090 [arXiv:2002.08387] [INSPIRE].

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

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

[31] 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].

[32] M.F. Paulos, J. Penedones, J. Toledo, B.C. van Rees and P. Vieira, The S-matrix bootstrap II: two dimensional amplitudes, JHEP 11 (2017) 143 [arXiv:1607.06110] [INSPIRE].

[33] D. Mazac and M.F. Paulos, The analytic functional bootstrap. Part I: 1D CFTs and 2D S-matrices, JHEP 02 (2019) 162 [arXiv:1803.10233] [INSPIRE].

[34] V.A. Kazakov, I.K. Kostov and N.A. Nekrasov, D particles, matrix integrals and KP hierarchy, Nucl. Phys. B 557 (1999) 413 [hep-th/9810035] [INSPIRE].

[35] V.A. Kazakov and P. Zinn-Justin, Two matrix model with ABAB interaction, Nucl. Phys. B 546 (1999) 647 [hep-th/9808043] [INSPIRE].

[36] I.K. Kostov, Exact solution of the six vertex model on a random lattice, Nucl. Phys. B 575 (2000) 513 [hep-th/9911023] [INSPIRE].

[37] P. Zinn-Justin, The Six vertex model on random lattices, EPL 50 (2000) 15 [cond-mat/9909250] [INSPIRE].

[38] R.G. Jha, Introduction to Monte Carlo for matrix models, SciPost Phys. Lect. Notes 46 (2022) 1 [arXiv:2111.02410] [INSPIRE].

[39] B. Eynard, Topological expansion for the 1-Hermitian matrix model correlation functions, JHEP 11 (2004) 031 [hep-th/0407261] [INSPIRE].

[40] P. Di Francesco, P.H. Ginsparg and J. Zinn-Justin, 2-D Gravity and random matrices, Phys. Rept. 254 (1995) 1 [hep-th/9306153] [INSPIRE].

[41] M. Reed and B. Simon, Fourier Analysis, Self-Adjointness. Volume II, Academic Press, Cambridge U.S.A. (1975) [ISBN: 9780125850025].

[42] V.A. Kazakov, The Appearance of Matter Fields from Quantum Fluctuations of 2D Gravity, Mod. Phys. Lett. A 4 (1989) 2125 [INSPIRE].

[43] M. Staudacher, Combinatorial solution of the two matrix model, Phys. Lett. B 305 (1993) 332 [hep-th/9301038] [INSPIRE].

[44] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, Cambridge U.K. (2004) [DOI].

[45] J.R. Hoppe, Quantum Theory of a Massless Relativistic Surface and a Two-Dimensional Bound State Problem, Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge U.S.A. (1982).
[46] S. El-Showk, M.F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin and A. Vichi, *Solving the 3D Ising Model with the Conformal Bootstrap*, Phys. Rev. D **86** (2012) 025022, [arXiv:1203.6064] [INSPIRE].

[47] Y. Shimamune, *On the Phase Structure of Large N Matrix Models and Gauge Models*, Phys. Lett. B **108** (1982) 407 [INSPIRE].

[48] M. Teper, *Large N*, PoS LATTICE2008 (2008) 022 [arXiv:0812.0085] [INSPIRE].

[49] E. Brézin and D.J. Gross, *The External Field Problem in the Large N Limit of QCD*, Phys. Lett. B **97** (1980) 120 [INSPIRE].

[50] I.K. Kostov, *Random Surfaces, Solvable Lattice Models And Discrete Quantum Gravity In Two-Dimension*, in XIX International Seminar on Theoretical Physics: nonperturbative Aspects of the Standard Model (GIFT Seminar), Saclay, France, June 1988 [INSPIRE].

[51] M. Yamashita, K. Fujisawa, M. Fukuda, K. Kobayashi, K. Nakata and M. Nakata, *Latest developments in the sdpa family for solving large-scale sdps*, in *Handbook on Semidefinite, Conic and Polynomial Optimization*, Springer US, Boston U.S.A (2012), pp. 687–713 [DOI].

[52] D. Simmons-Duffin, *A Semidefinite Program Solver for the Conformal Bootstrap*, JHEP **06** (2015) 174 [arXiv:1502.02033] [INSPIRE].

[53] B. Eynard, *Large N expansion of the 2 matrix model*, JHEP **01** (2003) 051 [hep-th/0210047] [INSPIRE].

[54] V.A. Kazakov and A. Marshakov, *Complex curve of the two matrix model and its tau function*, J. Phys. A **36** (2003) 3107 [hep-th/0211236] [INSPIRE].

[55] V.A. Kazakov and I.K. Kostov, *Instantons in noncritical strings from the two matrix model*, in *From Fields to Strings: Circumnavigating Theoretical Physics: Ian Kogan Memorial Collection*, M. Shifman et al. eds., World Scientific, Singapore (2005), pp. 1864–1894, DOI [hep-th/0403152] [INSPIRE].