Bootstrapping Matrix Quantum Mechanics

Recent work has shown that large $N$ (multi-) matrix integrals can be solved numerically by imposing positivity constraints on higher point correlation functions [1]. We have generalized this method to obtain the spectrum and simple expectation values of large $N$, gauged matrix quantum mechanics. In our approach, operator expectation values are related through conditions such as $\langle [H,{\mathcal{O}}] \rangle = \langle G {\mathcal{O}} \rangle = 0$. Here $H$ is the Hamiltonian, $G$ the generator of $SU(N)$ gauge transformations and ${\mathcal{O}}$ an arbitrary operator. Bounds on the energy and expectation values of short operators are then obtained from positivity constraints on the expectation values of certain longer operators. We first demonstrate how this method efficiently solves the conventional quantum anharmonic oscillator. We then reproduce the known solution of large $N$ single matrix quantum mechanics. Finally, we present new results on the ground state of large $N$ two matrix quantum mechanics.


Introduction
Large N matrices are at the heart of the holographic emergence of semiclassical, gravitating spacetime geometry [2]. Matrix quantum mechanics theories focus our attention on the essence of this phenomenon; unlike higher dimensional large N Yang-Mills theories, matrix quantum mechanics has no built in locality to start with, and thus space must emerge in its entirety. The simplest such theory is the single matrix quantum mechanics description of two dimensional string theory [3], while the richest are the maximally supersymmetric multi-matrix theories of BFSS [4] and BMN [5]. There are many theories in between, with varying numbers of matrices and degrees of supersymmetry [6]. Thus far, only the single matrix quantum mechanics has proved solvable at large N [7].
Nonzero temperature Monte Carlo studies of large N multi-matrix quantum mechanical systems have successfully captured aspects of a known dual spacetime in supersymmetric theories [8][9][10][11]. Substantial Monte Carlo studies have also been performed for nonzero temperature bosonic multi-matrix theories, e.g. [12,13]. However, recent work increasingly suggests that the quantum structure of holographic quantum states -revealed for instance in their entanglement [14][15][16][17] -plays a central role in the emergence of space. It therefore behooves us to find methods suitable for studying the zero temperature quantum states of multi-matrix quantum mechanics directly. Progress was made recently in this direction by using a neural network variational wavefunction [18]. Here we describe a different approach.
Our work is directly inspired by a recent beautiful paper by Lin [1]. That paper studied large N matrix integrals, which is an easier problem than large N quantum mechanics but shares important features. Lin showed how relatively simple positivity constraints on higher point correlation functions could be combined with the large N loop equations to efficiently produce strong numerical constraints on correlation functions of matrix integrals. In the following we will show how this methodology can be adapted to the quantum mechanical problem. The positivity constraints will be essentially the same, augmented to include momentum as well as position operators in the case of matrix theories, while certain operator identities in energy eigenstates (or more general Gibbs states) will play the role of the loop equations.
In §2 we consider a warm-up case of a conventional quantum mechanical anharmonic oscillator. The results are in Fig. 1, showing that the ground and first excited state energies E and expectation values x 2 can be strongly constrained with little numerical effort. This case is the closest to the matrix integral results of [1], because expectation values of operators depending on the momentum can be solved for explicitly and a recursion relation (9) can be found that determines all expectation values x t in terms of just x 2 and E.
In §3 we use the bootstrap methods to solve a large N one matrix quantum mechanics.
Here the momentum operators cannot be eliminated explicitly, and we do not use a closed form recursion relation for all expectation values. However, we find that the energy and expectation values of short operators can be efficiently constrained by applying positivity constraints to a matrix generated by all strings of operators of length ≤ L. There are 2 L such strings, leading to a matrix with 2 2L entries. The results in Fig. 2 show that the known analytic results for the ground state energy E 0 and expectation value tr X 2 can be reproduced to very high accuracy already with L = 3. We also constrain tr X 2 in excited states with energy E, which is continuous in the large N limit, as shown in Fig. 3.
In §4 we turn to a two matrix quantum mechanics. The methodology is the same as for the one matrix case, except that now there will be 4 2L matrix elements to consider. This model is not solvable, so we have corroborated our numerical results with a Born-Oppenheimer approximation, that gives rigorous lower and upper bounds on the ground state energy. The results in Fig. 4 show that the minimal energy allowed by the bootstrap equations at L = 4 is within the narrow region allowed by the Born-Oppenheimer bounds, and quite close to the lower bound, suggesting that the numerics is close to convergence onto the true ground state energy. That figure also shows the bootstrap results for the expectation values tr X 2 and tr[X, Y ] 2 .

The quantum anharmonic oscillator
The essence and effectiveness of our bootstrap method can be demonstrated in a toy example of a quantum anharmonic oscillator, with Hamiltonian Here [p, x] = −i. We will use positivity constraints on expectation values to obtain the results shown in Fig. 1 below for the energy E and expectation value x 2 of the ground state and first excited state. In later sections we will generalize this approach to the quantum mechanics of matrices.
The first step is to relate the expectation values of different operators. We will obtain the recursion relation (9) below. For expectation values in energy eigenstates, and for any operator O, For example, let O = xp. Equation (2) then gives the Virial theorem, This relation links the expectation values of the operators p 2 , x 2 and x 4 . Following from the Virial theorem, the energy of the eigenstate is More systematically, take O = x s and O = x t p in (2) for integers s, t ≥ 0. Commuting the operators x, p with the identity [p, x r ] = −irx r−1 gives the relations Eliminating the terms with a single p operator, we arrive at the relation A powerful identity in this single particle example is that which is a strengthened version of (2). We emphasize (2) instead of (7) because, as we will see later, (7) becomes less useful in the matrix case. Nonetheless, in the present anharmonic oscillator case, we can set O = x t−1 in (7) to obtain Plugging (8) into (6) gives a recursive relation between expectation values of powers of x: where E is given by (3). Also we know that x 0 = 1 and x t = 0 if t is odd, so all expectation values of x t can be computed from E and x 2 with (9).
With the recursion relation (9) at hand we move onto the second step. We wish to solve for E and x 2 , the only two unknown variables, by bootstrapping. This step works as in [1].  Upper plot: the allowed region for (E, x 2 ) near the ground state solution (marked by the red cross) for different sizes of the bootstrap matrix K = 7, 8, 9; lower plot: the allowed region near the first excited state.
The basic positivity constraint is that The bootstrap result is shown in Fig. 1. Even for moderate K the values of E and x 2 can be determined quite accurately. We also see that the region of allowed values splits into a discrete set of islands, corresponding to the discrete spectrum of the Hamiltonian. (The spectrum will be continuous later in the infinite-N matrix models we consider.) Both the spectrum and observables at different energies can be extracted in this way, but higher-order operators or higher-energy states will require more constraints to be computed accurately.

One matrix quantum mechanics
The bootstrap method discussed in the anharmonic oscillator example can be generalized to matrix quantum mechanics at N = ∞. To start, we consider the single-matrix quantum mechanics with Hamiltonian: where P and X are N -by-N Hermitian matrices with quantum commutators [P ij , X kl ] = −iδ il δ jk . The theory (11) can be mapped onto the quantum mechanics of N free fermions and is hence solvable [7]. In this section we will show how positivity constraints reproduce the known solution. The matching is shown in Figs. 2 and 3 below.
To apply the bootstrap ideas, we relate operator expectation values by symmetries. In Such a state ρ could be a pure energy eigenstate or a mixed thermal state, for example.
Similar to the single-particle case in the previous section, the Virial theorem is derived from (11) by choosing O = tr XP : The Hamiltonian (11) also has an SU (N )-adjoint symmetry generated by the matrix where the extra identity piece ensures that tr G = 0, with the operator ordering [X, P ] = XP − P X in (14). We are interested in gauged matrix quantum mechanics, so that physical states must be invariant under this symmetry. Thus in particular, The simplest implication of this constraint is as follows. From tr G = 0 we have and on the other hand because [H, tr X 2 ] = 0, From symmetries we have determined the expectation values of tr XP and tr P X.
Cyclicity of the trace gives another set of relations between operators. However, note that commuting quantum operators may be necessary in applying the cyclic formula. For example, as an operator equation, tr XP 3 = tr P 3 X + i tr I tr P 2 + i tr P tr P + i tr P 2 tr I.
We can further use large-N factorization to replace expectation values of multi-trace operators by products of single-trace operators. Therefore, to leading order in N → ∞, tr XP 3 = tr P 3 X + 2iN tr P 2 + i tr P tr P .
Equations (12), (15), cyclicity of the trace, and reality conditions O † = O * generate all relations between expectation values that we will use for the bootstrap.
As a mini-bootstrap example, consider trial operators I, X, X 2 and P . From the condition (10), the following bootstrap matrix should be positive semidefinite: Unlike in the single-particle case of the previous section, here we are considering operators built from both X and P . We will not be able to explicitly eliminate P in favor of the energy in the matrix case. This is because only tr P 2 appears in the Hamiltonian, but the matrix P can appear in many other combinations. We also note that the expectation value for an odd number of matrices vanishes. Positivity of (20) tells us that where equations (13) and (17) are used. The inequalities (21) are the bootstrap constraints on expectation values in this simple example. When g = 0, tr X 2 = 1 2 N 2 and tr X 4 = 1 2 N 3 , so the last inequality in (21) is saturated and the other two are not.
The bootstrap constraint will be stronger if we include more trial operators. We can do the following. Firstly, take all possible strings of X and P that have length ≤ L. Given these operators, write down the matrix analogous to (20). This matrix must be positive semidefinite. Secondly, we regard each entry in the bootstrap matrix as a variable (which is the expectation value of a single-trace operator with length ≤ 2L), and write down the equalities between them following from (12), (15) (taking O in those equations to be all strings of matrices of length ≤ 2L, any constraints generating operators outside of the trial subset are discarded), cyclicity of the trace (an example is shown in (19)), O † = O * and that the expectation value of an odd number of matrices vanishes. Thirdly, for numerical efficiency, although not strictly necessary, we solve the equalities between observables that are linear (mostly from (12) and (15); note that using cyclicity of the trace and large-N factorization might introduce terms quadratic in expectation values) to obtain a reduced set of linearly free variables. In this way we obtain a number of constraints (quadratic cyclicity and bootstrap positivity constraints) on a reduced set of variables.
Unlike in the single-particle case, we do not necessarily require that the state is an energy eigenstate and the energy E does not appear explicitly in the bootstrap constraints.
At infinite N the matrix quantum mechanics has a continuous spectrum and therefore we do not expect to see separated islands corresponding to excited states. The bootstrap constraints allow for a continuum of energies, as they must. We therefore proceed to use gradient descent to minimize the energy in the allowed region of expectation values. In this way we obtain a lower bound on the ground state energy of the theory. The result is a lower bound because for any number of constraints, certainly the true ground state energy is allowed, and hence the true ground state energy is above the minimum that we find. In Fig. 2 we observe that the lower bound is very close to the true ground state value, already for L = 3, and other observables, such as tr X 2 , are also solved accurately.
Compared to the single-particle example of the previous section, in the matrix quantum mechanics case we cannot determine all expectation values from a finite subset. However the bootstrap constraints seem still to be sufficiently restrictive to solve for few-matrix expectation values. The uncertainties on our result can be obtained by minimizing or maximizing the allowed expectation values while keeping the energy fixed -in this way we determine a possible window for the observable at any given energy. In Fig. 3 this window is seen to be small near the lowest energies, and shrinks as we increase the number of constraints. These facts are consistent with the bootstrap converging towards the correct ground state energy as L is increased, as evidenced already in Fig. 2.

Two matrix quantum mechanics
One matrix quantum mechanics are tractable analytically as one can always diagonalize the matrix. This is not the case for multi-matrix quantum mechanics, where the matrices are not simultaneously diagonalizable in typical states of interest. In this section we will illustrate how bootstrap methods can successfully be used for such theories, focussing on a relatively simple two-matrix quantum mechanics with a global O(2) symmetry (in addition to the large N gauge symmetry). The Hamiltonian is given by with X and Y being N -by-N Hermitian matrices, with conjugate momenta P X and P Y , and m and g coupling constants. This theory is not exactly solvable. An early discussion of the massless (m = 0) limit of the theory is [19]. More generally, by rescaling the matrices we see that physical quantities can only depend on the ratio g/m 2 .
We will impose rotational invariance to obtain even more relations between observables in the two-matrix case. We expect the ground state to be rotationally invariant. The Hamiltonian (22) has an SO(2) symmetry generated by Similar to previous discussions, for states ρ with [S, ρ] = 0, including eigenstates of S, Thus in the two matrix quantum mechanics, equations (12), (15), (24), cyclicity of the trace, and reality relations O † = O * will be used to generate all equations between expectation values that we will use for the bootstrap. The bootstrap then proceeds in exactly the same way as for the case of a single matrix. The results for the ground state energy, tr X 2 +tr Y 2 and tr[X, Y ] 2 are shown in figure 4.
In order to corroborate the accuracy of the L = 4 results, we obtain rigorous upper and lower bounds on the true ground state energy using a trial Born-Oppenheimer wavefunction.
We see in figure 4 that the L = 4 bootstrap results indeed lie within a narrow window allowed by these bounds. This suggests that the bootstrap results are close to convergence. We now briefly describe the trial wavefunction, with details given in the appendices A and B.
The SU (N ) gauge invariance allows us to diagonalize one of the two matrices, say X.
Let the eigenvalues be x i . The Hamiltonian for the entries y ij of the remaining matrix is a sum of harmonic oscillators, with frequencies ω 2 ij = m 2 + g 2 (x i − x j ) 2 . We can therefore write down a Born-Oppenheimer wavefunction in which these oscillators are placed in their ground state: This ansatz can be expected to describe the actual ground state wavefunction if the ω ij become sufficiently large, so that the off-diagonal components of Y are 'fast' compared to the eigenvalues x i . This will not be the case here. However, as we recall in the appendix, Born-Oppenheimer wavefunctions lead to both upper and lower bounds on the ground state energy. The upper bound follows from treating the wavefunction as a variational ansatz.
The lower bound is obtained by finding the ground state of the eigenvalues in an effective potential due to the zero point energy of the y ij oscillators. In Fig. 4 we see that the bounds following from the wavefunction (25) turn out to be remarkably tight. This provides a nontrivial check on the bootstrap results.
The wavefunction (25) is not rotationally invariant, but this is unimportant for bounding the ground state energy. Furthermore, rotational invariance is restored by acting on the wavefunction with the generator S. The advantage of the form (25) is that computing the energy reduces to a single-matrix large N eigenvalue problem, as we now describe.
Both the variational upper bound and the effective potential lower bound amount to minimizing an effective Hamiltonian for the eigenvalues. These effective Hamiltonians, H var and H BO respectively, are derived in appendix A. The quantum mechanics of eigenvalues can be solved using well-established collective field methods [20]. We give details in appendix B. In particular, as N → ∞ the eigenvalue distribution ρ(x) = i δ(x − x i ) must minimize certain functionals, so that the true ground state energy E 0 is bounded by In the appendices we show that with ω(x, y) = m 2 + g 2 (x − y) 2 , while It is straightforward to perform the minimizations in (26) numerically by discretizing ρ(x).
This leads to the results for the energy shown in Fig. 4.
From the results in Fig. 4 one can verify that the ratio N tr[X, Y ] 2 /(tr X 2 ) 2 tends to a nonzero constant at large N g 2 . This means that the matrices do not become commuting in this limit. This can be constrasted with the analogous two matrix integral, with no time, that does become commuting at large N g 2 [21]. This is consistent with the fact that the two matrix integral diverges in the massless limit [22,23], as the eigenvalues spread far apart along the classically flat directions of the potential due to commuting matrices, while the massless matrix quantum mechanics still has a discrete spectrum of normalizable states [24].
That is, the flat directions are lifted quantum mechanically.

Discussion
In summary, we have introduced a systematic numerical method to obtain energies and expectation values of large N matrix quantum mechanics states. The method involves establishing relationships between expectation values and then imposing positivity of a certain matrix of expectation values, in the spirit of [1]. In Fig. 2 we see that the known analytic results for one matrix large N quantum mechanics are readily reproduced. In Fig. 4 we have obtained new results for the ground state energy and expectation values of a two matrix large N quantum mechanics.
There are many natural extensions to our work. The extension to more matrices should be possible with increased computing power or perhaps by optimizing the algorithm. Looking at supersymmetric states in supersymmetric theories may allow for stronger relationships between expectation values, using the supersymmetry generators. Both more matrices and supersymmetry will of course be necessary to tackle the full blown BFSS and BMN theories.
Finally, extensions to Gibbs states (or, to high energy eigenstates) may allow nonzero temperature quantum physics to be accessed with our bootstrap methods. The main challenge here will be to understand the density of states. This could give an alternative probe of the thermal phase transitions studied via Monte Carlo in e.g. [12,13], as well as a new window onto black hole microstates.

Acknowledgements
This work arose from discussions with Edward Mazenc and Daniel Ranard, who also collaborated on the early stages of the project. JK is supported by the Simons Foundation. SAH is partially supported by DOE award de-sc0018134 and by a Simons Investigator award.

A Born-Oppenheimer wavefunction and effective Hamiltonian
This appendix gives details of computations involving a Born-Oppenheimer wavefunction for the two matrix quantum mechanics (22). The wavefunction that we are searching for is a complex function Ψ(X, Y ) of Hermitian matrices X and Y . The state should be SU (N ) gauge invariant and hence for any unitary matrix W ∈ SU (N ), It will be convenient to parametrize such a state with the following set of variables: a diagonal real matrix x i , a Hermitian matrix y ij and a unitary matrix U ∈ SU (N ), such that In these variables we can write down the following Born-Oppenheimer ansatz, in which the y ij oscillators are put in their ground state for a fixed configuration of eigenvalues x i : with ω 2 ij = m 2 + g 2 (x i − x j ) 2 . Equation (31) defines a gauge invariant wavefunction by specifying its values on the gauge slice where X is diagonal. However, we should check that (31) is well-defined because (30) does not uniquely determine x i and y ij as a function of X and Y . Indeed, there is a residual U (1) N −1 gauge symmetry after fixing X to be diagonal: if (31) is invariant under this residual gauge symmetry as well, Ψ(X, Y ) in (31) is well-defined.
To obtain a variational upper bound, we wish to find an effective Hamiltonian for the 'slow' x i degrees of freedom that calculates the expectation value of the full Hamiltonian (22) in the variational state (31). The expectation value of the Hamiltonian in the state Ψ consists of a kinetic part and a potential part: We discuss these in turn. The kinetic energy is Here ∂/∂X ij = 1 2 (∂/∂ReX ij − i∂/∂ImX ij ) are complex derivatives because the matrices are Hermitian. Because the kinetic energy operator is also gauge invariant, the integrand in (33) is constant along gauge orbits. So it suffices to evaluate it on the gauge slice where U in (30) is the identity. Then by the chain rule and (30), at U = I, and Because Ψ is gauge invariant as in (29), ∂Ψ/∂U = 0 so for i = j, Plug (34) and (36) into (33) and evaluate the y ij integrals in the state (31), where ∆ = i<j (x i − x j ) 2 is the usual Vandermonde determinant, with dXdY = ∆dx i dy ij .
The potential term on the gauge slice U = I is and thus Overall the effective variational Hamiltonian on x i , such that Ψ|H|Ψ = ψ|H var |ψ , is therefore The choice of gauge and the form of the ansatz (31) break rotational symmetry. We have done this because it has allowed the problem to be reduced to a single-matrix eigenvalue Hamiltonian (40), which we will be able to solve explicitly. It is possible to restore rotational symmetry by acting on the wavefunction with the generator of rotations. This will not change the energy of the variational state.
From the variational principle we know that the ground state energy of the reduced Hamiltonian (40) is an upper bound on the ground state energy of the original Hamiltonian (22). However, it is well-known that Born-Oppenheimer wavefunctions also give a lower bound on the ground state energy. In the present context (as we prove below) this means that if we drop the final term in (40), the ground state energy of the Born-Oppenheimer Hamiltonian is a lower bound on the ground state energy of (22).
A short proof of this fact is as follows: split the kinetic term into three parts H kin = where E BO (x i ) = N i,j=1 ω ij is the ground state energy of the harmonic oscillator Hamiltonian for y ij 's:

B Large N collective field solution
In this appendix we solve for the ground state energies of the effective eigenvalue Hamiltonians (40) and (41), using the large N collective field method. We thereby obtain an upper and a lower bound for the ground state energy of (22). As is well known, at large N the collective field of eigenvalues becomes classical. We can follow the established steps [20] to obtain the energy as a functional of this collective field. To obtain the Hamiltonian for ρ(x) we must relate the derivative ∂ x i to the conjugate collective variable π(x) = −iδ/δρ(x). The chain rule shows that Plugging these into (41) and defining where P denotes taking the principal value, one finds with We also used the fact that The Hamiltonian in (47) is not manifestly Hermitian. This can be cured by performing a canonical transformation that shifts π by iρ H , resulting in the Hamiltonian, With this Hamiltonian we can straightforwardly compute the ground state energy and certain observables in the ground state. At large N the eigenvalue distribution becomes classical and hence the momentum π(x) vanishes in the ground state. Therefore it is sufficient to minimize the potential energy functional. Using the identity (here π is the irrational number, not the conjugate momentum) this can be written as with Equation (52) is the expression given in (27) in the main text. It must be minimized subject to the normalization constraint dxρ(x) = N and the constraint that ρ(x) be pointwise non-negative. In the large N limit, this normalization combined with balancing the terms in the energy functional and taking the mass to be fixed at order one (recall that the mass can be removed by rescaling the matrices) requires the scaling This is the familiar large N scaling of these quantities. In particular the 't Hooft coupling λ = g 2 N is finite in this limit.
The minimization of (52) is straightforward to perform numerically, by discretizing the integral. With the numerical solution at hand one can evaluate the energy E BO of the state.
These results are shown in Fig. 4 in the main text.
Similarly we can minimize the effective variational Hamiltonian (40) to obtain an upper bound on E 0 . The steps are the same as above, and the functional to minimize is now E var [ρ] = dxρ(x) π 2 3 ρ(x) 2 + m 2 x 2 + dxdyρ(x)ρ(y)ω(x, y) + dxdydzρ(x)ρ(y)ρ(z) (ω(x, z) − ω(y, z)) 2 4ω(x, z)ω(y, z)(x − y) 2 . (55) As discussed in Appendix A, we expect that the true ground state energy E 0 is bounded above and below as We can verify explicitly that these inequalities are obeyed in perturbation theory in small with P (x) a polynomial (whose degree increases order by order in perturbation theory). At large N we obtain (with λ = N g 2 and m = 1) In these expressions we see that the Born-Oppenheimer results only start to differ from the full answer at order λ 2 and that the inequalities (56) are obeyed.
In Fig. 4 of the main text we see that for all couplings the L = 4 bootstrap results lie within a narrow range bounded by (56).
The expectation values tr X 2 and tr[X, Y ] 2 in the trial wavefunction (31) do not provide bounds in the way that the energy does, and therefore we have not included them in Fig. 4. For completeness we note that these expectation values can be computed from the minimizing numerical distribution ρ(x) as The wavefunction (31) is not rotationally symmetric and hence tr X 2 = tr Y 2 in general.