††institutetext: Department of Physics and Astronomy, University of Southern
California,
Los Angeles, CA 90089, U.S.A.

# Perturbative Unorientable JT Gravity and Matrix Models

Wasif Ahmed and  Ashton Lowenstein<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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.

##  1 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 Brezin:1977sv ; BREZIN1990144 ; DOUGLAS1990635
; PhysRevLett.64.127 ; GROSS1990333 ; HOOFT1974461 ; KAZAKOV1985295 ;
DAVID198545 ; BESSIS1980109 ; Kazakov:1989bc . 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\to\infty$ 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 Wigner .
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 Johnson:2022wsr ;
Johnson:2022hrj . 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 JACKIW1985343 ;
TEITELBOIM198341 ; Almheiri:2014cka . 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 Almheiri:2014cka
. In Saad:2019lba  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
Mertens_2021 . It was subsequently shown in Johnson:2020heh  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 Johnson:2019eik ;
Johnson:2020exp ; Johnson:2021tnl .

Generalizations of JT gravity (including orientability, defects,
supersymmetry, etc.) have been studied from both a field theory
Chamseddine:1991fg ; Forste:2017kwy ; Fan:2021wsb  and random matrix
perspective Stanford:2019vob ; Johnson:2021owr ; Johnson:2020heh . 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 forrester2010log ), the double scaling limits of
the more general Wigner-Dyson $\beta$-ensembles is not as investigated. Other
descriptions of unoriented gravity in two dimensions can be found in
JaumeGomis_2004 ; OrenBergman_2004 ; HARRIS1991685 . Since these models
describe unoriented minimal strings any effort in this direction provides nice
contact between math and physics.

In Stanford:2019vob  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 Saad:2019lba  they conclude that
the presence of even a single cross-cap can cause divergences in otherwise
well-behaved objects. See Stanford:2023dtm  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 $\beta=1$ ensembles with
arbitrary polynomial matrix potential. Results for the Gaussian case are known
already (see e.g. mehta2004random ; forrester2010log ) and progress has been
made in the direction of the string equations for more general models
Brezin:1990dk ; Carlisle_2005 .

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 Mahoux  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 Brezin:1990dk  the authors showed
how to double scale the derivative with respect to the matrix eigenvalue
$\lambda$

$\frac{d}{d\lambda}\xrightarrow{\text{double scale}}P_{k}.$ (1)

The operator $P_{k}$ is the Lax pair of a Sturm-Liouville operator
$\mathcal{H}=-\partial^{2}+u$, familiar in the KP/KdV hierarchy; in particular
it generates the $k^{\text{th}}$ KdV flow of the operator $\mathcal{H}$. By
imposing an auxiliary condition called the string equation one finds that
$P_{k}=\left(\mathcal{H}^{k-\frac{1}{2}}\right)_{+}$, 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 Adler_2000  the authors
derived an expression for the finite-$N$ eigenvalue density of a $\beta=1$
matrix model in terms of these polynomials. By adapting results from
Brezin:1990dk  we present the double scaling limit of the density (eqn. (57)).
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 Drinfeld:1984qv

$P_{k}=-\prod_{i=1}^{k}\left(\hbar\partial-\frac{1}{2}g_{i}\right)\hbar\partial\prod_{j=1}^{k}\left(\hbar\partial+\frac{1}{2}g_{k-j}\right).$
(2)

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 $\left(\mathcal{H}^{k-\frac{1}{2}}\right)_{+}$.

In analogy to the $\beta=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. (88)), 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 $\beta=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. In
section 6 we apply our results to obtain a formaluation of unoriented JT
gravity as an interpolation between an infinite number of minimal models. A
formal expression for the eigenvalue density is provided.

##  2 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 ginsparg1993lectures , and for more mathematical reviews
mehta2004random ; forrester2010log .

Given a set of $N\times 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

$Z=\int dM\,e^{-N\tr V(M)},$ (3)

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 $M\to
U^{\dagger}\operatorname{diag}(\lambda_{0},\cdots,\lambda_{N-1})U$, with $U$
in the appropriate symmetry group, we can rewrite the partition function after
integrating out $U$

$Z=C\int\prod_{i=0}^{N-1}d\lambda_{i}\,|\Delta(\lambda)|^{\beta}e^{-N\sum_{i=0}^{N-1}V(\lambda_{i})},$
(4)

where $\Delta(\lambda)=\prod_{i<j}(\lambda_{j}-\lambda_{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 $\beta$ is determined
by the symmetry group used to diagonalize $M$: for the group $U(N)$ $\beta=2$
and for $O(N)$ $\beta=1$.

The domain of integration for the eigenvalues is also affected by the value of
$\beta$. 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 $\Delta$. However, for even values of $\beta$ we can extend the
integration to all of $\mathbb{R}^{N}$ by dividing by a combinatorial factor
to take into account relabelling the indices anytime $\lambda_{j}-\lambda_{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 $\beta$ though, and must use alternate methods.

It can be shown that $\Delta$ is the determinant of a matrix consisting of
powers of the eigenvalues, $\Delta(\lambda)=\det\lambda_{i}^{j}$. By taking
linear combinations of the rows with real coefficients, we can introduce a
family of polynomials $p_{i}(x)$ so that $\Delta(\lambda)=\det
p_{j}(\lambda_{i})$. The polynomials will be normalized so that
$p_{i}(\lambda)=\lambda^{i}+\cdots.$ We stress here that the introduction of
these polynomials is possible for any value of $\beta$, 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

$\int_{-\infty}^{\infty}d\lambda
e^{-V(\lambda)}p_{i}(\lambda)p_{j}(\lambda)=h_{i}\delta_{ij}.$ (5)

This is particularly helpful when the values of $\lambda$ can be expanded to
all of $\mathbb{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

$\psi_{i}(\lambda)=h_{i}^{-1/2}e^{-V(\lambda)/2}\,p_{i}(\lambda).$ (6)

The choice of name has a two-fold motivation. In the Gaussian matrix model,
the orthogonal polynomials are the Hermite polynomials and the functions
$\psi_{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 $\psi_{i}$ is the same as the index on
$\lambda_{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
$\psi_{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^{\text{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 level111Note that since we begin
counting at $i=0$, the $N$th level is at $i=N-1$., and all lower values the
Fermi sea.

The oscillator wavefunctions provide an orthonormal basis for the Hilbert
space $L^{2}(\mathbb{R})$ with the standard inner product. Denote the basis as
$\\{\ket{\psi_{i}}\\}_{i=0}^{\infty}$. Operators on $L^{2}$ can thus be
decomposed in terms of these functions in the usual way

$\mathcal{O}_{ij}\equiv\matrixelement{\psi_{j}}{\mathcal{O}}{\psi_{i}}=\int_{-\infty}^{\infty}dx\,\psi_{j}(x)\hat{\mathcal{O}}\psi_{i}(x).$
(7)

By $\hat{\mathcal{O}}$ we denote the action of the operator on the actual
function $\psi_{i}(x)$.

The polynomials $p_{i}$ satisfy the following recursion relation

$\lambda p_{i}(\lambda)=p_{i+1}(\lambda)+R_{i}p_{i-1}(\lambda),$ (8)

with $R_{i}=h_{i}/h_{i-1}$. The corresponding recursion relation for the
oscillator wavefunctions is

$\lambda\psi_{i}(\lambda)=\sqrt{R_{i+1}}\psi_{i+1}(\lambda)+\sqrt{R_{i}}\psi_{i-1}(\lambda).$
(9)

The self-adjoint operator $\Lambda$ that implements this transformation on
$L^{2}$ has matrix elements

$\Lambda_{ij}=\sqrt{R_{i+1}}\delta_{j,i+1}+\sqrt{R_{i}}\delta_{j,i-1}.$ (10)

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|\leq 2k-1$, are determined by the
recursion coefficients $R_{i}$, and are anti-symmetric because $L$ is anti-
self-adjoint. The operators $\Lambda$ and $L$ satisfy the canonical
commutation relation $[\Lambda,L]\propto 1$.

Since $V$ is an even function it can be shown that each polynomial $p_{i}$ has
definite parity, given by $p_{i}(-\lambda)=(-1)^{i}p_{i}(\lambda)$. This
induces the $\mathbb{Z}_{2}$-grading
$L^{2}=\text{span}\\{\ket{\psi_{2\alpha}}\\}\oplus\text{span}\\{\ket{\psi_{2\beta+1}}\\}$.
Thus we can conveniently decompose our operators into $2\times 2$ block form.
For example, both $\Lambda$ and $L$ are odd under the grading, so we designate

$\displaystyle\Lambda$ $\displaystyle=\left(\begin{matrix}0&\ell\\\
\ell&0\end{matrix}\right),$ (11) $\displaystyle L$
$\displaystyle=\left(\begin{matrix}0&-c^{\dagger}\\\ c&0\end{matrix}\right).$

Objects that are labelled $0,\dots,N-1$ will carry a latin index $i,j$ (but
note that $k$ is reserved for the order of the matrix potential). When
considering $\mathbb{Z}_{2}$-grading we will decompose some of these objects
into even and odd parts, and it will be convenient to write $i=2\alpha$ or
$j=2\gamma+1$, respectively. Objects that have definite $\mathbb{Z}_{2}$
parity will thus carry a Greek index $\alpha,\gamma$ (but $\beta$ 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_{\gamma\alpha}=-L_{2\alpha,2\gamma+1}$. Since the matrix
representing $L$ has $2k-1$ non-zero off-diagonals, 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 $\psi_{2\alpha}$ can be written

$\psi_{2\alpha}^{\prime}=-\sum_{\gamma=0}^{\frac{N}{2}+k-1}c_{\gamma\alpha}\psi_{2\gamma+1}.$
(12)

For this range of $\gamma$ in the sum, some of the matrix elements
$c_{\alpha\gamma}$ will be zero depending on the value of $\alpha$, but these
bounds of summation are guaranteed to produce all necessary non-zero terms.

The double scaling limit BREZIN1990144 ; DOUGLAS1990635 ; GROSS1990333 ;
KAZAKOV1985295  is a combination of taking the size of the matrix $N\to\infty$
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
$\hbar$. The index $i$ is replaced by a continuous variable
$x\in(-\infty,\mu]$, where $\mu$ is the Fermi level. The large-$N$ limit
confines the spectrum to a single interval222This is the one-cut case., and in
the double scaling limit we zoom in on one of the edges of the inetrval,
replacing $\lambda$ with a new variabel $E\in[0,\infty)$. The oscillator
wavefunctions become functions of $\psi(x,E)$, the recursion coefficients
$R_{i}$ are replaced by a function $u(x)$, and the operator $\Lambda$ that
implements the recursion becomes a Schrodinger Hamiltonian
$\mathcal{H}=-\hbar^{2}\partial_{x}^{2}+u(x)$. All together, the recursion
relation becomes the Schrodinger equation

$\mathcal{H}\psi(x,E)=E\psi(x,E).$ (13)

The derivative operator $L=\frac{d}{d\lambda}$ is replaced by a differential
operator with respect to $x$, called $P_{k}$. The unique differential operator
satisfying $[P_{k},\mathcal{H}]\propto 1$ is given by
$P_{k}\propto\left(\mathcal{H}^{k-\frac{1}{2}}\right)_{+}$, where the $+$
subscript denotes keeping only non-negative powers of the derivative
$\partial_{x}$.

##  3 The $\beta=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
$\beta=1$ $O(N)$ symmetry class. The matrix $H$ is in the bi-fundamental
representation $\mathbf{N}\otimes\mathbf{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

$N\tr(H^{2})=\frac{N}{2}\sum_{i,j,m,n=0}^{N-1}H_{im}\Big{(}\delta_{in}\delta_{mj}+\delta_{ij}\delta_{mn}\Big{)}H_{jn}.$
(14)

This implies that the propagator is
$\frac{2}{N}(\delta_{in}\delta_{mj}+\delta_{ij}\delta_{mn})$ for $\beta=1$
ensembles. The second term comes from the fact that the indices are
indistinguishable for our matrices: $H_{ij}=H_{ji}$.

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 surfaces333This is, of
course, up to equivalences between even numbers of crosscaps and handles.. 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^{\tilde{F}}$ can be
expanded as

$\tilde{F}=\sum_{g=0}^{\infty}\sum_{c=0}^{\infty}N^{2-2g-c}\tilde{F}_{g,c}.$
(15)

The number $\chi(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 quantity444We remove the tilde to
distinguish between the two scaling regimes. $F_{g,c}$ represents the
contribution to the free energy of that topology.

Since $\beta=1$ is odd, we must take care to specify a domain of integration
in $Z$ that allows us to remove the absolute value on $\Delta$. We write

$Z=\int_{\mathcal{I}}\prod_{i=0}^{N-1}d\lambda_{i}\Delta(\lambda)e^{-N\sum_{j=0}^{N-1}V(\lambda_{j})},$
(16)

where $\mathcal{I}$ denotes the ordering
$-\infty<\lambda_{N-1}\leq\lambda_{N-2}\leq\cdots\leq\lambda_{1}\leq\lambda_{0}<\infty$.
General $n$-point functions of the eigenvalues can be computed using an idea
pioneered by Dyson doi:10.1063/1.1703773  and applied to the GOE by Mehta and
Mahoux Mahoux , 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

$G_{N}(u,v)=\left\langle{\prod_{\alpha=0}^{\frac{N}{2}-1}u(\lambda_{2\alpha})\prod_{\beta=0}^{\frac{N}{2}-1}v(\lambda_{2\beta+1})}\right\rangle.$
(17)

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.

###  3.1 Orthogonal Polynomial Approach

First, introduce a set of orthogonal polynomials $p_{i}(\lambda)$ satisfying

$\int_{-\infty}^{\infty}d\lambda\,e^{-2NV(\lambda)}p_{i}(\lambda)p_{j}(\lambda)=h_{i}\delta_{ij}.$
(18)

Notice that the integration measure includes a different power of $e^{-NV}$
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

$\psi_{i}(\lambda)=\frac{1}{\sqrt{h_{i}}}p_{i}(\lambda)e^{-NV(\lambda)}.$ (19)

Introduce the functions

$F_{i}(v;\lambda_{2\alpha})=\int_{-\infty}^{\lambda_{2\alpha}}d\lambda_{2\alpha+1}v(\lambda_{2\alpha+1})\psi_{i}(\lambda_{2\alpha+1}).$
(20)

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

$\displaystyle
G_{N}(u,v)=C_{N}\int_{-\infty}^{\infty}d\lambda_{0}e^{-NV(\lambda_{0})}u(\lambda_{0})\int_{-\infty}^{\lambda_{0}}d\lambda_{2}e^{-NV(\lambda_{2})}u(\lambda_{2})\cdots$
(21)
$\displaystyle\cdots\int_{-\infty}^{\lambda_{N-4}}d\lambda_{N-2}e^{-NV(\lambda_{N-2})}u(\lambda_{N-2})\det(m(v;\lambda_{2\alpha})),$

where $C_{N}$ is a numerical factor depending on the normalizations $h_{i}$
and the matrix $m$ is given by

$\displaystyle m_{i,2\alpha}$ $\displaystyle=p_{i}(\lambda_{2\alpha}),$ (22)
$\displaystyle m_{i,2\alpha+1}$ $\displaystyle=F_{i}(v;\lambda_{2\alpha}).$

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 $(-\infty,\infty)$
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

$G_{N}(u,v)=C^{\prime}_{N}\sum_{\sigma\in
S_{N}}\prod_{\alpha=0}^{\frac{N}{2}-1}J_{\sigma(2\alpha),\sigma(2\alpha+1)}(u,v),$
(23)

where

$J_{ij}(u,v)=\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy\,u(x)v(y)\psi_{[i}(x)\psi_{j]}(y),$
(24)

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\times 2$ block form. Begin by setting
$J_{2\alpha,2\gamma}(u,v)=f_{\alpha\gamma}(u,v)$. Next, define

$\displaystyle g_{\alpha\gamma}(u,v)$
$\displaystyle=-\sum_{\sigma=0}^{\frac{N}{2}+k-1}c_{\sigma\gamma}J_{2\alpha,2\sigma+1}(u,v)$
(25)
$\displaystyle=\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy\,u(x)v(y)\Big{(}\psi_{2\alpha}(x)\psi^{\prime}_{2\gamma}(y)-\psi_{2\alpha}(y)\psi_{2\gamma}^{\prime}(x)\Big{)}.$

Notice that the values of the dummy index $\sigma$ go outside of the “$N\times
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

$\displaystyle\mu_{\alpha\gamma}(u,v)$
$\displaystyle=\sum_{\mu,\nu=0}^{\frac{N}{2}+k-1}c_{\mu\alpha}J_{2\mu+1,2\nu+1}c_{\nu\gamma}$
(26)
$\displaystyle=\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy\,u(x)v(y)\Big{(}\psi^{\prime}_{2\alpha}(x)\psi^{\prime}_{2\gamma}(y)$
$\displaystyle\quad-\psi^{\prime}_{2\alpha}(y)\psi^{\prime}_{2\gamma}(x)\Big{)}.$

After some rearranging to ensure antisymmetry, we find the Pfaffian can be
written

$\text{Pf}(J)\to\text{det}^{1/2}\left(\begin{matrix}f_{\alpha\gamma}&g_{\alpha\gamma}\\\
-g_{\gamma\alpha}&\mu_{\alpha\gamma}\end{matrix}\right),$ (27)

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

$\rho_{n}(x_{0},\dots,x_{n-1})=\frac{N!}{(N-n)!}\int_{-\infty}^{\infty}d\lambda_{n}\int_{-\infty}^{\infty}d\lambda_{N-1}\,\prod_{i<j}|\lambda_{i}-\lambda_{j}|e^{-N\sum_{k=0}^{N-1}V(\lambda_{k})}.$
(28)

These correlators are related to $G_{N}$ via mehta2004random

$G_{N}(1+a,1+a)=\sum_{n=1}^{N}\frac{1}{n!}\int_{-\infty}^{\infty}\prod_{i=0}^{n-1}a(x_{i})dx_{i}\rho_{n}(x_{0},\dots,x_{n-1})$
(29)

In particular, this means that to compute the eigenvalue density
$\rho_{1}(\lambda)\equiv\tilde{\rho}(\lambda)$, we take one functional
derivative of $G_{N}$ and set the source to $0$. Mehta shows that for the GOE

$\tilde{\rho}(\lambda)=\sum_{\alpha=0}^{\frac{N}{2}-1}\left[\psi_{2\alpha}(\lambda)^{2}-\psi^{\prime}_{2\alpha}(\lambda)\int_{0}^{\lambda}d\lambda^{\prime}\psi_{2\alpha}(\lambda^{\prime})\right].$
(30)

The result (30) nearly holds for $k>1$. Denoting the density for the $k$th
multicritical model by $\tilde{\rho}^{(k)}$, the result is really

$\tilde{\rho}^{(k)}=\tilde{\rho}^{(1)}+\text{ surface terms},$ (31)

where the oscillator functions used in $\tilde{\rho}$ are $k$-dependent as
well. The additional terms, called surface terms here, are localized around
$\alpha=N/2$ and arise from the fact that we included terms outside of the
Fermi sea to obtain $\psi_{2\alpha}^{\prime}$ when rewriting the matrix $J$.
The number of extra terms is set by $k$ in the upper limit of summation in
(12). The exact form of these extra terms is straightforward, but reveals
little useful information even in the large-$N$ limit.

###  3.2 Skew Orthogonal Polynomial Approach

While the orthogonal polynomial approach bears some fruit when applied to the
$\beta=1$ models, there is another set of functions more naturally suited to
solving the problem. We begin again with (16) and compute the correlator
$G_{N}(u,u)$, setting $u=v$ for convenience. We still arrive at the expression

$\displaystyle G_{N}(u,v)=C_{N}$
$\displaystyle\int_{-\infty}^{\infty}d\lambda_{0}e^{-NV(\lambda_{0})}u(\lambda_{0})\int_{-\infty}^{\lambda_{0}}d\lambda_{2}e^{-NV(\lambda_{2})}u(\lambda_{2})\cdots$
(32)
$\displaystyle\cdots\int_{-\infty}^{\lambda_{N-4}}d\lambda_{N-2}e^{-NV(\lambda_{N-2})}u(\lambda_{N-2})\det(m(v;\lambda_{2\alpha})),$

where now we define

$\displaystyle m_{i,2\alpha}$ $\displaystyle=q_{i}(\lambda_{2\alpha}),$ (33)
$\displaystyle m_{i,2\alpha+1}$
$\displaystyle=\tilde{F}_{i}(u;\lambda_{2\alpha})=\int_{-\infty}^{\lambda_{2\alpha}}d\lambda^{\prime}e^{-NV(\lambda)}q_{i}(\lambda^{\prime}).$

The polynomials $q_{i}(\lambda)$ are monic and chosen to obey the following
relation Brezin:1990dk ; Adler_2000 ; Mahoux

$\langle
q_{2\alpha},q_{2\gamma+1}\rangle_{s}\equiv\iint_{-\infty}^{\infty}d\lambda
d\lambda^{\prime}\,e^{-NV(\lambda)-NV(\lambda^{\prime})}\varepsilon(\lambda-\lambda^{\prime})q_{i}(\lambda)q_{j}(\lambda^{\prime})=-r_{\alpha}\delta_{\alpha\gamma},$
(34)

where
$\varepsilon(\lambda-\lambda^{\prime})=\text{sgn}(\lambda-\lambda^{\prime})$,
and with all other $\langle q_{i},q_{j}\rangle_{s}=0$. The bilinear form
$\langle\cdot,\cdot\rangle_{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 $\tilde{J}$ given by

$\tilde{J}_{ij}=\langle uq_{i},uq_{j}\rangle_{s}.$ (35)

If we set $u=1$, this computes the full matrix integral.

In analogy to the oscillator wavefunctions $\psi_{i}$ defined above, introduce
the skew oscillator wavefunctions

$\zeta_{i}(\lambda)=r_{i}^{-1/2}e^{-NV(\lambda)}q_{i}(\lambda).$ (36)

Notice that the matrix $\tilde{J}$ has the $2\times 2$ block structure

$\tilde{J}[u]=\left(\begin{matrix}f_{\alpha\gamma}[u]&g_{\alpha\gamma}[u]\\\
-g_{\gamma\alpha}[u]&\mu_{\alpha\gamma}[u]\end{matrix}\right),$ (37)

where

$\displaystyle f_{\alpha\gamma}[u]$
$\displaystyle\equiv\tilde{J}_{2\alpha,2\gamma}[u],$ (38) $\displaystyle
g_{\alpha\gamma}[u]$ $\displaystyle\equiv\tilde{J}_{2\alpha,2\gamma+1}[u],$
$\displaystyle\mu_{\alpha\gamma}[u]$
$\displaystyle\equiv\tilde{J}_{2\alpha+1,2\gamma+1}[u].$

The functions $f,g$ and $\mu$ 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\times 2$
structure of $\tilde{J}$ can be written (schematically) as

$\tilde{J}=\left(\begin{matrix}\epsilon f&1+\epsilon\nu\\\
-(1+\epsilon\nu)&\epsilon\mu\end{matrix}\right),$ (39)

with $\epsilon$ small and related to the function $a$. The pfaffian
pf$\tilde{J}$ admits a straightforward expansion in $\epsilon$ (see Appendix
A.7 in mehta2004random ). The density can once again be extracted from
$G_{N}(1+a)$ by taking a functional derivative with respect to $a$. The result
is Adler_2000 ; Mahoux

$\tilde{\rho}(\lambda)=\sum_{\alpha=0}^{\frac{N}{2}-1}\Big{(}\zeta_{2\alpha}(\lambda)[\hat{\varepsilon}\cdot\zeta_{2\alpha+1}](\lambda)-\zeta_{2\alpha+1}(\lambda)[\hat{\varepsilon}\cdot\zeta_{2\alpha}](\lambda)\Big{)}.$
(40)

This expression for the eigenvalue density is exact for any value of $N$. This
is in contrast to (30), 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 Adler_2000 . 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

$q_{i}=\sum_{j=0}^{i}O_{ij}p_{j}.$ (41)

Representing the skew oscillator wavefunctions with the kets
$\ket{\zeta_{i}}$, the defining relation for the polynomials $q_{i}$ is then
given by

$\langle
q_{i},q_{j}\rangle_{s}=2\sqrt{r_{i}r_{j}}\matrixelement{\zeta_{i}}{\hat{\varepsilon}}{\zeta_{j}},$

where $\hat{\varepsilon}$ is the integral operator with kernel
$\frac{1}{2}\varepsilon(\lambda-\lambda^{\prime})$ 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^{-2NV}$. By
taking into account the normalizations, the skew oscillator wavefunctions can
be expressed in terms of the oscillator wavefunctions in a way related ro (41)

$\zeta_{i}=\sum_{j=0}^{i}\tilde{O}_{ij}\psi_{j},\quad\tilde{O}_{ij}=\sqrt{h_{j}/r_{i}}O_{ij}.$
(42)

The skew oscillator wavefunctions have the same $\mathbb{Z}_{2}$-grading as
the oscillator wavefunctions, so the matrix $\tilde{O}$ must be even under the
$\mathbb{Z}_{2}$-grading. Define the matrices $a,b$ so that

$\tilde{O}=\left(\begin{matrix}a&0\\\ 0&b\end{matrix}\right).$

In Brezin:1990dk  they show that $a,b$ are related to the derivative matrix
$c$ by

$c=b^{T}a.$ (43)

Hence the eigenvalue density is written (schematically) as

$\tilde{\rho}\sim(a\cdot\psi_{\text{even}})\cdot(\hat{\varepsilon}\cdot
b\cdot\psi_{\text{odd}})-(b\cdot\psi_{\text{odd}})\cdot(\hat{\varepsilon}\cdot
b\cdot\psi_{\text{even}}).$ (44)

To conclude, it should be noted that $\tilde{\rho}$ is the diagonal part of a
larger function Adler_2000 . Define

$s(\lambda,\lambda^{\prime})=\sum_{\alpha=0}^{\frac{N}{2}-1}\Big{(}\zeta_{2\alpha}(\lambda)[\hat{\varepsilon}\cdot\zeta_{2\alpha+1}](\lambda^{\prime})-\zeta_{2\alpha+1}(\lambda^{\prime})[\hat{\varepsilon}\cdot\zeta_{2\alpha}](\lambda)\Big{)}.$
(45)

The analog of the self-reproducing kernel familiar from studies of $\beta=2$
theories is given by

$f_{1}(\lambda,\lambda^{\prime})=\left(\begin{matrix}s(\lambda,\lambda^{\prime})&I(\lambda,\lambda^{\prime})\\\
\partial_{\lambda}s(\lambda,\lambda^{\prime})&s(\lambda,\lambda^{\prime})\end{matrix}\right),$
(46)

where

$I(\lambda,\lambda^{\prime})=\frac{1}{2}\int_{-\infty}^{\infty}s(\lambda,z)\varepsilon(z-\lambda^{\prime})dz-\frac{1}{2}\varepsilon(\lambda-\lambda^{\prime}).$
(47)

The eigenvalue density is given by

$\tilde{\rho}(\lambda)=\lim_{\lambda^{\prime}\to\lambda}\text{qdet}f_{1}(\lambda,\lambda^{\prime}),$
(48)

where qdet denotes the quaternionic determinant (see Mahoux , for example).
The result follows from the fact that $I(\lambda,\lambda)=0$.

###  3.3 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
$\beta=1,2,4$ theories, can be determined by using known facts about the
Hermite polynomials. See Forste:2017kwy  for a review. In particular, the
double scaled eigenvalue density for the GOE can be expressed in terms of the
Airy function

$\rho_{GOE}(E)=\rho^{(2)}(E)+\frac{1}{2}\text{Ai}(-E)\left(1-\int_{-\infty}^{E}\text{Ai}(-E^{\prime})dE^{\prime}\right),$
(49)

where $\rho^{(2)}$ 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 Brezin:1990dk . In the normalization used here, they are

$\displaystyle a$ $\displaystyle\to
T_{k}=\hbar^{k-1}\partial^{k-1}+\sum_{i=0}^{k-2}\hbar^{i}g_{i}\partial^{i},$
(50) $\displaystyle b$ $\displaystyle\to
S_{k}=\hbar^{k}\partial^{k}+\sum_{i=0}^{k-1}\hbar^{i}h_{i}\partial^{i}.$

The double scaled version of (43) is $S_{k}^{\dagger}T_{k}=P_{k}$. The authors
of Brezin:1990dk  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$).

##  4 New Results in Matrix Models

###  4.1 Double Scaled Density

We begin by double scaling the modified version of the finite-$N$ GOE density
given in (31). First, in the strictly large-$N$ limit we change the index $i$
into a continuous variable $X\in[0,1]$. Then introduce a small parameter
$\delta$ that picks out the scaling parts

$\displaystyle X$ $\displaystyle\to\delta^{4}x,\quad\lambda\to\delta^{2}E,$
(51) $\displaystyle\psi_{i}(\lambda)$
$\displaystyle\to\psi(x,E),\quad\frac{d}{d\lambda}\to\delta^{-2}P_{k}.$

In particular in the large-$N$ limit, the eigenvalue distribution for these
models settles into a finite window $[-\lambda_{c},\lambda_{c}]$ which, due to
the symmetry of the theory, is symmetric about the origin. The full scaling
ansatz for the eigenvalue is $\lambda=\lambda_{c}-\delta^{2}E$; using this as
a coordinate transformation implies that integrals over $\lambda$ have the
scaling part

$\int d\lambda\to-\delta^{2}\int dE,$ (52)

where the bounds of integration are updated on a case-by-case basis.

By changing the large-$N$ integration over $\alpha/N$ to $i/N$ we introduce a
factor of $\frac{1}{2}$. The double scaled density is then

$\rho(E)=\frac{1}{2}\int_{-\infty}^{0}\left[\psi(x,E)^{2}-P_{k}\psi(x,E)\int_{E}^{\infty}\psi(x,E^{\prime})dE^{\prime}\right]dx+\mathcal{F}_{k}[\psi](0,E).$
(53)

The double scaling limit of the surface terms from (31) are denoted by
$\mathcal{F}_{k}[\psi]$. Since before double scaling there are
$\mathcal{O}(1)$ surface terms around $\alpha=N/2$, $\mathcal{F}_{k}[\psi]$
should involve derivatives of $\psi$ 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 $\beta=1$ density written in terms of the skew
oscillator wavefunctions (40) we need to discuss the double scaling limit of
the operator $\hat{\varepsilon}$. Consider the application of
$\hat{\varepsilon}$ to functions with definite parity. The results are

$\displaystyle\text{even}:$ $\displaystyle\quad\hat{\varepsilon}\cdot
f(\lambda)=\int_{0}^{\lambda}f(\lambda^{\prime})d\lambda^{\prime},$ (54)
$\displaystyle\text{odd}:$ $\displaystyle\quad\hat{\varepsilon}\cdot
f(\lambda)=\int_{\lambda}^{\lambda_{c}}f(\lambda^{\prime})d\lambda^{\prime},$

where we’ve assumed that $f$ only has support on $[-\lambda_{c},\lambda_{c}]$.
By the scaling ansatz for the $\lambda$ integral,

$\displaystyle\text{even}:$ $\displaystyle\quad\hat{\varepsilon}\cdot
f(\lambda)\xrightarrow{\text{double
scale}}\delta^{2}\int_{E}^{\infty}f(E^{\prime})dE^{\prime},$ (55)
$\displaystyle\text{odd}:$ $\displaystyle\quad\hat{\varepsilon}\cdot
f(\lambda)\xrightarrow{\text{double
scale}}-\delta^{2}\int_{0}^{E}f(E^{\prime})dE^{\prime}.$

The formula for the finite-$N$ density in (40) involves the skew
wavefunctions, which are related to the oscillator wavefunctions via the
transformations $a$ and $b$ as depicted in (44). 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 $\zeta$ and
$\psi$ is expressed as555The powers of $\delta$ appearing in these
relationships may not actually both be $-1$. What is important is that the sum
of the two powers is $-2$, as it is for the operator $P_{k}$.

$\displaystyle\text{even}:$
$\displaystyle\quad\zeta_{2\alpha}(\lambda)=[a\cdot\psi]_{2\alpha}(\lambda)\xrightarrow{\text{double
scale}}\delta^{-1}\hat{T}_{k}\psi(x,E)$ (56) $\displaystyle\text{odd}:$
$\displaystyle\quad\zeta_{2\alpha+1}(\lambda)=[b\cdot\psi]_{2\alpha+1}(\lambda)\xrightarrow{\text{double
scale}}\delta^{-1}\hat{S}_{k}\psi(x,E)$

Putting together these pieces, we find the expression for the double scaled
$\beta=1$ density

$\rho(E)=\frac{1}{2\hbar}\int_{-\infty}^{0}dx\Bigg{[}T_{k}\psi(x,E)\left(\int_{0}^{E}S_{k}\psi(x,E^{\prime})dE^{\prime}\right)+S_{k}\psi(x,E)\left(\int_{E}^{\infty}T_{k}\psi(x,E^{\prime})dE^{\prime}\right)\Bigg{]}.$
(57)

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^{\prime}$
integral, we would pick up a contribution at $x=0$ and get the term

$\psi(x,E)\int_{0}^{E}T^{\dagger}_{k}S_{k}\psi(x,E^{\prime})dE^{\prime},$ (58)

where the adjoint is taken in $L^{2}(\mathbb{R})$. Now, the specific
combination $T_{k}^{\dagger}S_{k}$ is proportional to $\hat{P}_{k}$, which is
canonically conjugate to the Hamiltonian $\mathcal{H}$ for which $\psi(x,E)$
is an eigenfunction. Thus $P_{k}$ has a representation
$P_{k}\sim\partial_{E}$, and

$\psi(x,E)\int_{0}^{E}T^{\dagger}_{k}S_{k}\psi(x,E^{\prime})dE^{\prime}\propto\psi(x,E)\Big{(}\psi(x,E)-\psi(x,0)\Big{)},$
(59)

which produces the $\psi^{2}$ term familiar from the study of $\beta=2$
theories.

Consider the second term in the density. Once again we integrate by parts to
shift the location of $\hat{T}_{k}$. This picks up more terms at $x=0$ and
gives

$T_{k}^{\dagger}S_{k}\psi(x,E)\int_{E}^{\infty}\psi(x,E^{\prime})dE^{\prime}\propto
P_{k}\psi(x,E)\int_{E}^{\infty}\psi(x,E^{\prime})dE^{\prime}.$ (60)

Being more careful about signs, we conclude

$\rho(E)=\frac{1}{2}\int_{-\infty}^{0}dx\,\Bigg{[}\psi(x,E)^{2}-P_{k}\psi(x,E)\int_{E}^{\infty}\psi(x,E^{\prime})dE^{\prime}\Bigg{]}+\text{
surface terms},$ (61)

where the surface terms are generated by moving the operator $\hat{T}_{k}$
around.

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

$\rho_{\text{GOE}}(E)=\frac{1}{2}\int_{-\infty}^{0}dx\,\Bigg{[}\psi(x,E)^{2}-P_{k}\psi(x,E)\int_{E}^{\infty}\psi(x,E^{\prime})dE^{\prime}\Bigg{]}.$
(62)

We confirmed numerically that (62) matches the previously known result.

###  4.2 New String Equations

It is evident from (57) that we have to determine the operators $T_{k},S_{k}$,
and there is a well defined algorithm to determine them outlined in
Brezin:1990dk . For the $k^{\text{th}}$ minimal model we can write the
following factorization equation

$P_{k}=-\prod_{i=1}^{k-1}\left(\hbar\partial-\frac{1}{2}g_{i}\right)\hbar\partial\prod_{j=1}^{k-1}\left(\hbar\partial+\frac{1}{2}g_{k-j}\right).$
(63)

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}^{\dagger}T_{k}$ with $T_{k},S_{k}$ previously
defined in eqn. (50). 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-self-adjoint and setting some
integration constants to 0 (see Brezin:1990dk ) one can reduce the number of
functions so that $T_{k},S_{k}$ depend on the same family of functions. The
new factorization in eqn. (63) amounts to setting

$\displaystyle S_{k}$
$\displaystyle=\prod_{j=1}^{k-1}\left(\hbar^{2}\partial^{2}+\frac{\hbar}{2}g_{k-j}\partial+\frac{\hbar}{2}g_{k-j}^{\prime}\right),$
(64) $\displaystyle T_{k}$
$\displaystyle=\prod_{j=1}^{k-1}\left(\hbar\partial+\frac{1}{2}g_{j}\right).$

Expanding both sides of (63) and equating powers of $\partial$ 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 $\partial^{2k-3}$ is always666This exact form is
technically only true for a single multicritical model. See the following
discussion of massive interpolations for more.Brezin:1990dk

$\left(k-\frac{1}{2}\right)u=\sum_{j=1}^{k-1}\left[-\hbar
jg_{j}^{\prime}+\frac{1}{4}g_{j}^{2}\right].$ (65)

Moreover the equation determined by $\partial^{2k-4}$ is simply the derivative
of the above equation. This implies that $g_{0}\sim\sqrt{u_{0}}$, even for
more complicated massive interpolations.

The preceding algorithm can be adapted to consider an interpolation up to the
$k^{\text{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=\sum_{i}t_{i}P_{i}$ as a differential operator is set by $k$.
Since (65) is determined by $\partial^{2k-3}$, the highest order term with a
nontrivial coefficient, the only two operators that can contribute to that
equation are $P_{2k-1}$ and $P_{2k-3}$. Therefore the proper way to modify
(65) is

$\left(k-\frac{1}{2}\right)t_{k}u+t_{k-1}=t_{k}\sum_{j=1}^{k-1}\left[-\hbar
jg_{j}^{\prime}+\frac{1}{4}g_{j}^{2}\right].$ (66)

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 $\partial^{2k-5}$, and therefore will explicitly include the
coupling constants $t_{k},t_{k-1}$, and $t_{k-2}$.

##  5 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.

###  5.1 The GOE

We can use the fact that the wavefunctions $\psi(x,E)$ and the eigenvalue
density $\rho$ 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 is777Strictly speaking this is
true for the GUE. The Schrodinger equation for the GOE is related to this by a
rescaling of the energy. See the appendix for discussion.

$\hbar^{2}\partial_{x}^{2}\psi+(x+E)\psi=0.$ (67)

Taking the standard WKB form of the wavefunction, we find the two
possibilities

$\psi(x,E)=\frac{1}{\sqrt{2\pi\hbar}}(x+E)^{-1/4}\left[1-\left(\frac{i5\hbar}{48(x+E)^{3/2}}+\frac{5\hbar^{2}}{64(x+E)^{3}}\right)\right]e^{\pm
i\frac{2}{3\hbar}(x+E)^{3/2}}+\cdots.$ (68)

Of course, the differential equation can be solved exactly in terms of the
Airy function, whose Taylor series expansion in small-$\hbar$ is a linear
combination of the two perturbative solutions obtained above

$\displaystyle\hbar^{-2/3}\text{Ai}\left(-\hbar^{-2/3}(x+E)\right)$
$\displaystyle\approx\frac{1}{\sqrt{\pi\hbar}}(x+E)^{-1/4}e^{-\frac{5\hbar^{2}}{64}(x+E)^{-3}}$
(69)
$\displaystyle\times\cos\left[\frac{2}{3\hbar}(x+E)^{3/2}-\frac{5\hbar}{48}(x+E)^{-3/2}-\frac{\pi}{4}\right].$

The WKB expansion is inherently an expansion in small $\hbar$. However, the
non-perturbative oscillating piece has to be dealth with since it does not
simplify under the assumption that $\hbar\ll 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 $\psi$ 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,

$\displaystyle T\propto 1,$ (70) $\displaystyle
S\propto\frac{\partial}{\partial x}.$

Using the GOE wavefunctions $\psi(x,2\sqrt{2}E)$ one finds

$\rho(E)=\frac{2^{3/4}\sqrt{E}}{\pi\hbar}+\frac{7\hbar}{2^{3/4}256\pi
E^{5/2}}+\cdots.$ (71)

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 $\mathcal{O}(\hbar^{1})$ term has a different coefficient than the
$\beta=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 $\mathcal{O}(\hbar^{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 $\mathcal{O}(\hbar^{0})$ integral, but it actually contributes at
$\mathcal{O}(\hbar^{1})$. As a corollary of this, there are no terms in the
perturbative series at $\mathcal{O}(\hbar^{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.

###  5.2 The $k=1$ to $k=2$ Interpolation

Take an interpolation between $k=1$ and $k=2$ with coupling constant $t$
(we’ll set $t_{2}=1$ and $t_{1}=t$). The momentum operator for this
interpolation is given by

$P=\hbar^{3}\partial^{3}+\left(\frac{3}{2}u+t\right)\hbar\partial+\frac{3}{4}\hbar
u^{\prime}.$ (72)

There is only one independent function $g$ in the factorization of $P$, so

$\displaystyle S_{2}$
$\displaystyle=\hbar^{2}\partial^{2}+\frac{1}{2}g\hbar\partial+\frac{1}{2}\hbar
g^{\prime},$ (73) $\displaystyle T_{2}$
$\displaystyle=\hbar\partial+\frac{1}{2}g.$

The independent equation of motion for $g$ is

$\frac{3}{2}u+t=-\hbar g^{\prime}+\frac{1}{4}g^{2}.$ (74)

Each function will have a perturbative expansion in $\hbar$888The subscripts
here are used to denote the order of the perturbative expansion of the single
function $g$, whereas before the subscript enumerated a family of different
functions in the factorization of $P$.

$\displaystyle u$ $\displaystyle=u_{0}+\hbar^{2}u_{2}+\hbar^{4}u_{4}+\cdots$
(75) $\displaystyle g$ $\displaystyle=g_{0}+\hbar
g_{1}+\hbar^{2}g_{2}+\cdots.$

The first orders in the expansion of $g$ are given by

$\displaystyle g_{0}$ $\displaystyle=\pm\sqrt{4t+6u_{0}},$ (76) $\displaystyle
g_{1}$ $\displaystyle=\frac{3u^{\prime}_{0}}{2t+3u_{0}},$ $\displaystyle
g_{2}$
$\displaystyle=\frac{3\left(8t^{2}u_{2}-8tu_{0}^{\prime\prime}+24tu_{0}u_{2}-12u_{0}u_{0}^{\prime\prime}+15u_{0}^{\prime
2}+18u_{0}^{2}u_{2}(x)\right)}{2\sqrt{2}(2t+3u_{0})^{5/2}}.$

It is easy to show that as $t\to 0$, the perturbative expansion smoothly
transitions to the single $k=2$ case.

###  5.3 Pure Unoriented Gravity

Consider just the $k=2$ minimal model. The leading contribution to the density
of states for the $k^{\text{th}}$ minimal model is proportional to
$E^{k-\frac{1}{2}}$, so we expect to see $E^{3/2}$ behavior. The density
admits a topological expansion in powers of $\hbar$,

$\rho(E)=\sum_{g,c=0}^{\infty}\hbar^{2g+c-1}\rho_{g,c}(E),$ (77)

where $g$ is the number of handles of the surface and $c$ is the number of
crosscaps. By convention we define $\rho_{0,0}\equiv\rho_{0}$.

By keeping the leading order contributions from $T_{2},S_{2}$ we find the disk
density

$\rho_{0}(E)=\int_{-\infty}^{0}dx\left[\frac{1}{4}\partial(g_{0}\psi)\int_{E}^{\infty}(g_{0}\psi)dE^{\prime}-\frac{1}{4}(g_{0}\psi)\int_{0}^{E}\partial(g_{0}\psi)dE^{\prime}\right].$
(78)

Isolating the first subleading contributions gives the crosscap density

$\displaystyle\rho_{0,1}(E)$
$\displaystyle=\int_{-\infty}^{0}dx\Bigg{[}\frac{1}{4}\partial(g_{1}\psi)\int_{E}^{\infty}(g_{0}\psi)dE^{\prime}+\frac{1}{4}\partial(g_{0}\psi)\int_{E}^{\infty}(g_{1}\psi)dE^{\prime}$
(79)
$\displaystyle-\frac{1}{4}(g_{1}\psi)\int_{0}^{E}\partial(g_{0}\psi)dE^{\prime}-\frac{1}{4}(g_{0}\psi)\int_{0}^{E}\partial(g_{1}\psi)dE^{\prime}$
$\displaystyle+\frac{1}{2}\partial(g_{0}\psi)\int_{E}^{\infty}(\partial\psi)dE^{\prime}+\frac{1}{2}\partial^{2}\psi\int_{E}^{\infty}(g_{0}\psi)dE^{\prime}$
$\displaystyle-\frac{1}{2}(\partial\psi)\int_{0}^{E}\partial(g_{0}\psi)dE^{\prime}-\frac{1}{2}(g_{0}\psi)\int_{0}^{E}(\partial^{2}\psi)dE^{\prime}\Bigg{]}.$

The functions $g_{0}$ and $g_{1}$ are obtained from (76) by setting $t=0$.

Both (78) and (79) contain more information than is desired, unless we compute
the wavefunctions perturbatively using WKB analysis. Such an analysis yields

$\psi(x,E)=\frac{1}{\sqrt{\pi\hbar}}\left(E-\sqrt{-x}\right)^{-1/4}\cos\left[\frac{4}{15\hbar}\sqrt{E-\sqrt{-x}}\left(2E^{2}+E\sqrt{-x}+3x\right)\right].$
(80)

Unfortunately it is difficult to maintain neither analytical nor numerical
control over the calculation of the density using the WKB wavefunctions. As a
last resort we appeal to the power of the non-perturbative numerical framework
described in Johnson:2021tnl .

The function $u$ satisfies the differential equation

$u\mathcal{R}^{2}-\frac{\hbar^{2}}{2}\mathcal{R}^{\prime\prime}+\frac{\hbar^{2}}{4}(\mathcal{R}^{\prime})^{2}=0,$
(81)

where $\mathcal{R}=-\frac{\hbar^{2}}{3}u^{\prime\prime}+u^{2}+x$. Equation
(81) can be solved numerically by putting the system in a box and imposing the
boundary conditions

$u(x\to-\infty)=\sqrt{-x},\quad u(x\to+\infty)=-\frac{1}{4x^{2}}.$ (82)

Then the Schrodinger equation for this potential is discretized and solved
using the methods of Numerov . In this work, a spatial grid size of $\Delta
x=.03$ is used. The presence of the walls needed to obtain numerical solutions
constrains the maximum trusted energy. The average spacing between energy
levels is $\Delta E\approx.02$.

Using the numerically determined wavefunctions produces the disk level density
in fig. (1).

Figure 1: $\rho_{0}$ vs E curve. We can see that $\rho_{0}$ perfectly wraps
around $E^{\frac{3}{2}}$ in this energy window.

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 Johnson:2022wsr . However, here they are
meaningless because the formalism is only meant to be trusted at the level of
perturbation theory.

If we plot $\rho_{0}$ up to higher energies we see a deviation of the curve
from $E^{\frac{3}{2}}$ as shown in figure (2).

Figure 2: $\rho_{0}$ vs E curve. We can see the deviation at higher energies.

The deviation is there because we have an integral from $E$ to $\infty$ in
(78), and this upper limit was replaced by the highest energy eigenvalue we
had in our numerical solution. In this plot, that cutoff was roughly
$E_{\text{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
(3).

(a) $\rho_{0}$ with different cutoffs.

(b) Zoomed in

Figure 3: 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^{\frac{3}{2}}$ curve.

We can similarly plot $\rho_{0,1}$ using equation (79). To compare the values
of $\rho_{0}$ and $\rho_{0,1}$, we plot them together in figure (4). We can
see that $\rho_{0}$ is more dominant, and $\rho_{1}$ is just like a damped
oscillation around the $E$ axis, as it should be for a perturbative
correction. The sum $\rho_{0}+\rho_{0,1}$ still oscillates with the
approximate $E^{3/2}$ behavior of the leading contribution.

(a) $\rho_{0}$ and $\rho_{0,1}$ separately

(b) $\rho_{0}+\rho_{0,1}$

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

The amplitude of $\rho_{0,1}$ is decreasing as $E\to\infty$. The topological
expansion of the density in (77) is equivalently an expansion in large $E$, so
at large energy we expect to recover the disk contribution.

##  6 Unoriented JT Gravity

As shown in Johnson:2020heh ; Johnson:2021owr ; Johnson:2021tnl ;
Johnson:2022wsr , JT gravity is perturbatively defined as a massive
interpolation between an infinite number of multicritical models. The full
equation of motion obeyed by $u$, called the string equation, is given by

$\mathcal{R}\equiv\sum_{k=1}^{\infty}t_{k}R_{k}[u]+x=0,$ (83)

where the couplings $t_{k}$ are determined by the leading density of states of
the gravity theory

$t_{k}=\frac{\pi^{2k-2}}{2k!(k-1)!},$ (84)

and $R_{k}$ is the $k$th Gel’fand-Dikii polynomial.

As is clear from the formalism developed above, the function $u$ is seeded
from the $\beta=2$ theory into the $\beta=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=\sum_{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

$\mathcal{P}=\sum_{k=1}^{\infty}\frac{\pi^{2k-2}}{2k!(k-1)!}P_{k}.$ (85)

We can, in principle, define operators $\mathcal{T},\mathcal{S}$ such that
$\mathcal{P}=\mathcal{S}^{\dagger}\mathcal{T}$ as before, and which map the
wavefunctions into the skew wavefunctions. Schematically, the factorization
looks like

$\mathcal{P}=-\prod_{i=1}^{\infty}\left(\hbar\partial-\frac{1}{2}g_{i}\right)\hbar\partial\prod_{j=\infty}^{1}\left(\hbar\partial+\frac{1}{2}g_{j}\right).$
(86)

The flipped bounds on the right-most product indicate that the functions
$g_{i}$ appear in reverse order in that part

$\mathcal{P}=-\left(\hbar\partial-\frac{1}{2}g_{1}\right)\left(\hbar\partial-\frac{1}{2}g_{2}\right)\cdots\hbar\partial\cdots\left(\hbar\partial+\frac{1}{2}g_{2}\right)\left(\hbar\partial+\frac{1}{2}g_{1}\right).$
(87)

In the same manner as before, we find enough differential equations by
comparing coefficients of $\partial$ to determine the $g_{i}$.

Given the extension of the operator $\mathcal{P}$ to the infinite
interpolation that defines JT gravity and its formal factorization in terms of
operators $\mathcal{T}$ and $\mathcal{S}$, we expect that (57) remains true
for unoriented JT gravity, giving the density

$\rho(E)=\frac{1}{2\hbar}\int_{-\infty}^{0}dx\Bigg{[}\mathcal{T}\psi(x,E)\left(\int_{0}^{E}\mathcal{S}\psi(x,E^{\prime})dE^{\prime}\right)-\mathcal{S}\psi(x,E)\left(\int_{E}^{\infty}\mathcal{T}\psi(x,E^{\prime})dE^{\prime}\right)\Bigg{]},$
(88)

to all orders in perturbation theory.

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 $\mathcal{P}$ is infinite order, and so too are
$\mathcal{T}$ and $\mathcal{S}$. The solution is to notice that $t_{k}$ gets
smaller and smaller for higher $k$ and $\lim_{k\rightarrow\infty}t_{k}=0$. In
Johnson:2020exp  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_{\text{max}}$ by finding the function $u$
satisfying999This is only true perturbatively. There are difficulties
associated to restricting oneself to perturbation theory, as discussed earlier
in the $k=2$ example. The truncation performed in Johnson:2020exp  was for the
non-perturbative formulation of JT gravity.

$\sum_{k=1}^{k_{\text{max}}}t_{k}R_{k}+x=0,$

where $t_{k}$ are the same coupling constants for JT gravity (84). 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 $\psi$. Then
determine the interpolated operator
$\mathcal{P}\approx\sum_{k=1}^{k_{\text{max}}}P_{k}$ in terms of $u$, and use
that to find the $k_{\text{max}}-1$ functions $g_{i}$. Then build the
truncated operators $\mathcal{T}$ and $\mathcal{S}$ out of the $g_{i}$, and
finally use the formula for the density (88).

##  7 Discussion

The purpose of this paper was to further explore the study of double scaled
$\beta=1$ Wigner-Dyson models initiated in papers like Brezin:1990dk ;
HARRIS1990384 . Recent advancements in two-dimensional 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 Stanford:2019vob  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 $\beta=1$ models, by unifying the more popular approach in the
physics community at the time – orthogonal polynomials – with the more popular
approach in the integrable systems at the time – skey 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 $\rho$. Further, the numerics displayed above had
$\hbar=1$, which is not particularly small. There are two logistical problems
associated to decreasing $\hbar$: 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 non-perturbative 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 $\hbar$.

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 quarks101010This is actually a problem in Zee’s textbook zee2010quantum ..
It would be interesting to explore this further.

Given that the minimal model construction in $\beta=2$ theories is so
intimately connected to the KdV hierarchy, there must be a similar story
behind the scenes of the $\beta=1$ construction above. Two promising avenues
to pursue are Drinfeld-Sokolov systems Drinfeld:1984qv , due to the
fundamental dependence of the double scaled $\beta=1$ theory on factorizing a
Lax operator, and the Pfaff lattice of Adler, et al adler1999pfaff , due to
its fundamental connection with the $\beta=1$ matrix integral at finite $N$.

Finally, one should expect that a non-perturbative treatment of the $\beta=1$
Wigner-Dyson models is possible. The $\beta=2$ Wigner-Dyson models find as a
possible non-perturbative completion the Altand-Zirnbauer (AZ)
$(\alpha,\beta)=(1+2\Gamma,2)$ models, where the parameter $\Gamma$ can be
interpreted as the number of background branes in the theory Dalley:1991vr ;
Johnson:2019eik ; Johnson:2020exp ; Johnson:2021tnl  . The logical prediction
is that the AZ $(\alpha,\beta)=(1+\gamma,1)$ models, for constant $\gamma$,
will provide a non-perturbative completion of the models considered here.
Demonstrating the connection between the $\beta=2$ and
$(\alpha,\beta)=(1+2\Gamma,2)$ models requires the use of the KdV integrable
structure of the theory, so it is possible that a non-perturbative completion
of the $\beta=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.

###### Acknowledgements.

We are grateful to Clifford V. Johnson for many helpful conversations and his
guidance in the preperation of this manuscript. Wasif Ahmed thanks Sushmit
Hossain for helping initially with MATLAB. This work is supported in part by
the US Department of Energy under grant DE-SC0011687.

##  Appendix A The GOE

###  A.1 Saddle point analysis

For any $\beta$ and matrix potential $V$, the matrix integral of the Dyson
$\beta$-ensemble is

$\tilde{Z}=\int\prod_{i}d\lambda_{i}\exp\left\\{-N\sum_{i}V(\lambda_{i})+\frac{\beta}{2}\sum_{i\neq
j}\log|\lambda_{i}-\lambda_{j}|\right\\}.$ (89)

Introducing the continuous spectral density $\tilde{\rho}_{0}$ in the usual
fashion and passing to the large-$N$ limit, the exponential in the integrand
can be written

$-N^{2}\left[\int\tilde{\rho}(\lambda)V(\lambda)d\lambda-\frac{\beta}{2}\iint\tilde{\rho}(\lambda)\tilde{\rho}(\mu)\log|\lambda-\mu|d\lambda
d\mu\right].$ (90)

The large-$N$ saddle point solution results from the condition

$V^{\prime}(\lambda)=\beta
P\int_{-a}^{a}\frac{\tilde{\rho}(\mu)}{\lambda-\mu}d\mu.$ (91)

Next, introduce the analytic function

$F(\lambda)=\int\frac{\tilde{\rho}(\mu)}{\lambda-\mu},$ (92)

which, via application of the Sokhotski–Plemelj formula, has the discontinuity

$F(\lambda\pm i\epsilon)=\frac{V^{\prime}(\lambda)}{\beta}\mp
i\pi\tilde{\rho}(\lambda).$ (93)

A solution for $\tilde{\rho}$ is found by using the ansatz

$F(\lambda)=\frac{V^{\prime}(\lambda)}{\beta}-Q(\lambda){\pi}\sqrt{\lambda^{2}-a^{2}},$
(94)

where the polynomial $Q$ and the endpoints $\pm a$ are determined by requiring
that $F\to 1/\lambda$ as $|\lambda|\to\infty$. The density is then given by

$\tilde{\rho}(\lambda)=\frac{Q(\lambda)}{\pi}\sqrt{a^{2}-\lambda^{2}}.$ (95)

For the Gaussian case $V(\lambda)=\alpha\lambda^{2}$ one finds

$Q(\lambda)=\frac{2\alpha}{\beta}\quad\&\quad a^{2}=\frac{\beta}{\alpha},$
(96)

so that the density is

$\tilde{\rho}(\lambda)=\frac{2\alpha}{\pi\beta}\sqrt{\frac{\beta}{\alpha}-\lambda^{2}}=\frac{2}{\pi
a^{2}}\sqrt{a^{2}-\lambda^{2}}.$ (97)

Given the canonical normalization $\alpha=\frac{1}{2}$, we have the following
results for $\beta=2,1$ (respectively)

$\displaystyle\tilde{\rho}(\lambda)$
$\displaystyle=\frac{\sqrt{4-\lambda^{2}}}{2\pi},$ (98)
$\displaystyle\tilde{\rho}(\lambda)$
$\displaystyle=\frac{\sqrt{2-\lambda^{2}}}{\pi}.$

Upon performing double scaling with the relation
$\lambda=\lambda_{c}^{(\beta)}-\delta^{2}E$, with $\lambda_{c}^{(2)}=2$ and
$\lambda_{c}^{(1)}=\sqrt{2}$, we have

$\displaystyle\rho(E)$ $\displaystyle=\frac{\sqrt{E}}{\pi\hbar},$ (99)
$\displaystyle\rho(E)$ $\displaystyle=\frac{2^{3/4}\sqrt{E}}{\pi\hbar},$

for $\beta=2,1$ respectively111111Some compensating factors of $N$ need to be
included to extract the scaling part of the density. These factors of $N$
combine with the left over power of $\delta$ to form the renormalized
parameter $\hbar$..

Notice that the $\beta=1$ density has the same functional form as the
$\beta=2$ density. Indeed, if we plug $\beta=2$ and $\alpha=1$ into eq. (97)
we obtain the lower result in eq. (99). This is neither a surprise nor an
accident. Of course we expect both ensembles to have the same leading
dependence on $E$ (i.e. the Wigner semi-circle), but the exact relationship
between the normalizations can be anticipated by looking at orthogonal
polynomials.

###  A.2 The crosscap resolvent and matrix model recursion

For general $\beta$, the recursion relation determining the resolvent at genus
$g=0,\frac{1}{2},1,\cdots,$ and number of boundaries $n$ is Stanford:2019vob

$\sqrt{\sigma(\lambda)}R_{g,n}(\lambda_{1},\cdots,\lambda_{n})=\frac{1}{2\pi
i}\int_{\mathcal{C}}\frac{d\lambda}{\lambda-\lambda_{1}}F_{g,n}(\lambda,I)\frac{\sqrt{\sigma(\lambda)}}{2y(\lambda)},$
(100)

where $\sigma(\lambda)=\lambda^{2}-a^{2}$, $\mathcal{C}$ is a contour
surrounding the cut $[-a,a]$, $y(\lambda)$ is the spectral curve, $I$ denotes
the set $\\{\lambda_{2},\cdots,\lambda_{n}\\}$, and

$\displaystyle F_{g,n}(\lambda,I)$
$\displaystyle=\left(1-\frac{2}{\beta}\right)\partial_{\lambda}R_{g-\frac{1}{2},n}(\lambda,I)+R_{g-1,n+1}(\lambda,\lambda,I)$
(101)
$\displaystyle\quad+\sum_{\text{stable}}R_{h,1+|J|}(\lambda,J)R_{h^{\prime},1+|J^{\prime}|}(\lambda,J^{\prime})$
$\displaystyle\quad+2\sum_{k=1}^{n}\left(R_{0,2}(\lambda,\lambda_{k})+\frac{1}{\beta}\frac{1}{(\lambda-\lambda_{k})^{2}}\right)R_{g,|I|-1}(\lambda,I\setminus\lambda_{k}).$

The subscript ‘stable’ on the sum means we ignore terms where one of the
factors is $R_{0,1}(\lambda)$ or $R_{0,2}(\lambda,\lambda_{k})$.

The only term relevant to the crosscap (with one boundary) is the first one.
We have

$F_{\frac{1}{2},1}(\lambda)=-\partial_{\lambda}y(\lambda).$ (102)

The spectral curve is proportional to the leading eigenvalue density

$y(\lambda)=-i\pi\tilde{\rho}(\lambda).$ (103)

For the GOE, the spectral curve $y$ is also proportional to $\sqrt{\sigma}$,
as can be seen in the saddle point analysis above. Before proceeding
explicitly, we perform double scaling $\lambda\to-a+E$, after which the
spectral curve is $y(E)\propto\sqrt{-E}$ which is pure imaginary on the
positive real axis. In terms of the double scaled spectral curve, the factor
$F$ in the integrand is

$F_{\frac{1}{2},1}(E)\propto\frac{1}{\sqrt{-E}}.$ (104)

The contour $\mathcal{C}$ is formally a loop around the positive $E$-axis
where the branch cut is. Depending on whether or not the integrand is single-
valued, we will either be able to take the residue at the origin or use a
standard keyhole contour. As an integral over a real variable, the domain of
integration is the nonnegative real line, $\mathbb{R}^{+}$.

The double scaled resolvent is given by

$\sqrt{-E_{1}}R_{\frac{1}{2},1}(E_{1})\propto\int_{C}\frac{dE}{E-E_{1}}\frac{1}{\sqrt{-E}}.$
(105)

The integrand has a simple pole at $E=E_{1}$ and a branch cut from $E=0$ to
$E=\infty$. Since it is not single valued, the integral cannot be written as a
residue at the origin.

One can derive a general formula for the integral

$I(f;a)=\int_{0}^{\infty}dxf(x)x^{a},$ (106)

where $a>-1$ is a non-integer real number and $f$ is a function with simple
poles on the positive real axis that satisfies

$|f(R)|\leq\frac{M}{R^{p}},$ (107)

with $p>a+1$ and $R\gg 1$121212This is enough (along with $a>-1$) to ensure
that the integral vanishes on the circles around the origin and infinity. (see
for example hayek2000advanced ). Do this by passing to the complex plane with
the complex coordinate $E$ and considering a generalization of the keyhole
contour. Beginning with the standard keyhole contour, we make indentations in
the parts parallel to the positive real axis at the poles of the function $f$.
The result is

$I(f;a)=-\frac{\pi e^{-i\pi a}}{\sin(\pi a)}\sum_{j}r_{j}+\pi\text{cotan}(\pi
a)\sum_{j}r^{\ast}_{j},$ (108)

where $r_{j}$ is the residue of $f(x)x^{a}$ at the $j$th pole in
$\mathbb{C}\setminus\mathbb{R}^{+}$ and $r^{\ast}_{j}$ is the residue at the
$j$th pole on $\mathbb{R}^{+}$ evaluated above the branch cut.

The crosscap resolvent is clearly proportional to the integral
$I\left(f;-\frac{1}{2}\right)$ with the function $f(E)=1/(E-E_{1})$. For $R\gg
1$ and $R\gg E_{1}$ we have $|f(R)|\sim 1/R$, and $1>-\frac{1}{2}+1$.
Therefore we can apply the above formula. Since we take
$E_{1}\in\mathbb{R}^{+}$, there are no residues in the rest of the complex
plane. But, since cotan$(-\pi/2)=0$, the second term in $I$ is identically 0.
Therefore we conclude that

$R_{\frac{1}{2},1}(E)=0.$ (109)

To see this another way, introduce the uniformizing coordinate $z$ via
$z^{2}=-E$. Then

$z_{1}R_{\frac{1}{2},1}(z_{1})\propto\int_{C}\frac{dz}{z^{2}-z_{1}^{2}}.$
(110)

The factor of $1/z$ in the integrand is cancelled by the Jacobian of the
change of variables. If after the mapping the contour $\mathcal{C}$ contains
both $\pm z_{1}$, then we find

$R_{\frac{1}{2},1}(z)=0,$ (111)

because the residues contribute with opposite signs. If the contour contains
neither of $\pm z_{1}$, then the resolvent is identically 0. In both cases,
the crosscap density is 0.

## References

  * (1) E. Brezin, C. Itzykson, G. Parisi and J.B. Zuber, _Planar Diagrams_ , _Commun. Math. Phys._ 59 (1978) 35.
  * (2) E. Brézin and V. Kazakov, _Exactly solvable field theories of closed strings_ , _Physics Letters B_ 236 (1990) 144.
  * (3) M.R. Douglas and S.H. Shenker, _Strings in less than one dimension_ , _Nuclear Physics B_ 335 (1990) 635.
  * (4) D.J. Gross and A.A. Migdal, _Nonperturbative two-dimensional quantum gravity_ , _Phys. Rev. Lett._ 64 (1990) 127.
  * (5) D.J. Gross and A.A. Migdal, _A nonperturbative treatment of two-dimensional quantum gravity_ , _Nuclear Physics B_ 340 (1990) 333.
  * (6) G. Hooft, _A planar diagram theory for strong interactions_ , _Nuclear Physics B_ 72 (1974) 461.
  * (7) V. Kazakov, I. Kostov and A. Migdal, _Critical properties of randomly triangulated planar random surfaces_ , _Physics Letters B_ 157 (1985) 295.
  * (8) F. David, _Planar diagrams, two-dimensional lattice gravity and surface models_ , _Nuclear Physics B_ 257 (1985) 45.
  * (9) D. Bessis, C. Itzykson and J. Zuber, _Quantum field theory techniques in graphical enumeration_ , _Advances in Applied Mathematics_ 1 (1980) 109.
  * (10) V.A. Kazakov, _The Appearance of Matter Fields from Quantum Fluctuations of 2D Gravity_ , _Mod. Phys. Lett. A_ 4 (1989) 2125.
  * (11) E.P. Wigner, _Characteristic vectors of bordered matrices with infinite dimensions_ , _Annals of Mathematics_ 62 (1955) 548.
  * (12) C.V. Johnson, _The Microstate Physics of JT Gravity and Supergravity_ , 2201.11942.
  * (13) C.V. Johnson, _Wigner meets ’t Hooft near the black hole horizon_ , 2206.03509.
  * (14) R. Jackiw, _Lower dimensional gravity_ , _Nuclear Physics B_ 252 (1985) 343.
  * (15) C. Teitelboim, _Gravitation and hamiltonian structure in two spacetime dimensions_ , _Physics Letters B_ 126 (1983) 41.
  * (16) A. Almheiri and J. Polchinski, _Models of AdS 2 backreaction and holography_, _JHEP_ 11 (2015) 014 [1402.6334].
  * (17) P. Saad, S.H. Shenker and D. Stanford, _JT gravity as a matrix integral_ , 1903.11115.
  * (18) T.G. Mertens and G.J. Turiaci, _Liouville quantum gravity — holography, JT and matrices_ , _Journal of High Energy Physics_ 2021 (2021) .
  * (19) C.V. Johnson, _Jackiw-Teitelboim supergravity, minimal strings, and matrix models_ , _Phys. Rev. D_ 103 (2021) 046012 [2005.01893].
  * (20) C.V. Johnson, _Nonperturbative Jackiw-Teitelboim gravity_ , _Phys. Rev. D_ 101 (2020) 106023 [1912.03637].
  * (21) C.V. Johnson, _Explorations of nonperturbative Jackiw-Teitelboim gravity and supergravity_ , _Phys. Rev. D_ 103 (2021) 046013 [2006.10959].
  * (22) C.V. Johnson, _Consistency Conditions for Non-Perturbative Completions of JT Gravity_ , 2112.00766.
  * (23) A.H. Chamseddine, _Superstrings in arbitrary dimensions_ , _Phys. Lett. B_ 258 (1991) 97.
  * (24) S. Forste and I. Golla, _Nearly AdS 2 sugra and the super-Schwarzian_, _Phys. Lett. B_ 771 (2017) 157 [1703.10969].
  * (25) Y. Fan and T.G. Mertens, _Supergroup Structure of Jackiw-Teitelboim Supergravity_ , 2106.09353.
  * (26) D. Stanford and E. Witten, _JT gravity and the ensembles of random matrix theory_ , _Adv. Theor. Math. Phys._ 24 (2020) 1475 [1907.03363].
  * (27) C.V. Johnson, F. Rosso and A. Svesko, _Jackiw-Teitelboim supergravity as a double-cut matrix model_ , _Phys. Rev. D_ 104 (2021) 086019 [2102.02227].
  * (28) P.J. Forrester, _Log-gases and random matrices (lms-34)_ , in _Log-Gases and Random Matrices (LMS-34)_ , Princeton University Press (2010).
  * (29) J. Gomis and A. Kapustin, _Two-dimensional unoriented strings and matrix models_ , _Journal of High Energy Physics_ 2004 (2004) 002.
  * (30) O. Bergman and S. Hirano, _The cap in the hat: unoriented 2d strings and matrix(-vector) models_ , _Journal of High Energy Physics_ 2004 (2004) 043.
  * (31) G.R. Harris, _Loops and unoriented strings_ , _Nuclear Physics B_ 356 (1991) 685.
  * (32) D. Stanford, _A Mirzakhani recursion for non-orientable surfaces_ , 2303.04049.
  * (33) M.L. Mehta, _Random matrices_ , Elsevier (2004).
  * (34) E. Brezin and H. Neuberger, _Multicritical points of unoriented random surfaces_ , _Nucl. Phys. B_ 350 (1991) 513.
  * (35) J.E. Carlisle and C.V. Johnson, _Unoriented minimal type 0 strings_ , _Journal of High Energy Physics_ 2004 (2005) 073.
  * (36) G. Mahoux and M.L. Mehta, _A method of integration over matrix variables: Iv_ , _J. Phys. I France_ 1 (1991) 1093.
  * (37) M. Adler, P.J. Forrester, T. Nagao and P. van Moerbeke, _Classical Skew Orthogonal Polynomials and Random Matrices_ , _Journal of Statistical Physics_ 99 (2000) 141.
  * (38) V.G. Drinfeld and V.V. Sokolov, _Lie algebras and equations of Korteweg-de Vries type_ , _J. Sov. Math._ 30 (1984) 1975.
  * (39) P. Ginsparg and G. Moore, _Lectures on 2d gravity and 2d string theory (tasi 1992)_ , hep-th/9304011.
  * (40) F.J. Dyson, _Statistical theory of the energy levels of complex systems. i_ , _Journal of Mathematical Physics_ 3 (1962) 140.
  * (41) M. Pillai, J. Goglio and T.G. Walker, _Matrix Numerov method for solving Schrödinger’s equation_ , _American Journal of Physics_ 80 (2012) 1017.
  * (42) G.R. Harris and E.J. Martinec, _Unoriented strings and matrix ensembles_ , _Physics Letters B_ 245 (1990) 384.
  * (43) A. Zee, _Quantum Field Theory in a Nutshell: (Second Edition)_ , In a Nutshell, Princeton University Press (2010).
  * (44) M. Adler, E. Horozov and P. van Moerbeke, _The pfaff lattice and skew-orthogonal polynomials_ , 1999.
  * (45) S. Dalley, C.V. Johnson and T.R. Morris, _Nonperturbative two-dimensional quantum gravity_ , _Nucl. Phys. B_ 368 (1992) 655.
  * (46) S. Hayek, _Advanced Mathematical Methods in Science and Engineering_ , Taylor & Francis (2000).