Large N optimization for multi-matrix systems

In this work we revisit the problem of solving multi-matrix systems through numerical large N methods. The framework is a collective, loop space representation which provides a constrained optimization problem, addressed through master-field minimization. This scheme applies both to multi-matrix integrals (c = 0 systems) and multi-matrix quantum mechanics (c = 1 systems). The complete fluctuation spectrum is also computable in the above scheme, and is of immediate physical relevance in the later case. The complexity (and the growth of degrees of freedom) at large N have stymied earlier attempts and in the present work we present significant improvements in this regard. The (constrained) minimization and spectrum calculations are easily achieved with close to 104 variables, giving solution to Migdal-Makeenko, and collective field equations. Considering the large number of dynamical (loop) variables and the extreme nonlinearity of the problem, high precision is obtained when confronted with solvable cases. Through numerical results presented, we prove that our scheme solves, by numerical loop space methods, the general two matrix model problem.


Introduction
Large N multi-matrix problems are at the center of many theories of current interest, involving membranes [1], reduced super Yang-Mills theories [2][3][4][5][6], field theory of critical and noncritical strings [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22], phase transitions and black holes [23][24][25][26][27][28][29][30][31], and M-atrix theory [32]. At the same time, apart from very special cases, these systems are not solvable due to the fact that they are highly nonlinear. At finite N one has numerical Monte-Carlo methods [33], which have provided definite and most relevant results [32] with increase of simulations towards large N . The limiting theory, infinite N , features a rapid growth of JHEP01(2022)168 the degrees of freedom, represented by (Wilson) loop variables, which are physical, gauge invariant collective variables for the description of matrix and non-Abelian gauge theories. As is well known, (Wilson) loops become independent degrees of freedom at infinite N , and the exact theory is governed by non-linear Schwinger-Dyson (Migdal-Makeenko) equations [34], or alternatively in the collective field theory representation, in terms of an effective action [35] and/or a collective Hamiltonian [36]. The later provides a unified approach in which 1/N = G appears as a coupling constant, and has been used through the years for non perturbative studies. For coupled systems (of matrices) at quantum level, only simple systems are solvable. However, numerical approaches were developed in [37,38]. In these previous studies the nature of the (infinite N ) planar solution was understood as related to a constrained minimization problem, where a significant role is played by a set of inequalities associated with invariants (loops) in the collective description. An effective scheme for dealing with this constrained minimization was identified, through the use of master-field variables [38]. These methods were also seen to apply at sub-leading order in N and in particular, to consideration of the spectrum [39,40].
Recently, there is a renewed interest in large N optimization [41], with studies [42][43][44] that overlap with earlier work, and which re-discover the importance of loop space inequality constraints.
Due to potential high relevance in problems of emergent geometry, thermalization and black hole formation we revisit the earlier collective field constrained minimization and numerical master field methods, with interest in increasing the numbers of degrees of freedom, and the potential for high precision results. These are developed in the present work. For concreteness, and with purpose of making comparisons (with analytical results, when possible) we mostly deal with systems of two hermitian matrices, but the methods are seen to apply for any number, with arbitrary single or multi trace interactions.
The content of this paper is as follows: in the overview section 2 we provide a short summary of the collective large N method [38]. In particular, the infinite N constrained minimization scheme is summarized, featuring the associated complete set of (loop space) inequalities. We then give a summary of the unconstrained master variable method which applies for the large N ground state and also for spectrum studies [39,40]. We explain the study of multi-matrix integrals through the Hamiltonian description [45] given by Fokker-Planck. In section 3 we describe the numerical methods used for constrained minimization. In section 4 we give a list of models studied and present the numerical solution of the associated matrix integrals. In section 5 we present results, giving numerical solutions for one-and two-matrix large N quantum Hamiltonian problems. At the large N level these essentially correspond to giving a (numerical) solution of the fully nonlinear Schwinger-Dyson equations with approximately 10 4 loop variables. Correspondingly methods (and evaluation) of large N ground state energies, gaps and low lying spectra are also given. Conclusions and future applications to problems of interest are commented in section 6.

JHEP01(2022)168 2 Overview
In this article we will study and develop numerical techniques for solving the large N multimatrix theories. We will consider two classes of the multi-matrix systems. The first one is the multi-matrix integral (c = 0 systems) whose partition function reads where S(M ) is a multi trace action. The second one is matrix quantum mechanics (MQM) (c = 1 systems), whose dynamics is given by a Hamiltonian of the form where Π l is the canonical conjugate of M l . The large N expansion in the collective field formulation is developed after a change of variables from the original matrix valued variables to invariant variables, which we refer to as "loops". This terminology has its origins in lattice gauge theory, where the basic field degrees of freedom are unitary matrices U l , one for each link in the lattice. In that case the invariant variables, Wilson loops φ(C), are obtained by taking the trace of an ordered product of unitary matrices, one for each link of a closed path C. It is common to continue to refer to invariant variables as "loops" even for theories of hermitian matrices M l where invariants are given by the trace of products of the matrices. As an example, for the case of two hermitian matrices we have In this case the invariant φ(C) is not labeled by a closed path, but rather by specifying a word C in the alphabet of the matrices. The word specifies the order in which matrices are multiplied before tracing. The invariant loop variables are then described by all of the words with cycling identification. For example, Tr(M 1 M 1 M 2 ) is equal to Tr(M 1 M 2 M 1 ) and Tr(M 2 M 1 M 1 ) due to the cyclicity of trace, hence they all refer to the same invariant loop variable. It is in this sense that we use the loop terminology. For the purpose of counting the number of loops, it is useful to enumerate the loops with a permutation σ.
The loop corresponding to permutation σ is denoted φ σ . Consider loops built as single trace products of n M 1 s and m M 2 s. The permutation σ ∈ S n+m specifies the loop φ σ as follows The cycle structure of the permutation translates into the trace structure of φ σ . For a single trace σ must be an n + m cycle. This parametrization is not unique as distinct permutations do not necessarily define distinct loops. This follows by noting that In general, two loops φ σ 1 and φ σ 2 are the same if for some τ ∈ S n × S m . Removing this redundancy, we are left with a complete set of (infinitely many) loops. To develop some intuition for this description, note that if τ belongs to the cyclic group Z n+m generated by the n + m cycle given by (123 · · · n + m), then the equality (2.6) expresses nothing but the cyclicity of the trace. In general, we need to divide by more than just cyclicity and (2.6) is a convenient way to correctly account for all redundancies. A relevant class of observables for the c = 0 systems are provided by correlation functions of the loops These expectation values can be determined through equations of motion for the loop expectation values. The loop equation is derived as a Schwinger-Dyson equation for the loops. In the case of a theory of unitary matrices the equation of motion for the Wilson loops are known as the Migdal-Makeenko loop equations [34]. For the case of hermitian matrices, the loop equations follow by inserting the matrix derivative under the integral: The loop equation is a quadratic equation in the large N limit which takes the form The integer p(C; C 1 , C 2 ) specifies the number of ways in which a loop C can be split into loops C 1 and C 2 , while the integer j(C, s; C ) specifies the number of ways loops C and s can be joined to produce C . Both will be described in detail below. Solving these nonlinear coupled equations is highly non-trivial. Using collective field theory one obtains an effective potential which, when minimized, gives the large N solution to the Schwinger-Dyson equations. As for the c = 1 MQM systems, changing to invariant variables we obtain a collective potential V col , which when minimized, again determines the large N expectation values of loops. These will be discussed in detail below.

Collective (loop space) representation
The collective representation of multi-matrix systems can be given both at the action level for matrix integrals (c = 0), and at the Hamiltonian level for coupled quantum mechanical (c = 1) systems. These are closely related, so in formulating numerical methods one can work in the Hamiltonian framework [36]. To proceed, let us define the adjoint of the loop variable. For a loop φ(C), we define its adjoint as φ(C), whereC is the reverse of the word C. For example, for a word C = abaab, where a = M 1 and b = M 2 , its reverse is C = baaba. In hermitian matrix models, this corresponds to taking a hermitian conjugate of the matrix products. As such the adjoint φ(C) is simply the complex conjugate of φ(C) in hermitian matrix systems: The "bar" symbol on the right hand side of the above formula denotes complex conjugate. We note in the hermitian one-matrix models case, the adjoint of a loop is itself since all loops are real valued. The use the loop adjoints makes the hermiticity of the collective representation manifest as we will see below.
Matrix integrals. Let us start with the collective representation of the multi-matrix integral problem with any number d of matrices which at large N is efficiently described using an effective action. This is obtained by a change of integration variables (from matrices to loops) resulting in a large N collective action S eff = S − ln J. The main ingredient in this effective description is the Jacobian J and its form is in general specified by the collective formalism. Variation of this collective action correctly produces the Schwinger-Dyson (SD) equations (2.9). The derivative of the collective action [35] gives The functional ω(C) stands for representing the splitting of contour C into sub-contours (C 1 , C 2 ). 1 The split occurs at the pinched link, so in general there are several distinct ways for the process, which necessitates the sum. The integer p(C; C 1 , C 2 ) counts the number of ways loop C can be partitioned, by the splitting operation, into loops C 1 and C 2 . In a similar way one has The generators of the Lie algebra of U(N ) are normalized such that α t α ab t α a b = 1 2 δ ab δ a b . For hermitian matrix systems,Ê α l →Ê ab

JHEP01(2022)168
for the opposite operation of joining contours. We note the use of loop adjoint in the definition. The integer j(C 1 , C 2 ; C) counts the number of ways in which C 1 and C 2 can be joined to produce C. A computer algorithm was developed to generate loops and also these loop processes which we call loop algebra.
For numerical minimization it is useful to follow [45], where it was established (through a stochastic quantization) that the above SD problem can be represented through (Fokker-Planck type) the quantum mechanical Hamiltonian: H = K + V eff . The solution of nonlinear SD equations is then obtained by minimization of an effective potential of the following form , (2.16) where the "bar" symbols again denote complex conjugates. This effective potential also gives the leading large N configuration of the bosonic sector of multi matrix Marinari-Parisi [46] type models. Expansion around the stationary point leads to equations for small fluctuations and a systematic 1/N expansion scheme.

Matrix quantum mechanics.
We now describe the collective field formulation of the large N MQM in detail. Considering a general multi-matrix quantum mechanics problem, in the operator formalism, one has a transition to the collective description by performing a change to curvilinear (loop space) variables (with a Jacobian J) which induces the collective Hamiltonian, taking the form and π(C) representing the conjugates to the loops φ(C). Here V [φ] is the original potential written in terms of loops. The collective Hamiltonian is manifestly hermitian, and is equivalent to the original MQM Hamiltonian (2.2). We have used the notation V col to denote the collective potential obtained for MQM systems, which plays an analogue role of the effective potential V eff for matrix integral problems, as described above.
Notation emphasis. To proceed, let us emphasize that throughout we will use V eff to denote the effective potential associated with matrix integrals (c = 0 systems), and V col to denote the collective potential associated with MQM (c = 1 systems).
Large N background and fluctuation spectrum. In the next subsection we explain that the N → ∞ limit is obtained by minimizing the potential V eff or V col and that this minimization is subject to a sequence of inequalities which constrain the range of the loop variables. The relevance of constrained minimization becomes even more fundamental when

JHEP01(2022)168
one proceeds to study the spectrum at large N . Ordinarily this would be given directly in loop space by expanding about the stationary field 20) and reducing the collective Hamiltonian to a quadratic, small fluctuation Hamiltonian ) . (2.22) This H (2) , and its diagonalization provides the spectrum at large N . Essentially, it is determined by eigenvalues of the matrix: As we will explain next, these are to be found subject to obeying a set of (loop) space inequalities that are central for reaching the correct minima.

Loop space inequalities and constrained minimization
Positivity of the loop joining matrix Ω. In the collective Hamiltonian description, the N → ∞ limit (and the sum of planar diagrams) is given by the semiclassical approximation. The problem is therefore to solve for the static stationary configuration denoted φ 0 (C), which minimizes the potential V col . This would be generally given by the equation However, as was understood earlier [37,38], in the minimization a role is also played by a sequence of inequalities (analogous to Schwarz inequalities) which constrain the range of loop variables. Such a sequence of inequalities is generic, and will be present for any set of variables representing invariants. In the collective description the inequalities are directly visible and can be generally specified and given in terms of the loop space matrix Ω. This matrix participates in the kinetic term of the Hamiltonian and defines the loop space symplectic form. As such the matrix Ω must be positive semi-definite. Indeed from the definition of the matrix Ω one has the fact that it can be written as i.e. it is explicitly positive semi-definite. Indeed, in our case, defining 2

JHEP01(2022)168
we obtain Positivity constraints. The minimization therefore must be done subject to the positivity condition of the loop joining matrix Ω. In general the positive semi-definiteness of the loop valued matrix Ω gives the complete set of inequalities, which schematically reads We can write the complete set of these generalized loop space (Schwarz) inequalities as follows. A convenient basis (and explicit form) turns out to be given by the sequence of the sub-determinants det k (Ω) ≥ 0. The sub-determinant is defined by (repeated indices are summed) 28) where N Ω is the dimension of the loop space joining matrix Ω and The expression in the last line of (2.28) makes the positivity of the sub-determinant manifest. An alternative formula for the sub-determinant is provided by where χ (1 l ) (·) is a Schur polynomial and (1 l ) is the Young diagram with a single column and l rows. This formula is simple with the normalization chosen in (2.28). The subdeterminant basis is particularly useful when some of the constraints are saturated, in which case certain eigenvalues of Ω vanish. If p eigenvalues vanish, then there are p independent vanishing sub-determinants

JHEP01(2022)168
Once expressed in terms of loops, this positivity condition implies highly non-trivial inequality constraints among the φ(C)'s. The positivity inequalities constrain the eigenvalues of the loop space matrix Ω, leading to a constrained minimization of V eff or V col . 3 The stationary minimum solution is generally characterized by saturation of a certain number of inequalities, whereby the loop space matrix Ω develops zero eigenvalues.

Master field
A most complete way to deal with the constrained optimization at large N is through variables that automatically assure the positivity condition of the loop space matrix Ω. These are the master field variables, or simply the original matrix valued variables of the system. In terms of such variables the expectation values of loops are determined by a direct stationary-point equation of the loop space effective (c = 0) or collective (c = 1) Here the notation V eff/col represents either V eff or V col . This saddle-point equation in terms of the master variables is correct in all phases of the theory, both the weak-and strongcoupling. Let us also mention how the correct loop-space equation follows from the master equation. Multiplying and summing over the appropriate factors, we obtain This fact that the saddle-point equation in terms of the master variable produces the correct loop-space equation can be taken as another argument for the existence of the master field. Denote the set of master variables and their complex conjugates by φ α andφ α , respectively. To see that the positivity of Ω is assured when working in terms of the master variables, note that we can write Ω as The set {φ α ,φ α } is always assumed to be at least as large as the set of invariants {φ(C)}.
We can think of these variables as the original variables of the theory, in which case they transform non-trivially under the existing internal symmetries and form a larger set than that of the invariants that can be obtained from them. Nevertheless, we have also other situations in mind: one may truncate the effective potential, in which case the number of φ α variables must be larger than the number of loop variables included in the potential. In general, the set {φ α } is at least as large as the set of invariants. This is the case in one-matrix models, for instance, where one identifies the φ α 's with the matrix eigenvalues.

JHEP01(2022)168
The important characteristic of these variables then is that they solve the positivity constraint explicitly. Therefore, in the large-N limit, the ground state configuration is determined by the saddle-point condition It is clear that it is useful to reduce the problem to unconstrained minimization. In [37,38] we have presented and tested numerically such a framework developing a hybrid loop space + master field approach. The use of master variables [38] is simple. In terms of the original variables the inequalities are automatically obeyed so we can think of changing back to these master variables but keeping the collective loop space Hamiltonian. This is because it is the loop space Hamiltonian that generates the 1/N expansion. Now the ground state and the fluctuation spectrum will be given by unconstrained minimization, satisfying This represents an implicit equation for the master field φ α since it enters V eff/col through the loop variables φ(C{φ α }). As such it is not very useful at the analytic level but it can be easily implemented numerically [37]. We will consider systems of two hermitian matrices in this article. Systems of hermitian matrices are always in the phase in which some of the constraints are saturated. Consequently for these systems the use of master variables is of particular importance. The loop invariants consist of single trace products of these matrices. As a result, one of the matrices can be chosen to be diagonal, and the other will be parametrized in the Lie algebra of the unitary group. The master variables are then real, which is something we use in the following. 4 The master eigenvalue equations for the fluctuation spectrum are similarly obtained in this hybrid scheme. Denote in our Hamiltonian formulation the master field solution by φ 0 α . We are then led to appropriate spectrum eigenvalue equations through a shift in the loop space collective Hamiltonian. One essentially considers a 'canonical' transformation [40] from loops to master fields to obtain the following quadratic Hamiltonian Note that the term linear in η vanishes as a result of equation (2.35). The mass matrix M αβ is then to be computed from

JHEP01(2022)168
Here the repeated indices C and C are summed, which in principle range from 1 to infinity. To compute it numerically we must perform the truncation discussed below. The fluctuation spectrum is obtained by solving for the square root of the nonzero eigenvalues of the mass matrix. In section 5 we will show in detail that this is equivalent to the direct computation of fluctuation spectrum in loop space.
We will see that in numerical evaluations of multi-matrix problems, the hybrid loop space + master field formalism turns out to be advantageous, converging rapidly and giving excellent agreement already at small loop and color cutoffs. With interest in precision optimization we will furthermore test the method for larger sizes and number of minimization variables (∼ 9 × 10 3 ).

Loop truncation
A numerical implementation of the large N minimization of the collective potential necessarily involves a truncation of the infinite dimensional loop space. We will truncate the loop space by restricting to loops that limit the number of matrices appearing in the trace to be smaller than some fixed cut off L max . This loop truncation is reminiscent of level truncation in string field theory [51]. To gather some insight into this truncation, in this section we consider the counting of loops as a function of the cut off. A convenient approach towards this counting uses permutations to enumerate loops.
Loop counting. Consider the hermitian two-matrix model. The truncation of loop space that we employ restricts n+m ≤ L max a fixed maximum loop length, where n and m denote the number of M 1 and M 2 in a loop, respectively. It is interesting to count the number of loops as a function of L max . This counting is a useful input to the numerical implementation since it specifies how many variables appear in the numerical minimization. There is an extremely rapid growth of the number of loops with increasing L max . The counting problem can be solved with a standard application of Polya theory [52]. We start by introducing the single letter partition function, which for two matrices is The single trace partition function, which counts the number of loops, is now given by The sum is over all integers n. At each n there is a second sum over d which runs over the divisors of n, i.e. all the integers that can be divided into n without remainder. The function ϕ(d) is the Euler totient function. The degree L contribution to F (x, y) counts single trace operators constructed from L matrices. The coefficient N n,m of the monomial JHEP01(2022)168 of degree n in x and degree m in y counts the number of loops that can be constructed using n M 1 s and m M 2 s. The first few terms of F (x, y) are To interpret this answer, note that for example, the contribution 2x 3 y 2 implies that there are two independent loops that can be constructed using 3M 1 s and 2M 2 s. These two loops are To explore how the total number of loops grows, it is useful to consider the blind partition function, obtained by setting y = α = x. The coefficient of α n counts the total number of loops constructed using n matrices. The blind partition function is F (α, α) = 2α + 3α 2 + 4α 3 + 6α 4 + 8α 5 + 14α 6 + 20α 7 + 36α 8 + 60α 9 + 108α 10 + 188α 11 + 352α 12 + 632α 13 + 1182α 14 + 2192α 15 + 4116α 16 + 7712α 17 + +14602α 18 demonstrating an extremely rapid growth in the number of invariants (loops).
Loop space truncation. Our numerical implementation of loop space dynamics truncates to the subspace of invariants, given by all loops with L max = 2l − 2 matrices or less in the trace. In this scheme Ω is an N Ω × N Ω matrix, where N Ω is the number of loops with l matrices or less in the trace. Ω itself depends on a total of N Loops , which is the number of loops with 2l − 2 matrices or less. For this reason, our minimization scheme minimizes with respect to N Loops independent variables. The values of L max , N Ω and N Loops for values 3 ≤ l ≤ 10 are given in table 1 below. For the minimization, we parametrize M 1 and M 2 using N + N 2 real valued master variables. These are the N eigenvalues for M 1 and the N 2 matrix elements in M 2 . For a given l we must choose the number of colors N large enough that N 2 + N ≥ N Loops . To obtain the large N background, our numerical experiments focus on l = 9, so that we keep a total of 8923 loops. We choose N = 94 so that there are a total of 8930 master variables.

Spectrum calculation.
To obtain the fluctuation spectrum based on master variables, one needs to preserve the matrix structure of the loop derivatives (3.7) with respect to both M 1 and M 2 . We then truncate the mass matrix equation (2.39) as follows 4  9  15  6  15  37  8  23  93  10  37  261  12  57  801  14  93  2615  16  153  8923  18 261 31237 Here the "bar" symbols represent the complex conjugate. The mass matrix M αβ is a matrix of dimension obtained from the multiplication of matrices. It is essential that one sums over all N Loops in the equation above. 5 One obtains N Ω nonzero eigenvalues and (2N 2 ) − N Ω (numerically) zero eigenvalues. We now observe that is a (N Loops × N Loops ) matrix which, in loop space, would include all loops up to length 4l − 6. In practice, it is not feasible to obtain such Ω directly in loop space, given the size of the truncations considered in this article (e.g., for l = 9 , 4l − 6 = 30). However, it can be generated from the master variables. The nonzero eigenvalues of (3.5) can then be matched with those of the loop space spectrum matrix This is a N Loops × N Loops matrix, and is expressed explicitly in terms of loop variables. It has N Ω nonzero eigenvalues and N Loops − N Ω (numerically) zero eigenvalues. Throughout, we have checked that the nonzero eigenvalues of (3.5) and (3.6) are identical. 5 Indeed, it was shown in [40] that if the sum is restricted to N Ω loops, every non-zero eigenvalue of (2.23) is also an eigenvalue of (3.5), with N Loops replaced by N Ω . Except in the strong coupling phase of unitary matrix systems, the spectrum obtained simply on the basis of (2.23), with Ω 0 a (N Ω × N Ω ) matrix is not accurate.

Optimization procedure
A Python code was developed to obtain the computational results. The first step is to generate all the distinct single trace loops of a given length l. There are different ways of generating them, including using Mathematica and Polya theory. A simple procedure is to generate the list of C n l combinations for 0 ≤ n ≤ l, which index the position of (say) the matrix M 2 in the string (word) of M 1 and M 2 matrices, and then remove loops which are identical up to cyclic permutations. They are then indexed and stored as a list of arrays, eg. [1, 1, 1], [1, 1, 2], [1,2,2], [2,2,2], etc., and stacked for different lengths, to obtain the list of N Ω and N Loops loops described in section 3.1, reproducing table 1. The zeroth indexed element of the list is the empty array corresponding to φ(0) = Tr(I)/N = 1. It is fixed throughout.
The next step is to generate the loop joining matrix Ω Here we use little 'c' instead of capital 'C' to emphasize the loop truncation: the loop joining matrix Ω now is a finite matrix of dimension N Ω × N Ω instead of an infinite dimensional matrix. The code implements explicitly the first equality in the above equation, recalling that 6 Tr(M 1 g(· · · M 1 · · · M 2 · · · )) + · · · = g ji (· · · M 1 · · · M 2 · · · ) + · · · (3.7) (g ij is obtained by extracting M 1 from the loop). The joined loop φ(c ) is identified, and the nonzero joining coefficients j(c, c ; c ) are stored. It should be emphasized that, through a joining process higher loops are generated, and the N Ω × N Ω matrix Ω contains higher loops up to length N Loops . It is this full set of loops contained in and generated by Ω that will participate in the optimization process. A similar procedure is followed to generate ω defined through loop splitting ω(c) =

JHEP01(2022)168
Here t ab ij (a < b) is the set of real off-diagonal generators (σ 1 in the entries (ij) and (ji)). t ab ij (a > b) are the purely imaginary generators (σ 2 in the entries (ij) and (ji)), and t aa ij are the entries of a diagonal matrix. In order to extract the explicit dependence on the powers of N from the loops, they are defined as We used a standard minimize function from the scipy.optimize library. We used two methods, the BFGS and the CG methods. We found that the BFGS method is slightly faster, but for large loop truncations, the CG method is more stable. Both these methods require the evaluation of the gradient. This is achieved by calculating for each iteration the derivatives of the loops with respect to the master variables The initial master variables configuration consists of a randomly generated real vector and of a randomly generated real matrix. For Fokker-Planck (and the underlying c = 0) type systems, we have set as convergence criteria that the norm of the gradient vector becomes less than N (N + 1) 10 −16 . In other words, at convergence, a typical gradient vector element has norm of order 10 −16 . Convergence of the algorithm is remarkably stable, with the energy monotonically decreasing to zero in successive iterations. Depending on the size of the truncation, the energy at convergence is ∼ 10 −24 − 10 −31 . The norm of the gradient components typically range from ∼ 10 −15 −10 −20 . At convergence, the Schwinger-Dyson equations are satisfied to typical accuracy ∼ 10 −10 −10 −18 . With a 3.0 GHz Mac, the codes take from about a few seconds for N Loops = 37, about two hours for N Loops = 2615 and more than a day for N Loops = 8923.
For spectra of the MQM (c = 1) systems discussed in 5, we use N = 94 to generate large N background. The initial master variables configuration again consists of a randomly generated real vector and of a randomly generated real matrix. The convergence criteria is that the norm of the gradient vector becomes less than N (N + 1) 10 −16 . Again, for N Loops = 37, converge is achieved in hours, while for N Loops = 8923 convergence takes days.
For a given loop truncation size N Loops , starting with the lowest N satisfying N (N + 1) ≥ N Loops and increasing N , loop values are seen not to change much. Loop values quickly converge to their exact values (when known) as N Loops is increased. This is particularly so for the N Ω small loops. This is evidenced, for example, in the matrix integrals case, in table 8.

SD models
As a first application of the methods outlined above, we will study matrix integrals. For these models (in the decoupled case) there are exact analytic calculations which can be used to validate our numerical results. In all cases we are able to confirm that our numerical results are essentially exact for loops of lower length, with the accuracy falling off as we approach L max which is the maximal length of loops admitted in the minimization. We note that all these models have a Hamiltonian quantum mechanical interpretation (Fokker-Planck). We also consider a variety of problems with two matrices to demonstrate the methods. It will be clear that the extensions to more matrices proceed in the same way, without difficulty.
Numerically we are evaluating the integral (2.8) with the action S given by where the potential is given by When the coupling k = 0 we will refer to the system as a "single matrix system" and when k = 0 as a "two-matrix system". Note however, that for both types of systems we evaluate mixed loops φ(C) obtained by tracing products involving both matrices M 1 and M 2 , so that even when k = 0 the problem is still a multi matrix problem.

Single matrix systems
Free theory. The potential of the zero dimensional model is so that we have set g 3 = g 4 = 0. We also set k = 0. The effective potential is Quartic theory. The potential of the corresponding zero dimensional model is so that we have set k = g 3 = 0. The effective potential is

JHEP01(2022)168
Cubic model. We use so that we have set k = g 4 = 0. The effective potential is

Two-matrix systems
Quadratic model. The potential of the zero dimensional model is so that we have set g 3 = g 4 = 0. In this case we keep k = 0. The effective potential is Coupled cubic two-matrix model. We keep k > 0 and we set so that we have set g 4 = 0. The effective potential is

Results
In tables 2-10 we present results for SD models obtained with l = 9 and N = 94. 7 The results show that for small loops we essentially obtain the exact answer. For larger loops (with more than l matrices in the trace, in the notation of section 3) the results are less accurate, but even for the longest loops accuracy is typically always better than 5%.   Table 10. The two lowest lying states for k = g 3 = 0 and g 4 = 10.

Matrix quantum mechanics spectrum
In this section, we consider multi-matrix quantum mechanical models, which includes the evaluation of (Wilson) loop expectation values at large N , ground state energies at large N , and the spectrum of fluctuations (which corresponds to N 0 , i.e. order 1). The coupled two-matrix Hamiltonian reads

JHEP01(2022)168
We note, that the SD models featured in the previous section are also of this form, but an effective potential containing double trace couplings. Consequently the numerical methods for evaluating the stationary points and the spectrum apply to both with no difference in degree of difficulty. Likewise, our methods can be applied to multi-matrix models with arbitrary number of matrices. We describe evaluation of the spectrum for the above case. For normalization purposes, we will also give numerical (and analytical) one-matrix results.

Fluctuations
The general strategy for the spectrum calculation of small fluctuations in loop space, which holds for all multi-matrix quantum mechanics in the limit of large N , has been described in previous sections, and is easily implementable on a computer. In this subsection we explain the subtly about truncation further. As in the optimization procedure, we truncate the infinite dimensional loop space. Let l denote the chosen loop length truncation, and N Ω denote the number of loops whose lengths are less or equal to l. Ω then is a N Ω by N Ω matrix which involves N Loops loops in total, and the loop of maximal length it contains is L = 2l − 2.
Sending Ω → N −2 Ω, we obtain the collective Hamiltonian Here is the original potential written in terms of loops. We then expand H col to order O(1), which is accomplished by shifting loop variables around the ground state and omitting the constant background terms. We use Capital C to emphasize that here the loop labels run from 1 to N Loops , instead of N Ω , because V col in equation (5.4) involves N Loops independent loops. Accordingly there are also N Loops canonical conjugates, which requires us to use the N Loops by N Loops dimensional loop joining matrix denoted as Ω.
The choice of the loop joining matrix is illustrated in figure 1, where the deep blue blocks are used for truncation in V col , and the deep plus light blue blocks, i.e. Ω, are used in the following H col . The Taylor expansion then gives is the Hessian matrix of V col at the ground state. We note that all elements in Ω are linear functions of loops. Their second derivatives therefore are 0, and hence do not contribute to V  To evaluate the spectrum, one can diagonalize the kinetic term in (5.6) and then solve the eigenvalues of the resulting mass matrix. As pointed out earlier, we see that using Ω also resolves the mismatch of the dimensions between Ω 0 and V (2) 0 . Since Ω 0 is positive definite, one can perform a canonical transformation The spectrum is then given, in terms of the nonzero eigenvalues of the spectrum matrix where eig n denotes the nth nonzero eigenvalue. Here n starts from 1 instead of 0, because the zero mode is excluded in the definition of Ω (or Ω). The spectrum is independent of truncation loop length l, as we will see in the following concrete examples. In principle, the size of the spectrum one can obtain is equal to the N Ω . Besides, there are precisely N Loops − N Ω zero eigenvalues of Ω 0 V (2) 0 , therefore the truncation scheme automatically projects out the higher modes. As l is increased, higher modes are included, and one is able to obtain higher frequencies.
Depending on the size of the truncation and particularly for multi-matrix systems, it is not feasible in general to obtain Ω 0 directly in loop space. But Ω 0 can be always be constructed with master variables. The spectrum equation (5.8) then takes the form For example, the spectrum that is presented in section 5.3 is obtained in the background of 801 (this is the number of loops for a cut off of L max = 2×7−2) loops, whose value is determined by the master field after minimization. This results in a spectrum eigenvalue matrix of size ≈ 10 3 × 10 3 . Note also that 401 947 loops are effectively included in the computation of Ω 0 . This is only possible through the use of master variables and a direct loop space evaluation would not be possible. It is visible that the spectrum calculation in JHEP01(2022)168 exact results 0 0. terms of master fields should give the same results, except for different numbers of zero eigenvalues of Ω 0 V 0 . We observe that in the free theory cases, interestingly, one can actually obtain exact results by working with a smaller V (2) 0 matrix, whose matrix indices range only from 1 to N Ω , and correspondingly with the smaller Ω 0 matrix. This has been verified in both one-and two-matrix quantum mechanics. This is no longer the case the moment coupling constants are switched on.

One-matrix example
We then proceed to employ our general strategy to compute the spectrum of the hermitian one-matrix quantum mechanics:  4) we can also obtain the spectrum using a computer. A Mathematica program was developed for the spectrum calculation. The one-matrix case has the simple feature that N Ω = l and also N Loops = L max , which provides us a canonical example to illustrate our general strategy. To calculate the spectrum using the general strategy at loop length truncation l = 6, for example, there are totally L max = 10 loops contained in V col . Ω 0 then is a 6 × 6 matrix, but Ω 0 and V   reduces to decoupled harmonic oscillators, therefore we have ε n = n, as is verified. For other couplings we present the first level frequency ε 1 in figure 2b. The exact and numerical results again agree very well, including the critical region g ∼ g c = −1/3 √ 2π.

Two-matrix example
By applying the same strategy, and using the numerical minimization results, we are also able to evaluate the spectrum of the hermitian two-matrix quantum mechanics (5.1) and (5.2), using a Mathematica program. At first we present in table 12 the comparison of the ground state energies, in which the two-matrix results are obtained from numerical minimization, and the one-matrix results are obtained from the analytical expression (B.5). The two-matrix model in c = 0 cases doubles the one-matrix model interaction, and one expects that the ground state energy should therefore be twice that of the corresponding one-matrix result. Table 12 exhibits excellent agreement, verifying the expectation.
As for the spectrum computation, the multi-matrix models lose the simple feature N Ω = l, and a rapid growth of N Ω has already been observed in the two-matrix case. Based on the l = 9 (l = 7) numerical minimization data for interacting (free) case, we present some spectrum results evaluated at truncation l = 4, 5, 6, 7 in tables 13, 14, 15, 16 and 17. Taking the l = 6 case for example, one has N Ω = 37, so that Ω in V eff is 37 × 37.  Table 13. Two-matrix spectrum with loop length truncation l = 4.   Table 14. Two-matrix spectrum with loop length truncation l = 5. On the other hand, the fact L = 10 requires us to use the N Loops = 261 dimensional Ω 0 and V (2) 0 , which implies a large amount (261 − 37 = 224) of the zero eigenvalues of Ω 0 V (2) 0 . This was verified by our Mathematica program. Comparing the two-matrix spectrum results, we see that as l is increased, low lying spectrum are stable and convergent, and higher level frequencies are obtained.  Table 16. Two-matrix spectrum with loop length truncation l = 7. Table 17. Two-matrix spectrum with loop length truncation l = 7 continued.

JHEP01(2022)168
Let us now examine in more detail the primary characteristics of the two-matrix spectra. As the results show, in free theory, k = g 4 = 0, the spectrum values coincide with the one-matrix case. Besides, a high degeneracy pattern is observed. Let N (n) denote the number of inequivalent loops of length n. The degeneracy at level n is precisely equal to N (n). For instance, at level 2 the degeneracy is 3, since there are 3 independent loops, namely Tr(M 1 M 1 ), Tr(M 1 M 2 ), and Tr(M 2 M 2 ). As we turn on the quartic coupling g 4 while fixing k = 0, the degeneracy is lifted slightly, and a different degeneracy pattern is obtained. Comparing with the one-matrix results table 11, one can obviously see that the two-matrix spectrum contains the corresponding one-matrix spectrum as a subset, and with a degeneracy number 2. These subsets are filled in light blue colors in each k = 0 and g = 0 column of the two-matrix spectrum results. When turning on coupling k, the degeneracy patterns are nearly destroyed. The lifting of the degeneracy due to k and g 4 is presented in figure 3 based on data in table 15. We conjecture that these properties are universal for multi-matrix models.

JHEP01(2022)168
son loop background, it also establishes in concrete terms the existence of the master field. This existence. was questioned through the years in various works. Constrained minimization, which is accomplished with fairly large number of (loop) variables (∼ 10 4 ) is seen to give essentially exact results for large N expectation values, ground state energy and low lying spectra. This size can be further increased for added precision. The formulation, and methods developed are such that there is no difference in having larger number of matrix variables. The large N dynamics is evaluated in loop space, which is parametrized the same way (and easily computer augmented). Regarding further works, and applications of the methods developed in this work, one can contemplate many physical problems requiring understanding through matrix models range from confinement [53] to cosmology [54]. More specifically, extensions to supersymmetric versions of the multi-matrix QM can be considered. The four matrix BMN quantum mechanics (at large N ) is clearly accessible by this scheme, and we plan to present results in future work. Considering the two matrix case, whose solution is accomplished in this work, an interesting application to entanglement [55], thermodynamics and the matrix thermofield double (TFD) states. In general this requires a study of QM defined on the Schwinger-Keldysh contour [56,57]. The TFD state and the corresponding wave functional of the O(N ) vector theory based on the collective field formulation was recently studied in [58]. Adjusting our optimization to matrix systems appears possible. An approximate version of the TFD involves coupling the system through single trace interaction [59][60][61], this corresponds to interaction term in QM that we studied, therefore these models provide a possibility to simulate temperature in the ground state. Likewise will be a fuller exploration of the phase structure of this theory.

Acknowledgments
The work of RdMK, AJ and JPR is partially funded by a Simons Foundation Grant, Award ID 509116. The work of RdMK is also supported by the South African Research Chairs initiative of the Department of Science and Technology and the National Research Foundation.

A Analytic results for matrix integrals
For the cases with k = 0 we can use the techniques [62] to get values of loops that the numerics must reproduce. When k = 0 but g 3 = g 4 = 0 we are left with a quadratic integral which is easily performed.

A.1 Quartic model
For the matrix integral with action The solutions to these equations are ugly, so we will plug in definite values of g 3 . Note that we must take g 2 3 < 1

B Analytical results of the one-matrix quantum mechanics
We consider the hermitian one-matrix quantum mechanics with a quartic interaction This model is dual to the D = 1 string theory in the double scaling limit. The model was first solved in [19] and its spectrum was solved in [11] using a field theoretic approach. We briefly summarize the analytical results derived there.
Eigenvalue distribution. In the large N limit the eigenvalue distribution of the ground state [36,62] is where x ∈ [−Λ, Λ], and Λ is determined by the constraint The constraint is saturated at the critical value g c = −1/3 √ 2π.
Loop values. The loops at the ground state are (B.5) In the free case we have E (0) gs = 1/2. Figure 4 shows the results obtained from the above analytic expression, as well as the collective potential V eff values, demonstrating an excellent agreement.
Spectrum. In the large N limit the spectrum of small fluctuations is (B.6) The integral can be evaluated analytically where K is the complete elliptic integral of the first kind. The general feature of (B.7) is that the level n frequency is proportional to n, no matter how g 4 varies. One thus expects that the spectrum fits a straight line when plotted as a function of level number.  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.