Perturbative Unorientable JT Gravity and Matrix Models

We consider an orthogonal polynomial formulation of the double scaling limit of multicritical matrix models in the $\beta=1$ Dyson-Wigner class. They capture the physics of 2D quantum gravity coupled to minimal matter on unorientable surfaces, otherwise called unoriented minimal strings. We derive a formula for the density of states valid to all orders in perturbation theory. We show how to define an interpolation between the multicritical models and that a certain interpolation among an infinite number of them provides an alternative definition of unoriented JT gravity. We discuss the strengths and weaknesses of our formulation.


Introduction
The double scaling limit of random matrix models can be used to formulate many two dimensional (2D) Euclidean quantum gravity theories.The success of the double scaling limit in describing smooth macroscopic geometries (i.e.gravity) is most easily understood in the context of the 't Hooft large-N limit of matrix/gauge theories [1][2][3][4][5][6][7][8][9][10].In a random matrix model, the ribbon diagram expansion of the partition function is interpreted as literally enumerating tessellations of surfaces of increasingly more complicated topology, with 1/N (where N is the size of the matrix) determining the topological expansion parameter.Double scaling the matrix model then amounts to sending N → ∞ while simultaneously shrinking the characteristic area a of a face in a given tessellation.The result is a sum over smooth closed surfaces indexed by their Euler number.
Random matrix theory has a rich history in both physics and math.Its original use by Wigner to describe nuclear physics embodies its statistical nature, where the fundamental object (the matrix) is drawn from a randomly distributed ensemble with probabilities being determined by the matrix potential [11].While the interpretation in the spirit of 't Hooft connects more directly with geometry and gravity, it appears the Wignerian approach to random matrix theory may shed light on non-perturbative aspects of gravity [12,13].Nevertheless, the interplay between gravity and random matrix theory also provides a way to expand the math literature on the topic.
A particular 2D theory of interest lately, Jackiw-Teitelboim (JT) gravity, has shed a lot of light on properties of quantum gravity [14][15][16].Classical solutions of JT gravity on a closed manifold describe rigid hyperbolic spaces.It is more interesting to consider the theory on a manifold with at least one boundary, in which case the boundary dynamics are described by the Schwarzian theory [16].In [17] Saad, Shenker, and Stanford showed that perturbation theory in JT gravity is captured by the double scaling limit of a certain matrix model, which can be cast in terms of a limit of minimal string models [18].It was subsequently shown in [19] that the matrix model describing JT gravity perturbation theory is a special interpolation between an infinite number of multicritical matrix models.Those models have reliable non-perturbative completions, making it possible to fully define JT gravity (and supergravity) in a non-perturbative way [20][21][22].
Generalizations of JT gravity (including orientability, defects, supersymmetry, etc.) have been studied from both a field theory [23][24][25] and random matrix perspective [19,26,27].Considering unoriented spacetimes in Euclidean quantum gravity is interesting for several reasons.First, an unoriented bulk spacetime is appropriate for describing a boundary theory with time reversal symmetry in the context of the AdS/CFT correspondence.Second, although the Gaussian Orthogonal Ensemble (GOE) and Gaussian Symplectic Ensemble (GSE) have been studied extensively in the math literature (see for example [28]), the double scaling limits of the more general Wigner-Dyson β-ensembles is not as investigated.Other descriptions of unoriented gravity in two dimensions can be found in [29][30][31].Since these models describe unoriented minimal strings any effort in this direction provides nice contact between math and physics.
In [26] Stanford and Witten discuss unoriented JT gravity and its connection to random matrix theory.They develop perturbation theory using the loop equation formalism, which focuses on the resolvent operator of the random matrix.By applying techniques developed in [17] they conclude that the presence of even a single cross-cap can cause divergences in otherwise well-behaved objects.See [32] for recent progress relating to the loop equation formulation of unoriented JT gravity.An effort is made in this work to define perturbation theory from a different perspective, in the hopes that a more matrix model-focused approach may resolve the apparent issues with the loop equation formalism.This alternate point of view is met with mixed results, as it turns out to be difficult to extract results from the formalism even for the simplest non-trivial minimal model.
This paper details two methods by which one can analytically compute the large-N and double scaled eigenvalue densities for β = 1 ensembles with arbitrary polynomial matrix potential.Results for the Gaussian case are known already (see e.g.[28,33]) and progress has been made in the direction of the string equations for more general models [34,35].
The first method we discuss employs intrinsic properties of the aforementioned orthogonal polynomials to compute the n-point correlators of the eigenvalues.We argue that the result for the GOE found in [36] is nearly valid for arbitrary matrix potential, missing only surface term corrections.This formula for the density directly involves the orthogonal polynomials, so the pre-existing double scaling technology can be applied with the inclusion of some extra ingredients.In particular, in [34] the authors showed how to double scale the derivative with respect to the matrix eigenvalue λ d dλ The operator P k is the Lax pair of a Sturm-Liouville operator H = −∂ 2 + u, familiar in the KP/KdV hierarchy; in particular it generates the k th KdV flow of the operator H.By imposing an auxiliary condition called the string equation one finds that , making P k a (2k − 1)-order differential operator with polynomials in u and its derivatives as coefficients.
The second method we discuss utilizes a new set of polynomials.This set of functions has a sort of orthogonality relation with respect to a skew-symmetric bilinear form called a "skew inner product."For this reason the functions are called skew orthogonal polynomials.In [37] the authors derived an expression for the finite-N eigenvalue density of a β = 1 matrix model in terms of these polynomials.By adapting results from [34] we present the double scaling limit of the density (eqn.(4.7)).The result involves factorizing P k in terms of a new set of functions g i , making contact with the KdV hierarchy's relationship to more general Drinfeld-Sokolov hierarchies [38] The functions g i are used to define differential operators T k and S k that facilitate a transformation between the skew orthogonal polynomials and a set of associated orthogonal polynomials.Equations of motion for the g i are determined in terms of u by comparing the factorization to In analogy to the β = 2 models, a method for implementing an interpolation between minimal models is proposed.The Lax operator P of a such an interpolation is a simple linear combination of the individual Lax operators.The functions g i do not add linearly in the interpolation, but instead have to be determined by factorizing the entire Lax operator P .Utilizing this, we obtain a formal expression for the eigenvalue density of unoriented JT gravity (eqn.(6.6)), interpreted as an interpolation between an infinite number of these unoriented minimal models.
The rest of this paper is organized as follows.In section 2 we review the conventional orthogonal polynomial approach to matrix models.In section 3 we review known results in the finite-N study of β = 1 Wigner-Dyson models and discuss generalizations of the GOE results to more general matrix potentials.In section 4 we implement the double scaling limit and obtain a formula for the eigenvalue density, valid to all orders in perturbation theory.The equations of motion for the new functions g i , called new string equations, are discussed.In section 5 we perform some example calculations in the first few multicritical models to demonstrate the use of our results.In particular we show how to perturbatively determine the functions g i in terms of the function u for individual multicritical models as well as massive interpolations between them.We then use the WKB approximation to compute perturbative contributions to the double scaled density.A peculiarity involving the GOE is investigated along these lines.Finally, we apply consider the simplest nontrivial unoriented minimal string, which provides a way to test the interpolation formulae.In section 6 we apply our results to obtain a formluation of unoriented JT gravity as an interpolation between an infinite number of minimal models.A formal expression for the eigenvalue density is provided.

Matrix Models and Orthogonal Polynomials
We begin by reviewing some basics about random matrix theory (RMT) and establishing some conventions.For a review in the context of two-dimensional gravity see [39], and for more mathematical reviews [28,33].
Given a set of N × N matrices {M }, its corresponding random matrix probability density is determined by a function V called the potential and is given by e −N tr V (M ) dM , where dM refers to the flat measure on the space of N 2 variables.We will take V to be an even polynomial of rank 2k, and always assume that N is even.The elements of M are random variables, distributed according to this probability density.
An object of central importance in our endeavor is the matrix integral Z given by with respect to which scalar functions of the random matrix, called observables, can be computed as expectation values.It is common for the set of matrices in consideration to have a symmetry, for example being hermitian or real symmetric.In these cases there is a sort of gauge redundancy in Z that can be dealt with by diagonalizing M .Performing the change of variables U , with U in the appropriate symmetry group, we can rewrite the partition function after integrating out U where ∆(λ) = i<j (λ j − λ i ) is the Vandermonde determinant.The numerical prefactor C is related to the volume of the symmetry group and so we will drop it.The power β is determined by the symmetry group used to diagonalize M : for the group U (N ) β = 2 and for O(N ) β = 1.The domain of integration for the eigenvalues is also affected by the value of β.Since N is finite, we can assume that the eigenvalues are ordered.The matrix integral can be defined with a domain of integration that respects the ordering in order to not over count, and to eliminate the absolute value around ∆.However, for even values of β we can extend the integration to all of R N by dividing by a combinatorial factor to take into account relabelling the indices anytime λ j − λ i changes sign.Having each eigenvalue take any real value is convenient, as will be explored shortly.We are not so fortunate when considering odd values of β though, and must use alternate methods.
It can be shown that ∆ is the determinant of a matrix consisting of powers of the eigenvalues, ∆(λ) = det λ j i .By taking linear combinations of the rows with real coefficients, we can introduce a family of polynomials p i (x) so that ∆(λ) = det p j (λ i ).The polynomials will be normalized so that p i (λ) = λ i + • • • .We stress here that the introduction of these polynomials is possible for any value of β, and solely depends on the presence of the Vandermonde determinant.
A convenient choice for the polynomials p i is to make them orthogonal with respect to a power of the measure e −V , for example This is particularly helpful when the values of λ can be expanded to all of R. Note that the family of these polynomials is infinite, which will facilitate the large-N limit to be taken later.Define the oscillator wavefunctions The choice of name has a two-fold motivation.In the Gaussian matrix model, the orthogonal polynomials are the Hermite polynomials and the functions ψ i are identical to the normalized wavefunctions of the harmonic oscillator.Additionally, we have the Fermi gas, or log-gas, perspective that arises from interpreting the Vandermonde determinant as a logarithmic interaction potential.The eigenvalues behave like fermions and in a Fock space representation each one is represented by a sort of oscillator.
While formally the index i on ψ i is the same as the index on λ i , we know that even for a finite number of eigenvalues there are infinitely many wavefunctions.So, consider a many-body system comprised of N identical fermions.The wavefunction of a single fermion is represented by ψ i , where the index i indicates the energy level occupied by that particular particle.Since the fermions are identical only a single particle can occupy the i th energy level.By virtue of the fact that the Vandermonde determinant in the matrix model partition function can be written in terms of the orthogonal polynomials, so too can it be written in terms of the first N oscillator wavefunctions.Hence the matrix model takes into account the first N energy levels of the fermions.In a fermionic many-body system, the energy at which all lower energy levels are occupied is called the Fermi level, and the underlying energy levels the Fermi sea.As such, the index value i = N − 1 is called the Fermi level1 , and all lower values the Fermi sea.
The oscillator wavefunctions provide an orthonormal basis for the Hilbert space L 2 (R) with the standard inner product.Denote the basis as {|ψ i ⟩} ∞ i=0 .Operators on L 2 can thus be decomposed in terms of these functions in the usual way (2.5) By Ô we denote the action of the operator on the actual function ψ i (x).
The polynomials p i satisfy the following recursion relation with The corresponding recursion relation for the oscillator wavefunctions is The self-adjoint operator Λ that implements this transformation on L 2 has matrix elements A second operator of interest implements the derivative, which we denote by L in the finite-N regime.One can show that the matrix elements L ij are nonzero in general only when |i−j| ≤ 2k−1, are determined by the recursion coefficients R i , and are anti-symmetric because L is anti-self-adjoint.The operators Λ and L satisfy the canonical commutation relation [Λ, L] ∝ 1.
Since V is an even function it can be shown that each polynomial p i has definite parity, given by p i (−λ) = (−1) i p i (λ).This induces the Z 2 -grading L 2 = span{|ψ 2α ⟩} ⊕ span{|ψ 2β+1 ⟩}.Thus we can conveniently decompose our operators into 2 × 2 block form.For example, both Λ and L are odd under the grading, so we designate (2.9) Objects that are labelled 0, . . ., N −1 will carry a Latin index i, j (but note that k is reserved for the order of the matrix potential).When considering Z 2 -grading we will decompose some of these objects into even and odd parts, and it will be convenient to write i = 2α or j = 2γ + 1, respectively.Objects that have definite Z 2 parity will thus carry a Greek index α, γ (but β is reserved to denote the type of matrix ensemble).For example, take the operator L. By definition the operator c has matrix elements given by c γα = −L 2α,2γ+1 .Since the matrix representing L has 2k − 1 non-zero offdiagonals, the matrix representing c will have 2k − 1 non-zero off-diagonals (and the main diagonal will be nonzero as well, which is not true for L).The action of the derivative on the even wavefunctions ψ 2α can be written For this range of γ in the sum, some of the matrix elements c αγ will be zero depending on the value of α, but these bounds of summation are guaranteed to produce all necessary non-zero terms.
The double scaling limit [2,3,5,7] is a combination of taking the size of the matrix N → ∞ while also focusing the objects of the theory to neighborhoods around specific values.In this limit the small parameter 1/N has a renormalized form ℏ. The index i is replaced by a continuous variable x ∈ (−∞, µ], where µ is the Fermi level.The large-N limit confines the spectrum to a single interval2 , and in the double scaling limit we zoom in on one of the edges of the ineterval, replacing λ with a new variable E ∈ [0, ∞).The oscillator wavefunctions become functions of ψ(x, E), the recursion coefficients R i are replaced by a function u(x), and the operator Λ that implements the recursion becomes a Schrodinger Hamiltonian H = −ℏ 2 ∂ 2 x + u(x).All together, the recursion relation becomes the Schrodinger equation ( The derivative operator L = d dλ is replaced by a differential operator with respect to x, called P k .In the Lax formalism of integrable systems, the operators (H, P k ) form the Lax pair of the k th model in the KdV hierarchy, where P k is determined by requiring that the commutator takes a special value 3 .The unique differential operator satisfying that condition is given by P , where the + subscript denotes keeping only non-negative powers of the derivative ∂ x .See Appendix B for more details about pseudodifferential operators.
If one is interested in solving the Schrodinger equation, the function u must be determined first.In Hermitian matrix models, the equation of motion is precisely the canonical commutation relation [P k , H] = 1, which is the double scaled version of the relation satisfied by Λ and L mentioned previously.The outcome is expressed in terms of the Gelfand-Dikii differential polynomials R k .The equation of motion, called a string equation for historical reasons, is (2.12) Two important examples of Gelfand-Dikii polynomials that will be used later are More general models can be defined by introducing coupling constants t k and considering the equation of motion Such a model is called an interpolation between the individual multicritical models.This can be obtained from the canonical commutation relation by defining the generalized operator P = k t k P k .

The β = 1 Wigner-Dyson ensembles
To construct an enumeration of unoriented surfaces using a matrix integral, we restrict ourselves to real symmetric matrices H, which falls under the β = 1 O(N ) symmetry class.The matrix H is in the bi-fundamental representation N⊗N of O(N ).The two fundamental representations are identical, making the indices of H indistinguishable.As such, it will be more convenient not to keep track of upper and lower indices.
To motivate considering this theory as unoriented, we will briefly touch on the 't Hooft diagrams of the theory.The Gaussian action can be cast in propagator language This implies that the propagator is 2 N (δ in δ mj + δ ij δ mn ) for β = 1 ensembles.The second term comes from the fact that the indices are indistinguishable for our matrices: When one adds higher order terms to the potential and performs a diagrammatic expansion of Z, the presence of the new term in the propagator adds twists to the diagrams, which can only be drawn on unoriented surfaces4 .In the large-N limit, observables in the theory will have a topological 1/N expansion.For instance, the free energy defined by Z = e F can be expanded as The number χ(g, c) = 2 − 2g − c is the Euler characteristic of the diagram, corresponding to a tessellation of a surface with g handles and c crosscaps.In the double scaling limit, the quantity5 F g,c represents the contribution to the free energy of that topology.Since β = 1 is odd, we must take care to specify a domain of integration in Z that allows us to remove the absolute value on ∆.We write where General n-point functions of the eigenvalues can be computed using an idea pioneered by Dyson [40] and applied to the GOE by Mehta and Mahoux [36], called integration over alternate variables.The analysis begins by introducing two arbitrary functions u and v, and considering the following expectation value with respect to the matrix ensemble Further progress will be made by judiciously introducing a set of polynomials whose properties simplify the calculation.There are two natural choices one can make, which are the subjects of the following subsections.

Orthogonal Polynomial Approach
First, introduce a set of orthogonal polynomials Notice that the integration measure includes a different power of e −N V than used previously.Unless stated otherwise, the polynomials p i will be orthogonal with respect to this new measure.Consequently, the oscillator wavefunctions will be defined as Introduce the functions By absorbing the factors of v into the determinant and taking linear combinations of the columns of the matrix, it can be shown that G N is given by where C N is a numerical factor depending on the normalizations h i and the matrix m is given by The integrand is now symmetric under the swap of any two variables (due to the invariance of the determinant under swapping an even number of rows and columns), so we can enlarge the domain of integration to (−∞, ∞) for each eigenvalue at the cost of a numerical prefactor, which we disregard.Subsequently, the integrals over each eigenvalue are decoupled and G N can be written where and the brackets denote anti-symmetrization.This determines the correlator G N as the Pfaffian of J.
By taking linear combinations of the rows and columns of J, we can put it in a 2 × 2 block form.Begin by setting J 2α,2γ (u, v) = f αγ (u, v).Next, define (3.12) Notice that the values of the dummy index σ go outside of the "N × N block" in index space in which the matrix model technically lives.The impact of this on the eigenvalue density will be discussed shortly.Finally define After some rearranging to ensure antisymmetry, we find the Pfaffian can be written and hence the correlator G N can be expressed in terms of this new determinant.If u, v are both even functions, the diagonal blocks of this new matrix are zero and the Pfaffian reduces to det(g).We now use the two functions u, v as sources and treat G N as a generating function for the eigenvalue correlators These correlators are related to G N via [33] G In particular, this means that to compute the eigenvalue density ρ 1 (λ) ≡ ρ(λ), we take one functional derivative of G N and set the source to 0. Mehta shows that for the GOE The result (3.17) nearly holds for k > 1. Denoting the density for the kth multicritical model by ρ(k) , the result is really ρ(k) = ρ(1) + surface terms, (3.18)where the oscillator functions used in ρ are k-dependent as well.The additional terms, called surface terms here, are localized around α = N/2 and arise from the fact that we included terms outside of the Fermi sea to obtain ψ ′ 2α when rewriting the matrix J.The number of extra terms is set by k in the upper limit of summation in (2.10).The surface terms, denoted by K(λ), have the form In the simplest case k = 2, It is not possible to determine any of the quantities in this formula in closed form for finite N .

Skew Orthogonal Polynomial Approach
While the orthogonal polynomial approach bears some fruit when applied to the β = 1 models, there is another set of functions more naturally suited to solving the problem.We begin again with (3.3) and compute the correlator G N (u, u), setting u = v for convenience.We still arrive at the expression where now we define The polynomials q i (λ) are monic and chosen to obey the following relation [34,36,37] where ε(λ − λ ′ ) = sgn(λ − λ ′ ), and with all other ⟨q i , q j ⟩ s = 0.The bilinear form ⟨•, •⟩ s is skew-symmetric, and hence is referred to as a skew inner product.The polynomials q i are called skew orthogonal polynomials.The correlator G N (u, u) is once again given by the pfaffian of a matrix J given by Jij = ⟨uq i , uq j ⟩ s . (3.24) If we set u = 1, this computes the full matrix integral.
In analogy to the oscillator wavefunctions ψ i defined above, introduce the skew oscillator wavefunctions Notice that the matrix J has the 2 × 2 block structure where (3.27) The functions f, g and µ referenced here are not the same as the functions with the same names used in the orthogonal polynomial approach.Consider the function u = 1 + a where a is small.Then the 2 × 2 structure of J can be written (schematically) as with ϵ small and related to the function a.The pfaffian pf J admits a straightforward expansion in ϵ (see Appendix A.7 in [33]).The density can once again be extracted from G N (1 + a) by taking a functional derivative with respect to a.The result is [36,37] This expression for the eigenvalue density is exact for any value of N .This is in contrast to (3.17), which required the inclusion of k-dependent surface terms to be correct.The construction of skew orthogonal polynomials and their properties are explored in detail in [37].However, the ultimate goal of this work is to reach the double scaling limit, and there does not appear to be a process through which one can directly do this to the skew orthogonal polynomials.This difficulty is circumvented by once again appealing to orthogonal polynomials.
The orthogonal polynomials p i introduced above form a basis for the ring of polynomials.The skew orthogonal polynomials q i can be expanded in terms of the p i : there exists a lower triangular matrix O with 1's on the diagonal such that (3.30) Representing the skew oscillator wavefunctions with the kets |ζ i ⟩, the defining relation for the polynomials q i is then given by where ε is the integral operator with kernel 1 2 ε(λ − λ ′ ) and the L 2 inner product is used.An implication of this is that the skew inner product can be implemented by an operator acting on L 2 with the measure e −2N V .By taking into account the normalizations, the skew oscillator wavefunctions can be expressed in terms of the oscillator wavefunctions in a way related to (3.30) (3.31) The skew oscillator wavefunctions have the same Z 2 -grading as the oscillator wavefunctions, so the matrix Õ must be even under the Z 2 -grading.Define the matrices a, b so that Õ = a 0 0 b .
In [34] they show that a, b are related to the derivative matrix c by Hence the eigenvalue density is written (schematically) as To conclude, it should be noted that ρ is the diagonal part of a larger function [37].Define The analog of the self-reproducing kernel familiar from studies of β = 2 theories is given by where The eigenvalue density is given by ρ(λ) = lim where qdet denotes the quaternionic determinant (see [36], for example).The result follows from the fact that I(λ, λ) = 0.

Previous Double Scaling Results
The term "double scaling limit" used in this paper corresponds to the "soft edge" of an eigenvalue distribution in the math literature.The soft edge limits of numerous matrix models with Gaussian potentials, including the β = 1, 2, 4 theories, can be determined by using known facts about the Hermite polynomials.See [24] for a review.In particular, the double scaled eigenvalue density for the GOE can be expressed in terms of the Airy function where ρ Airy is the eigenvalue density for the double scaled GUE (also known as the Airy model).
The double scaling limits of the operators a, b introduced above were computed in [34].In the normalization used here, they are The double scaled version of (3.32) is S † k T k = P k .The authors of [34] outline an algorithm to determine the coefficients in T k , S k which we will describe below.Using this, they compute the free energy of a pure unoriented gravity theory (corresponding to the single minimal model k = 2).

Double Scaled Density
We begin by double scaling the modified version of the finite-N GOE density given in (3.18).First, in the strictly large-N limit we change the index i into a continuous variable X ∈ [0, 1].Then introduce a small parameter δ that picks out the scaling parts In particular in the large-N limit, the eigenvalue distribution for these models settles into a finite window [−λ c , λ c ] which, due to the symmetry of the theory, is symmetric about the origin.The full scaling ansatz for the eigenvalue is λ = λ c − δ 2 E; using this as a coordinate transformation implies that integrals over λ have the scaling part where the bounds of integration are updated on a case-by-case basis.
By changing the large-N integration over α/N to i/N we introduce a factor of 1 2 .The double scaled density is then The double scaling limit of the surface terms from (3.18) are denoted by F k [ψ].Since before double scaling there are O(1) surface terms around α = N/2, F k [ψ] should involve derivatives of ψ evaluated at the Fermi surface x = 0.For larger k, there are more surface terms before double scaling and hence higher order derivatives in the double scaling limit.
In order to double scale the β = 1 density written in terms of the skew oscillator wavefunctions (3.29) we need to discuss the double scaling limit of the operator ε.Consider the application of ε to functions with definite parity.The results are where we've assumed that f only has support on [−λ c , λ c ].By the scaling ansatz for the λ integral, The formula for the finite-N density in (3.29) involves the skew wavefunctions, which are related to the oscillator wavefunctions via the transformations a and b as depicted in (3.33).The operators a, b double scale to differential operators T k , S k , which are order−(k− 1) and k respectively.The double scaled limit of the relationship between ζ and ψ is expressed as6 even : Putting together these pieces, we find the expression for the double scaled β = 1 density (4.7)This expression is exact and does not ignore any surface terms.The normalization is chosen to extract the finite piece after double scaling.
The full expression for the density gives us the opportunity to roughly determine what the surface terms were in the orthogonal polynomial approach.Consider the first term in the above formula.If we were to integrate by parts with respect to x in order to shift T k inside the E ′ integral, we would pick up a contribution at x = 0 and get the term where the adjoint is taken in L 2 (R).Now, the specific combination T † k S k is proportional to Pk , which is canonically conjugate to the Hamiltonian H for which ψ(x, E) is an eigenfunction.Thus P k has a representation P k ∼ ∂ E , and which produces the ψ 2 term familiar from the study of β = 2 theories.
Consider the second term in the density.Once again we integrate by parts to shift the location of Tk .This picks up more terms at x = 0 and gives Being more careful about signs, we conclude where the surface terms are generated by moving the operator T k around.In the case k = 2, One thing that is obscured by this analysis is what happens in the GOE.For that theory, where k = 1, the operator T 1 is a constant, meaning it costs nothing to move it around inside the integral.Thus there are no surface terms in that case We can confirm numerically that (4.13) matches the previously known result given in eqn.(3.38) and the plots are shown in (1).

New String Equations
It is evident from (4.7) that we have to determine the operators T k , S k , and there is a well defined algorithm to determine them outlined in [34].For the k th minimal model we can write the following factorization equation 7 The operator P k is (2k−1)-order differential operator where the coefficients are functions of u and its derivatives.Notice how this differs from writing P k = S † k T k with T k , S k previously defined in eqn.(3.39).There, the two individual operators were defined in terms of separate families of functions.which we called g i , h i .By imposing that the factorization be anti-selfadjoint and setting some integration constants to 0 (see [34]) one can reduce the number of functions so that T k , S k depend on the same family of functions.The new factorization 7 Technically we will be working with the operator −iP k , if . This factor of i can be restored by modifying the factorization in terms of T k and S k , but the functions gi would remain unchanged, and so would all of our other results.The details can be found in Appendix B. in eqn.(4.14) amounts to setting Expanding both sides of (4.14) and equating powers of ∂ we get (2k −1) coupled differential equations for (k − 1) g i functions.These equations will not all be independent, and we can find (k − 1) independent differential equations.While we do not have exact forms for each equation obtained using the algorithm outlined above for arbitrary k, it is straightforward to show that the equation determined by ∂ 2k−3 is always8 [34] (see Appendix B for details) Moreover the equation determined by ∂ 2k−4 is simply the derivative of the above equation.This implies that g 0 ∼ √ u 0 , even for more complicated massive interpolations.
The preceding algorithm can be adapted to consider an interpolation up to the k th model with coupling constants t i .Since the order of P i as a differential operator increases with i, the total order of the superposition P = i t i P i as a differential operator is set by k.Since (4.16) is determined by ∂ 2k−3 , the highest order term with a nontrivial coefficient, the only two operators that can contribute to that equation are P k and P k−1 .Therefore the proper way to modify (4.16) is (4.17) The right-hand side is unchanged because the total differential operator being factorized has order 2k − 1.The other independent differential equations are determined by looking at the coefficients of lower order derivatives.The next one comes from ∂ 2k−5 , and therefore will explicitly include the coupling constants t k , t k−1 , and t k−2 .

Topological and Pure Unoriented Gravity
The simplest minimal model, k = 1, is dual to topological gravity.Although it is not a genuine theory of surfaces it is interesting to study in its own right, and provides the opportunity to do many analytical calculations.The first nontrivial minimal model, k = 2, is dual to pure gravity.

The GOE
We can use the fact that the wavefunctions ψ(x, E) and the eigenvalue density ρ are both known exactly in terms of the Airy function to test some of the framework above.Note that the GOE does not require an extra function g to compute the density.
The Schrodinger equation for the k = 1 model is9 Taking the standard WKB form of the wavefunction, we find the two possibilities (5.2) Of course, the differential equation can be solved exactly in terms of the Airy function, whose Taylor series expansion in small-ℏ is a linear combination of the two perturbative solutions obtained above ( The WKB expansion is inherently an expansion in small ℏ.However, the non-perturbative oscillating piece has to be dealt with since it does not simplify under the assumption that ℏ ≪ 1.We circumvent that with the following trick.The formula for the density is ambiguous as it does not allow for the possibility that the wavefunction ψ is complex.The formula's derivation from the matrix model does not fix this ambiguity, so it may be done ad hoc.The simplest option is to replace one wavefunction in each term with its conjugate.Taking the permutations of which wavefunctions get replaced ensures the outcome is real.Upon doing this, the non-perturbative pieces generally cancel out.
For the Gaussian theory, T ∝ 1, S ∝ ∂ ∂x . ( Using the GOE wavefunctions ψ(x, 2 √ 2E) one finds We note three things about the outcome.First, the disk level density has the normalization one would expect from the saddle point analysis (see appendix).Second, the O(ℏ 1 ) term has a different coefficient than the β = 2 density, but the same functional dependence on E. The two are not related to each other by the simple rescaling of E -this is not an issue because we do not expect the theories to be related in a simple way past the leading order in perturbation theory.Third, the O(ℏ 0 ) term has a coefficient of 0 according to the calculation.If this were truly a (perturbative) topological expansion in a theory of unoriented surfaces, then the surface with one boundary and one crosscap (i.e. the Mobius strip) would appear at this order.There is a non-vanishing, non-perturbative contribution to the O(ℏ 0 ) integral, but it actually contributes at O(ℏ 1 ).As a corollary of this, there are no terms in the perturbative series at O(ℏ 2n ).As we show in the appendix, the absence of the crosscap term is not unique to the WKB approach, but is also seen in the loop equation formalism.

The
Take an interpolation between k = 1 and k = 2 with coupling constants t 1 and t 2 .The momentum operator for this interpolation is given by There is only one independent function g in the factorization of P , so The independent equation of motion for g is (5.8) Each function will have a perturbative expansion in ℏ10 The first orders in the expansion of g are given by (5.10) It is easy to show that as t 1 → 0, the perturbative expansion smoothly transitions to the single k = 2 case.
In the following two subsections we will use these results to study the k = 2 theory, which is dual to pure unoriented gravity, and the (2, 3) unoriented minimal string.

Pure Unoriented Gravity
Consider just the k = 2 minimal model.The leading contribution to the density of states for the k th minimal model is proportional to E k− 1 2 , so we expect to see E 3/2 behavior.The density admits a topological expansion in powers of ℏ, where g is the number of handles of the surface and c is the number of crosscaps.By convention we define ρ 0,0 ≡ ρ 0 .By keeping the leading order contributions from T 2 , S 2 we find the disk density Isolating the first subleading contributions gives the crosscap density ( The functions g 0 and g 1 are obtained from (5.10) by setting t = 0.Both (5.12) and (5.13) contain more information than is desired, unless we compute the wavefunctions perturbatively using WKB analysis.The leading order contribution to u, found by solving R 2 + x = 0 (see (2.13) for the definition of R 2 ), is u 0 = √ −x.The WKB analysis yields Unfortunately it is difficult to maintain neither analytical nor numerical control over the calculation of the density using the WKB wavefunctions.In particular, the combination of the complicated E-dependence of the WKB wavefunction and the presence of the E integrals in the formula for the density renders analytical computation very difficult.Ordinarily in a Hermitian matrix model, one would use an averaging argument to get rid of the oscillatory part of the WKB wavefunction to compute perturbative contributions to the density.However in this case the E-integrals interact with this oscillatory functions and impact the counting of powers of ℏ in a nontrivial way, as is evident from performing the analogous calculation in the GOE.
As a last resort we appeal to the power of the non-perturbative numerical framework described in [22], which we summarize here.It has been known since the early 90s that ensembles of positive Hermitian matrices provide a viable nonperturbative completion for models of two-dimensional quantum gravity (see [41] for example).The string equation for u in such a model is modified to be where R = k t k R k + x.Equation (5.15), with the coupling constants specified, can be solved numerically by putting the system in a box and imposing the boundary conditions (5.16) Then the Schrodinger equation for this potential is discretized and solved using the methods of [42].
In this work, a spatial grid size of ∆x = .03is used.The presence of the walls needed to obtain numerical solutions constrains the maximum trusted energy.The average spacing between energy levels is ∆E ≈ .02.
The numerically determined wavefunctions are inserted into eq.( 5.12) to produce the disk level density in fig.(2).The curve is oscillating because the wavefunctions are non- perturbative.In a more complete non-perturbative framework the oscillations reveal the discrete nature of the quantum system [12].However, here they are meaningless because the formalism is only meant to be trusted at the level of perturbation theory.
If we plot ρ 0 up to higher energies we see a deviation of the curve from E  In this plot, that cutoff was roughly E max = 14.This argument can be solidified by checking the plot of the function with different cutoffs, i.e. upper limits, as shown in figure (4).
We can similarly plot ρ 0,1 by using the numerically determined wavefunctions in eq.(5.13).To compare the values of ρ 0 and ρ 0,1 , we plot them together in figure (5).We can see that ρ 0 is more dominant, and ρ 1 is just like a damped oscillation around the E axis, as it should be for a perturbative correction.The sum ρ 0 + ρ 0,1 still oscillates with the approximate E 3/2 behavior of the leading contribution.
The amplitude of ρ 0,1 is decreasing as E → ∞.The topological expansion of the density in (5.11) is equivalently an expansion in large E, so at large energy we expect to recover the disk contribution.The (2, 2p−1) oriented minimal string is described by a matrix model with the coupling constants [18] (5.17) In particular, for the choice p = 2 the nonzero couplings are t 1 = 1 and t 2 = 4π 2 /9.The leading order contribution to u is found to be (5.18) It is important to consider this model for two reasons.First, it represents a proof-of-concept for the interpolation technology described here, that we apply to unoriented JT gravity in the following section, because the leading eigenvalue density can be computed analytically.Second, it is well known that JT gravity is in some sense a p → ∞ limit of the (2, 2p − 1) minimal string (see [18], for example).Thus studying the first nontrivial interpolation in this family is a step toward the full theory.
It would be harder than the pure unoriented gravity case to maintain control over the WKB wavefunction in an analytical calculation, so we once again resort to the nonperturbative framework for numerics.Using the numerically determined wavefunctions produces the disk level density shown in fig.(6).The expected analytical result, given by is shown for comparison.
The displayed energy window is smaller in this case than the preceding subsection due to the numerical procedure.The cutoff on the energy in this data is E max ≈ 8.The first perturbative correction can also be calculated using (5.13).The results are displayed in fig.(7).

Unoriented JT Gravity
As shown in [12,19,22,27], JT gravity is perturbatively defined as a massive interpolation between an infinite number of multicritical models.The full equation of motion obeyed is where the couplings t k are determined by the leading density of states of the gravity theory and R k is the kth Gelfand-Dikii polynomial.
As is clear from the formalism developed above, the function u is seeded from the β = 2 theory into the β = 1 theory to determine the functions g i .Further, the communication between the two theories is facilitated by the factorization of the differential operator P , the Lax pair of the Schrodinger Hamiltonian.The operator P is a linear combination of P = k t k P k for an interpolation with coupling constants t k .Since the coupling constants are determined by the leading density of states, which is the same in the oriented and unoriented theories (up to normalization), the full operator for unoriented JT gravity will be We can, in principle, define operators T , S such that P = S † T as before, and which map the wavefunctions into the skew wavefunctions.Schematically, the factorization looks like The flipped bounds on the right-most product indicate that the functions g i appear in reverse order in that part In the same manner as before, we find enough differential equations by comparing coefficients of ∂ to determine the g i .
Given the extension of the operator P to the infinite interpolation that defines JT gravity and its formal factorization in terms of operators T and S, we predict that (4.7) remains true for unoriented JT gravity, giving the density ) to all orders in perturbation theory.This formula is both self-consistent and consistent with the infinite interpolation used to define JT gravity in terms of multicritical matrix models.
This theory and orientable JT gravity share the issue of having formally infinite derivatives in the definition.The string equation for u is infinite order, the operator P is infinite order, and so too are T and S. The solution is to notice that t k gets smaller and smaller for higher k and lim k→∞ t k = 0.In [21] it was noted that a truncation can be made to get reliable numerical results depending on the maximum energy one is interested in probing while computing the eigenvalue density.Practically, one studies the finite interpolation up to some k max by finding the function u satisfying 11 kmax k=1 t k R k + x = 0, where t k are the same coupling constants for JT gravity (6.2).Then the Schrodinger equation for u can be solved for the wavefunctions.
We believe that, in principle at least, the same can be done in unoriented JT gravity.The algorithm for computing the density is as follows.First, perform the truncated analysis to find u and the wavefunctions ψ.Then determine the interpolated operator P ≈ kmax k=1 P k in terms of u, and use that to find the k max − 1 functions g i .Then build the truncated operators T and S out of the g i , and finally use the formula for the density (6.6).

Discussion
The purpose of this paper was to further explore the study of double scaled β = 1 Wigner-Dyson models initiated in papers like [34,43].Recent advancements in twodimensional quantum gravity have brought to the forefront statistical properties of the theories, with a large part being played by the eigenvalue density.One model, JT gravity, has been of particular importance, and the last several years have seen the development of a cottage industry involving generalizations, deformations, and refinements of that theory.
Stanford and Witten [26] have explored many such extensions of JT gravity from the matrix model perspective, including its definition on unoriented surfaces.Meanwhile, Johnson and others have developed a complementary way of defining JT gravity and its generalizations via a certain combination of simpler matrix models.It is in the spirit of the latter that this work has proceeded.
We have shown that the modern perspective on double scaled matrix models, e.g. a focus on the eigenvalue density, is compatible with the results obtained in the 90s for β = 1 models, by unifying the more popular approach in the physics community at the timeorthogonal polynomials -with the more popular approach in the integrable systems at the time -skew orthogonal polynomials.
While it is possible to formally define the objects necessary to compute the eigenvalue density, we found that it is difficult to extract results for specific models.We believe the primary source of difficulty is the fact that the eigenvalue density is non-local in energy space via the presence of the energy integrals.Their mere presence is enough to render the WKB approach to perturbation theory intractable.Despite appealing to the (usually) powerful numerical approach, we found that many efforts to reduce error in one place, e.g. the fidelity of the wavefunctions, resulted in increased error elsewhere, e.g. the energy integrals in ρ.Further, the numerics displayed above had ℏ = 1, which is not particularly small.There are two logistical problems associated to decreasing ℏ: first, the string equation for u is actually harder to solve; second, the wavefunctions oscillate more rapidly, making the diagonalization problem in solving the Schrodinger equation more resource-intensive.
Despite the numerical difficulties, the leading order contribution to the density in the k = 2 model has the desired behavior for small E -albeit with some fictitious nonperturbative oscillation.The crosscap contribution to the density behaves in accordance with what one would expect of a pertubative correction with a relatively large value of ℏ.
As we have pointed out, analysis of the double scaled GOE yields a puzzle: it is missing unoriented contributions.Although it is not a genuine theory of surfaces, one can show that before double scaling there are non-vanishing unoriented contributions to the eigenvalue density when the theory is coupled to quarks 12 .It would be interesting to explore this further.
Given that the minimal model construction in β = 2 theories is so intimately connected to the KdV hierarchy, there must be a similar story behind the scenes of the β = 1 construction above.Two promising avenues to pursue are Drinfeld-Sokolov systems [38], due to the fundamental dependence of the double scaled β = 1 theory on factorizing a Lax operator, and the Pfaff lattice of Adler, et al [45], due to its fundamental connection with the β = 1 matrix integral at finite N .
Finally, one should expect that a non-perturbative treatment of the β = 1 Wigner-Dyson models is possible.The β = 2 Wigner-Dyson models find as a possible nonperturbative completion the Altand-Zirnbauer (AZ) (α, β) = (1 + 2Γ, 2) models, where the parameter Γ can be interpreted as the number of background branes in the theory [20][21][22]41] .The logical prediction is that the AZ (α, β) = (1 + γ, 1) models, for constant γ, will provide a non-perturbative completion of the models considered here.Demonstrating the connection between the β = 2 and (α, β) = (1 + 2Γ, 2) models requires the use of the KdV integrable structure of the theory, so it is possible that a non-perturbative completion of the β = 1 models will require a more careful treatment of their integrable structure as well.It would be interesting if such a non-perturbative completion of the unoriented theory could be used to study O-planes in two dimensions.

B.3 Factorization and New String Equations
Consider the factorization of P k in eq.(4.14), displayed again here for convenience The first condition on the functions g j will be found by comparing the coefficients of ∂ 2k−3 on both sides.In the preceding subsection we determined this coefficient in P for an arbitrary superposition.
There are two types of term that show up in the coefficient of ∂ 2k−3 on the right hand side of (B.10).The first comes from commuting all derivatives past functions.Since none of the functions will be differentiated, the mixed terms involving g i g j with i ̸ = j will cancel by symmetry.This leaves The other type of term comes from letting one of the derivatives in the factorization act on one of the g i .The possible outcomes can be visualized by breaking the factorization up into the operators T † and S as in eq.

Figure 1 :
Figure 1: Our result from (4.13) exactly overlaps with the standard result (3.38), as can be seen from (c).Two plots are indistinguishable.

Figure 2 : 3 2
Figure 2: ρ 0 vs E curve.We can see that ρ 0 perfectly wraps around E 3 2 in this energy window.

Figure 3 :
Figure 3: ρ 0 vs E curve.We can see the deviation at higher energies.

Figure 4 :
Figure 4: For lower energies all the curves with different cutoffs merge together.The curve deviates earlier as the cutoff is decreased.If we zoom in, we can see that for a particular energy the curve with a lower cutoff is more off from the E 3 2 curve.

1 Figure 5 :
Figure 5: The first unoriented contribution to the eigenvalue density is small compared to the disk level contribution, even with ℏ = 1.The combined contribution to the total density still has fictitious oscillations.

Figure 6 :
Figure 6: ρ 0 vs E curve.We can see that ρ 0 perfectly wraps around the expected result in this energy window.

1 Figure 7 :
Figure 7: The first unoriented contribution to the eigenvalue density in the interpolated matrix model.
(4.15), corresponding to lumping the terms of the form (∂ − g) into T † and the ones of the form ∂ (∂ + g) into S.The function g j can be differentiated once k − j times in T † and j + k times in S, since the terms in S can also be differentiated by the terms in T † .Taking into account the minus signs and factors of 2j = (−1)k−1 it k ℏ 2k−1 ∂ 2k−1 + (−1) k−1 it k 2k−3 ∂ 2k−3 + • • • .(B.13)