Notes on matrix models (matrix musings)

We begin by considering the following family of integrals

$$\begin{equation}{\mathcal{I}}_{N}={\int }_{\mathbb{R}}\mathrm{d}x\enspace {\mathrm{e}}^{-Nf\left(x\right)},\end{equation} \tag{ 2.1 }$$

where N is a positive integer which will eventually be taken to be large, and f(x) is a real valued function. In the limit where N becomes large, the exponential causes the integrand to peak sharply at the minima of the function f(x), while all other values are suppressed. Of all extrema, the integral ${\mathcal{I}}_{N}$ will be dominated by the one which minimises f(x) as N becomes large, and we denote this by x_e. Expanding f(x) in a Taylor series around x_e

$$\begin{equation}f\left(x\right)=f\left({x}_{e}\right)+\frac{{f}^{{\prime\prime}}\left({x}_{e}\right)}{2!}\enspace \delta {x}^{2}+\dots ,\quad \delta x\equiv x-{x}_{e},\end{equation} \tag{ 2.2 }$$

one can approximate ${\mathcal{I}}_{N}$ as

$$\begin{equation}{\mathcal{I}}_{N}\approx {\mathrm{e}}^{-Nf\left({x}_{e}\right)}{\int }_{\mathbb{R}}\mathrm{d}\delta x\enspace {\mathrm{e}}^{-\frac{N}{2}{f}^{{\prime\prime}}\left({x}_{e}\right)\delta {x}^{2}}\left(1-\frac{N}{3!}{f}^{\left(3\right)}\left({x}_{e}\right)\enspace \delta {x}^{3}+\dots \right).\end{equation} \tag{ 2.3 }$$

Upon redefining $\delta x=\delta \xi /\sqrt{N}$, we see that all higher powers of δξ are suppressed at large N, and we can compute the sub-leading correction to ${\mathcal{I}}_{N}$ in the large N expansion:

$$\begin{equation}{\mathcal{I}}_{N}\approx \sqrt{\frac{2\pi }{N{f}^{{\prime\prime}}\left({x}_{e}\right)}}\enspace {\mathrm{e}}^{-Nf\left({x}_{e}\right)}\left(1+\mathcal{O}\left(1/N\right)\right).\end{equation} \tag{ 2.4 }$$

If our function f(x) has multiple minima, we must scan for the global one since the others will give contributions that are exponentially suppressed. If there are degenerate minima, on which f(x) takes the same value, we must include them all to get a good approximation. Finally, if there are an infinite number of minima care must be exercised in our approximation scheme.

Example. As a simple example we consider the integral

$$\begin{equation}{\mathcal{J}}_{N}={\int }_{0}^{\infty }\mathrm{d}x\enspace {x}^{N}\enspace {\mathrm{e}}^{-Nx}={\int }_{0}^{\infty }\mathrm{d}x\enspace {\mathrm{e}}^{-N\left(x-\mathrm{log}x\right)}.\end{equation} \tag{ 2.5 }$$

This integral can be solved exactly by performing N partial integrations yielding (N − 1)!/N^N. Moreover, we note that taking f(x) = x − log x the integral is of the form (2.1). The extremum of f(x) is given by x_e = 1, with f(x_e) = 1. Expanding f(x) around x_e in a Taylor series up to second order and inserting it into (2.1) we get

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\enspace {\mathcal{J}}_{N}\approx \sqrt{\frac{2\pi }{N}}\enspace {\mathrm{e}}^{-N}\left(1+\mathcal{O}\left(1/N\right)\right).\end{equation} \tag{ 2.6 }$$

The same expression can be obtained by expanding the exact result ${\mathcal{J}}_{N}=\left(N-1\right)!/{N}^{N}$ using the Stirling approximation for the factorial.

2.1.1. Complex saddles and contour deformations

As a next step we consider the saddle point approximation for integrals with N variables ${x}_{I}\in \mathbb{R}$ containing a vector index I = 1, 2, ..., N. Consider the family of integrals

$$\begin{equation}{\mathcal{V}}_{N}={\int }_{{\mathbb{R}}^{N}}\prod _{I=1}^{N}\mathrm{d}{x}_{I}\enspace {\mathrm{e}}^{-Nf\left({x}_{I}{x}_{I}/N\right)},\end{equation} \tag{ 2.7 }$$

where we use the double index notation ${x}_{I}{x}_{I}\equiv {\sum }_{I=1}^{N}{x}_{I}{x}_{I}$ here and throughout. The integrals ${\mathcal{V}}_{N}$ are invariant under O(N) rotations of the x_I. By changing to spherical coordinates

$$\begin{equation}{x}_{I}=\sqrt{N}R\enspace {{\Omega}}_{I},\qquad \sum _{I=1}^{N}{{\Omega}}_{I}^{2}=1,\end{equation} \tag{ 2.8 }$$

we reduce the integral to the one-dimensional case we already worked out in the previous section. Importantly, this change of variables produces a nontrivial Jacobian leading to the following expression

$$\begin{equation}{\mathcal{V}}_{N}=\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}^{N-1}\enspace {N}^{N/2}\enspace {\int }_{0}^{\infty }\mathrm{d}R\enspace {\mathrm{e}}^{\left(N-1\right)\mathrm{log}R-Nf\left({R}^{2}\right)}.\end{equation} \tag{ 2.9 }$$

Happily, in the large N limit the Jacobian of the coordinate change competes at the same order with the original integrand in (2.7). Thus, we can employ the saddle point approximation in its simplest form to estimate ${\mathcal{V}}_{N}$ at large N.

Example. As a straightforward example we can approximate the volume of the N-sphere S^N. To this end, we take f(x_Ix_I) = x_Ix_I in (2.7) and compute the integral

$$\begin{equation}{\mathcal{V}}_{N}={\int }_{{\mathbb{R}}^{N}}\prod _{I=1}^{N}\mathrm{d}{x}_{I}\enspace {\mathrm{e}}^{-{x}_{I}{x}_{I}}\end{equation} \tag{ 2.10 }$$

in two different ways. An explicit calculation using Gaussian integrals tells us that the value of (2.10) gives ${\mathcal{V}}_{N}={\pi }^{N/2}$. On the other hand, we can change to spherical coordinates and perform a saddle point approximation. Taking care of the Jacobian we rewrite (2.10) as

$$\begin{equation}{\mathcal{V}}_{N}=\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}^{N-1}\enspace {N}^{N/2}\enspace {\int }_{0}^{\infty }\mathrm{d}R\enspace {\mathrm{e}}^{-NF\left(R\right)},\end{equation} \tag{ 2.11 }$$

with F(R) ≡ R² − (1 − 1/N)log R which at large N can be approximated as F(R) ≈ R² − log R. The dominant extremum of F(R) is given by ${R}_{e}=1/\sqrt{2}$ and so we obtain for the leading term of the integral

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\enspace {\mathcal{V}}_{N}=\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}^{N-1}\sqrt{\frac{\pi }{8N}}{\left(\frac{N}{2}\right)}^{N/2}\enspace {\mathrm{e}}^{-N/2}.\end{equation} \tag{ 2.12 }$$

Setting this equal to (2.10) we find an expression for the volume of S^N in the limit N ≫ 1

$$\begin{equation}\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}^{N}\approx \sqrt{N}{\left(\frac{2\pi e}{N}\right)}^{N/2}.\end{equation} \tag{ 2.13 }$$

This agrees with the large N limit of the volume of an N-sphere with known exact expression vol S^N−1 = Nπ^N/2/Γ(N/2 + 1).

We would now like to consider the vector-like integrals (2.7) from a slightly different perspective that is motivated by perturbative expansions often employed in the context of quantum field theory. We will do so by studying a specific example given by

$$\begin{equation}{\mathcal{Z}}_{N}\left(\alpha \right)={\int }_{{\mathbb{R}}^{N}}\prod _{I=1}^{N}\mathrm{d}{x}_{I}\enspace {\mathrm{e}}^{-N\left({x}_{I}{x}_{I}/N+\alpha {\left({x}_{I}{x}_{I}/N\right)}^{2}\right)},\end{equation} \tag{ 2.14 }$$

where $\alpha \in \mathbb{R}$ will be taken to be a small parameter. For ${\mathcal{Z}}_{N}\left(\alpha \right)$ to be well-defined we should further take α to be positive, but we shall see shortly that in the large N limit, it may be sensible to allow α to also take small negative values. We will first consider the integral (2.14) in the small α limit, for which we can Taylor expand the exponential and calculate the correction terms using Wick contractions. To this end, it is convenient to introduce the 'propagator' ⟨x_Ix_J⟩ given by

$$\begin{equation}\langle {x}_{I}{x}_{J}\rangle ={\mathcal{Z}}_{N}^{-1}\left(0\right){\int }_{{\mathbb{R}}^{N}}\prod _{K=1}^{N}\mathrm{d}{x}_{K}\enspace {x}_{I}{x}_{J}\enspace {\mathrm{e}}^{-{x}_{I}{x}_{I}}=\frac{1}{2}\enspace {\delta }_{IJ}.\end{equation} \tag{ 2.15 }$$

The corresponding graphical representation for the propagator and the quartic vertex in the integrand of (2.15) are dispayed in figure 1.

Zoom In Zoom Out Reset image size

Using the propagator and the quartic vertex we can form closed 'bubble' diagrams which compute perturbative contributions to ${\mathcal{Z}}_{N}\left(\alpha \right)$. The set of connected bubble diagrams is generated by ${\mathcal{F}}_{N}\left(\alpha \right)\equiv -\mathrm{log}\left({\mathcal{Z}}_{N}\left(\alpha \right)/{\mathcal{Z}}_{N}\left(0\right)\right)$. To leading order, we can approximate ${\mathcal{Z}}_{N}\left(\alpha \right)$ by

$$\begin{equation}{\mathcal{Z}}_{N}\left(\alpha \right)={\int }_{{\mathbb{R}}^{N}}\prod _{K=1}^{N}\enspace \mathrm{d}{x}_{K}\enspace {\mathrm{e}}^{-{x}_{I}{x}_{I}}\left(1-\alpha \enspace N{\left(\frac{{x}_{I}{x}_{I}}{N}\right)}^{2}+\cdots \enspace \right).\end{equation} \tag{ 2.16 }$$

One can evaluate the perturbative terms in the above integral by standard Gaussian integration. The first order correction in the small α expansion yields

$$\begin{equation}\alpha \enspace N\enspace {\int }_{{\mathbb{R}}^{N}}\prod _{K=1}^{N}\enspace \mathrm{d}{x}_{K}\enspace {\mathrm{e}}^{-{x}_{I}{x}_{I}}{\left(\frac{{x}_{I}{x}_{I}}{N}\right)}^{2}=\frac{\alpha }{4N}\sum _{I,J}\left({\delta }_{II}{\delta }_{JJ}+2{\delta }_{IJ}{\delta }_{IJ}\right)=\frac{\alpha }{4}\left(N+2\right).\end{equation} \tag{ 2.17 }$$

Diagrammatically, the above integral corresponds to a bubble diagram where we have closed the propagators emanating from the quartic vertex. Thus, to linear order in α we have

$$\begin{equation}{\mathrm{lim}}_{\alpha \to 0}\enspace {\mathcal{F}}_{N}\left(\alpha \right)=\frac{\alpha }{4}\left(N+2\right)+\mathcal{O}\left({\alpha }^{2}\right).\end{equation} \tag{ 2.18 }$$

As we proceed with higher order corrections, we encounter a variety of diagrams. For any given power of α a certain class of diagrams, known as cactus diagram, will dominate at large N. In the figure below, we demonstrate several diagrams of which the first three belong to the cactus family (figure 2).

Zoom In Zoom Out Reset image size

The fact that at large N only a subclass of diagrams survives leads to a dramatic simplification of the theory, and to the hope that a resummation can be performed. This resummation is precisely the saddle-point approximation of ${\mathcal{Z}}_{N}\left(\alpha \right)$. Let us perform the saddle point analysis to confirm the $\mathcal{O}\left(\alpha \right)$ correction we calculated perturbatively. To do so, let us go to spherical coordinates and perform the same steps as in (2.9). We find

$$\begin{equation}{\mathcal{Z}}_{N}\left(\alpha \right)=\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}^{N-1}{N}^{N/2}{\int }_{0}^{\infty }\mathrm{d}R\enspace {\mathrm{e}}^{-N{F}_{\alpha }\left(R\right)},\end{equation} \tag{ 2.19 }$$

where we have defined ${F}_{\alpha }\left(R\right)\equiv \left({R}^{2}+\alpha {R}^{4}-\left(1-1/N\right)\mathrm{log}\enspace R\right)$. To leading order in the large N limit, the extrema of F_α(R) are given by

$$\begin{equation}{R}_{{\pm}}^{2}=\frac{1}{4\alpha }\left(-1{\pm}\sqrt{1+4\alpha }\right).\end{equation} \tag{ 2.20 }$$

Of the two saddles, only ${R}_{+}^{2}$ lies on the original contour of integration and has a well-defined limit as α tends to zero. The other saddle ${R}_{-}^{2}$ goes to infinity in the limit of vanishing α and hence cannot be relevant to compare to the perturbative analysis. Using the saddle point approximation, we are thus led to the leading N approximation of ${\mathcal{F}}_{N}\left(\alpha \right)$:

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\enspace {\mathcal{F}}_{N}\left(\alpha \right)=\frac{N}{8\alpha }\left(-1+\sqrt{1+4\alpha }-\alpha \left(2-4\enspace \mathrm{log}\enspace 2\right)-4\alpha \enspace \mathrm{log}\left(\frac{-1+\sqrt{1+4\alpha }}{\alpha }\enspace \right)\right)+\mathcal{O}\left(1\right).\end{equation} \tag{ 2.21 }$$

The above result is valid to all orders in the small α expansion. Expanding (2.21) to leading order in α one can readily recover our perturbative correction (2.17). It is remarkable that the resummation of all the diagrams was encoded in the solution of the simple algebraic equation ∂_RF_N(R) = 0. Moreover (2.21) reveals that the expression gives a sensible result for Re α ⩾ −1/4 which includes a negative interval. To continue beyond the critical value α_c = −1/4 one has to specify a branch cut due to the non-analytic behaviour. That negative values of α might be sensible in the large N limit, even though the finite N integral is ill-defined for negative α, comes about precisely because at large N most of the contribution to the integral is localised near the critical value of R₊. Fortuitously, the integral (2.19) can be evaluated exactly for any value of N. One finds

$$\begin{equation}\mathrm{exp}\left(-{\mathcal{F}}_{N}\left(\alpha \right)\right)={\left(\frac{N}{4\alpha }\right)}^{\frac{N}{4}}U\left(\frac{N}{4},\frac{1}{2},\frac{N}{4\alpha }\right),\end{equation} \tag{ 2.22 }$$

where U(a, b, z) is Tricomi's confluent hypergeometric function. We have

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\enspace {\mathrm{lim}}_{\alpha \to 0}{\left(\frac{N}{4\alpha }\right)}^{\frac{N}{4}}U\left(\frac{N}{4},\frac{1}{2},\frac{N}{4\alpha }\right)=1-\frac{\alpha }{4}\left(N+2\right)+\mathcal{O}\left({\alpha }^{2}\right).\end{equation} \tag{ 2.23 }$$

This agrees again with the exponential of (2.18). In relation to the critical value R₊, it is worth mentioning a quantity that will become of increasing relevance throughout our discussion. Recall that the generating function of connected bubble diagrams admits a large N expansion of the form

$$\begin{equation}{\mathcal{F}}_{N}\left(\alpha \right)=N{\mathcal{F}}^{\left(0\right)}\left(\alpha \right)+{\mathcal{F}}^{\left(1\right)}\left(\alpha \right)+\frac{1}{N}\enspace {\mathcal{F}}^{\left(2\right)}\left(\alpha \right)+\cdots .\end{equation} \tag{ 2.24 }$$

Moreover, each of the functions ${\mathcal{F}}^{\left(k\right)}\left(\alpha \right)$ themselves admits a Taylor expansion near α = 0. For instance

$$\begin{equation}{\mathcal{F}}^{\left(0\right)}\left(\alpha \right)=\sum _{n=0}^{\infty }{f}_{n}^{\left(0\right)}{\alpha }^{n}.\end{equation} \tag{ 2.25 }$$

The power of α counts the number of bubbles in a given cactus diagram. Consequently, we can define the moments of the 'bubble number' as

$$\begin{equation}\langle {n}^{L}\rangle ={\partial }_{\mathrm{log}\alpha }^{\left(L\right)}\mathrm{log}\enspace {\mathcal{F}}^{\left(0\right)}\left(\alpha \right).\end{equation} \tag{ 2.26 }$$

Near the critical value α_c, and for L ⩾ 2, the above moments will diverge. For instance $\langle {n}^{2}\rangle \sim {\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{-1/2}$. This suggests that in the limit α → α_c the contribution of cacti with arbitrarily large numbers of bubbles becomes increasingly important. Moreover, the behaviour near α_c is somewhat reminiscent of critical phenomena.

2.4.1. Order of limits and analytic continuation

Every Hermitian matrix can be diagonalised using a unitary matrix $\mathcal{U}$ as $M=\mathcal{U}{D}_{M}{\mathcal{U}}^{{\dagger}}$ with D_M ≡ diag(λ₁, λ₂, ..., λ_N) a real diagonal matrix. Consequently, the exponent of the integrand in (3.1) depends purely on the N eigenvalues ${\lambda }_{I}\in \mathbb{R}$ of M. Thus, it is natural to consider the problem using a new set of variables given by an N × N unitary matrix $\mathcal{U}\in U\left(N\right)/U{\left(1\right)}^{N}$, and N real variables λ_I, such that $M=\mathcal{U}{D}_{M}{\mathcal{U}}^{{\dagger}}$. Notice that $\mathcal{U}$ resides in the coset space U(N)/U(1)^N, rather than the full U(N), because those unitary matrices where each diagonal entry is multiplied by a pure phase are redundant. One must also compute the Jacobian associated to the change of variables $M=\mathcal{U}{D}_{M}{\mathcal{U}}^{{\dagger}}$. To this end, it is useful to write $\mathcal{U}={\text{e}}^{\mathrm{i}\mathcal{L}}$, with $\mathcal{L}$ an N × N Hermitian generator of U(N)/U(1)^N. The infinitesimal line element on the space of Hermitian matrices M is

$$\begin{equation}\mathrm{d}{s}^{2}=\mathrm{T}\mathrm{r}\enspace \mathrm{d}M\mathrm{d}{M}^{{\dagger}}=\sum _{I}\mathrm{d}{\lambda }_{I}^{2}+2\sum _{I{< }J}{\left({\lambda }_{I}-{\lambda }_{J}\right)}^{2}\vert \mathrm{d}{\mathcal{L}}_{IJ}{\vert }^{2}.\end{equation} \tag{ 3.3 }$$

The Jabobian is the determinant of the above line element. It follows that

$$\begin{equation}\left[\mathrm{d}M\right]=\frac{1}{\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}_{N}}\prod _{I}\mathrm{d}{\lambda }_{I}\prod _{I\ne J}\vert {\lambda }_{I}-{\lambda }_{J}\vert \left[\mathrm{d}\mathcal{L}\right],\end{equation} \tag{ 3.4 }$$

where $\left[\mathrm{d}\mathcal{L}\right]$ is the volume element of the coset space U(N)/U(1)^N and we have further divided by the volume of the permutation group S_N to avoid multiple counting of eigenvalue matrices D_M with permuted elements.

Example. As an example, let us consider the case N = 2. We can parametrise the unitary matrix $\mathcal{U}$ as

$$\begin{equation}\mathcal{U}=\left(\begin{matrix}\hfill \mathrm{sin}\enspace \frac{\theta }{2}\hfill & \hfill {\mathrm{e}}^{\mathrm{i}\varphi }\enspace \mathrm{cos}\enspace \frac{\theta }{2}\hfill \\ \hfill -{\mathrm{e}}^{-\mathrm{i}\varphi }\enspace \mathrm{cos}\enspace \frac{\theta }{2}\hfill & \hfill \mathrm{sin}\enspace \frac{\theta }{2}\hfill \end{matrix}\right),\end{equation} \tag{ 3.5 }$$

where θ ∈ [0, π) and φ ∈ (0, 2π] are coordinates on U(2)/U(1)² whose Haar metric is the round metric on S² with radius squared 1/2. Using that the relation $\mathcal{U}={\text{e}}^{\mathrm{i}\mathcal{L}}$ implies

$$\begin{equation}\mathrm{d}\mathcal{L}=-\mathrm{i}\enspace {\mathcal{U}}^{-1}\mathrm{d}\mathcal{U}.\end{equation} \tag{ 3.6 }$$

The line element (3.3) can be written out explicitly:

$$\begin{equation}\mathrm{d}{s}^{2}=\mathrm{d}{\lambda }_{1}^{2}+\mathrm{d}{\lambda }_{2}^{2}+\frac{{\left({\lambda }_{1}-{\lambda }_{2}\right)}^{2}}{2}\left(\mathrm{d}{\theta }^{2}+{\mathrm{sin}}^{2}\enspace \theta \mathrm{d}{\varphi }^{2}\right).\end{equation} \tag{ 3.7 }$$

Defining $r=\left({\lambda }_{1}-{\lambda }_{2}\right)/\sqrt{2}$ and $z=\left({\lambda }_{1}+{\lambda }_{2}\right)/\sqrt{2}$ the above becomes

$$\begin{equation}\mathrm{d}{s}^{2}=\mathrm{d}{z}^{2}+\mathrm{d}{r}^{2}+{r}^{2}\left(\mathrm{d}{\theta }^{2}+{\mathrm{sin}}^{2}\enspace \theta \mathrm{d}{\varphi }^{2}\right).\end{equation} \tag{ 3.8 }$$

We observe that the above metric is indeed the flat metric on ${\mathbb{R}}^{4}$ in cylindrical coordinates up to the fact that $r\in \mathbb{R}$ as opposed to r ⩾ 0. However, the map r → −r is in fact a permutation of the eigenvalues λ₁ ↔ λ₂. This is a double counting in our configuration space, since conjugation with respect to $\mathcal{U}$ in (3.5) with θ = 0 exchanges λ₁ and λ₂. Thus, after dividing by the volume of the two-dimensional permutation group the volume element becomes

$$\begin{equation}\left[\mathrm{d}M\right]={r}^{2}\enspace \mathrm{d}r\enspace \mathrm{d}z\enspace \mathrm{d}{{\Omega}}_{2},\quad r{ >}0,\end{equation} \tag{ 3.9 }$$

where dΩ₂ is the standard volume element on the unit sphere S². The right hand side is indeed the flat metric on ${\mathbb{R}}^{4}$ expressed in cylindrical coordinates.

3.1.1. The Vandermonde contribution

The product appearing in the volume element (3.4) is the determinant of the Vandermonde matrix

$$\begin{equation}{\mathbf{V}}_{N}\equiv \left(\begin{matrix}\hfill 1\hfill & \hfill {\lambda }_{1}\hfill & \hfill {\lambda }_{1}^{2}\hfill & \hfill \cdots \hfill & \hfill {\lambda }_{1}^{N-1}\hfill \\ \hfill 1\hfill & \hfill {\lambda }_{2}\hfill & \hfill {\lambda }_{2}^{2}\hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill \hfill \\ \hfill 1\hfill & \hfill {\lambda }_{N}\hfill & \hfill {\lambda }_{N}^{2}\hfill & \hfill \dots \hfill & \hfill {\lambda }_{N}^{N-1}\hfill \end{matrix}\right).\end{equation} \tag{ 3.10 }$$

Since it will appear several times throughout our discussion, we introduce the following notation for the determinant of the Vandermonde matrix

$$\begin{equation}{{\Delta}}_{N}\left(\lambda \right)\equiv \prod _{I{< }J}\left({\lambda }_{I}-{\lambda }_{J}\right).\end{equation} \tag{ 3.11 }$$

We can bring our original integral (3.1) to the following form

$$\begin{equation}{\mathcal{M}}_{N}=\mathrm{v}\mathrm{o}\mathrm{l}\enspace \enspace \frac{U\left(N\right)}{U{\left(1\right)}^{N}{\times}{S}_{N}}\enspace {\int }_{{\mathbb{R}}^{N}}\prod _{I=1}^{N}\mathrm{d}{\lambda }_{I}\enspace {\mathrm{e}}^{-{N}^{2}S\left[\lambda \right]},\end{equation} \tag{ 3.12 }$$

where we have defined the multivariable function

$$\begin{equation}S\left[\lambda \right]\equiv \frac{1}{N}\sum _{I=1}^{N}V\left({\lambda }_{I}\right)-\frac{1}{{N}^{2}}\sum _{J\ne I}\mathrm{log}\vert {\lambda }_{I}-{\lambda }_{J}\vert .\end{equation} \tag{ 3.13 }$$

As for the vector integral case, the change of coordinates adds a nontrivial Jacobian that will compete with the original integrand. Moreover, the Jacobian causes an effective 'repulsive' pressure against the eigenvalues all lying on top of each other. The repulsive effect is the matrix analogue of the logarithmic term encountered in (2.9) for the vector integrals.

Before proceeding to solve (3.1) in the large N limit, let us summarise what we have achieved so far. By a unitary transformation of the Hermitian matrix M we reduced the matrix integral from a theory with N² independent degrees of freedom to a theory described by the N real eigenvalues of the Hemitian matrix M. The Jacobian arising from this change of coordinates gives rise to the square of the Vandermonde determinant Δ_N(λ)². This contribution competes with V(λ) such that the eigenvalues are distributed around the mininum of V(λ).

In order to solve the matrix integrals at large N, we will once again employ the saddle point approximation. The saddle point equations for S[λ] in (3.13) are given by

$$\begin{equation}{V}^{\prime }\left({\lambda }_{I}\right)=\frac{2}{N}\sum _{J\ne I}{\left({\lambda }_{I}-{\lambda }_{J}\right)}^{-1}.\end{equation} \tag{ 3.14 }$$

To solve these equations, it is convenient to introduce the following normalised eigenvalue distribution

$$\begin{equation}\rho \left(\lambda \right)=\frac{1}{N}\sum _{I=1}^{N}\delta \left(\lambda -{\lambda }_{I}\right),\qquad {\int }_{-a}^{a}\mathrm{d}\lambda \enspace \rho \left(\lambda \right)=1.\end{equation} \tag{ 3.15 }$$

In the large N limit, and under the assumption that the range [−a, a] is compact and does not grow with N, we can approximate ρ(λ) by a continuous, non-negative function. For simplicity, we further assume that the range of ρ(λ) is symmetric about the origin. This is not at great cost since we can always shift all the λ by a constant to accommodate this. We can then rewrite (3.14) as the following integral equation

$$\begin{equation}{V}^{\prime }\left(\lambda \right)=2\enspace {\int }_{-a}^{a} -\enspace \mathrm{d}\mu \enspace \frac{\rho \left(\mu \right)}{\lambda -\mu },\end{equation} \tag{ 3.16 }$$

where it is understood that we are taking the principal value for the integral. Our goal is to find the eigenvalue distribution ρ(λ) which solves the above equation. To do so, it is instructive to introduce the resolvent [39]

$$\begin{equation}{\mathrm{R}}_{N}\left(z\right)\equiv \frac{1}{N}\enspace \mathrm{T}\mathrm{r}{\left(z\enspace {\mathbb{I}}_{N}-M\right)}^{-1}=\frac{1}{N}\sum _{I=1}^{N}\frac{1}{z-{\lambda }_{I}},\quad z\in \mathbb{C}/\left\{{\lambda }_{I}\right\}.\end{equation} \tag{ 3.17 }$$

Sending N → ∞ the sum can be replaced by an integral where each eigenvalue is weighted by its average density

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\enspace {\mathrm{R}}_{N}\left(z\right)\equiv \mathrm{R}\left(z\right)={\int }_{-a}^{a}\mathrm{d}\mu \enspace \frac{\rho \left(\mu \right)}{z-\mu }.\end{equation} \tag{ 3.18 }$$

By evaluating the resolvent close to the real axis z = x ± i we find

$$\begin{equation}\mathrm{R}\left(x+\mathrm{i}{\epsilon}\right)={\int }_{-a}^{a}\mathrm{d}\mu \enspace \frac{\rho \left(\mu \right)}{x-\mu +\mathrm{i}{\epsilon}}={\int }_{-a}^{a}\mathrm{d}\mu \enspace \frac{\rho \left(\mu \right)\left(x-\mu \right)}{{\left(x-\mu \right)}^{2}+{{\epsilon}}^{2}}-\mathrm{i}{\int }_{-a}^{a}\mathrm{d}\mu \enspace \frac{\rho \left(\mu \right){\epsilon}}{{\left(x-\mu \right)}^{2}+{{\epsilon}}^{2}},\end{equation} \tag{ 3.19 }$$

where in the limit → 0 the first integral is a principal value integral, while we evaluate the second integral using delta function identities. Finally, we arrive at

$$\begin{equation}\mathrm{R}\left(x{\pm}\mathrm{i}{\epsilon}\right)={\int }_{-a}^{a} -\enspace \mathrm{d}\mu \enspace \frac{\rho \left(\mu \right)}{\left(x-\mu \right)}\mp \mathrm{i}\pi \rho \left(x\right),\end{equation} \tag{ 3.20 }$$

where we have employed the Sokhotski–Plemelj theorem. Using (3.14) and (3.15)

$$\begin{equation}\mathrm{R}\left(x{\pm}\mathrm{i}{\epsilon}\right)=\frac{1}{2}\enspace {V}^{\prime }\left(x\right)\mp \mathrm{i}\pi \rho \left(x\right).\end{equation} \tag{ 3.21 }$$

Using (3.18) and the above expressions we obtain the following properties the resolvent must satisfy:

$$\begin{equation}{\mathrm{r}\mathrm{e}\mathrm{s}}_{a}:\enspace {\mathrm{lim}}_{z\to \infty }\enspace \mathrm{R}\left(z\right)=\frac{1}{z},\end{equation} \tag{ 3.22 }$$

$$\begin{equation}{\mathrm{r}\mathrm{e}\mathrm{s}}_{b}:\enspace \rho \left(x\right)=\frac{1}{2\pi \mathrm{i}}\left(\mathrm{R}\left(x-\mathrm{i}{\epsilon}\right)-\mathrm{R}\left(x+\mathrm{i}{\epsilon}\right)\right),\quad x\in \mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\left(\rho \right),\end{equation} \tag{ 3.23 }$$

$$\begin{equation}{\mathrm{r}\mathrm{e}\mathrm{s}}_{c}:\enspace {V}^{\prime }\left(x\right)=\mathrm{R}\left(x+\mathrm{i}{\epsilon}\right)+\mathrm{R}\left(x-\mathrm{i}{\epsilon}\right),\quad x\in \mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\left(\rho \right).\end{equation} \tag{ 3.24 }$$

The condition res_b specifies the jump of the resolvent across the branch cut x ∈ [−a, a]. Elsewhere, the resolvent is analytic. Finally, res_c fixes the real part of the resolvent to the derivative of the potential. Thus, the class of problems we are trying to solve is a Riemann–Hilbert type problem.

In the large N limit the exponent of our original integral (3.12) takes the form

$$\begin{equation}S\left[\rho \left(\lambda \right)\right]={\int }_{-a}^{a}\mathrm{d}\lambda \enspace \rho \left(\lambda \right)V\left(\lambda \right)-{\int }_{-a}^{a}\mathrm{d}\lambda \enspace \rho \left(\lambda \right)\enspace {\int }_{-a}^{a}\mathrm{d}\mu \enspace \rho \left(\mu \right)\mathrm{log}\vert \lambda -\mu \vert .\end{equation} \tag{ 3.25 }$$

Using the saddle point equations, we can further simplify

$$\begin{equation}S\left[{\rho }_{\text{ext}}\left(\lambda \right)\right]=\frac{1}{2}{\int }_{-a}^{a}\mathrm{d}\lambda \enspace {\rho }_{\text{ext}}\left(\lambda \right)V\left(\lambda \right)-{\int }_{-a}^{a}\mathrm{d}\lambda \enspace {\rho }_{\text{ext}}\left(\lambda \right)\enspace \mathrm{log}\vert \lambda \vert ,\end{equation} \tag{ 3.26 }$$

where, for simplicity, we have further assumed that V(0) = 0 and ρ_ext(λ) is that solution of (3.16) which minimises S[ρ(λ)] in (3.25). Thus, our original integral (3.12) is approximated by the expression

$$\begin{equation}{\mathcal{M}}_{N}\approx \mathrm{v}\mathrm{o}\mathrm{l}\enspace U\left(N\right)\enspace {\mathrm{e}}^{-{N}^{2}S\left[{\rho }_{\text{ext}}\left(\lambda \right)\right]},\end{equation} \tag{ 3.27 }$$

where we have dropped the volume of U(1)^N × S_N since it does not compete with terms that grow exponentially in N².

Gaussian example. We now consider the simplest concrete example. Let us take $V\left(M\right)=\frac{1}{2}{M}^{2}$ such that

$$\begin{equation}{\mathcal{M}}_{N}={\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]\enspace {\mathrm{e}}^{-\frac{N}{2}\mathrm{T}\mathrm{r}\enspace {M}^{2}}={\left(\frac{2\pi }{N}\right)}^{{N}^{2}/2},\end{equation} \tag{ 3.28 }$$

where we obtained the exact result using Gaussian integrals. We will now approximate (3.28) using the saddle point approximation. As was already outlined, this involves finding the resolvent R(z). Let us consider the following ansatz

$$\begin{equation}\mathrm{R}\left(z\right)=\frac{1}{2}{V}^{\prime }\left(z\right)-P\left(z\right)\sqrt{{z}^{2}-{a}^{2}},\end{equation} \tag{ 3.29 }$$

where P(z) is a polynomial in z. The degree of P(z) must be chosen such that it cancels all positive powers of R(z) to guarantee res_a. For V(z) = z²/2 we find that P(z) must be constant. Moreover, imposing the remaining conditions res_b,c on the resolvent one obtains

$$\begin{equation}\mathrm{R}\left(z\right)=\frac{1}{2}z-\frac{1}{2}\sqrt{{z}^{2}-4}\quad \text{and}\quad {\rho }_{\text{ext}}\left(x\right)=\frac{1}{2\pi }\sqrt{4-{x}^{2}},\quad x\in \left[-2,2\right].\end{equation} \tag{ 3.30 }$$

The eigenvalue distribution (3.30) is known as Wigner's semicircle law and appears in a wide range of physical and mathematical examples. Notice that it is connected and has compact support. Using the eigenvalue distribution we obtain the large N approximation

$$\begin{equation}{\mathcal{M}}_{N}\approx \mathrm{v}\mathrm{o}\mathrm{l}\enspace U\left(N\right)\enspace {\mathrm{e}}^{-\frac{3}{4}{N}^{2}}.\end{equation} \tag{ 3.31 }$$

Setting this expression equal to the exact result in (3.28), we can provide an approximation for the volume of the unitary group in the limit of large rank

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\enspace \mathrm{log}\enspace \mathrm{v}\mathrm{o}\mathrm{l}\enspace U\left(N\right)\approx \frac{{N}^{2}}{2}\left(\mathrm{log}\enspace \frac{2\pi }{N}+\frac{3}{2}\right)+\mathcal{O}\left(N\right).\end{equation} \tag{ 3.32 }$$

We can readily verify that the above expression agrees with the large N approximation of the exact expression

$$\begin{equation}\mathrm{v}\mathrm{o}\mathrm{l}\enspace U\left(N\right)=\frac{{\left(2\pi \right)}^{N\left(N+1\right)/2}}{G\left(N+1\right)},\end{equation} \tag{ 3.33 }$$

where G(N + 1) is Barnes G-function and we recall that at large N

$$\begin{equation}\mathrm{log}\enspace G\left(N+1\right)=-\frac{3{N}^{2}}{4}+\frac{{N}^{2}}{2}\enspace \mathrm{log}\enspace N+\frac{N}{2}\enspace \mathrm{log}\enspace 2\pi -\frac{1}{12}\enspace \mathrm{log}\enspace N+\cdots .\end{equation} \tag{ 3.34 }$$

Quartic example. Having done the Gaussian case, we now move on to a slightly more involved example. We will consider the integral

$$\begin{equation}{\mathcal{Z}}_{N}\left(\alpha \right)={\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]\enspace {\mathrm{e}}^{-N\enspace \mathrm{T}\mathrm{r}{V}_{\alpha }\left(M\right)},\quad {V}_{\alpha }\left(M\right)=\frac{1}{2}{M}^{2}+\alpha {M}^{4}.\end{equation} \tag{ 3.35 }$$

Using the ansatz (3.29) and following the approach delineated for the Gaussian case, we can write down the resolvent with P(z) a quadratic polynomial

$$\begin{equation}\mathrm{R}\left(z\right)=\frac{1}{2}z+2\alpha {z}^{3}-\left({p}_{1}+{p}_{2}z+{p}_{3}{z}^{2}\right)\sqrt{{z}^{2}-{a}^{2}}.\end{equation} \tag{ 3.36 }$$

The parameters p₁, p₂ and p₃ are fixed by res_a.

$$\begin{equation}{p}_{1}=\frac{1}{2}+\alpha {a}^{2},\qquad {p}_{2}=0,\qquad {p}_{3}=2\alpha .\end{equation} \tag{ 3.37 }$$

The parameter a solves the equation

$$\begin{equation}3\alpha \enspace {a}^{4}+{a}^{2}-4=0{\Rightarrow}{a}_{{\pm}}^{2}=-\frac{1}{6\alpha }{\pm}\frac{1}{6\alpha }\sqrt{1+48\alpha }.\end{equation} \tag{ 3.38 }$$

We are interested in the solution that is continuously connected to the Gaussian solution as α → 0, which is a₊. Notice that the saddle point solution (3.38) exhibits non-analytic behaviour as α approaches α_c = −1/48, somewhat reminiscent of what we saw in our vector integral example (2.20). If we wish to continue past α = α_c, the argument of the square root becomes negative and we have to specify a branch cut. Our expression for the resolvent and eigenvalue density are given by

$$\begin{equation}\mathrm{R}\left(z\right)=\frac{1}{2}z+2\alpha {z}^{3}-\left(\frac{1}{2}+\alpha \enspace {a}_{+}^{2}+2\alpha {z}^{2}\right)\sqrt{{z}^{2}-{a}_{+}^{2}},\end{equation} \tag{ 3.39 }$$

$$\begin{equation}{\rho }_{\text{ext}}\left(\lambda \right)=\frac{1}{\pi }\left(\frac{1}{2}+\alpha \enspace {a}_{+}^{2}+2\alpha {\lambda }^{2}\right)\sqrt{{a}_{+}^{2}-{\lambda }^{2}},\end{equation} \tag{ 3.40 }$$

with a₊ given by (3.38) (figure 3).

Zoom In Zoom Out Reset image size

To approximate our original integral (3.35) at large N, we must evaluate

$$\begin{equation}S\left[{\rho }_{\text{ext}}\left(\lambda \right)\right]={\int }_{0}^{{a}_{+}}\mathrm{d}\lambda \enspace {\rho }_{\text{ext}}\left(\lambda \right)\left(\frac{1}{2}{\lambda }^{2}+\alpha {\lambda }^{4}\right)-2{\int }_{0}^{{a}_{+}}\mathrm{d}\lambda \enspace {\rho }_{\text{ext}}\left(\lambda \right)\mathrm{log}\vert \lambda \vert \end{equation} \tag{ 3.41 }$$

on the saddle point solution ρ_ext(λ). In the above expression we also used that our integrands are symmetric around the origin to adjust the boundaries of the integrals. Evaluating this integral at the saddle ρ_ext(λ) we find

$$\begin{equation}S\left[{\rho }_{\text{ext}}\left(\lambda \right)\right]-S{\left[{\rho }_{\text{ext}}\left(\lambda \right)\right]}_{\alpha =0}=-\frac{3}{8}-\frac{1}{384}\left({a}_{+}^{2}-40\right){a}_{+}^{2}+\mathrm{log}\enspace \frac{2}{{a}_{+}},\end{equation} \tag{ 3.42 }$$

where we used (3.38) to write the above as a function of a₊. Substituting a₊, a small α expansion reveals the following [39–41]

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\enspace \frac{{\mathcal{F}}_{N}\left(\alpha \right)}{{N}^{2}}=2\alpha -18{\alpha }^{2}+288{\alpha }^{3}-6048{\alpha }^{4}+\frac{746\enspace 496}{5}{\alpha }^{5}+\mathcal{O}\left({\alpha }^{6}\right),\end{equation} \tag{ 3.43 }$$

where we have defined ${\mathcal{F}}_{N}\left(\alpha \right)\equiv -\mathrm{log}\enspace {\mathcal{Z}}_{N}\left(\alpha \right)/{\mathcal{Z}}_{N}\left(0\right)$, analogously to the vector integral case. More generally, and somewhat similar to the vector integral case (2.24), the matrix integrals allow for an expansion in large N²

$$\begin{equation}{\mathcal{F}}_{N}\left(\alpha \right)={N}^{2}{\mathcal{F}}^{\left(0\right)}\left(\alpha \right)+{\mathcal{F}}^{\left(1\right)}\left(\alpha \right)+\frac{1}{{N}^{2}}{\mathcal{F}}^{\left(2\right)}\left(\alpha \right)+\cdots \enspace ,\end{equation} \tag{ 3.44 }$$

where each of the ${\mathcal{F}}^{\left(n\right)}\left(\alpha \right)$ itself admits a power series expansion around α = 0. A careful examination of (3.42) and (3.43) reveals [39]

$$\begin{equation}{\mathcal{F}}^{\left(0\right)}\left(\alpha \right)=-\sum _{n=1}^{\infty }{\left(-12\alpha \right)}^{n}\frac{\left(2n-1\right)!}{n!\left(n+2\right)!}\enspace =2{\alpha \enspace }_{3}{F}_{2}\left(1,1,\frac{3}{2};2,4;-48\alpha \right).\end{equation} \tag{ 3.45 }$$

We can expand the above near the critical value α_c = −1/48 to find the leading non-analytic behaviour (which we denote with a subscript n.a.)

$$\begin{equation}{\mathrm{lim}}_{\alpha \to {\alpha }_{\mathrm{c}}}\enspace {\partial }_{\alpha }^{\left(3\right)}\enspace {\mathcal{F}}_{\text{n.a.}}^{\left(0\right)}\left(\alpha \right)=4608\sqrt{3}{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{-1/2}+\cdots .\end{equation} \tag{ 3.46 }$$

We will analyse the terms ${\mathcal{F}}^{\left(h\right)}\left(\alpha \right)$ in (3.44) in the next section. In particular we will provide detailed evidence that near the critical value α_c the leading non-analytic behaviour goes as

$$\begin{equation}{\mathrm{lim}}_{\alpha \to {\alpha }_{\mathrm{c}}}\enspace {\mathcal{F}}_{\text{n.a.}}^{\left(h\right)}\left(\alpha \right)\sim {\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{5{\chi }_{h}/4},\quad h\in \mathbb{N},\end{equation} \tag{ 3.47 }$$

where a logarithmic behaviour is understood for h = 1 and χ_h = 2 − 2h.

The planar solution for the eigenvalue distribution stemming from a general V(λ) can be written down rather concisely. We discuss here the case for which the eigenvalue distribution is connected and has compact support on a single real interval λ ∈ [b, a]. In that case the resolvent is given by [48]

$$\begin{equation}\text{R}\left(z\right)={\frac{1}{2}\oint }_{\mathcal{C}}\frac{\mathrm{d}u}{2\pi \mathrm{i}}\frac{{V}^{\prime }\left(u\right)}{z-u}\sqrt{\frac{\left(z-a\right)\left(z-b\right)}{\left(u-a\right)\left(u-b\right)}},\end{equation} \tag{ 3.48 }$$

where the contour $\mathcal{C}$ goes around the branch cut z ∈ [b, a]. The end points a and b follow from the conditions (3.22) and (3.24) imposed on the resolvent. Given R(z), we can extract the eigenvalue density from (3.23).

It is interesting to note that there exist a variety of polynomials which upon tuning the coefficients give rise to non-analytic behaviour different from (3.47). These matrix integrals are known as multicritical matrix integrals and were originally introduced in [43]. For V(M) an even polynomial of degree 2m we denote the matrix integral as an mth multicritical model. For instance, using the above techniques, the following polynomial [43, 45]

$$\begin{equation}V\left(M\right)=\frac{1}{\gamma }\left(\frac{1}{2}{M}^{2}-\frac{1}{3}{M}^{4}+\frac{4}{45}{M}^{6}\right),\quad \gamma =\frac{1}{12}-\frac{1}{12}{\left(1-{a}^{2}\right)}^{3},\end{equation} \tag{ 3.49 }$$

with the eigenvalues of M distributed in the range [−a, a] gives rise to a non-analytic behaviour of the form

$$\begin{equation}{\mathrm{lim}}_{\gamma \to {\gamma }_{\mathrm{c}}}\enspace {\partial }_{\gamma }^{\left(3\right)}{\mathcal{F}}_{\text{n.a.}}^{\left(0\right)}\left(\gamma \right)\sim {\left(\gamma -{\gamma }_{\mathrm{c}}\right)}^{-2/3},\end{equation} \tag{ 3.50 }$$

with γ_c = 1/12. Note that we have determined γ in (3.49) by combining (3.48) with res_a (3.22). We note that the branch cut of the density is still determined by a square root (figure 4):

$$\begin{equation}{\rho }_{\text{ext}}\left(\lambda \right)=\frac{1}{\gamma \enspace \pi }\left(\frac{1}{2}-\frac{1}{3}{a}^{2}+\frac{1}{10}{a}^{4}+\frac{2}{15}\left({a}^{2}-5\right){\lambda }^{2}+\frac{4}{15}{\lambda }^{4}\right)\sqrt{{a}^{2}-{\lambda }^{2}},\quad \lambda \in \left[-a,a\right].\end{equation} \tag{ 3.51 }$$

Zoom In Zoom Out Reset image size

A natural set of integrals living in the same class as (3.1) are the following

$$\begin{equation}{\mathcal{M}}_{N}^{\left({k}_{1},{k}_{2},\dots ,{k}_{n}\right)}\equiv \frac{1}{{\mathcal{M}}_{N}}\enspace {\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]{\mathrm{e}}^{-N\mathrm{T}\mathrm{r}V\left(M\right)}\prod _{i=1}^{n}\frac{1}{N}\enspace \mathrm{T}\mathrm{r}\enspace {M}^{{k}_{i}}=\frac{1}{{N}^{n}}\left\langle \prod _{i=1}^{n}\mathrm{T}\mathrm{r}\enspace {M}^{{k}_{i}}\right\rangle .\end{equation} \tag{ 3.52 }$$

These integrals also preserve the U(N) symmetry. To leading order in the large N limit we can use the eigenvalue density to estimate the ${\mathcal{M}}_{N}^{\left({k}_{1},{k}_{2},\dots ,{k}_{n}\right)}$. For example

$$\begin{equation}{\mathcal{M}}_{N}^{\left({k}_{1},{k}_{2},\dots ,{k}_{n}\right)}=\prod _{i=1}^{n}{\int }_{-a}^{a}\mathrm{d}\lambda \enspace {\rho }_{\text{ext}}\left(\lambda \right)\enspace {\lambda }^{{k}_{i}}+\mathcal{O}\left(1/{N}^{2}\right).\end{equation} \tag{ 3.53 }$$

So long as neither n nor k_i scale with N we see that to leading order in the large N limit we cannot distinguish ${\mathcal{M}}_{N}^{\left({k}_{1},\dots ,{k}_{n}\right)}$ from ${\mathcal{M}}_{N}^{\left({k}_{1}\right)}{\mathcal{M}}_{N}^{\left({k}_{2}\right)}\dots {\mathcal{M}}_{N}^{\left({k}_{n}\right)}$. This is a characteristic property of large N systems known as large N factorisation [42]. From this perspective, Tr M^k/N can be viewed as a collection of weakly correlated quantities, the correlation strength going as some inverse power of N. The phenomenon of large N factorisation thus provides us with a novel type of perturbative expansion for a certain class of matrix functions.

As an example, let us use large N factorisation to recover the equations governing the resolvent at large N. The starting point is the observation that the integral ${\mathcal{M}}_{N}$ should be left unchanged under a change in variables M → F(M) where F is a matrix-valued function. For instance, we can consider the family of functions F_ℓ(M) = M + g e^ℓM where g is taken to be parametrically small and ℓ is a real number. Invariance of ${\mathcal{M}}_{N}$ under these transformations imposes a set of interesting constraints [48]. For parametrically small g one finds

$$\begin{equation}\frac{\ell }{{N}^{2}}\sum _{k=0}^{\infty }\frac{{\ell }^{k}}{k!}\sum _{j=0}^{k}\left\langle \mathrm{T}\mathrm{r}\enspace {M}^{j}\enspace \mathrm{T}\mathrm{r}\enspace {M}^{k-j}\right\rangle ={\hat{G}}_{\ell }^{\left(V\right)}\enspace \frac{1}{N}\left\langle \mathrm{T}\mathrm{r}\enspace {\mathrm{e}}^{\ell M}\enspace \right\rangle .\end{equation} \tag{ 3.54 }$$

The left-hand side stems from the variation of the measure, whereas the right hand side stems from the variation of the exponent. ${\hat{G}}_{\ell }^{\left(V\right)}$ is a differential operator acting on functions of ℓ. If V(M) admits a power series, $V\left(M\right)={\sum }_{n=0}^{\infty }{\alpha }_{n}{M}^{n}$, we have

$$\begin{equation}{\hat{G}}_{\ell }^{\left(V\right)}=\sum _{n=1}^{\infty }n\enspace {\alpha }_{n}\enspace {\partial }_{\ell }^{\left(n-1\right)}.\end{equation} \tag{ 3.55 }$$

To leading order in the large N expansion, we can use large N factorisation to break the left-hand side of (3.54) into a sum over a product of traces. The leading large N expression is then given by

$$\begin{equation}{\int }_{0}^{\ell }\mathrm{d}u\enspace \left\langle {W}_{u}\right\rangle \left\langle {W}_{\ell -u}\right\rangle ={\hat{G}}_{\ell }^{\left(V\right)}\left\langle {W}_{\ell }\right\rangle .\end{equation} \tag{ 3.56 }$$

The above equation is a type of equation known as a loop equation [46–49, 172], and the object

$$\begin{equation}{W}_{\ell }\equiv \frac{1}{N}\enspace \mathrm{T}\mathrm{r}\enspace {\mathrm{e}}^{\ell M}\end{equation} \tag{ 3.57 }$$

is known as a macroscopic loop operator. The resolvent R_N(z) defined in (3.17) can be obtained from the loop operator by the following Laplace transform

$$\begin{equation}{\text{R}}_{N}\left(z\right)={\int }_{0}^{\infty }\mathrm{d}\ell \enspace {\mathrm{e}}^{-\ell z}\enspace {W}_{\ell }\enspace {\Rightarrow}\enspace \text{R}\left(z\right)={\mathrm{lim}}_{N\to \infty }\enspace {\int }_{0}^{\infty }\mathrm{d}\ell \enspace {\mathrm{e}}^{-\ell z}\left\langle {W}_{\ell }\right\rangle .\end{equation} \tag{ 3.58 }$$

We can take the Laplace transform of (3.56) with respect to ℓ and find

$$\begin{equation}\text{R}{\left(z\right)}^{2}={V}^{\prime }\left(z\right)\mathrm{R}\left(z\right)-\mathcal{P}\left(z\right),\qquad \mathcal{P}\left(z\right)=-{\int }_{0}^{\infty }\mathrm{d}\ell \enspace {\hat{G}}_{\ell }^{\left(V\right)}\left\langle {W}_{\ell }\enspace {\mathrm{e}}^{-\ell z}\right\rangle .\end{equation} \tag{ 3.59 }$$

Solving the quadratic equation leads to

$$\begin{equation}\mathrm{R}\left(z\right)=\frac{1}{2}\enspace {V}^{\prime }\left(z\right)-\frac{1}{2}\sqrt{{V}^{\prime }{\left(z\right)}^{2}-4\mathcal{P}\left(z\right)},\end{equation} \tag{ 3.60 }$$

where we choose the negative root to satisfy res_a (3.22). Equation (3.60) constitutes a derivation of (3.48).

Examples. We present the loop operators for the three polynomials discussed in the beginning of this section. For the Gaussian case we can use the density in (3.30) to calculate the loop operator

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\left\langle {W}_{\ell }\right\rangle =\frac{1}{\ell }\enspace {I}_{1}\left(2\ell \right),\end{equation} \tag{ 3.61 }$$

where I_n(z) is the modified Bessel function of the first kind. The differential operator is ${\hat{G}}_{\ell }^{\left(V\right)}={\partial }_{\ell }$ and it confirms (3.56). In the large z limit we obtain

$$\begin{equation}\mathcal{P}\left(z\right)=-{\int }_{0}^{\infty }\mathrm{d}\ell \enspace {\partial }_{\ell }\left(\frac{1}{\ell }\enspace {I}_{1}\left(2\ell \right)\enspace {\mathrm{e}}^{-\ell z}\right)=1,\end{equation} \tag{ 3.62 }$$

such that (3.60) agrees perfectly with the resolvent in (3.30). For the quartic polynomial we find ${\hat{G}}_{\ell }^{\left({V}_{\alpha }\right)}={\partial }_{\ell }+4\alpha {\partial }_{\ell }^{3}$ and using the density in (3.39) we obtain

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\left\langle {W}_{\ell }\right\rangle =\frac{{a}_{+}\left(1+6\enspace {a}_{+}^{2}\alpha \right)}{2\ell }\enspace {I}_{1}\left({a}_{+}\enspace \ell \right)-\frac{6\enspace {a}_{+}^{2}\alpha }{{\ell }^{2}}\enspace {I}_{2}\left({a}_{+}\enspace \ell \right),\end{equation} \tag{ 3.63 }$$

where a₊ has been defined in (3.38). Combining W_ℓ with ${\hat{G}}_{\ell }^{\left({V}_{\alpha }\right)}$ we easily confirm (3.56) and combining it with (3.60) we confirm the polynomial in (3.39). Finally, for the multicritical polynomial introduced in (3.49) we obtain

$$\begin{equation}{\mathrm{lim}}_{N\to \infty }\left\langle {W}_{\ell }\right\rangle =\frac{a}{2\gamma {\ell }^{4}}\left(8{a}^{2}\ell +{\left(1-{a}^{2}\right)}^{2}{\ell }^{3}\right){I}_{1}\left(a\ell \right)-\frac{2{a}^{2}}{\gamma {\ell }^{4}}\left(8-\left(1-{a}^{2}\right){\ell }^{2}\right){I}_{2}\left(a\ell \right),\end{equation} \tag{ 3.64 }$$

and the differential operator ${\hat{G}}_{\ell }^{\left(V\right)}$ is given by

$$\begin{equation}{\hat{G}}_{\ell }^{\left(V\right)}=\frac{1}{\gamma }\left({\partial }_{\ell }-\frac{4}{3}\enspace {\partial }_{\ell }^{3}+\frac{4}{45}\enspace {\partial }_{\ell }^{5}\right),\end{equation} \tag{ 3.65 }$$

with γ defined in (3.49). Using (3.64), it is straightforward to verify (3.56).

We now proceed to study how the perturbative expansion of a matrix type integral is organised in the large N limit. We will uncover a remarkable connection, originally observed by 't Hooft [7], to Riemann surfaces.

Let us focus on matrix integrals with the following structure

$$\begin{equation}{\mathcal{M}}_{N}={\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]\enspace {\mathrm{e}}^{-N\enspace \mathrm{T}\mathrm{r}\enspace V\left(M\right)},\end{equation} \tag{ 3.66 }$$

where V(M) is an arbitrary polynomial in M containing parameters α which will be taken to be small. Though the lessons will be general, for the sake of concreteness we will consider again a purely quartic example V_α(M) (3.35). Since we will be interested in a perturbative expansion in small α, it is convenient to write down the propagator

$$\begin{equation}\langle {M}_{IK}{M}_{JL}\rangle ={\mathcal{Z}}_{N}^{-1}\left(0\right)\enspace {\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]\enspace {\mathrm{e}}^{-\frac{N}{2}\mathrm{T}\mathrm{r}\enspace {M}^{2}}{M}_{IJ}{M}_{KL}=\frac{1}{N}{\delta }_{IL}{\delta }_{KJ}.\end{equation} \tag{ 3.67 }$$

For small coupling α we expand the exponential (3.35) leading to

$$\begin{equation}{\mathcal{Z}}_{N}\left(\alpha \right)={\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]\enspace {\mathrm{e}}^{-\frac{N}{2}\mathrm{T}\mathrm{r}\enspace {M}^{2}}\sum _{k\in \mathbb{N}}{\left(-1\right)}^{k}\frac{{\left(N\alpha \right)}^{k}}{k!}\enspace {\left(\mathrm{T}\mathrm{r}\enspace {M}^{4}\right)}^{k}.\end{equation} \tag{ 3.68 }$$

The graphical representation of the propagator and the quartic vertex is given in figure 5:

Zoom In Zoom Out Reset image size

Using standard Gaussian integration, we can compute the various terms in the small α expansion. For example, to linear order in α the large N integral (3.35) becomes

$$\begin{equation}{\mathcal{F}}_{N}\left(\alpha \right)\equiv -\mathrm{log}\enspace \frac{{\mathcal{Z}}_{N}\left(\alpha \right)}{{\mathcal{Z}}_{N}\left(0\right)}=\left(2{N}^{2}+1\right)\alpha +\mathcal{O}\left({\alpha }^{2}\right).\end{equation} \tag{ 3.69 }$$

Each of the terms in the α expansion has a diagrammatic representation. We give some examples in the figure above. Due to their double line representation, these diagrams are often called ribbon diagrams (figure 6).

Zoom In Zoom Out Reset image size

From the small α expansion it is clear that the diagrams contributing to the perturbative expansion of a matrix integral can scale with different powers of N. In (3.67) we observed that the propagator scales with inverse power of N. Moreover, each loop contributes a factor of N, and each vertex adds a factor αN. This implies that a diagram with L loops, V vertices, and P propagators scales as α^VN^V+L−P. Making the identification [7]

$$\begin{equation}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{p}\mathrm{s}\enspace L\enspace \hat{=}\enspace \mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}\mathrm{s}\enspace F,\qquad \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{a}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{s}\enspace P\enspace \hat{=}\enspace \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}\mathrm{s}\enspace E,\qquad \mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}\mathrm{s}\enspace V\enspace \hat{=}\enspace \mathrm{v}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}\mathrm{s}\enspace V,\end{equation} \tag{ 3.70 }$$

one can identify the power of N associated to a particular ribbon diagram with the Euler characteristic χ_h of the two-dimensional compact surface Σ_h of genus h it can be drawn on. Since χ_h = 2 − 2h, at least pictorially, there seems to be a natural genus expansion in the large N limit for the perturbative expansion of our matrix integral ${\mathcal{Z}}_{N}\left(\alpha \right)$. We leave it to the reader to verify that this conclusion will not be affected by allowing V_α(M) to be an arbitrary polynomial, so long as the couplings are scaled appropriately as we take N large.

To recapitulate: the large N organisation of our perturbative expansion can be expressed as

$$\begin{equation}{\mathcal{F}}_{N}\left(\alpha \right)=\sum _{h=0}^{\infty }{\mathrm{e}}^{{\chi }_{h}\mathrm{log}N}\enspace {\mathcal{F}}^{\left(h\right)}\left(\alpha \right),\end{equation} \tag{ 3.71 }$$

where h labels the genus of the Riemann surface Σ_h, χ_h its Euler characteristic, and each ${\mathcal{F}}_{h}\left(\alpha \right)$ is a sum of connected diagrams that can be drawn on a surface of genus h

$$\begin{equation}{\mathcal{F}}^{\left(h\right)}\left(\alpha \right)=\sum _{V=0}^{\infty }{f}_{V}^{\left(h\right)}{\alpha }^{V}.\end{equation} \tag{ 3.72 }$$

The leading contribution in the large N expansion is known as the planar contribution, and sub-leading contributions in the large N expansion are known as non-planar contributions. We argued that there exists a natural map between fat graphs and the polygonisation of Σ_h. For the quartic theory, the surfaces of the graph dual to the polygonisation of Σ_h are squares (figure 7).

Zoom In Zoom Out Reset image size

3.5.1. Riemann surfaces with boundaries

In addition to ribbon diagrams corresponding to discretised Riemann surfaces without boundary, we can also examine diagrams living on Riemann surfaces with boundaries. We start by coupling N dimensional complex Grassmann valued vectors ${\overline{\xi }}^{\left(s\right)}$ and ξ^(s) to the matrix integral (3.1)

$$\begin{equation}{\mathcal{K}}_{N}\left(\alpha ;z\right)\equiv \frac{1}{{\mathcal{M}}_{N}\left(\alpha \right)}{\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]\enspace \int \prod _{s=1}^{{N}_{f}}\prod _{I=1}^{N}\mathrm{d}{\overline{\xi }}_{I}^{\left(s\right)}\mathrm{d}{\xi }_{I}^{\left(s\right)}\enspace {\mathrm{e}}^{-N\mathrm{T}\mathrm{r}\enspace {V}_{\alpha }\left(M\right)-{\overline{\xi }}^{\left(s\right)}\left(z\enspace {\mathbb{I}}_{N}-M\right){\xi }^{\left(s\right)}},\quad z\in \mathbb{C}.\end{equation} \tag{ 3.73 }$$

The index s can be viewed as a flavour index ranging over N_f 'flavours'. For concreteness, although not necessary, we choose the quartic polynomial (3.35). The ribbon diagrams in figure 5 get enhanced by graphs with a single line (figure 8).

Zoom In Zoom Out Reset image size

Integrating out ${\overline{\xi }}^{\left(s\right)}$ and ξ^(s) in (3.73) we obtain

$$\begin{equation}{\mathcal{K}}_{N}\left(\alpha ;z\right)=\frac{1}{{\mathcal{M}}_{N}\left(\alpha \right)}{\int }_{{\mathbb{R}}^{{N}^{2}}}\left[\mathrm{d}M\right]{\text{e}}^{-N\enspace \mathrm{T}\mathrm{r}{V}_{\alpha }\left(M\right)+{N}_{f}\enspace \mathcal{W}\left(z\right)}=1+\enspace \left\langle {N}_{f}\enspace \mathcal{W}\left(z\right)\right\rangle +\frac{1}{2!}\left\langle {N}_{f}^{2}\enspace \mathcal{W}{\left(z\right)}^{2}\right\rangle \enspace +\cdots \enspace ,\end{equation} \tag{ 3.74 }$$

where we defined

$$\begin{equation}\mathcal{W}\left(z\right)\equiv \mathrm{T}\mathrm{r}\enspace \mathrm{log}\left(z\enspace {\mathbb{I}}_{N}-M\right).\end{equation} \tag{ 3.75 }$$

It is worth noting that ${\partial }_{z}\mathcal{W}\left(z\right)$ is equal to the resolvent introduced in (3.17). Further defining ${\mathcal{B}}_{N}\left(\alpha ;z\right)\equiv \mathrm{log}\enspace {\mathcal{K}}_{N}\left(\alpha ;z\right)$ we obtain the large N expansion

$$\begin{equation}{\mathcal{B}}_{N}\left(\alpha ;z\right)=\sum _{h=0}^{\infty }\sum _{b=1}^{\infty }\enspace {\mathrm{e}}^{{\chi }_{h,b}\mathrm{log}N}{N}_{f}^{b}\enspace {\mathcal{B}}^{\left(h,b\right)}\left(\alpha ;z\right),\quad {\chi }_{h,b}\equiv 2-2h-b.\end{equation} \tag{ 3.76 }$$

Here, ${\mathcal{B}}_{N}\left(\alpha ;z\right)$ encodes the sum of connected diagrams corresponding to discretised Riemann surfaces with h holes and b boundaries. In the figure below we display an example of a ribbon diagram and its dual lattice for h = 0 and b = 1 (figure 9).

Zoom In Zoom Out Reset image size

Each term ${\mathcal{B}}^{\left(h,b\right)}\left(\alpha ;z\right)$ in (3.76) is itself a sum of connected diagrams with fixed h and b:

$$\begin{equation}{\mathcal{B}}^{\left(h,b\right)}\left(\alpha ;z\right)=\sum _{V=0}^{\infty }\sum _{B=1}^{\infty }{f}_{V,B}^{\left(h,b\right)}{\alpha }^{V}{z}^{B}.\end{equation} \tag{ 3.77 }$$

The parameter z is the boundary analog of α, while N_f is the boundary analog of N. Keeping track of the powers of N_f in (3.76) allows one to distinguish Riemann surfaces which would be indistinguishable purely from their power of N (figure 10).

Zoom In Zoom Out Reset image size

3.5.2. Continuum limit and double scaling

In analogy to the bubble number we defined in (2.26) for cactus diagrams, we can define the 'vertex number' for the polygonisation of Σ_h as

$$\begin{equation}\langle {n}_{h}\rangle ={\partial }_{\mathrm{log}\alpha }\mathrm{log}\enspace {\mathcal{F}}^{\left(h\right)}\left(\alpha \right).\end{equation} \tag{ 3.78 }$$

For instance, given (3.45) we can compute ⟨n₀⟩ near α_c to find $\langle {n}_{0}\rangle \sim {\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{-1}$. More generally, granting (3.47), we find $\langle {n}_{h}\rangle \sim {\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{-1}$ for all h. Consequently, as α approaches α_c the number of vertices diverges suggesting the possibility for a continuum limit of the discretised surfaces Σ_h.

There is another limit of interest we can consider, known as the double scaling limit. The limit consists of simultaneously taking N large as well as α to it is critical value α_c, while keeping a particular combination α and N fixed. Recall (3.47), namely that the non-analytic part of ${\mathcal{F}}^{\left(h\right)}\left(\alpha \right)$ near α_c scales as

$$\begin{equation}{\mathrm{lim}}_{\alpha \to {\alpha }_{\mathrm{c}}}\enspace {\mathcal{F}}_{\text{n.a.}}^{\left(h\right)}\left(\alpha \right)={f}_{h}{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{5{\chi }_{h}/4},\quad h\in \mathbb{N},\end{equation} \tag{ 3.79 }$$

with f_h being some proportionality constant and h = 1 is understood to be logarithmic. It follows that the large N expansion near α = α_c is approximately

$$\begin{equation}{\mathcal{F}}_{N}\left(\alpha \right)\approx {N}^{2}{\mathcal{F}}^{\left(0\right)}\left({\alpha }_{\mathrm{c}}\right)+\sum _{h=0}^{\infty }{f}_{h}\enspace {\mathrm{e}}^{{\chi }_{h}\mathrm{log}N}{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{5{\chi }_{h}/4}.\end{equation} \tag{ 3.80 }$$

The above expression suggests the introduction of a new parameter

$$\begin{equation}{\kappa }^{-1}\equiv N{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{5/4}.\end{equation} \tag{ 3.81 }$$

The double scaling limit takes N → ∞ and α → α_c while keeping κ fixed. In this way, we see how the continuum limit and the large N limit work harmoniously in producing a new function of potential interest

$$\begin{equation}\mathcal{F}\left(\kappa \right)\equiv {\mathrm{lim}}_{N\to \infty ,\enspace \enspace \kappa \enspace \enspace \text{fixed}}\left({\mathcal{F}}_{N}\left(\alpha \right)-{N}^{2}{\mathcal{F}}^{\left(0\right)}\left(\alpha \right)\right)=\sum _{h=0}^{\infty }{f}_{h}\enspace {\kappa }^{-{\chi }_{h}}.\end{equation} \tag{ 3.82 }$$

In the next subsection, by introducing a new method for dealing with the matrix integrals, we will produce a non-linear differential equation whose solution encodes $\mathcal{F}\left(\kappa \right)$.

An interesting system which naturally exhibits a genus expansion is the perturbative expansion of worldsheet string theory. At this stage there is no immediate reason for the two systems to be related, but one may speculate so. We shall see in section 8 that such speculations are materialised in a concrete and elegant sense.

Two polynomials are said to be orthogonal with respect to a weight function w(x) if they satisfy

$$\begin{equation}{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}_{a}:\int \mathrm{d}x\enspace w\left(x\right)\enspace {p}_{n}\left(x\right){p}_{m}\left(x\right)={h}_{n}{\delta }_{m,n}\end{equation} \tag{ 4.1 }$$

for some h_n. A familiar example is given by the Hermite polynomials

$$\begin{equation}{H}_{n}\left(x\right)={\left(-1\right)}^{n}\enspace {\mathrm{e}}^{{x}^{2}}\frac{{\mathrm{d}}^{n}}{\mathrm{d}{x}^{n}}\enspace {\mathrm{e}}^{-{x}^{2}},\quad n=1,2,\dots \end{equation} \tag{ 4.2 }$$

which satisfy

$$\begin{equation}{\int }_{\mathbb{R}}\mathrm{d}x\enspace {\mathrm{e}}^{-{x}^{2}}\enspace {H}_{n}\left(x\right){H}_{m}\left(x\right)={2}^{n}\enspace \sqrt{\pi }\enspace n!\enspace {\delta }_{mn}.\end{equation} \tag{ 4.3 }$$

In addition to (4.1), orthogonal polynomials satisfy the three-term recurrence relation

$$\begin{align}\hfill {\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}_{b}:& \enspace x\enspace {p}_{n}\left(x\right)={A}_{n}\enspace {p}_{n}\left(x\right)+{S}_{n}\enspace {p}_{n+1}\left(x\right)+{R}_{n}\enspace {p}_{n-1}\left(x\right)\quad \text{for}\enspace n{ >}0,\hfill \\ \hfill & x\enspace {p}_{0}\left(x\right)={A}_{0}\enspace {p}_{0}\left(x\right)+{S}_{0}\enspace {p}_{1}\left(x\right),\hfill \end{align} \tag{ 4.4 }$$

where A_n, S_n, and R_n are some real constants. The Hermite polynomials (4.2) satisfy

$$\begin{equation}x\enspace {H}_{n}=\frac{1}{2}\enspace {H}_{n+1}+n\enspace {H}_{n-1},\end{equation} \tag{ 4.5 }$$

and thus for n ⩾ 0 we have A_n = 0, S_n = 1/2, and R_n = n. Other examples of orthogonal polynomials include the Legendre polynomials, Laguerre polynomials, and Chebychev polynomials.

In what follows we will focus on monic polynomials

$$\begin{equation}{P}_{n}\left(\lambda \right)\equiv {\lambda }^{n}+\sum _{j=0}^{n-1}{a}_{j}{\lambda }^{j},\quad n=0,\dots ,N-1.\end{equation} \tag{ 4.6 }$$

It is worth mentioning that when the measure in (4.1) is even under λ → −λ, the monic polynomials P_n(λ) transform as P_n(−λ) = (−1)ⁿP_n(λ). This implies that only even powers of λ appear in P_2n(λ) whereas only odd powers of λ appear in P_2n−1(λ). Moreover, invariance under λ → −λ implies that the coefficient A_n in the three-term recurrence relation vanishes, while being monic implies S_n = 1.

4.1.1. Orthogonal polynomials for matrix integrals

We will now use the properties of monic orthogonal polynomials to solve matrix integrals. In order to do so, we observe that the Vandermonde matrix can be expressed in the following way

$$\begin{equation}{\mathbf{V}}_{N}=\left(\begin{matrix}\hfill {P}_{0}\left({\lambda }_{1}\right)\hfill & \hfill {P}_{1}\left({\lambda }_{1}\right)\hfill & \hfill \dots \hfill & \hfill {P}_{N-1}\left({\lambda }_{1}\right)\hfill \\ \hfill {P}_{0}\left({\lambda }_{2}\right)\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {\vdots}\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {P}_{0}\left({\lambda }_{N}\right)\hfill & \hfill {P}_{1}\left({\lambda }_{N}\right)\hfill & \hfill \dots \hfill & \hfill {P}_{N-1}\left({\lambda }_{N}\right)\hfill \end{matrix}\right).\end{equation} \tag{ 4.7 }$$

The above expression is related to (3.10) by a similarity transformation. The determinant Δ_N(λ) depends only on the leading degree of each polynomial. Using the Leibniz formula for determinants, we can recast Δ_N(λ) as

$$\begin{equation}{{\Delta}}_{N}\left(\lambda \right)=\sum _{\sigma \in {S}_{N}}\mathrm{sgn}\left(\sigma \right)\prod _{I=1}^{N}{P}_{\sigma \left(I\right)-1}\left({\lambda }_{I}\right).\end{equation} \tag{ 4.8 }$$

Consequently, we can re-express the matrix integral (3.1) as

$$\begin{equation}\text{vol}\frac{U{\left(1\right)}^{N}{\times}{S}_{N}}{U\left(N\right)}\enspace {\mathcal{M}}_{N}=\sum _{\sigma ,\tau \in {S}_{N}}\mathrm{sgn}\left(\sigma \right)\enspace \mathrm{sgn}\left(\tau \right)\prod _{I=1}^{N}{\int }_{\mathbb{R}}\mathrm{d}\mu \left({\lambda }_{I}\right){P}_{\sigma \left(I\right)-1}\left({\lambda }_{I}\right){P}_{\tau \left(I\right)-1}\left({\lambda }_{I}\right)=N!\prod _{I=0}^{N-1}{h}_{I},\end{equation} \tag{ 4.9 }$$

where for notational convenience we have defined dμ(λ) ≡ dλ e^−NV(λ). Using the three-term recurrence relation we can relate h_n to R_n in (4.4) and rewrite (4.9) entirely in terms of R_n, which are themselves determined by the weight function w(λ) = e^−NV(λ). It will also prove convenient to express h_n in the following way

$$\begin{equation}{h}_{n}={\int }_{\mathbb{R}}\mathrm{d}\mu \left(\lambda \right)\lambda {P}_{n-1}\left(\lambda \right){P}_{n}\left(\lambda \right).\end{equation} \tag{ 4.10 }$$

An additional relation that will be useful is

$$\begin{equation}n{h}_{n}={\int }_{\mathbb{R}}\mathrm{d}\mu \left(\lambda \right){P}_{n}^{\prime }\left(\lambda \right)\enspace \lambda {P}_{n}\left(\lambda \right)=N{R}_{n}{\int }_{\mathbb{R}}\mathrm{d}\mu \left(\lambda \right){V}^{\prime }\left(\lambda \right){P}_{n}\left(\lambda \right){P}_{n-1}\left(\lambda \right).\end{equation} \tag{ 4.11 }$$

For example, when combined with (4.10) the above relation informs us that for a Gaussian potential V(λ) = λ²/2 we have ${h}_{0}=\sqrt{2\pi /N}$ and R_n = n/N. Hence the orthogonal polynomials are the Hermite polynomials for V(λ) = λ²/2. From the above expressions, we can recursively calculate the sub-leading coefficients a_J in (4.6). Applying ortho_b to (4.10) leads to the recursion h_n = R_nh_n−1, which implies that (4.9) can be rewritten as

$$\begin{equation}{\mathcal{M}}_{N}=\text{vol}\;\frac{U\left(N\right)}{U{\left(1\right)}^{N}{\times}{S}_{N}}{\times}N!\enspace {h}_{0}^{N}{R}_{1}^{N-1}{R}_{2}^{N-2}\dots {R}_{N-1}.\end{equation} \tag{ 4.12 }$$

As a simple check of the above expression we note that for α = 0

$$\begin{equation}{\mathcal{M}}_{N}=\text{vol}\;\frac{U\left(N\right)}{U{\left(1\right)}^{N}{\times}{S}_{N}}{\times}N!{\left(\frac{2\pi }{N}\right)}^{N/2}\prod _{n=1}^{N-1}{\left(\frac{n}{N}\right)}^{N-n}={\left(\frac{2\pi }{N}\right)}^{{N}^{2}/2},\end{equation} \tag{ 4.13 }$$

where we have used (3.33), vol S_N = N!, and vol U(1) = 2π.

Finally, we obtain a relation between the quantities defined and ${\mathcal{F}}_{N}\left(\alpha \right)$:

$$\begin{equation}\frac{1}{{N}^{2}}{\mathcal{F}}_{N}\left(\alpha \right)=-\frac{1}{N}\enspace \mathrm{log}\enspace \frac{{h}_{0}\left(\alpha \right)}{{h}_{0}\left(0\right)}-\frac{1}{N}\sum _{n=1}^{N-1}\left(1-\frac{n}{N}\right)\mathrm{log}\enspace \frac{{R}_{n}\left(\alpha \right)}{{R}_{n}\left(0\right)}.\end{equation} \tag{ 4.14 }$$

The argument of h_n and R_n indicates their dependence on the coefficient α parametrising a non-Gaussian piece of V(λ). We now delve into a detailed example.

We are now ready to see how the orthogonal polynomials can be used to go beyond the planar approximation of the large N expansion of ${\mathcal{F}}_{N}\left(\alpha \right)$. For concreteness we take V_α(λ) = λ²/2 + αλ⁴. From (4.11) we find

$$\begin{equation}\frac{n}{N}={R}_{n}\left(\alpha \right)\left(1+4\alpha \left({R}_{n+1}\left(\alpha \right)+{R}_{n}\left(\alpha \right)+{R}_{n-1}\left(\alpha \right)\right)\right),\end{equation} \tag{ 4.15 }$$

where the term multiplying 4α follows from the quartic interaction term in the potential and the relation

$$\begin{equation}{\int }_{\mathbb{R}}\mathrm{d}\mu \left(\lambda \right){\lambda }^{3}{P}_{n-1}\left(\lambda \right){P}_{n}\left(\lambda \right)=\left({R}_{n+1}\left(\alpha \right)+{R}_{n}\left(\alpha \right)+{R}_{n-1}\left(\alpha \right)\right){h}_{n}\left(\alpha \right),\end{equation} \tag{ 4.16 }$$

which is obtained using the three-term-recurrence relation (4.4). In appendix A we provide a graphical procedure to obtain relations similar to (4.16) also for higher order potentials. We can use equation (4.15) to study the large N limit. Let us define the variables ɛ ≡ 1/N and x ≡ nɛ. In the large N limit, x is well approximated by a continuous parameter. In view of this, it is convenient to set r(x, α) ≡ R_n(α). We note that r(x, α) is also a function of N, but we suppress this dependence for notational simplicity. We can rewrite (4.15) as

$$\begin{equation}x=r\left(x,\alpha \right)+4\alpha r\left(x,\alpha \right)\left[r\left(x+\varepsilon ,\alpha \right)+r\left(x,\alpha \right)+r\left(x-\varepsilon ,\alpha \right)\right].\end{equation} \tag{ 4.17 }$$

It follows from (4.17) that r(x, α) is symmetric under ɛ ↔ −ɛ and we can expand it in even powers of ɛ

$$\begin{equation}r\left(x,\alpha \right)={r}_{0}\left(x,\alpha \right)+{\varepsilon }^{2}\enspace {r}_{2}\left(x,\alpha \right)+{\varepsilon }^{4}\enspace {r}_{4}\left(x,\alpha \right)+\cdots \enspace .\end{equation} \tag{ 4.18 }$$

Collecting terms with equal powers of ɛ we obtain the expression:

$$\begin{equation}r\left(x+\varepsilon ,\alpha \right)+r\left(x-\varepsilon ,\alpha \right)=2\sum _{n=0}^{\infty }{\varepsilon }^{2n}\sum _{k+p=n}\frac{1}{\left(2p\right)!}\frac{{\mathrm{d}}^{2p}}{\mathrm{d}{x}^{2p}}\enspace {r}_{2k}\left(x,\alpha \right).\end{equation} \tag{ 4.19 }$$

Re-inserting this into (4.17) and comparing powers of ɛ we conclude

$$\begin{equation}x\enspace {\delta }_{s,0}={r}_{2s}\left(x,\alpha \right)+4\alpha \sum _{m+n=s}{r}_{2m}\left(x,\alpha \right)\left[{r}_{2n}\left(x,\alpha \right)+2\sum _{k+p=n}\frac{1}{\left(2p\right)!}\frac{{\mathrm{d}}^{2p}}{\mathrm{d}{x}^{2p}}\enspace {r}_{2k}\left(x,\alpha \right)\right].\end{equation} \tag{ 4.20 }$$

For the cases s = 0 and s = 1 we readily find

$$\begin{equation}{r}_{0}\left(x,\alpha \right)=\frac{-1+\sqrt{1+48\alpha x}}{24\alpha }\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad {r}_{2}\left(x,\alpha \right)=\frac{96{\alpha }^{2}{r}_{0}\left(x,\alpha \right)}{{\left(1+48\alpha x\right)}^{2}}.\end{equation} \tag{ 4.21 }$$

Our final ingredient will be the Euler–Maclaurin formula

$$\begin{equation}\frac{1}{N}\sum _{n=1}^{N}f\left(\frac{n}{N}\right)={\left.{\int }_{0}^{1}\mathrm{d}xf\left(x\right)+\frac{1}{2N}f\left(x\right)\right\vert }_{0}^{1}+{\left.\sum _{n=1}^{p-1}\frac{{B}_{2n}}{\left(2n\right)!}\frac{1}{{N}^{2n}}f{\left(x\right)}^{\left(2n-1\right)}\right\vert }_{0}^{1}+{\mathcal{R}}_{N}.\end{equation} \tag{ 4.22 }$$

In the above, f(x) is a 2p times continuously differentiable function, ${\mathcal{R}}_{N}$ is a remainder term scaling as $\mathcal{O}\left(1/{N}^{2p+1}\right)$, and the B_2n denote the Bernoulli numbers. Applying the Euler–Maclaurin formula to

$$\begin{equation}f\left(x\right)=\left(1-x\right)\mathrm{log}\enspace \frac{r\left(x,\alpha \right)}{x},\end{equation} \tag{ 4.23 }$$

and expanding (4.14) in inverse powers of N, we find

$$\begin{align}\hfill \frac{1}{{N}^{2}}{\mathcal{F}}_{N}\left(\alpha \right)=& -{\int }_{0}^{1}\mathrm{d}x\left(1-x\right)\mathrm{log}\enspace \frac{r\left(x,\alpha \right)}{x}-\frac{1}{N}\enspace \mathrm{log}\enspace \frac{{h}_{0}\left(\alpha \right)}{{h}_{0}\left(0\right)}+\frac{1}{2N}\enspace {\mathrm{lim}}_{x\to 0}\enspace \mathrm{log}\enspace \frac{r\left(x,\alpha \right)}{x}\hfill \\ \hfill & {\left.-\frac{1}{12{N}^{2}}{\left(\left(1-x\right)\mathrm{log}\enspace \frac{r\left(x,\alpha \right)}{x}\right)}^{\left(1\right)}\right\vert }_{0}^{1}\hfill \end{align} \tag{ 4.24 }$$

up to order $\mathcal{O}\left(1/{N}^{4}\right)$ corrections. To obtain h₀(α) we simply evaluate

$$\begin{equation}{h}_{0}\left(\alpha \right)={\int }_{\mathbb{R}}\mathrm{d}\mu \left(\lambda \right)=\frac{{\mathrm{e}}^{\frac{N}{32\alpha }}}{2\sqrt{2}}\enspace \sqrt{\frac{1}{\alpha }}\enspace {K}_{\frac{1}{4}}\left(\frac{N}{32\alpha }\right)=\sqrt{\frac{2\pi }{N}}\left(1-\frac{3}{N}\enspace \alpha +\frac{105}{2{N}^{2}}\enspace {\alpha }^{2}+\cdots \enspace \right),\end{equation} \tag{ 4.25 }$$

where K_n(x) is the modified Bessel function of the second kind. Expanding all three terms in (4.18) up to powers of order $\mathcal{O}\left(1/{N}^{2}\right)$ we find [41]

$$\begin{align}\hfill \frac{1}{{N}^{2}}{\mathcal{F}}_{N}\left(\alpha \right)=& -{\int }_{0}^{1}\mathrm{d}x\enspace \left(1-x\right)\mathrm{log}\enspace \frac{{r}_{0}\left(x,\alpha \right)}{x}\hfill \\ \hfill & \left.-\frac{1}{{N}^{2}}{\left[{\int }_{0}^{1}\mathrm{d}x\enspace \left(1-x\right)\enspace \frac{{r}_{2}\left(x,\alpha \right)}{{r}_{0}\left(x,\alpha \right)}+\frac{1}{12}{\left[\left(1-x\right)\mathrm{log}\enspace \frac{{r}_{0}\left(x,\alpha \right)}{x}\right]}^{\left(1\right)}\right\vert }_{0}^{1}-3\enspace \alpha \right].\hfill \end{align} \tag{ 4.26 }$$

Using (4.21) both integrals in (4.26) can be evaluated analytically. The expression obtained for ${\mathcal{F}}^{\left(0\right)}\left(\alpha \right)$ agrees with the re-summed expression in (3.45) and near α_c = −1/48 we recover

$$\begin{equation}{\mathrm{lim}}_{\alpha \to {\alpha }_{\mathrm{c}}}\enspace {\partial }_{\alpha }^{\left(3\right)}\enspace {\mathcal{F}}_{\text{n.a.}}^{\left(0\right)}\left(\alpha \right)\sim {\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{-1/2}+\cdots \enspace .\end{equation} \tag{ 4.27 }$$

For ${\mathcal{F}}^{\left(1\right)}\left(\alpha \right)$ we obtain

$$\begin{equation}{\mathcal{F}}^{\left(1\right)}\left(\alpha \right)=\frac{1}{24}\left(\mathrm{log}\enspace 4+\mathrm{log}\left(1+48\alpha \right)-2\enspace \mathrm{log}\left(1+\sqrt{1+48\alpha }\right)\right).\end{equation} \tag{ 4.28 }$$

Interestingly, in the case of ${\mathcal{F}}^{\left(1\right)}\left(\alpha \right)$, the non-analyticity near α_c comes in the form of a logarithm. Finally, a small α expansion leads to

$$\begin{equation}{\mathcal{F}}^{\left(0\right)}\left(\alpha \right)=2\alpha -18{\alpha }^{2}+288{\alpha }^{3}-6048{\alpha }^{4}+\frac{746\enspace 496}{5}{\alpha }^{5}-4105\enspace 728{\alpha }^{6}+\mathcal{O}\left({\alpha }^{7}\right),\end{equation} \tag{ 4.29 }$$

$$\begin{equation}{\mathcal{F}}^{\left(1\right)}\left(\alpha \right)=\alpha -30{\alpha }^{2}+1056{\alpha }^{3}-40\enspace 176{\alpha }^{4}+\frac{8004\enspace 096}{5}{\alpha }^{5}-65\enspace 774\enspace 592{\alpha }^{6}+\mathcal{O}\left({\alpha }^{7}\right).\end{equation} \tag{ 4.30 }$$

In addition to ${\mathcal{F}}^{\left(0\right)}\left(\alpha \right)$ and ${\mathcal{F}}^{\left(1\right)}\left(\alpha \right)$, we present ${\mathcal{F}}^{\left(h\right)}\left(\alpha \right)$ for h = 2, 3, and 4 in appendix C. Their non-analytic behaviour near α_c agrees with the general form (3.79).

The orthogonal polynomial method allows us to systematically compute non-planar contributions. Near the critical coupling, however, one might expect that only a small piece of the detailed functions ${\mathcal{F}}^{\left(h\right)}\left(\alpha \right)$ should matter. To this end, let us revisit equation (4.17). Recall that ɛ is parameterically small at large N, admitting an expansion of the type (4.18). We also note from our expressions (4.21) and (C.1) that if we also take α → α_c, expressions are dominated by the region near x = 1. To render the expressions finite near x = 1 we must keep ${\kappa }^{-1}={\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{5/4}N$ fixed as we take N → ∞. This is the double scaling limit we encountered earlier. Thus, we are prompted to study the full equation (4.17) in the double scaling limit. If we take

$$\begin{equation}x=1+\left(\alpha -{\alpha }_{\mathrm{c}}\right)z,\qquad r\left(x,\alpha \right)={r}_{0}\left(1,{\alpha }_{\mathrm{c}}\right)+{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{1/2}\delta r\left(z\right),\qquad \varepsilon ={\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{5/4}\kappa ,\end{equation} \tag{ 4.31 }$$

we can readily show that in the limit α → α_c, and recalling that α_c = −1/48, equation (4.17) implies

$$\begin{equation}\frac{1}{4}\delta r{\left(z\right)}^{2}+\frac{{\kappa }^{2}}{6}\delta {r}^{{\prime\prime}}\left(z\right)+z=0.\end{equation} \tag{ 4.32 }$$

Thus, at least in the double scaling limit, the full genus expansion (3.82) is reduced to solving the above Painlevé I equation [65–68]. This is a remarkable simplification of the original problem.

It is worth mentioning some features of the Painlevé I equation. For instance, the Painlevé I equation remains unchanged under the rescaling

$$\begin{equation}z\to \lambda z,\qquad \delta r\left(z\right)\to {\lambda }^{1/2}\delta r\left(z\right),\qquad \kappa \to {\lambda }^{5/4}\kappa .\end{equation} \tag{ 4.33 }$$

Given the above scaling symmetry, it only makes sense to build perturbative expansions out of a scale invariant quantity such as κ(−z)^−5/4 rather than, say, small κ. The first few terms of this expansion are given by

$$\begin{equation}\delta {r}_{{\pm}}\left(z\right)={\left(-z\right)}^{1/2}\left({\pm}2+\frac{1}{12}\frac{{\kappa }^{2}}{{\left(-z\right)}^{5/2}}\mp \frac{49}{576}\frac{{\kappa }^{4}}{{\left(-z\right)}^{5}}+\frac{1225}{3456}\frac{{\kappa }^{6}}{{\left(-z\right)}^{15/2}}+\sum _{n=4}^{\infty }{a}_{n}^{\left({\pm}\right)}\frac{{\kappa }^{2n}}{{\left(-z\right)}^{5n/2}}\right).\end{equation} \tag{ 4.34 }$$

Expanding further, one observes that the coefficients of the negative branch δr₋(z) are all positive while those of the positive branch δr₊(z) alternate. To specify a solution of the Painlevé I equation, we must pick a set of initial conditions for δr(z) and δr'(z) at some value of z. Upon fixing half of the initial data, most choices for the remaining data lead to the oscillatory branch δr₊(z) [71] (figures 11 and 12).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Of the two branches (4.34), only the negative branch δr₋(z) connects to the double scaling limit stemming from the matrix integral. From now on we focus on this branch. The coefficients for δr₋(z) satisfy the recursion relation

$$\begin{equation}{a}_{n+1}=\frac{25{n}^{2}-1}{24}\enspace {a}_{n}+\frac{1}{4}\sum _{m=1}^{n}{a}_{m}{a}_{n+1-m},\qquad {a}_{0}=-2,\quad n{\geqslant}0.\end{equation} \tag{ 4.35 }$$

For large n the coefficients above grow as $\sim 5/\left(4\sqrt{6}\right)\enspace {\Gamma}\left(2n-1/2\right)$ [63]. This implies that the branch reproducing the perturbative expansion (4.34) is not Borel summable [66]. In fact, there is a family of solutions to (4.32) containing non-perturbative terms [60, 63] of the type

$$\begin{equation}{\epsilon}\left(z\right)=c\enspace {\left(\frac{-z}{{\kappa }^{4/5}}\right)}^{-1/8}{\times}{\mathrm{e}}^{-\frac{4}{5\kappa }\sqrt{6}\enspace {\left(-z\right)}^{5/4}}.\end{equation} \tag{ 4.36 }$$

This family admits the same perturbative expansion (4.34) at large (−z)^5/4/κ. The parameter c is an undetermined integration constant. The form (4.36) is obtained by considering a small deviation δr(z) = δr₀(z) + (z) and solving for (z) in a WKB approximation. To linear order in (z) we have

$$\begin{equation}\frac{1}{4}\delta {r}_{0}{\left(z\right)}^{2}+\frac{{\kappa }^{2}}{6}\delta {r}_{\text{0}}^{{\prime\prime}}\left(z\right)+z=0,\qquad \frac{{\kappa }^{2}}{6}{{\epsilon}}^{{\prime\prime}}\left(z\right)+\frac{1}{2}\delta {r}_{0}\left(z\right){\epsilon}\left(z\right)=0.\end{equation} \tag{ 4.37 }$$

For completeness, we mention that one can also construct a perturbative expansion for solutions near a double pole at z = z₀

$$\begin{equation}\delta r\left(z\right)=-\frac{4}{{\left(z-{z}_{0}\right)}^{2}}+\frac{3{z}_{0}}{5}{\left(z-{z}_{0}\right)}^{2}+{\left(z-{z}_{0}\right)}^{3}+{a}_{4}{\left(z-{z}_{0}\right)}^{4}-\frac{3{z}_{0}^{2}}{100}{\left(z-{z}_{0}\right)}^{6}+\cdots \enspace ,\end{equation} \tag{ 4.38 }$$

where we have set κ = 1. Both the location of the pole, z₀, as well as the value of a₄ are free parameters.

We now relate δr₋(z) to the double scaling limit expression $\mathcal{F}\left(\kappa \right)$ introduced in (3.82). To do so, we note that in the double scaling limit, the expression (4.24) is dominated by

$$\begin{equation}{\mathcal{F}}_{N}\left(\alpha \right)\approx -{N}^{2}\enspace {\int }_{0}^{1}\mathrm{d}x\left(1-x\right)\mathrm{log}\enspace \frac{r\left(x,\alpha \right)}{x}+\cdots \enspace .\end{equation} \tag{ 4.39 }$$

Combining the above expression with (4.31) we obtain

$$\begin{align}\hfill \mathcal{F}\left(\kappa \right)& \equiv {\mathrm{lim}}_{\text{d.s.l.}}\left({\mathcal{F}}_{N}\left(\alpha \right)-{N}^{2}{\mathcal{F}}^{\left(0\right)}\left(\alpha \right)\right)\hfill \\ \hfill & =\frac{12\enspace 288\sqrt{3}}{5{\kappa }^{2}}+\frac{1}{2{\kappa }^{2}}\enspace {\mathrm{lim}}_{\varepsilon \to {0}^{+}}{\int }_{-{\left(\kappa /\varepsilon \right)}^{4/5}}^{0}\mathrm{d}zz\left(\delta {r}_{-}\left(z\right)+2{\left(-z\right)}^{1/2}\right).\hfill \end{align} \tag{ 4.40 }$$

It from (3.45) that ${\mathcal{F}}^{\left(0\right)}\left({\alpha }_{\mathrm{c}}\right)=\left(7/24-\mathrm{log}\enspace 2/2\right)$. We now recall the perturbative expansion (4.34) which is reliable for (κ/ɛ)^4/5 ⩾ −z ≫ κ^4/5. Since the first term in the expansion (4.34) is κ independent and grows at large z it is convenient to subtract it from δr₋(z), as we have done in (4.40). It is of interest to see how the perturbative expansion (4.34) is encoded in $\mathcal{F}\left(\kappa \right)$. We must recall that there is in fact no actual small κ expansion since the scaling symmetry (4.33) always allows us to rescale κ. Instead, it is convenient to consider a slightly generalised integral given by

$$\begin{equation}\mathcal{F}\left(\kappa ;\eta \right)\equiv \frac{1}{2{\kappa }^{2}}\enspace {\mathrm{lim}}_{\varepsilon \to {0}^{+}}{\int }_{-{\left(\kappa /\varepsilon \right)}^{4/5}}^{-\eta }\mathrm{d}z\enspace \left(z+\eta \right)\left(\delta {r}_{-}\left(z\right)+2{\left(-z\right)}^{1/2}\right).\end{equation} \tag{ 4.41 }$$

We have introduced a cutoff at some small value −η ≳ κ^4/5 to avoid entering a region where the perturbative expansion (4.34) is invalidated. We note that ${\partial }_{\eta }^{2}\mathcal{F}\left(\kappa ;\eta \right)=-\left(\delta {r}_{-}\left(-\eta \right)+2{\eta }^{1/2}\right)/2{\kappa }^{2}$. Evaluating $\mathcal{F}\left(\kappa ;\eta \right)$ in a small κ/η^5/4 expansion we find

$$\begin{equation}\mathcal{F}\left(\kappa ;\eta \right)=\frac{1}{24}\left(\mathrm{log}\left(\eta {\left(\frac{\varepsilon }{\kappa }\right)}^{4/5}\right)+1\right)-\frac{7{\kappa }^{2}}{1440\enspace {\eta }^{5/2}}-\frac{245{\kappa }^{4}}{41\enspace 472{\eta }^{5}}+\cdots \enspace .\end{equation} \tag{ 4.42 }$$

We recognise the coefficient of the logarithm from our previous expression (4.28).

4.3.1. Double scaling limit for multicritical models

One can perform a similar analysis for the multicritical model introduced in (3.49). Using the methods discussed in appendix A it is straightforward to derive the equation

$$\begin{equation}{r}_{0}\left(x,\gamma \right)-4{r}_{0}{\left(x,\gamma \right)}^{2}+\frac{16}{3}{r}_{0}{\left(x,\gamma \right)}^{3}=\gamma \enspace x.\end{equation} \tag{ 4.43 }$$

The function r₀(x, γ) encodes the planar contribution. Solving for r₀(x) and selecting the branch relevant to the matrix integral, we obtain

$$\begin{equation}{r}_{0}\left(x,\gamma \right)=\frac{1}{4}\left(1+{\left(-1+12\gamma x\right)}^{1/3}\right).\end{equation} \tag{ 4.44 }$$

From the above find that γ_c = 1/12, consistent with our previous discussion. We now define

$$\begin{equation}z\equiv {N}^{6/7}\left(1-\frac{\gamma }{{\gamma }_{\mathrm{c}}}x\right),\qquad r\left(x,\gamma \right)\equiv {r}_{0}\left(1,{\gamma }_{\mathrm{c}}\right)\left(1-{N}^{-2/7}\delta r\left(z\right)\right),\qquad \varepsilon \equiv 1/N.\end{equation} \tag{ 4.45 }$$

The double scaling limit is given by taking N large while keeping ${\vartheta }^{-1}\equiv {\left(\gamma -{\gamma }_{\mathrm{c}}\right)}^{7/6}N$ fixed. The equation governing δr(z) in the double scaling limit reads

$$\begin{equation}\delta r{\left(z\right)}^{3}-\frac{1}{2}\delta {r}^{\prime }{\left(z\right)}^{2}-\delta r\left(z\right)\delta {r}^{{\prime\prime}}\left(z\right)+\frac{1}{10}\delta {r}^{\left(4\right)}\left(z\right)-z=0.\end{equation} \tag{ 4.46 }$$

We can solve the above equation in a large z expansion

$$\begin{equation}\delta r\left(z\right)={z}^{1/3}\left(1+\sum _{n=1}^{\infty }{b}_{k}{z}^{-7n/3}\right),\end{equation} \tag{ 4.47 }$$

with the first few terms given by

$$\begin{align}\hfill \delta r\left(z\right)={z}^{1/3}& \left(1-\frac{1}{18}{z}^{-7/3}-\frac{7}{108}{z}^{-14/3}-\frac{4199}{17\enspace 496}{z}^{-21/3}-\frac{409\enspace 297}{262\enspace 440}{z}^{-28/3}\right.\hfill \\ \hfill & \left.-\frac{101\enspace 108\enspace 329}{9447\enspace 840}{z}^{-35/3}+\frac{25\enspace 947\enspace 984\enspace 239}{191\enspace 318\enspace 760}{z}^{-42/3}+\cdots \enspace \right).\hfill \end{align} \tag{ 4.48 }$$

Expanding further the coefficients of the above expansion are found to be alternating. Contrary to the coefficients (4.35) of the Painlevé I equation, the coefficients (4.47) of the perturbative expansion of (4.46) are Borel summable [63]. Related to this, ${\mathcal{M}}_{N}$ (3.1) for the polynomial V(M) (3.49) is well defined at finite N for all positive γ.

It is also interesting to assess whether the above series can admit non-perturbative corrections. Expanding about the δr(z) solution given by the large z expansion (4.48), we find that to linear order the small deviation (z) must satisfy

$$\begin{equation}\frac{1}{10}{{\epsilon}}^{\left(4\right)}\left(z\right)-{z}^{1/3}{{\epsilon}}^{{\prime\prime}}\left(z\right)-\frac{1}{3}{z}^{-2/3}{{\epsilon}}^{\prime }\left(z\right)+3{z}^{2/3}{\epsilon}\left(z\right)=0.\end{equation} \tag{ 4.49 }$$

Performing a WKB analysis we obtain ${\epsilon}\left(z\right)\approx {c}_{{\pm}}\enspace {z}^{-1/4}\enspace {\mathrm{e}}^{-\frac{6}{7}\sqrt{5{\pm}\mathrm{i}\sqrt{5}}\enspace {z}^{7/6}}$. It is interesting to note that the exponent is now complex. Interestingly, the polynomial governing the eigenvalues, V(λ) given in (3.49), has no maxima. Rather it has critical points at the complex values ${\lambda }_{{\pm}}^{2}=\frac{1}{4}\left(5{\pm}\mathrm{i}\sqrt{5}\right)$. Further multicritical models and their properties are presented in appendix B.

As a final note, we discuss a connection between the non-perturbative corrections of the Painlevé I equation and the original matrix integral (3.1). To do so, we consider configurations for which one (or multiple) of the eigenvalues is sitting outside of the densely filled eigenvalue distribution. Such a configuration will give a subdominant contribution to the matrix integral whose size we will now compute.

For simplicity, we focus on the case where a single eigenvalue λ ≡ λ_N is separated from the remaining N − 1 eigenvalues. The matrix integral ${\mathcal{M}}_{N}$ can be written in the following form

$$\begin{equation}\text{vol}\;\frac{U{\left(1\right)}^{N}{\times}{S}_{N}}{U\left(N\right)}{\mathcal{M}}_{N}={\int }_{\mathbb{R}}\mathrm{d}\lambda \enspace {\mathrm{e}}^{-NV\left(\lambda \right)}{\int }_{{\mathbb{R}}^{N-1}}\prod _{I=1}^{N-1}\mathrm{d}{\mu }_{I}{{\Delta}}_{N-1}{\left(\mu \right)}^{2}\enspace {\mathrm{e}}^{-N\enspace \sum _{I=1}^{N-1}V\left({\mu }_{I}\right)}\prod _{I=1}^{N-1}{\left(\lambda -{\mu }_{I}\right)}^{2},\end{equation} \tag{ 4.50 }$$

where the expectation value is defined in (3.52). We define the effective polynomial

$$\begin{equation}{V}_{\text{eff}}\left(\lambda \right)\equiv V\left(\lambda \right)-\frac{1}{N}\enspace \left\langle \mathrm{T}\mathrm{r}\enspace \mathrm{log}{\left(\lambda \enspace {\mathbb{I}}_{N-1}-{M}_{N-1}\right)}^{2}\right\rangle ,\end{equation} \tag{ 4.51 }$$

which is the leading order correction of the polynomial V(λ) in the large N limit. Explicitly, assuming that ${\mathcal{M}}_{N}$ admits a single cut eigenvalue density ρ(λ) symmetric about the origin, we find

$$\begin{equation}{V}_{\text{eff}}\left(\lambda \right)=V\left(\lambda \right)-{\int }_{-a}^{a}\mathrm{d}\mu \enspace \rho \left(\mu \right)\mathrm{log}{\left(\lambda -\mu \right)}^{2},\end{equation} \tag{ 4.52 }$$

to leading order in the large N limit. It follows from (3.16) that within the range λ ∈ (−a, a), V_eff(λ) must be a constant V₀ given by

$$\begin{equation}{V}_{0}\equiv V\left(0\right)-{\int }_{-a}^{a}\mathrm{d}\mu \enspace \rho \left(\mu \right)\mathrm{log}\enspace {\mu }^{2}.\end{equation} \tag{ 4.53 }$$

The effective polynomial ceases to be constant for λ² > a². To leading order in the large N limit, the piece of the integral pertinent to the separated eigenvalue is given by

$$\begin{equation}{\mathcal{I}}_{N}\equiv {\int }_{\mathbb{R}}\mathrm{d}\lambda \enspace {\mathrm{e}}^{-N{V}_{\text{eff}}\left(\lambda \right)}=2a\enspace {\mathrm{e}}^{-N{V}_{0}}\left(1+\frac{1}{2a}{\int }_{\lambda \notin \left[-a,a\right]}\mathrm{d}\lambda \enspace {\mathrm{e}}^{-N\left({V}_{\text{eff}}\left(\lambda \right)-{V}_{0}\right)}\right).\end{equation} \tag{ 4.54 }$$

The critical points of V_eff(λ) for λ² > a² can be found by solving

$$\begin{equation}{V}_{\text{eff}}^{\prime }\left(\lambda \right)=\mathrm{sign}\enspace \lambda \enspace \text{Re}\sqrt{{V}^{\prime }{\left(\lambda \right)}^{2}-4\mathcal{P}\left(\lambda \right)}=0,\end{equation} \tag{ 4.55 }$$

where $\mathcal{P}\left(\lambda \right)$ is defined in (3.60).

At large N we can apply the saddle point approximation to the integral in (4.54), which requires evaluating its exponent at its critical points. We are interested in contributions that stem from a single eigenvalue sitting at a critical point outside the interval containing the dense set of remaining N − 1 eigenvalues. These contributions can lead to corrections which are exponentially suppressed in N. Consequently they do not contribute to the perturbative 1/N expansion—they are non-perturbative. In certain contexts, such suppressed configurations bear the name instantons.

Example. As a concrete example, we consider the quartic polynomial (3.35). Taking α ∈ (α_c, 0), we find the following two maxima for the effective polynomial:

$$\begin{equation}{\lambda }_{{\pm}}={\pm}\sqrt{-\frac{1}{6\alpha }-\frac{\sqrt{1+48\alpha }}{12\alpha }}.\end{equation} \tag{ 4.56 }$$

For λ ∈ (−a₊, a₊), with a₊ given in (3.38) and the eigenvalue density in (3.39), we find

$$\begin{equation}{V}_{0}=-\frac{1}{48\alpha }+\frac{\sqrt{1+48\alpha }+24\alpha \left(1+\mathrm{log}\enspace 576-2\enspace \mathrm{log}\enspace \frac{\sqrt{1+48\alpha }-1}{\alpha }\right)}{48\alpha }.\end{equation} \tag{ 4.57 }$$

Evaluating V_eff(λ) at λ = λ_± one finds

$$\begin{equation}{V}_{\text{eff}}\left({\lambda }_{{\pm}}\right)-{V}_{0}=\sqrt{3}\frac{\sqrt{2+\sqrt{y}}}{\left(1-y\right)}\enspace {y}^{1/4}-2\enspace \mathrm{R}\mathrm{e}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\frac{\sqrt{2+\sqrt{y}}}{\sqrt{3}{y}^{1/4}}\right),\end{equation} \tag{ 4.58 }$$

where we have defined y ≡ 1 − α/α_c. A small y expansion reveals

$$\begin{equation}{V}_{\text{eff}}\left({\lambda }_{{\pm}}\right)-{V}_{0}=\frac{4}{5}\sqrt{6}\enspace {y}^{5/4}+\frac{2}{7}\sqrt{\frac{3}{2}}\enspace {y}^{7/4}+\cdots \enspace .\end{equation} \tag{ 4.59 }$$

Interestingly, in the double scaling limit, where N is taken to infinity while keeping Ny^5/4 fixed, the leading term in the above expansion gives a contribution to (4.54) that remains finite [63]. The correction is reminiscent of the non-perturbative correction (4.36) found in our discussion of the Painlevé I equation.

Following an approach inspired by that pursued in section 3, we parameterise A and B in terms of their eigenvalues

$$\begin{equation}A={\mathcal{U}}_{A}{D}_{A}\enspace {\mathcal{U}}_{A}^{{\dagger}},\qquad {D}_{A}=\mathrm{diag}\left({a}_{1},\dots ,{a}_{N}\right)\end{equation} \tag{ 5.2 }$$

and similarly for B. Recalling the discussion in subsection 3.1, the matrices ${\mathcal{U}}_{A}$ and ${\mathcal{U}}_{B}$ are general elements of U(N)/U(1)^N. Integrating over the angular variables in (5.1) leads to

$$\begin{equation}{\mathcal{T}}_{N}\left(c\right)=\mathrm{v}\mathrm{o}\mathrm{l}\enspace \frac{U\left(N\right)}{U{\left(1\right)}^{N}{\times}{S}_{N}}{\int }_{{\mathbb{R}}^{2N}}\prod _{I=1}^{N}\mathrm{d}{a}_{I}\enspace \mathrm{d}{b}_{I}\enspace {{\Delta}}_{N}{\left(a\right)}^{2}{{\Delta}}_{N}{\left(b\right)}^{2}\enspace {\mathrm{e}}^{-N\enspace \sum _{I=1}^{N}\left(V\left({a}_{I}\right)+V\left({b}_{I}\right)\right)}I\left(a,b,c\right),\end{equation} \tag{ 5.3 }$$

where I(a, b, c) is known as the Harish–Chandra–Itzykson–Zuber (HCIZ) integral

$$\begin{equation}I\left(a,b,c\right)\equiv \frac{1}{\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}_{N}}\int \left[\mathrm{d}\mathcal{L}\right]\enspace {\mathrm{e}}^{2cN\enspace \mathrm{T}\mathrm{r}\enspace {D}_{A}\mathcal{U}{D}_{B}{\mathcal{U}}^{{\dagger}}},\qquad \mathcal{U}\equiv {\mathcal{U}}_{A}^{{\dagger}}{\mathcal{U}}_{B}.\end{equation} \tag{ 5.4 }$$

The integral I(a, b, c) evaluates to [86]

$$\begin{equation}I\left(a,b,c\right)=\mathrm{v}\mathrm{o}\mathrm{l}\enspace \frac{U\left(N\right)}{U{\left(1\right)}^{N}{\times}{S}_{N}}\frac{G\left(N+1\right)}{{\left(2cN\right)}^{N\left(N-1\right)/2}}\frac{\mathrm{det}{\left({\mathrm{e}}^{2cN{a}_{I}{b}_{J}}\right)}_{1{\leqslant}I,J{\leqslant}N}}{{{\Delta}}_{N}\left(a\right){{\Delta}}_{N}\left(b\right)},\end{equation} \tag{ 5.5 }$$

where $G\left(N+1\right)={\prod }_{K=1}^{N-1}K!$ denotes Barnes G-function. Finally, using the Leibniz formula for the determinant we obtain for a measure $\mathrm{d}\tilde {\mu }\left(a,b\right)$ completely anti-symmetric under a_I ↔ a_J and b_I ↔ b_J

$$\begin{equation}\int \mathrm{d}\tilde {\mu }\left(a,b\right)\mathrm{det}{\left({\mathrm{e}}^{2cN{a}_{I}{b}_{J}}\right)}_{1{\leqslant}I,J{\leqslant}N}=\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}_{N}\int \mathrm{d}\tilde {\mu }\left(a,b\right)\enspace {\mathrm{e}}^{2cN\sum _{I=1}^{N}{a}_{I}{b}_{I}}.\end{equation} \tag{ 5.6 }$$

We can thus write (5.3) as

$$\begin{equation}{\mathcal{T}}_{N}\left(c\right)={\left(\mathrm{v}\mathrm{o}\mathrm{l}\enspace \frac{U\left(N\right)}{U{\left(1\right)}^{N}}\right)}^{2}\frac{G\left(N+1\right)}{\mathrm{v}\mathrm{o}\mathrm{l}\enspace {S}_{N}}\int \prod _{I=1}^{N}\mathrm{d}{a}_{I}\mathrm{d}{b}_{I}\enspace \frac{{{\Delta}}_{N}\left(a\right){{\Delta}}_{N}\left(b\right)}{{\left(2cN\right)}^{N\left(N-1\right)/2}}\enspace {\mathrm{e}}^{-N\enspace \sum _{I}\left(V\left({a}_{I}\right)+V\left({b}_{I}\right)-2\enspace c\enspace {a}_{I}{b}_{I}\right)},\end{equation} \tag{ 5.7 }$$

The above expression resembles the single matrix expression (3.12). However, it has some differences such as the Vandermonde determinant for each collection of eigenvalues appearing with a single power.

HCIZ integral for N = 2. As a simple example, we consider the case N = 2. We note that $\left[\mathrm{d}\mathcal{L}\right]$ is the line element on the sphere

$$\begin{equation}\mathrm{d}{s}^{2}=\frac{1}{2}\left(\mathrm{d}{\theta }^{2}+\mathrm{sin}\enspace {\theta }^{2}\mathrm{d}{\phi }^{2}\right).\end{equation} \tag{ 5.8 }$$

The HCIZ integral for N = 2 is thus equal to

$$\begin{align}\hfill I\left(a,b,c\right)& =\frac{1}{4}{\int }_{0}^{2\pi }\mathrm{d}\phi {\int }_{0}^{\pi }\mathrm{d}\theta \enspace \mathrm{sin}\enspace \theta \enspace {\mathrm{e}}^{4c\enspace \left({\mathrm{sin}}^{2}\frac{\theta }{2}\left({a}_{1}{b}_{1}+{a}_{2}{b}_{2}\right)+{\mathrm{cos}}^{2}\frac{\theta }{2}\left({a}_{1}{b}_{2}+{a}_{2}{b}_{1}\right)\right)}\hfill \\ \hfill & =\frac{\pi }{4c}\enspace \frac{1}{{{\Delta}}_{2}\left(a\right){{\Delta}}_{2}\left(b\right)}\enspace \mathrm{det}\left(\begin{matrix}\hfill {\mathrm{e}}^{4c\enspace {a}_{1}{b}_{1}}\hfill & \hfill {\mathrm{e}}^{4c\enspace {a}_{1}{b}_{2}}\hfill \\ \hfill {\mathrm{e}}^{4c\enspace {a}_{2}{b}_{1}}\hfill & \hfill {\mathrm{e}}^{4c\enspace {a}_{2}{b}_{2}}\hfill \end{matrix}\right),\hfill \end{align} \tag{ 5.9 }$$

confirming (5.5) for N = 2.

Having reached the expression (5.7), we can proceed to use the techniques developed in the previous to solve the integral at large N. One approach is to introduce an eigenvalue density for the a_I and b_I eigenvalues and study a generalisation of the loop equations introduced in section 3.4 for multi-matrix models [55, 56]. Here, we will consider the orthogonal polynomial approach. We take V(a) and V(b) to be the quartic polynomial previously considered in the single matrix case (3.35). Thus, we would like to solve

$$\begin{equation}{\mathcal{T}}_{N}\left(\alpha ,c\right)={\mathcal{C}}_{N}{\int }_{{\mathbb{R}}^{2N}}\prod _{I=1}^{N}\mathrm{d}{a}_{I}\enspace \mathrm{d}{b}_{I}\enspace {{\Delta}}_{N}\left(a\right){{\Delta}}_{N}\left(b\right)\enspace {\mathrm{e}}^{-N\sum _{I=1}^{N}W\left({a}_{I},{b}_{I}\right)},\end{equation} \tag{ 5.10 }$$

with

$$\begin{equation}W\left(a,b\right)=\frac{1}{2}\left({a}^{2}+{b}^{2}\right)+\alpha \left({a}^{4}+{b}^{4}\right)-2c\enspace a\enspace b,\quad {c}^{2}\in \left(0,1/4\right).\end{equation} \tag{ 5.11 }$$

We have absorbed the volume pre-factors into an N-dependent constant ${\mathcal{C}}_{N}$. It is worth reiterating that the Vandermonde for each set of eigenvalues only appears linearly in (5.10).

At this stage we introduce a set of monic polynomials Q_n(x), with $x\in \mathbb{R}$, subject to the condition Q_n(−x) = (−1)ⁿQ_n(x). A priori we should have introduced two different sets of monic polynomials for the eigenvalues of A and B respectively. However due to the symmetry of W(a, b) under interchanging these eigenvalues we can choose the same set of monic polynomials. Orthogonality of the Q_n(x) is expressed as

$$\begin{equation}{\int }_{{\mathbb{R}}^{2}}\mathrm{d}\mu \left(a,b\right)\enspace {Q}_{n}\left(a\right){Q}_{m}\left(b\right)={k}_{n}\left(\alpha ,c\right){\delta }_{mn},\qquad \mathrm{d}\mu \left(a,b\right)\equiv \mathrm{d}a\enspace \mathrm{d}b\enspace {\mathrm{e}}^{-NW\left(a,b\right)},\end{equation} \tag{ 5.12 }$$

and the three-term recurrence relation is given by

$$\begin{equation}x\enspace {Q}_{n}\left(x\right)={Q}_{n+1}\left(x\right)+{R}_{n}\enspace {Q}_{n-1}\left(x\right)+{S}_{n}\enspace {Q}_{n-3}\left(x\right).\end{equation} \tag{ 5.13 }$$

Notice that for the two-matrix model we have to introduce an additional variable S_n. We can only choose a polynomial that has the same parity as xQ_n(x) under x ↔ −x since we have an even measure. Thus, we exclude Q_n(x) and Q_n−2(x). Following the same steps leading to (4.12) we rewrite (5.10) as

$$\begin{equation}{\mathcal{T}}_{N}\left(\alpha ,c\right)={\mathcal{C}}_{N}N!\prod _{n=0}^{N-1}{k}_{n}\left(\alpha ,c\right).\end{equation} \tag{ 5.14 }$$

Combining orthogonality and the three-term recurrence relation we obtain, using an integration by parts, three independent equations for the three unknowns k_n, R_n and S_n

$$\begin{equation}{\int }_{{\mathbb{R}}^{2}}\mathrm{d}\mu \left(a,b\right)\enspace {Q}_{n-1}\left(a\right){Q}_{n}\left(b\right){\partial }_{a}W\left(a,b\right)=0,\end{equation} \tag{ 5.15 }$$

$$\begin{equation}{\int }_{{\mathbb{R}}^{2}}\mathrm{d}\mu \left(a,b\right)\enspace {Q}_{n}\left(a\right){Q}_{n-1}\left(b\right){\partial }_{a}W\left(a,b\right)=\frac{n}{N}\enspace {k}_{n-1}\left(\alpha ,c\right),\end{equation} \tag{ 5.16 }$$

$$\begin{equation}{\int }_{{\mathbb{R}}^{2}}\mathrm{d}\mu \left(a,b\right)\enspace {Q}_{n-3}\left(a\right){Q}_{n}\left(b\right){\partial }_{a}W\left(a,b\right)=0.\end{equation} \tag{ 5.17 }$$

Solving the three integrals on the left-hand side is straightforward. The only slightly more difficult one is the second integral where we have to apply the three-term recurrence relation multiple times. It is now convenient to define z_n(α, c) ≡ 2αk_n(α, c)/(ck_n−1(α, c)) for α ≠ 0 and z_n(0, c) ≡ 2k_n(0, c)/(ck_n−1(0, c)) for α = 0. For α ≠ 0, the above three equations then imply

$$\begin{equation}{z}_{n}\left(\alpha ,c\right)=4\alpha \enspace {R}_{n}\left(\alpha ,c\right){\left(1+4\alpha \left({R}_{n+1}\left(\alpha ,c\right)+{R}_{n}\left(\alpha ,c\right)+{R}_{n-1}\left(\alpha ,c\right)\right)\right)}^{-1},\end{equation} \tag{ 5.18 }$$

$$\begin{align}\hfill {z}_{n}\left(\alpha ,c\right)& =\frac{\alpha }{{c}^{2}}\enspace {R}_{n}\left(\alpha ,c\right)\left(1+4\alpha \left({R}_{n+1}\left(\alpha ,c\right)+{R}_{n}\left(\alpha ,c\right)+{R}_{n-1}\left(\alpha ,c\right)\right)\right)\enspace \hfill \\ \hfill & \quad +\frac{4{\alpha }^{2}}{{c}^{2}}\left({S}_{n+2}\left(\alpha ,c\right)+{S}_{n+1}\left(\alpha ,c\right)+{S}_{n}\left(\alpha ,c\right)\right)-\frac{\alpha }{{c}^{2}}\frac{n}{N},\hfill \end{align} \tag{ 5.19 }$$

and

$$\begin{equation}{z}_{n}\left(\alpha ,c\right){z}_{n-1}\left(\alpha ,c\right){z}_{n-2}\left(\alpha ,c\right)=\frac{4{\alpha }^{2}}{{c}^{2}}{S}_{n}\left(\alpha ,c\right).\end{equation} \tag{ 5.20 }$$

For α = 0 the three equations become

$$\begin{equation}{z}_{n}\left(0,c\right)=4{R}_{n}\left(0,c\right),\qquad {c}^{2}{z}_{n}\left(0,c\right)={R}_{n}\left(0,c\right)-\frac{n}{N},\qquad {S}_{n}\left(0,c\right)=0.\end{equation} \tag{ 5.21 }$$

Upon introducing ɛ = 1/N and x = nɛ, at large N we can express z_n(α, c), R_n(α, c), and S_n(α, c) in terms of continuous functions

$$\begin{equation}{z}_{n}\left(\alpha ,c\right)\equiv z\left(x,\alpha ,c\right),\qquad {R}_{n}\left(\alpha ,c\right)\equiv r\left(x,\alpha ,c\right),\qquad {S}_{n+1}\left(\alpha ,c\right)\equiv s\left(x,\alpha ,c\right).\end{equation} \tag{ 5.22 }$$

To leading order in a small ɛ-expansion the expressions (5.18)–(5.20) are given by

$$\begin{equation}z\left(x,\alpha ,c\right)=4\alpha \enspace r\left(x,\alpha ,c\right){\left(1+12\alpha \enspace r\left(x,\alpha ,c\right)\right)}^{-1},\end{equation} \tag{ 5.23 }$$

$$\begin{equation}z\left(x,\alpha ,c\right)=\frac{\alpha }{{c}^{2}}\enspace r\left(x,\alpha ,c\right)\left(1+12\alpha \enspace r\left(x,\alpha ,c\right)\right)+\frac{12{\alpha }^{2}}{{c}^{2}}s\left(x,\alpha ,c\right)-\frac{\alpha x}{{c}^{2}},\end{equation} \tag{ 5.24 }$$

$$\begin{equation}z{\left(x,\alpha ,c\right)}^{3}=\frac{4{\alpha }^{2}}{{c}^{2}}\enspace s\left(x,\alpha ,c\right).\end{equation} \tag{ 5.25 }$$

Eliminating r(x, α, c) and s(x, α, c) in the above equations we find the following expressions

$$\begin{equation}\frac{\omega \left(z\left(x,\alpha ,c\right)\right)}{{\left(1-3z\left(x,\alpha ,c\right)\right)}^{2}}=0\quad \text{and}\quad z\left(x,0,c\right)=\frac{4x}{\left(1-4{c}^{2}\right)},\end{equation} \tag{ 5.26 }$$

where we have defined

$$\begin{equation}\omega \left(z\right)\equiv -108{c}^{2}{z}^{5}+72{c}^{2}{z}^{4}+24{c}^{2}{z}^{3}+3\left(12\alpha x-8{c}^{2}\right){z}^{2}-\left(24\alpha x-4{c}^{2}+1\right)z+4\alpha x.\end{equation} \tag{ 5.27 }$$

We have thus boiled down the large N approximation of our problem to solving a quintic polynomial. Generally speaking, quintic polynomials do not admit simple solutions. Nevertheless, for certain special values of α and c the quintic of interest (5.27) admits simple solutions. These turn out to encode a certain critical behaviour in the planar approximation. For instance, when

$$\begin{align}\hfill {\alpha }_{1}x& =\frac{1}{9}\left(\sqrt{2c}+3c-2{c}^{2}\right),\qquad {\alpha }_{2}x=-\frac{1}{9}\left(\sqrt{2c}-3c+2{c}^{2}\right),\qquad {\alpha }_{3}x=-\frac{1}{48}\left(1-\frac{32{c}^{2}}{3}\right),\hfill \\ \hfill {\alpha }_{4}x& =\frac{1}{9}\left(\sqrt{-2c}-3c-2{c}^{2}\right),\qquad {\alpha }_{5}x=-\frac{1}{9}\left(\sqrt{-2c}+3c+2{c}^{2}\right).\hfill \end{align} \tag{ 5.28 }$$

the discriminant of the quintic (5.27) vanishes. Consequently, ω(z) shares a zero with its derivative ω'(z). For α₂x and α₄x the discriminant of ω'(z), which now is a polynomial in c, vanishes provided c = ±1/8. A quartic polynomial with a vanishing discriminant contains a multiple root, which for the specific case of ω'(z) lies at z = −1/3. It can be easily checked that (5.27) at c = ±1/8 and αx = −5/288 can be written as

$$\begin{equation}\omega \left(z\right)=-\frac{3}{16}{\left(z+\frac{1}{3}\right)}^{3}\left(9{z}^{2}-15z+10\right).\end{equation} \tag{ 5.29 }$$

To identify the non-analytic structure of interest, it is enough to solve (5.26) with c = ±1/8, in a small $\delta =\left(\alpha x+5/288\right)$ expansion about the z = −1/3 solution. To leading order we find

$$\begin{equation}z\left(x,\alpha ,{\pm}1/8\right)=-\frac{1}{3}+\frac{1}{3}\enspace \sqrt[3]{\frac{5}{6}}\enspace \sqrt[3]{1+\frac{288}{5}\enspace \alpha x}+\cdots \enspace .\end{equation} \tag{ 5.30 }$$

The above expansion exhibits the leading non-analyticity we were after. Away from c = ±1/8 the quintic (5.27) has at most a second order zero leading to a different non-analytic structure more akin to that of the single matrix case.

Planar contribution. As for the study of single matrix integrals, we can apply the orthogonal polynomial technique to evaluate ${\mathcal{G}}_{N}\left(\alpha ,c\right)\equiv -\mathrm{log}\enspace {\mathcal{T}}_{N}\left(\alpha ,c\right)/{\mathcal{T}}_{N}\left(0,c\right)$ in a systematic large N expansion. Here ${\mathcal{T}}_{N}\left(\alpha ,c\right)$ is obtained from (5.14) in the same way that (4.12) was obtained for the single matrix integral:

$$\begin{equation}{\mathcal{T}}_{N}\left(\alpha ,c\right)={\mathcal{C}}_{N}N!{\left(\frac{c}{2\alpha }\right)}^{N\left(N-1\right)/2}{k}_{0}^{N}\enspace {z}_{1}^{N-1}{z}_{2}^{N-2}\dots {z}_{N-1}.\end{equation} \tag{ 5.31 }$$

Performing steps analogous to those leading to (4.14) and combining the resulting expression for ${\mathcal{G}}_{N}\left(\alpha ,c\right)$ with the Euler–Maclaurin formula (4.22) results in the expression

$$\begin{equation}\frac{1}{{N}^{2}}\enspace {\mathcal{G}}_{N}\left(\alpha ,c\right)=-{\int }_{0}^{1}\mathrm{d}x\enspace \left(1-x\right)\mathrm{log}\enspace \frac{z\left(x,\alpha ,c\right)}{z\left(x,0,c\right)}+\frac{1}{2}\enspace \mathrm{log}\enspace \alpha +\mathcal{O}\left({N}^{-2}\right).\end{equation} \tag{ 5.32 }$$

Evaluating (5.32) for c = ±1/8 by using the expansion (5.30), we conclude that the non-analytic behaviour is given by

$$\begin{equation}{\mathrm{lim}}_{\alpha \to {\alpha }_{\mathrm{c}}}\enspace {\partial }_{\alpha }^{\left(3\right)}{\mathcal{G}}_{\text{n.a.}}^{\left(0\right)}\left(\alpha ,{\pm}1/8\right)\sim {\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{-2/3},\end{equation} \tag{ 5.33 }$$

where we recall the critical value α_c = −5/288.

Non-planar contributions. One can continue along these lines and calculate non-planar contributions. As for the single matrix case, one can argue that the non-analytic structure of the non-planar contributions takes the form

$$\begin{equation}{\mathrm{lim}}_{\alpha \to {\alpha }_{\mathrm{c}}}\enspace {\mathcal{G}}_{\text{n.a.}}^{\left(h\right)}={g}_{h}{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{7{\chi }_{h}/6},\quad h\in \mathbb{N},\end{equation} \tag{ 5.34 }$$

where the h = 1 case is understood to be logarithmic. Near criticality, the essential difference [16, 44, 82] from the single matrix case (3.79) is that the critical exponent is now 7/6 instead of 5/4.

Double scaling limit. As for the single matrix mode, one can consider a double scaling limit in which the combination ${\kappa }^{-1}=N{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{7/6}$ is kept fixed as N tends to infinity. Though we do not provide the details, it is worth mentioning that the Painlevé I equation (4.32) uncovered in the double scaling limit of the single matrix model has an avatar for the two-matrix case under consideration. For the curious reader, we state that the analogous equation can be put in the following form [72, 73, 75]

$$\begin{equation}\delta r{\left(z\right)}^{3}-\delta r\left(z\right)\delta {r}^{{\prime\prime}}\left(z\right)-\frac{1}{2}\delta {r}^{\prime }{\left(z\right)}^{2}+\frac{2}{27}\enspace \delta {r}^{\left(4\right)}\left(z\right)-z=0.\end{equation} \tag{ 5.35 }$$

The assymptotic expansion is given by [63]

$$\begin{equation}\delta r\left(z\right)={z}^{1/3}\left(1+\sum _{n=1}^{\infty }{a}_{k}{z}^{-7n/3}\right),\end{equation} \tag{ 5.36 }$$

where the first few coefficients are given by

$$\begin{equation}\delta r\left(z\right)={z}^{1/3}\left(1-\frac{1}{18}\enspace {z}^{-7/3}-\frac{1925}{26\enspace 244}\enspace {z}^{-14/3}-\frac{509\enspace 575}{1417\enspace 176}\enspace {z}^{-21/3}-\frac{445\enspace 712\enspace 575}{114\enspace 791\enspace 256}\enspace {z}^{-28/3}+\cdots \enspace \right).\end{equation} \tag{ 5.37 }$$

As for the single-matrix case, the above expansion admits non-perturbative corrections of the type

$$\begin{equation}{\epsilon}\left(z\right)={c}_{1}\enspace {z}^{-1/4}\enspace {\mathrm{e}}^{-\frac{18}{7}\enspace {z}^{7/6}}+{c}_{2}\enspace {z}^{-1/4}\enspace {\mathrm{e}}^{-\frac{9\sqrt{2}}{7}\enspace {z}^{7/6}}.\end{equation} \tag{ 5.38 }$$

We note that the integration constants can not be fixed by the WKB analysis. The appearance of multiple solutions has been discussed in [64].

As a final note, we consider how the perturbative expansion in the 't Hooft limit is modified due to the appearance of a second matrix. Taking into account that the quadratic part of the exponent in the integral (5.1) contains a mixed term, we have the following set of propagators

$$\begin{equation}\langle {A}_{IK}{B}_{JL}\rangle =\frac{2c}{N\left(1-4{c}^{2}\right)}\enspace {\delta }_{IL}{\delta }_{JK},\qquad \langle {A}_{IK}{A}_{JL}\rangle =\langle {B}_{IK}{B}_{JL}\rangle =\frac{1}{N\left(1-4{c}^{2}\right)}\enspace {\delta }_{IL}{\delta }_{JK}.\end{equation} \tag{ 5.39 }$$

Thus, we have additional structure at each face indicating whether it has edges built from AA, BB or AB type propagators (figure 13).

Zoom In Zoom Out Reset image size

We might envision a continuum limit in which the vertices densely fill the discretised Riemann surface. In such a situation, the additional structure at the vertices may be viewed as an additional field taking two values living on each vertex of the continuous surface. More generally, we could have added an additional coupling weighting the two quartic vertices differently. This would generalise (5.11) to

$$\begin{equation}W\left(a,b\right)=\frac{1}{2}\left({a}^{2}+{b}^{2}\right)+\alpha \left(g\enspace {a}^{4}+{g}^{-1}\enspace {b}^{4}\right)-2c\enspace a\enspace b,\end{equation} \tag{ 5.40 }$$

where g > 0 parametrises the relative weight of the A and B vertices (figure 14).

Zoom In Zoom Out Reset image size

At this stage, we have no immediate reason to suspect that such a two-dimensional field theory is itself local. Nevertheless, and somewhat remarkably, we will discuss evidence for this in section 8 where we will also provide an interpretation for the three couplings in (5.40), namely α, c, and g.

Adding more matrices further decorates the Feynman diagrams. One might imagine a limit where we have an infinite chain of matrices such that they effectively carry a continuous label t. We will now consider precisely such a situation.

We begin by introducing the following class of classical actions

$$\begin{equation}{S}_{N}\left[M\left(t\right)\right]=N\enspace \mathrm{T}\mathrm{r}\int \mathrm{d}t\left(\frac{1}{2}\dot {M}{\left(t\right)}^{2}-V\left(M\left(t\right)\right)\right).\end{equation} \tag{ 6.1 }$$

Here, V(M(t)) is a polynomial potential governing the N × N Hermitian matrix valued path M(t), and the dot denotes a derivative with respect to time. We can view the above as a class of classical theories comprising N² degrees of freedom whose dynamical features are governed by the potential V(M(t)). Quantum mechanical transition amplitudes are given by the Feynman path integral

$$\begin{equation}{\mathcal{A}}_{N}\left({M}_{f},{M}_{i}\right)=\int \mathcal{D}M\left(t\right)\enspace {\mathrm{e}}^{\frac{i}{\hslash }{S}_{N}\left[M\left(t\right)\right]},\end{equation} \tag{ 6.2 }$$

with boundary conditions M(t_i) = M_i and M(t_f) = M_f, and our measure is now given by

$$\begin{equation}\mathcal{D}M\left(t\right)\equiv \prod _{t\in \mathbb{R}}\left[\mathrm{d}M\left(t\right)\right].\end{equation} \tag{ 6.3 }$$

The path integral (6.2) serves as the quantum mechanical generalisation of the matrix integrals previously explored.

As for the ordinary matrix integrals, we can exploit the global U(N) symmetry of (6.1) by parametrising M(t) as

$$\begin{equation}M\left(t\right)\to \mathcal{U}\left(t\right)\enspace {D}_{M}\left(t\right)\enspace \mathcal{U}{\left(t\right)}^{{\dagger}},\qquad {D}_{M}\left(t\right)=\mathrm{diag}\left({\lambda }_{1}\left(t\right),{\lambda }_{2}\left(t\right),\dots ,{\lambda }_{N}\left(t\right)\right).\end{equation} \tag{ 6.4 }$$

Here, $\mathcal{U}\left(t\right)$ is a time dependent element of U(N)/U(1)^N, and the λ_I(t) are the real valued time dependent eigenvalues of M(t). Under (6.4) the kinetic term in the action (6.1) is expressed as

$$\begin{equation}\mathrm{T}\mathrm{r}\enspace \dot {M}{\left(t\right)}^{2}=\mathrm{T}\mathrm{r}\enspace {\dot {D}}_{M}{\left(t\right)}^{2}+\mathrm{T}\mathrm{r}{\left[{D}_{M}\left(t\right),\dot {\mathcal{U}}\left(t\right)\mathcal{U}{\left(t\right)}^{{\dagger}}\right]}^{2}.\end{equation} \tag{ 6.5 }$$

We note that $\dot {\mathcal{U}}\left(t\right)\enspace \mathcal{U}{\left(t\right)}^{{\dagger}}$ is again a Hermitian matrix which can be expressed as

$$\begin{equation}\dot {\mathcal{U}}\left(t\right)\enspace \mathcal{U}{\left(t\right)}^{{\dagger}}=\frac{\mathrm{i}}{\sqrt{2}}\sum _{a=1}^{N\left(N-1\right)/2}\left({T}_{a}^{S}\enspace {\dot {\beta }}_{a}\left(t\right)+{T}_{a}^{A}\enspace {\dot {\gamma }}_{a}\left(t\right)\right),\end{equation} \tag{ 6.6 }$$

where ${T}_{a}^{S}$ and ${T}_{a}^{A}$ are a basis of real symmetric, and pure imaginary antisymmetric N × N generators of U(N)/U(1)^N, with normalisation $\mathrm{T}\mathrm{r}\enspace {T}_{a}^{{\epsilon}}{T}_{b}^{{{\epsilon}}^{\prime }}={\delta }_{ab}\enspace {\delta }^{{\epsilon}{{\epsilon}}^{\prime }}$ with , ' ∈ {S, A}. Thus,

$$\begin{equation}\mathrm{T}\mathrm{r}\enspace {\dot {M}}^{2}=\mathrm{T}\mathrm{r}\enspace {\dot {D}}_{M}^{2}+\frac{1}{2}\sum _{I{< }J}{\left({\lambda }_{I}-{\lambda }_{J}\right)}^{2}\left({\dot {\beta }}_{IJ}^{2}+{\dot {\gamma }}_{IJ}^{2}\right),\end{equation} \tag{ 6.7 }$$

where for notational convenience we suppressed the explicit time dependence, β_IJ are the matrix elements of ${\sum }_{a}{T}_{a}^{S}\enspace {\dot {\beta }}_{a}$ and analogously γ_IJ the matrix elements of ${\sum }_{a}\enspace {T}_{a}^{A}\enspace {\dot {\gamma }}_{a}$.

Turning to the Hamiltonian formalism, we calculate the canonical momenta for λ, β_IJ and γ_IJ respectively:

$$\begin{equation}{\pi }_{I}=N{\dot {\lambda }}_{I},\qquad {\pi }_{IJ}=\frac{N}{2}{\left({\lambda }_{I}-{\lambda }_{J}\right)}^{2}\enspace {\dot {\beta }}_{IJ},\qquad {\tilde {\pi }}_{IJ}=\frac{N}{2}\enspace {\left({\lambda }_{I}-{\lambda }_{J}\right)}^{2}\enspace {\dot {\gamma }}_{IJ}.\end{equation} \tag{ 6.8 }$$

From these, we can construct the classical Hamiltonian

$$\begin{equation}H=\sum _{I=1}^{N}\left(\frac{1}{2N}\enspace {\pi }_{I}^{2}+NV\left({\lambda }_{I}\right)\right)+\frac{1}{N}\sum _{I{< }J}\frac{{\pi }_{IJ}^{2}+{\tilde {\pi }}_{IJ}^{2}}{{\left({\lambda }_{I}-{\lambda }_{J}\right)}^{2}}.\end{equation} \tag{ 6.9 }$$

At the classical level, the lowest energy configuration lies at the minimum of V(λ_I) with all classical momenta vanishing. As we shall soon see, the quantum mechanical ground state differs considerably from its classical counterpart.

It is straightforward to promote (6.9) to a quantum Hamiltonian. Recall that in the eigenvalue basis (6.4), the N² degrees of freedom are living on a curved space. Moreover, this curved space is merely a coordinate transformation of the flat metric on ${\mathbb{R}}^{{N}^{2}}$. Consequently, the quantum kinetic term is nothing other than the Laplace–Beltrami differential operator associated to our curved space. The quantum version of the Hamiltonian in (6.9) is consequently given by

$$\begin{equation}\hat{H}=\sum _{I=1}^{N}\left(-\frac{1}{2N}\frac{1}{{{\Delta}}_{N}\left(\lambda \right)}\frac{{\partial }^{2}}{\partial {\lambda }_{I}^{2}}\enspace {{\Delta}}_{N}\left(\lambda \right)+NV\left({\lambda }_{I}\right)\right)+\frac{1}{N}\sum _{I{< }J}\frac{{\hat{\pi }}_{IJ}^{2}+{\hat{\tilde {\pi }}}_{IJ}^{2}}{{\left({\lambda }_{I}-{\lambda }_{J}\right)}^{2}},\end{equation} \tag{ 6.10 }$$

where we have explicitly written the eigenvalue dependence for the Laplacian, we are working in units for which ℏ = 1 and Δ_N(λ) is the determinant of the Vandermonde matrix introduced in (3.11). The angular part in (6.10) describes the motion on the compact coset space U(N)/U(1)^N.⁶ Due to the compactness of the coset space, the eigenfunctions of H can be annihilated by ${\hat{\pi }}_{IJ}$ and ${\hat{\tilde {\pi }}}_{IJ}$ while retaining normalisability. Thus, it is natural to study wavefunctions Ψ(λ_I) that are functions only of the eigenvalues, but not the coordinates on the compact coset space U(N)/U(1)^N. Acting on such wavefunctions, the Hamiltonian is reduced to

$$\begin{equation}\hat{H}=\sum _{I=1}^{N}\left(-\frac{1}{2N}\frac{1}{{{\Delta}}_{N}\left(\lambda \right)}\frac{{\partial }^{2}}{\partial {\lambda }_{I}^{2}}\enspace {{\Delta}}_{N}\left(\lambda \right)+NV\left({\lambda }_{I}\right)\right).\end{equation} \tag{ 6.11 }$$

A subsequent redefinition of the wavefunctions

$$\begin{equation}{\Psi}\left({\lambda }_{1},\dots ,{\lambda }_{N}\right)=\frac{\psi \left({\lambda }_{1},\dots ,{\lambda }_{N}\right)}{{{\Delta}}_{N}\left(\lambda \right)},\end{equation} \tag{ 6.12 }$$

leads to a Hamiltonian with standard kinetic term. Thus, we have reduced our eigenvalue problem to the following single particle problem

$$\begin{equation}\left(-\frac{1}{2N}\frac{{\partial }^{2}}{\partial {\lambda }^{2}}+NV\left(\lambda \right)\right){\psi }_{{\varepsilon }_{n}}\left(\lambda \right)=N{\varepsilon }_{n}{\psi }_{{\varepsilon }_{n}}\left(\lambda \right),\quad n=1,2,\dots ,N.\end{equation} \tag{ 6.13 }$$

Recall that Δ_N(λ) is antisymmetric under exchange of any two eigenvalues λ_I. It then follows from (6.12) that a general state Ψ invariant under the permutation subgroup S_N ⊂ U(N) enforces the rescaled wavefunction ψ to be an anti-symmetric function of the λ_I. In particular, the problem has been reduced to solving a system of N free particles subject to the Pauli exclusion principle [39]. It is already clear at this stage that the quantum mechanical state, particularly at large N, differs substantially from its classical counterpart. In other words, quantum effects play a dominant role in the description of the system.

The ground state of the matrix quantum mechanical systems is described by N free fermions, each governed by the potential V(λ). The fermions fill up the first N energy levels, all the way up to the Fermi level ɛ_F ≡ ɛ_N. We will be interested in small energy excitations or ripples above the filled Fermi sea. Rather generally, in the large N limit, the low energy physics near a filled Fermi sea has a linear dispersion relation in two-dimensions [97].

We now assume that the potential has a maximum at some value of λ. Recall that in our study of ordinary matrix integrals, we uncovered non-analytic behaviour for particular choices of the polynomial V(λ). For instance, for V_α(λ) = λ²/2 + αλ⁴ we uncovered non-analytic behaviour for a negative value of α = α_c < 0. For negative values of α, V_α(λ) contains a maximum. Thus, in tuning the value of α beyond α_c, the confining effect of the quadratic term in the Lagrangian will no longer be strong enough to compete with the repulsion of eigenvalues from both the Vandermonde and the quartic contribution to V_α(λ). Recall further that for the matrix integrals the non-analyticities were important in obtaining a continuum limit with a divergent number of vertices in the discretised Riemann surface.

It is natural, then, to ask what kind of non-analyticities might we expect from the free fermions stemming from quantum mechanical matrices? As a simple test observable, we can consider the Fermi energy for a potential given by

$$\begin{equation}V\left(\lambda \right)=\frac{1}{2}{\gamma }^{2}{\lambda }^{2}-\alpha {\lambda }^{4}.\end{equation} \tag{ 6.14 }$$

We have maxima at ${\lambda }_{{\pm}}={\pm}\gamma {\left(2\sqrt{\alpha }\right)}^{-1}$, such that V(λ_±) = γ⁴(16α)⁻¹. Given a fixed number, N, of fermions we can consider how the Fermi energy ɛ_F behaves as a function of the dimensionless parameter $\tilde {\alpha }\equiv \alpha \enspace {\gamma }^{-3}$. As we vary $\tilde {\alpha }$, the Fermi level will begin to approach the maximum of the potential. In the large N limit we can resort to a semi-classical approximation. This approximation can be cast in the form of a Bohr–Sommerfeld condition

$$\begin{equation}{\int }_{0}^{{\lambda }_{+}}\mathrm{d}\lambda \sqrt{{\varepsilon }_{\mathrm{F}}-V\left(\lambda \right)}=\frac{\pi }{2\sqrt{2}}\left(1+\frac{1}{2N}\right).\end{equation} \tag{ 6.15 }$$

Examining the above integral reveals that upon tuning ɛ_F to the maximum value of the potential gives rise to non-analytic behaviour. This stems from the part of the integral near the maxima. The same non-analyticity is captured by the following, simpler, integral which we evaluate in a small δɛ_F-expansion [88–91]:

$$\begin{equation}{\int }_{0}^{{\lambda }_{+}}\mathrm{d}\lambda \sqrt{\delta {\varepsilon }_{\mathrm{F}}-{V}^{{\prime\prime}}\left({\lambda }_{+}\right){\left(\lambda -{\lambda }_{+}\right)}^{2}/2}=-\frac{\delta {\varepsilon }_{\mathrm{F}}}{4\gamma }\enspace \mathrm{log}\enspace \tilde {\alpha }\frac{\delta {\varepsilon }_{\mathrm{F}}}{\gamma }+\text{analytic}.\end{equation} \tag{ 6.16 }$$

We should view −δɛ_F as the difference between the Fermi level and the maximum value of the potential. So long as the Fermi level is close to the top of the potential the detailed features of the potential will not affect the preceding analysis.

For the sake of generality, it is perhaps worth pointing out that tuning the potential to have a non-quadratic maximum would affect the non-analytic behaviour. For instance

$$\begin{equation}{\int }_{0}^{{\lambda }_{+}}\mathrm{d}\lambda \sqrt{\delta {\varepsilon }_{\mathrm{F}}+{\left(\lambda -{\lambda }_{+}\right)}^{m}}=\sqrt{\delta {\varepsilon }_{\mathrm{F}}}{\lambda }_{+}{{\times}}_{2}{F}_{1}\left(-\frac{1}{2};\frac{1}{m};1+\frac{1}{m};-\frac{{\left(-{\lambda }_{+}\right)}^{m}}{\delta {\varepsilon }_{\mathrm{F}}}\right),\end{equation} \tag{ 6.17 }$$

from which we extract a non-analytic behaviour of the type $\sim \delta {\varepsilon }_{\mathrm{F}}^{1/2+1/m}$ for m ≠ 2 [89, 91]. In what follows we focus on the m = 2 case.

A final (albeit somewhat a posteriori) motivation for assuming a local maximum is that we will eventually compare calculations for the matrix quantum mechanical theory to a certain worldsheet string theory. The natural observable from the worldsheet perspective is an S-matrix. In the presence of a local maximum, the low energy excitations naturally encode an S-matrix. This is most evident, if we consider the physics localised near the maximum, where the model reduces to N fermions in an inverse harmonic well. Let us introduce the following dimensionless quantities

$$\begin{equation}x={2}^{3/4}{N}^{1/2}\enspace \left(\lambda -{\lambda }_{+}\right)\enspace {\gamma }^{1/2},\qquad \nu =\frac{N}{{2}^{1/2}}\enspace \left(\frac{{\varepsilon }_{\mathrm{F}}+\delta \varepsilon -V\left({\lambda }_{+}\right)}{\gamma }\right),\end{equation} \tag{ 6.18 }$$

and consider excitations with energy near ɛ such that at large N we have δɛ/ɛ_F ≪ 1. The single particle Schrödinger equation (6.13) near the maximum of the potential now becomes

$$\begin{equation}\left(-{\partial }_{x}^{2}-\frac{{x}^{2}}{4}\right){\psi }_{\nu ,p}=\nu \enspace {\psi }_{\nu ,p}.\end{equation} \tag{ 6.19 }$$

Note further that the rescaling (6.18) involves N in such a way that the range of x and ν becomes effectively infinite as N tends to infinity. The parity index in ψ_ν,p, p ∈ {±}, reflects the behaviour of the wavefunction under the discrete symmetry x → −x, namely ψ_ν,p(−x) = pψ_ν,p(x).

When considering the scattering problem, we take the Fermi level ν_F of our new Schrödinger problem (6.19) to be negative, such that we can send in an incoming wave from negative x, let it bounce against the wall, and compute the reflected wavefunction. The Schrödinger equation (6.19) admits solutions whose explicit form is given by the parabolic cylinder functions D_a(z). Explicitly,

$$\begin{equation}{{\Psi}}_{\nu ,{\pm}}\left(x\right)={\mathrm{e}}^{-\mathrm{i}\pi \left(1/2-\mathrm{i}\nu \right)/4}\left({D}_{\mathrm{i}\left(\nu +\mathrm{i}/2\right)}\left({\mathrm{e}}^{\frac{3\pi \mathrm{i}}{4}}x\right){\pm}{D}_{\mathrm{i}\left(\nu +\mathrm{i}/2\right)}\left({\mathrm{e}}^{-\frac{\pi \mathrm{i}}{4}}x\right)\right).\end{equation} \tag{ 6.20 }$$

It will also prove convenient to express the above wavefunctions in a basis of incoming/outgoing waves near large and negative values of x. These can be expressed as follows

$$\begin{equation}{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right)=\frac{1}{2}\left({{\Psi}}_{\nu ,+}\left(x\right)+{{\Psi}}_{\nu ,-}\left(x\right)\right),\qquad {\psi }_{\nu }^{\left(\mathrm{i}\mathrm{n}\right)}\left(x\right)={\left({\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right)\right)}^{{\ast}}.\end{equation} \tag{ 6.21 }$$

The outgoing wavefunctions admit an asymptotic expansion of the following form

$$\begin{equation}{\mathrm{lim}}_{x\to -\infty }\enspace {\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right)={\mathrm{e}}^{\mathrm{i}{x}^{2}/4}{\left(-x\right)}^{-1/2+\mathrm{i}\nu }+\cdots \enspace .\end{equation} \tag{ 6.22 }$$

Notice that in addition to the energy, ν parameterises the wavelength of the spatial oscillations. It is easiest to see this by considering a coordinate transformation to the spatial coordinate x = −e^−u. The Wronskian for the above wavefunctions is given by

$$\begin{equation}W\left(\nu \right)={\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right){\partial }_{x}{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right)-{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right){\partial }_{x}{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right)=\mathrm{i}\frac{\sqrt{2\pi }}{{\Gamma}\left(1/2-\mathrm{i}\nu \right)}\enspace {\mathrm{e}}^{-\pi \nu /2}.\end{equation} \tag{ 6.23 }$$

As expected, W(ν) is independent of x. We can also express the parity eigenstates in terms of ψ^(out)(x) as follows

$$\begin{equation}{\psi }_{\nu ,+}\left(x\right)\equiv \frac{1}{\sqrt{2}}\enspace \frac{{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right)+{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right)}{W\left(\nu \right)},\qquad {\psi }_{\nu ,-}\left(x\right)\equiv \frac{1}{\sqrt{2}}\enspace \frac{{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right)-{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right)}{W{\left(\nu \right)}^{{\ast}}},\end{equation} \tag{ 6.24 }$$

where we have rescaled (6.20) with respect to the Wronskian.

By virtue of the discrete symmetry x ↔ −x, ${\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right)$ and its complex conjugate also describe solutions to (6.19). These are waves that are purely incoming or outgoing on the other side of the potential barrier.

So far we kept the particle number N fixed. It will be convenient to also allow the particle number to vary, thus constructing a second quantised picture of the quantum mechanical theory. To do so, we introduce lowering and raising operators, ${\hat{a}}_{\nu ,p}$ and ${\hat{a}}_{\nu ,p}^{{\dagger}}$ satisfying the standard anti-commutation relations

$$\begin{equation}\left\{{\hat{a}}_{\nu ,p},{\hat{a}}_{{\nu }^{\prime },{p}^{\prime }}\right\}=0,\qquad \left\{{\hat{a}}_{\nu ,p},{\hat{a}}_{{\nu }^{\prime },{p}^{\prime }}^{{\dagger}}\right\}=2\pi \enspace {\delta }_{\nu ,{\nu }^{\prime }}{\delta }_{p,{p}^{\prime }}.\end{equation} \tag{ 6.25 }$$

The ground state is filled until the Fermi level and we introduce the state |Fermi⟩ corresponding to a filled Fermi sea by

$$\begin{align}\hfill {\hat{a}}_{\nu ,p}\enspace \vert \text{Fermi}\rangle & =0,\quad \nu { >}{\nu }_{\mathrm{F}}\quad \text{and}\quad p\in \left\{{\pm}\right\},\hfill \\ \hfill {\hat{a}}_{\nu ,p}^{{\dagger}}\enspace \vert \text{Fermi}\rangle & =0,\quad \nu {< }{\nu }_{\mathrm{F}}\quad \text{and}\quad p\in \left\{{\pm}\right\},\hfill \end{align} \tag{ 6.26 }$$

where we recall that ν_F < 0. The above Fermi sea is filled on both sides of the potential.⁷ We can create multi-particle energy eigenstates on top of |Fermi⟩ by acting with ${\hat{a}}_{\nu ,p}^{{\dagger}}$ with ν > ν_F, and hole type excitations by acting with â_ν,p with ν < ν_F.

We can thus introduce field operators for the creation and annihilation of wavefunctions at x:

$$\begin{equation}{\hat{{\Psi}}}_{p}\left(t,x\right)={\int }_{\mathbb{R}}\frac{\mathrm{d}\nu }{2\pi }\enspace \enspace {\mathrm{e}}^{-\mathrm{i}\nu t}\enspace {\hat{a}}_{\nu ,p}\enspace {\psi }_{\nu ,p}^{{\ast}}\left(x\right),\qquad {\hat{{\Psi}}}_{p}^{{\dagger}}\left(t,x\right)={\int }_{\mathbb{R}}\frac{\mathrm{d}\nu }{2\pi }\enspace {\mathrm{e}}^{\mathrm{i}\nu t}\enspace {\hat{a}}_{\nu ,p}^{{\dagger}}\enspace {\psi }_{\nu ,p}\left(x\right),\end{equation} \tag{ 6.27 }$$

where the integration is over the continuous energy levels ν and ${\psi }_{\nu ,p}^{{\ast}}\left(x\right)$ and ψ_ν,p(x) are the parity eigenstates introduced in (6.24). In the expression above ∗ denotes the complex conjugation of $\mathbb{C}$-numbers, while † denotes the complex conjugation of quantum operators. Given the above fermionic operators, we can define a number density operator

$$\begin{equation}\hat{n}\left(t,x\right)={\hat{{\Psi}}}^{{\dagger}}\left(t,x\right)\hat{{\Psi}}\left(t,x\right).\end{equation} \tag{ 6.28 }$$

In the above, we have defined $\hat{{\Psi}}\left(t,x\right)\equiv {\hat{{\Psi}}}_{+}\left(t,x\right)+{\hat{{\Psi}}}_{-}\left(t,x\right)$ and its conjugate, satisfying

$$\begin{equation}\left\{\hat{{\Psi}}\left(t,x\right),{\hat{{\Psi}}}^{{\dagger}}\left(t,{x}^{\prime }\right)\right\}=2\enspace \delta \left(x-{x}^{\prime }\right).\end{equation} \tag{ 6.29 }$$

The operator ${\hat{{\Psi}}}^{{\dagger}}\left(t,x\right)$ creates states that can be perceived as single particle states on top of the filled Fermi sea when viewed near the left boundary at large negative values of x.

6.5.1. Feynman propagator

Given a second quantised theory we can construct various propagators for the fermions. For instance, the Feynman propagator is given by

$$\begin{equation}{S}_{\mathrm{F}}\left(t,x;{t}^{\prime },{x}^{\prime }\right)=-\mathrm{i}\enspace \langle \text{Fermi}\vert \enspace {T}_{\mathrm{F}}\enspace \hat{{\Psi}}\left(t,x\right){\hat{{\Psi}}}^{{\dagger}}\left({t}^{\prime },{x}^{\prime }\right)\vert \text{Fermi}\rangle ,\end{equation} \tag{ 6.30 }$$

where T_F denotes Fermionic time-ordering. It satisfies

$$\begin{equation}\left(-{\partial }_{x}^{2}-\frac{{x}^{2}}{4}-\mathrm{i}{\partial }_{t}\right){S}_{\mathrm{F}}\left(t,x;{t}^{\prime },{x}^{\prime }\right)=\delta \left(t-{t}^{\prime }\right)\delta \left(x-{x}^{\prime }\right).\end{equation} \tag{ 6.31 }$$

More explicitly, we can express S_F(t, x; t', x') as

$$\begin{equation}\mathrm{i}\enspace {S}_{\mathrm{F}}\left(t,x ;{t}^{\prime },{x}^{\prime }\right)={\Theta}\left(s\right){\int }_{{\nu }_{\mathrm{F}}}^{\infty }\frac{\mathrm{d}\nu }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}\nu s}\enspace {\psi }_{\nu ,p}^{{\ast}}\left(x\right){\psi }_{\nu ,p}\left({x}^{\prime }\right)-{\Theta}\left(-s\right){\int }_{-\infty }^{{\nu }_{\mathrm{F}}}\frac{\mathrm{d}\nu }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}\nu s}\enspace {\psi }_{\nu ,p}^{{\ast}}\left(x\right){\psi }_{\nu ,p}\left({x}^{\prime }\right),\end{equation} \tag{ 6.32 }$$

where we sum over the parity index, and for notational simplicity we have defined s ≡ t − t'. We can also write down the Fourier transform of the Fermion propagator with respect to s. For this we use the contour integral representation of the Heaviside function

$$\begin{equation}{\Theta}\left(s\right)={\mathrm{lim}}_{\varepsilon \to {0}^{+}}\enspace {\int }_{\mathbb{R}}\frac{\mathrm{d}{\Omega}}{2\pi \mathrm{i}}\frac{{\mathrm{e}}^{\mathrm{i}{\Omega}s}}{{\Omega}-\mathrm{i}\varepsilon }.\end{equation} \tag{ 6.33 }$$

The Fourier transform then reads

$$\begin{equation}{\tilde {S}}_{\mathrm{F}}\left(x,{x}^{\prime };\omega \right)=-{\mathrm{lim}}_{\varepsilon \to {0}^{+}}\enspace {\int }_{\mathbb{R}}\frac{\mathrm{d}\nu }{2\pi }\enspace {\psi }_{\nu ,p}^{{\ast}}\left(x\right){\psi }_{\nu ,p}\left({x}^{\prime }\right)\left(\enspace \frac{{\Theta}\left(\nu -{\nu }_{\mathrm{F}}\right)}{\nu -\omega -\mathrm{i}\varepsilon }+\frac{{\Theta}\left({\nu }_{\mathrm{F}}-\nu \right)}{\nu -\omega +\mathrm{i}\varepsilon }\enspace \right),\end{equation} \tag{ 6.34 }$$

and we use the following conventions for the Fourier transform:

$$\begin{equation}f\left(t\right)={\int }_{\mathbb{R}}\frac{\mathrm{d}\omega }{2\pi }\enspace {\mathrm{e}}^{-\mathrm{i}\omega t}\tilde {f}\left(\omega \right)\quad \text{and}\quad \tilde {f}\left(\omega \right)={\int }_{\mathbb{R}}\mathrm{d}t\enspace {\mathrm{e}}^{+\mathrm{i}\omega t}\enspace f\left(t\right).\end{equation} \tag{ 6.35 }$$

Complex conjugation at the level of the Fourier transform mirrors itself in x ↔ x' and a particle hole exchange. Finally, we note that for large negative values of x and x' the parity states behave as

$$\begin{align}\hfill & {\mathrm{lim}}_{x,{x}^{\prime }\to -\infty }\left({\psi }_{+,\nu }^{{\ast}}\left(x\right){\psi }_{+,\nu }\left({x}^{\prime }\right)+{\psi }_{-,\nu }^{{\ast}}\left(x\right){\psi }_{-,\nu }\left({x}^{\prime }\right)\right)\hfill \\ \hfill & \quad =\frac{1}{\sqrt{x{x}^{\prime }}}\left({\mathrm{e}}^{-\mathrm{i}\left({x}^{2}-{x}^{\prime 2}\right)/4-\mathrm{i}\nu \mathrm{log}x/{x}^{\prime }}+{R}_{\nu }^{{\ast}}\enspace {\mathrm{e}}^{-\mathrm{i}\left({x}^{2}+{x}^{\prime 2}\right)/4-\mathrm{i}\nu \mathrm{log}x{x}^{\prime }}\right)+\text{h.c.},\hfill \end{align} \tag{ 6.36 }$$

where we defined

$$\begin{equation}{R}_{\nu }\equiv -\mathrm{i}\enspace \frac{{\Gamma}\left(1/2-\mathrm{i}\nu \right)}{\sqrt{2\pi }}\enspace {\mathrm{e}}^{-\pi \nu /2}.\end{equation} \tag{ 6.37 }$$

We will provide a physical interpretation of the definition of R_ν in the next section. For large negative values of x and x', with x > x', we then find the asymptotic form for (6.36):

$$\begin{align}\hfill {\tilde {S}}_{\mathrm{F}}\left(x,{x}^{\prime } ;\omega \right)& \approx -\frac{\mathrm{i}}{\sqrt{x{x}^{\prime }}}\left({\mathrm{e}}^{\mathrm{i}\left({x}^{2}-{x}^{\prime 2}\right)/4+\mathrm{i}\omega \mathrm{log}x/{x}^{\prime }}+{R}_{\omega }\enspace {\mathrm{e}}^{\mathrm{i}\left({x}^{2}+{x}^{\prime 2}\right)/4+\mathrm{i}\omega \mathrm{log}x{x}^{\prime }}\right){\Theta}\left(\omega -{\nu }_{\mathrm{F}}\right)\hfill \\ \hfill & +\frac{\mathrm{i}}{\sqrt{x{x}^{\prime }}}\left({\mathrm{e}}^{-\mathrm{i}\left({x}^{2}-{x}^{\prime 2}\right)/4-\mathrm{i}\omega \mathrm{log}x/{x}^{\prime }}+{R}_{\omega }^{{\ast}}\enspace {\mathrm{e}}^{-\mathrm{i}\left({x}^{2}+{x}^{\prime 2}\right)/4-\mathrm{i}\omega \mathrm{log}x{x}^{\prime }}\right){\Theta}\left({\nu }_{\mathrm{F}}-\omega \right).\hfill \end{align} \tag{ 6.38 }$$

For x' > x, both large and negative we find

$$\begin{align}\hfill {\tilde {S}}_{\mathrm{F}}\left(x,{x}^{\prime } ;\omega \right)& \approx -\frac{\mathrm{i}}{\sqrt{x{x}^{\prime }}}\left({\mathrm{e}}^{-\mathrm{i}\left({x}^{2}-{x}^{\prime 2}\right)/4-\mathrm{i}\omega \mathrm{log}x/{x}^{\prime }}+{R}_{\omega }\enspace {\mathrm{e}}^{\mathrm{i}\left({x}^{2}+{x}^{\prime 2}\right)/4+\mathrm{i}\omega \mathrm{log}x{x}^{\prime }}\right){\Theta}\left(\omega -{\nu }_{\mathrm{F}}\right)\hfill \\ \hfill & +\frac{\mathrm{i}}{\sqrt{x{x}^{\prime }}}\left({\mathrm{e}}^{\mathrm{i}\left({x}^{2}-{x}^{\prime 2}\right)/4+\mathrm{i}\omega \mathrm{log}x/{x}^{\prime }}+{R}_{\omega }^{{\ast}}\enspace {\mathrm{e}}^{-\mathrm{i}\left({x}^{2}+{x}^{\prime 2}\right)/4-\mathrm{i}\omega \mathrm{log}x{x}^{\prime }}\right){\Theta}\left({\nu }_{\mathrm{F}}-\omega \right).\hfill \end{align} \tag{ 6.39 }$$

We will now put all these expressions to good use, as they have provided us the basic building blocks to construct an S-matrix for our quantum mechanical matrices.

Recall the Hamiltonian (6.19) we obtained at the end of the previous section

$$\begin{equation}\hat{H}=-{\partial }_{x}^{2}-\frac{{x}^{2}}{4}.\end{equation} \tag{ 7.1 }$$

A natural object to consider is given by the reflection coefficient, R_ν, of a wave of frequency ν coming in from the asymptotic region at large and negative values of x and reflecting back, as depicted in figure 15.

Zoom In Zoom Out Reset image size

As our boundary condition we impose that the incoming flux from large positive values of x vanishes. In this regard, it is convenient to note that ${\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right)$ is a wave that is purely outgoing for large positive values of x. Expanding for large negative values of x we find a superposition of an incoming and reflected wave

$$\begin{equation}\frac{{\psi }_{\nu }^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right)}{W\left(\nu \right)}\approx \frac{{\mathrm{e}}^{-\mathrm{i}{x}^{2}/4-\mathrm{i}\nu \mathrm{log}\left(-x\right)}}{\sqrt{-x}}-\mathrm{i}\enspace \frac{{\Gamma}\left(1/2-\mathrm{i}\nu \right)}{\sqrt{-2\pi x}}\enspace {\mathrm{e}}^{-\pi \nu /2}\enspace {\mathrm{e}}^{\mathrm{i}{x}^{2}/4+\mathrm{i}\nu \mathrm{log}\left(-x\right)}+\cdots \enspace ,\end{equation} \tag{ 7.2 }$$

where we have normalised the wavefunction such that the incoming (first) part has unit coefficient. To obtain the reflection and transmission coefficient we calculate the probability current

$$\begin{equation}\mathcal{J}\left(t,x\right)=-\mathrm{i}\left({\psi }^{{\ast}}{\partial }_{x}\psi -\psi {\partial }_{x}{\psi }^{{\ast}}\right),\qquad {\partial }_{t}\rho +{\partial }_{x}\mathcal{J}=0.\end{equation} \tag{ 7.3 }$$

From the probability current one can readily obtain the reflection and transmission coefficients

$$\begin{equation}\vert {R}_{\nu }{\vert }^{2}=\frac{\vert {\mathcal{J}}_{\mathrm{r}\mathrm{e}\mathrm{f}\vert }}{\vert {\mathcal{J}}_{\mathrm{i}\mathrm{n}}\vert }=\frac{1}{1+{\mathrm{e}}^{2\pi \nu }},\qquad \vert {T}_{\nu }{\vert }^{2}=\frac{\vert {\mathcal{J}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}\vert }{\vert {\mathcal{J}}_{\mathrm{i}\mathrm{n}}\vert }=\frac{1}{1+{\mathrm{e}}^{-2\pi \nu }}=1-\vert {R}_{\nu }{\vert }^{2},\end{equation} \tag{ 7.4 }$$

where we obtain ${\mathcal{J}}_{\mathrm{i}\mathrm{n}}$ from (7.3) by using the purely incoming part of (7.2) and similarly the purely reflected and transmitted wave for ${\mathcal{J}}_{\mathrm{r}\mathrm{e}\mathrm{f}}$ and ${\mathcal{J}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$ respectively. For negative values of ν, which is the situation we are most interested in, the transmission coefficient is exponentially suppressed. In fact, the expression (7.2) encodes more information. For instance, it encodes the reflection coefficient

$$\begin{equation}{R}_{\nu }=-\mathrm{i}\enspace \frac{{\Gamma}\left(1/2-\mathrm{i}\nu \right)}{\sqrt{2\pi }}\enspace {\mathrm{e}}^{-\pi \nu /2},\end{equation} \tag{ 7.5 }$$

which can be viewed as an S-matrix element for the scattering of a single wave against the potential barrier. Perturbative unitarity of the S-matrix is given by the fact that R_ν is a pure phase up to e^−πν/2 corrections. Using the Gamma function identity

$$\begin{equation}\vert {\Gamma}\left(1/2-\mathrm{i}\nu \right){\vert }^{2}=\frac{\pi }{\mathrm{cosh}\left(\pi \nu \right)},\end{equation} \tag{ 7.6 }$$

the absolute value of (7.5) reduces to (7.4).

The Green's function is given by evaluating the matrix elements of the resolvent

$$\begin{equation}G\left(x,{x}^{\prime } ;z\right)=\langle {x}^{\prime }\vert {\left(\hat{H}-z\right)}^{-1}\vert x\rangle ,\qquad {\left(G\left(x,{x}^{\prime } ;z\right)\right)}^{{\ast}}=G\left({x}^{\prime },x ;{z}^{{\ast}}\right),\end{equation} \tag{ 7.7 }$$

where z takes values across the whole complex plane, modulo the spectrum of $\hat{H}$. The Green's function satisfies the equation

$$\begin{equation}\left(-{\partial }_{x}^{2}-\frac{{x}^{2}}{4}-z\right)G\left(x,{x}^{\prime } ;z\right)=\delta \left(x-{x}^{\prime }\right).\end{equation} \tag{ 7.8 }$$

It is often convenient to express G(x, x'; z) in terms of eigenfunctions of $\hat{H}$ with complex eigenvalue z. Concretely,

$$\begin{equation}G\left(x,{x}^{\prime } ;z\right)=\frac{1}{W\left(z\right)}\left({\psi }_{z}^{R}\left(x\right){\psi }_{z}^{L}\left({x}^{\prime }\right){\Theta}\left({x}^{\prime }-x\right)+{\psi }_{z}^{R}\left({x}^{\prime }\right){\psi }_{z}^{L}\left(x\right){\Theta}\left(x-{x}^{\prime }\right)\right),\end{equation} \tag{ 7.9 }$$

where ${\psi }_{z}^{L}\left(x\right)$ decays at x = −∞, and ${\psi }_{z}^{R}\left(x\right)$ decays at x = ∞, and W(z) is the Wronskian (6.23). As an example, when z takes values in the upper half of the complex plane we can express ${\psi }_{z}^{L,R}\left(x\right)$ as

$$\begin{equation}{\psi }_{z}^{L}\left(x\right)={\psi }_{z}^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(x\right),\qquad {\psi }_{z}^{R}\left(x\right)={\psi }_{z}^{\left(\mathrm{o}\mathrm{u}\mathrm{t}\right)}\left(-x\right).\end{equation} \tag{ 7.10 }$$

For large negative and large positive values of x these wave functions scale as

$$\begin{equation}{\mathrm{lim}}_{x\to -\infty }\enspace {\psi }_{z}^{L}\left(x\right)=\frac{1}{\sqrt{-x}}\enspace {\mathrm{e}}^{\mathrm{i}{x}^{2}/4+\mathrm{i}z\mathrm{log}\left(-x\right)}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad {\mathrm{lim}}_{x\to \infty }\enspace {\psi }_{z}^{R}\left(x\right)=\frac{1}{\sqrt{x}}\enspace {\mathrm{e}}^{\mathrm{i}{x}^{2}/4+\mathrm{i}z\mathrm{log}x}.\end{equation} \tag{ 7.11 }$$

Using the expressions for the wavefunctions in (7.2), we can obtain an asymptotic expansion of G(x, x'; z) for large and negative values of x and x'. For z in the upper half plane we obtain

$$\begin{align}\hfill G\left(x,{x}^{\prime } ;z\right)\approx & -\frac{\mathrm{i}}{\sqrt{x{x}^{\prime }}}\left({\mathrm{e}}^{-\mathrm{i}\left({x}^{2}-{x}^{\prime 2}\right)/4-\mathrm{i}z\mathrm{log}x/{x}^{\prime }}+{R}_{z}\enspace {\mathrm{e}}^{\mathrm{i}\left({x}^{2}+{x}^{\prime 2}\right)/4+\mathrm{i}z\mathrm{log}x{x}^{\prime }}\right){\Theta}\left({x}^{\prime }-x\right)\hfill \\ \hfill & -\frac{\mathrm{i}}{\sqrt{x{x}^{\prime }}}\left({\mathrm{e}}^{\mathrm{i}\left({x}^{2}-{x}^{\prime 2}\right)/4+\mathrm{i}z\mathrm{log}x/{x}^{\prime }}+{R}_{z}\enspace {\mathrm{e}}^{\mathrm{i}\left({x}^{2}+{x}^{\prime 2}\right)/4+\mathrm{i}z\mathrm{log}x{x}^{\prime }}\right){\Theta}\left(x-{x}^{\prime }\right).\hfill \end{align} \tag{ 7.12 }$$

A similar expression can be derived for z in the lower half plane. The reflection coefficient R_ν given in (7.5) and the above Green's function are the basic building blocks for the scattering problem of interest—namely, the scattering of the probability density ρ(t, x) = ψ*(t, x)ψ(t, x) about the filled Fermi sea. We note that ρ(t, x) is a quantum mechanical analogue of the eigenvalue density ρ(λ) introduced in (3.15).

As a final remark, before embarking on the calculation of various scattering processes we would like to relate the Green's function introduced in (7.9) to the Feynman propagator (6.32) discussed in the previous section. These obey the same equation. Comparing the first line in (6.38) with the lower line of the Green's function in (7.12), we note that the two expressions take the same form. The difference lies in the Heaviside function appearing in ${\tilde {S}}_{\mathrm{F}}\left(x,{x}^{\prime } ;\omega \right)$. This is a simple manifestation of the fact that G(x, x'; z) stems from a first quantised setup, whereas ${\tilde {S}}_{\mathrm{F}}\left(x,{x}^{\prime } ;\omega \right)$ stems from a 'second quantised' picture allowing for particle creation and annihilation.

To ensure a large range of admissible scattering frequencies we take ${\nu }_{\mathrm{F}}^{2}\gg 1$. This condition will also ensure that the transmitted wave is exponentially suppressed, such that we can aptly refer to scattering processes leaking across the barrier as non-perturbative phenomena. The physical problem we are after consists of sending in an incoming density wave ρ_ω(x) at some incoming frequency ω near the spatial boundary at large and negative values of x, and measuring the scattered wave at some outgoing frequency ω', again at large negative values of x. Since we are scattering a composite object, the elements of our perturbative S-matrix will be given by certain convolutions of the reflection coefficients R_ν. In what follows, we will often express results in terms of the spatial coordinate

$$\begin{equation}u=-\mathrm{log}\left(-x\right),\quad x{< }0.\end{equation} \tag{ 7.13 }$$

Notice that for large negative values of x, u is also large and negative and u increases monotonically with x.

To compute scattering processes for multiple incoming and outgoing waves⁸ it is convenient to resort to the multi-particle picture described in section 6.5. We begin by considering the n-point function

$$\begin{equation}\mathcal{Q}\left(\left\{{t}_{k},{x}_{k}\right\}\right)=\langle \text{Fermi}\vert {T}_{\mathrm{B}}\prod _{k=1}^{n}\hat{n}\left({t}_{k},{x}_{k}\right)\vert \text{Fermi}\rangle ,\end{equation} \tag{ 7.14 }$$

where now T_B indicates bosonic time ordering and we defined $\hat{n}\left({t}_{k},{x}_{k}\right)$ in (6.28). Given that our theory is non-interacting, we can calculate $\mathcal{Q}$ via multiple Wick contractions using (6.32). We note that the above correlation function is invariant under (t_i, x_i) ↔ (t_j, x_j).

General approach. Before delving into concrete examples, it is worth explaining the general strategy. Starting with the quantum correlation function (7.14) we perform all possible Wick contractions. Keeping the fully connected part we obtain

$$\begin{equation}{\mathcal{Q}}_{\mathrm{c}}\left(\left\{{t}_{k},{x}_{k}\right\}\right)=-\frac{{\mathrm{i}}^{n}}{n}\sum _{\sigma \in \mathrm{P}\mathrm{e}\mathrm{r}\mathrm{m}\left(n\right)}\enspace \prod _{k=1}^{n}{S}_{\mathrm{F}}\left({t}_{{\sigma }_{k}},{x}_{{\sigma }_{k}};{t}_{{\sigma }_{k+1}},{x}_{{\sigma }_{k+1}}\right),\end{equation} \tag{ 7.15 }$$

where Perm(n) denotes the set of order n! permuting the elements in {1, 2, ..., n}. We also compute indices modulo n (e.g. σ_n+1 ≡ σ₁ and so on). The subscript c is a reminder that we are considering the fully connected part. Notice that (7.15) is invariant under (t_i, x_i) ↔ (t_j, x_j). Going to Fourier space the above leads to an expression manifestly invariant under (ω_i, x_i) ↔ (ω_j, x_j). From this expression we extract the fully connected part ${\tilde {\mathcal{Q}}}_{\mathrm{c}}\left(\left\{{\omega }_{i},{x}_{i}\right\}\right)$ and obtain the one-to-n scattering amplitude S(ω₁, ..., ω_n−1|ω_n) defined as the part of the correlation function taking the form

$$\begin{equation}\delta \left({\omega }_{1}+\cdots +{\omega }_{n}\right)\enspace \mathrm{exp}\left(\mathrm{i}{\omega }_{n}{u}_{n}-\mathrm{i}\sum _{k=1}^{n-1}{\omega }_{k}{u}_{k}\right)\enspace S\left({\omega }_{1},\dots ,{\omega }_{n-1}\vert {\omega }_{n}\right).\end{equation} \tag{ 7.16 }$$

Outgoing frequencies ω_k are positive for k = 1, ..., n − 1 and as a result, the frequency of the incoming wave ω_n must be negative. Since the original correlation function (7.14) is invariant under (t_i, x_i) ↔ (t_j, x_j) we impose a specific ordering of the spatial coordinates, namely x₁ > x₂ >⋯> x_n, to avoid overcounting. For m > 1 we can similarly define the m-to-n amplitude. We now proceed to study some concrete examples.

One-to-one scattering. As a warm up exercise, we compute the one-to-one scattering problem in the second quantised picture. We wish to calculate the Fourier transform of the connected part of

$$\begin{equation}\mathcal{Q}\left(\left\{{t}_{i},{x}_{i}\right\}\right)=\langle \text{Fermi}\vert \enspace {T}_{\mathrm{B}}\enspace \hat{n}\left({t}_{1},{x}_{1}\right)\enspace \hat{n}\left({t}_{2},{x}_{2}\right)\vert \text{Fermi}\rangle ,\end{equation} \tag{ 7.17 }$$

where now T_B corresponds to bosonic time ordering. Keeping track of the various Wick contractions and using the Fermion propagator (6.32), we find

$$\begin{equation}{\mathcal{Q}}_{\mathrm{c}}\left({t}_{1},{x}_{1} ;{t}_{2},{x}_{2}\right)=\frac{1}{2}\sum _{\sigma \in \mathrm{P}\mathrm{e}\mathrm{r}\mathrm{m}\left(2\right)}\enspace \prod _{k=1}^{2}{S}_{\mathrm{F}}\left({t}_{{\sigma }_{k}},{x}_{{\sigma }_{k}};{t}_{{\sigma }_{k+1}},{x}_{{\sigma }_{k+1}}\right).\end{equation} \tag{ 7.18 }$$

In appendix D we provide some steps leading to (7.18). Going to Fourier space we obtain

$$\begin{equation}{\tilde {\mathcal{Q}}}_{\mathrm{c}}\left({\omega }_{1},{x}_{1} ;{\omega }_{2},{x}_{2}\right)={\int }_{{\mathbb{R}}^{2}}\mathrm{d}{t}_{1}\mathrm{d}{t}_{2}\enspace {\mathrm{e}}^{+\mathrm{i}{\omega }_{1}{t}_{1}+\mathrm{i}{\omega }_{2}{t}_{2}}{\mathcal{Q}}_{\mathrm{c}}\left({t}_{1},{x}_{1};{t}_{2},{x}_{2}\right).\end{equation} \tag{ 7.19 }$$

Performing the integrals over t₁ and t₂, we are led to the following expression

$$\begin{align}\hfill {\tilde {\mathcal{Q}}}_{\mathrm{c}}\left({\omega }_{1},{x}_{1} ;{\omega }_{2},{x}_{2}\right)& =\frac{2\pi }{2}\delta \left({\omega }_{1}+{\omega }_{2}\right)\hfill \\ \hfill & \quad {\times}{\int }_{\mathbb{R}}\frac{\mathrm{d}\upsilon }{2\pi }\left({\tilde {S}}_{\mathrm{F}}\left({x}_{1},{x}_{2} ;\upsilon +{\omega }_{1}\right){\tilde {S}}_{\mathrm{F}}\left({x}_{2},{x}_{1} ;\upsilon \right)\right.\hfill \\ \hfill & \left.+{\tilde {S}}_{\mathrm{F}}\left({x}_{2},{x}_{1} ;\upsilon +{\omega }_{2}\right){\tilde {S}}_{\mathrm{F}}\left({x}_{1},{x}_{2} ;\upsilon \right)\right).\hfill \end{align} \tag{ 7.20 }$$

We can expand the above expression at large and negative values of x₁ and x₂. There are various pieces in this regime. Of these, we would like to keep those terms corresponding to particles that have been scattered by the barrier, and moreover do not exhibit large oscillations of the type $\sim {\mathrm{e}}^{{\pm}\mathrm{i}{x}^{2}/4}$. Expressing the result in terms of the spatial coordinate u = −log(−x) we find

$$\begin{equation}{\tilde {\mathcal{Q}}}_{\mathrm{c}}\left({\omega }_{1},{u}_{1} ;{\omega }_{2},{u}_{2}\right)=2\pi \delta \left({\omega }_{1}+{\omega }_{2}\right)\left({\omega }_{2}\enspace \frac{{\mathrm{e}}^{\mathrm{i}{\omega }_{1}{u}_{1}-\mathrm{i}{\omega }_{2}{u}_{2}}}{2\pi }+\frac{{\mathrm{e}}^{-\mathrm{i}{\omega }_{1}{u}_{1}+\mathrm{i}{\omega }_{2}{u}_{2}}}{2\pi }\enspace S\left({\omega }_{1}\vert {\omega }_{2}\right)\right).\end{equation} \tag{ 7.21 }$$

and S(ω₁|ω₂) is defined as

$$\begin{equation}S\left({\omega }_{1}\vert {\omega }_{2}\right)\equiv {\mathrm{e}}^{-\mathrm{i}{\omega }_{2}\mathrm{log}\left(-{\nu }_{\mathrm{F}}\right)}\enspace {\int }_{0}^{{\omega }_{1}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{2}}.\end{equation} \tag{ 7.22 }$$

The phase ${\mathrm{e}}^{-\mathrm{i}{\omega }_{2}\mathrm{log}\left(-{\nu }_{\mathrm{F}}\right)}$ subtracts the phase appearing in the perturbative expansion of the reflection coefficient (E.1). (One can also remove this phase by a redefinition of the u-coordinate.) Notice that the frequency of R_ν is always above the Fermi level since it corresponds to particle scattering, whereas the frequency of ${R}_{\nu }^{{\ast}}$ is always below the Fermi level since it corresponds to the scattering of a hole. Performing a large ν_F expansion, we obtain

$$\begin{equation}S\left({\omega }_{1}\vert {\omega }_{2}\right)={\omega }_{1}+\frac{1}{24{\nu }_{\mathrm{F}}^{2}}\left(\mathrm{i}{\omega }_{1}^{2}-{\omega }_{1}^{4}\left({\omega }_{1}-2\mathrm{i}\right)\right)+\mathcal{O}\left({\nu }_{\mathrm{F}}^{-4}\right).\end{equation} \tag{ 7.23 }$$

It is worth noting that the leading piece of S(ω₁|ω₂), including the energy conserving delta-function, is the one-to-one scattering amplitude we would have obtained for a theory of massless fields in two-spacetime dimensions in a theory subject to Poincaré invariance. The terms subleading in 1/ν_F are corrections which indicate any such Poincaré invariance is ultimately broken. We will comment on this from the perspective of a continuum description in the next section. Finally dimensional analysis may raise concerns about the above expression. We must remember, however, that in going to the Hamiltonian (7.1) we have set several dimensionful parameters to unity. Upon restoring them, we can also restore our faith in dimensional analysis.

One-to-two scattering. We now move on to the case of one-to-two scattering. This scattering amplitude encodes the type of interactions present in the theory. The correlation function of interest is

$$\begin{equation}\mathcal{Q}\left(\left\{{t}_{i},{x}_{i}\right\}\right)=\langle \text{Fermi}\vert \enspace {T}_{\mathrm{B}}\enspace \hat{n}\left({t}_{1},{x}_{1}\right)\hat{n}\left({t}_{2},{x}_{2}\right)\hat{n}\left({t}_{3},{x}_{3}\right)\vert \text{Fermi}\rangle ,\end{equation} \tag{ 7.24 }$$

The resulting expression for the connected part reads

$$\begin{equation}{\mathcal{Q}}_{\mathrm{c}}\left(\left\{{t}_{i},{x}_{i}\right\}\right)=\frac{\mathrm{i}}{3}\sum _{\sigma \in \mathrm{P}\mathrm{e}\mathrm{r}\mathrm{m}\left(3\right)}\enspace \prod _{k=1}^{3}{S}_{\mathrm{F}}\left({t}_{{\sigma }_{k}},{x}_{{\sigma }_{k}};{t}_{{\sigma }_{k+1}},{x}_{{\sigma }_{k+1}}\right).\end{equation} \tag{ 7.25 }$$

Fourier transforming with respect to the time coordinates, we obtain

$$\begin{equation}{\tilde {\mathcal{Q}}}_{\mathrm{c}}\left({\omega }_{1},{x}_{1} ;{\omega }_{2},{x}_{2} ;{\omega }_{3},{x}_{3}\right)={\int }_{{\mathbb{R}}^{3}}\mathrm{d}{t}_{1}\mathrm{d}{t}_{2}\mathrm{d}{t}_{3}\enspace {\mathrm{e}}^{\mathrm{i}{\omega }_{1}{t}_{1}+\mathrm{i}{\omega }_{2}{t}_{2}+\mathrm{i}{\omega }_{3}{t}_{3}}{\mathcal{Q}}_{\mathrm{c}}\left({t}_{1},{x}_{1} ;{t}_{2},{x}_{2} ;{t}_{3},{x}_{3}\right).\end{equation} \tag{ 7.26 }$$

Due to time-translation invariance, we obtain an overall delta function imposing the conservation of energy. Combining (7.16) with (7.26) we obtain

$$\begin{equation}S\left({\omega }_{1},{\omega }_{2}\vert \enspace {\omega }_{3}\right)\equiv {\mathrm{e}}^{-\mathrm{i}{\omega }_{3}\mathrm{log}\left(-{\nu }_{\mathrm{F}}\right)}\enspace \left[{\int }_{0}^{{\omega }_{1}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{3}}-{\int }_{{\omega }_{2}}^{{\omega }_{1}+{\omega }_{2}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{3}}\enspace \right].\end{equation} \tag{ 7.27 }$$

We can express the above expression in a large ν_F expansion using (E.1)

$$\begin{align}\hfill S\left({\omega }_{1},{\omega }_{2}\vert \enspace {\omega }_{3}\right)& =\frac{\mathrm{i}}{{\nu }_{\mathrm{F}}}\enspace {\omega }_{1}{\omega }_{2}{\omega }_{3}-\frac{\mathrm{i}}{24{\nu }_{\mathrm{F}}^{3}}{\omega }_{1}{\omega }_{2}{\omega }_{3}\left({\omega }_{3}+i\right)\left({\omega }_{3}+2i\right)\hfill \\ \hfill & {\times}\left({\omega }_{1}\left({\omega }_{1}-i\right)+{\omega }_{2}\left({\omega }_{2}-i\right)+1\right)+\mathcal{O}\left({\nu }_{\mathrm{F}}^{-5}\right).\hfill \end{align} \tag{ 7.28 }$$

Note that the alternating signs arise from the alternating signs in (6.38) and we take the frequencies of the scattered waves ω₁, ω₂ > 0.

One-to-three scattering. In direct analogy to the one-to-one and one-to-two case, we can work out the one-to-three scattering amplitude. After obtaining the possible Wick contractions of

$$\begin{equation}\mathcal{Q}\left(\left\{{t}_{k},{x}_{k}\right\}\right)=\langle \text{Fermi}\vert \enspace {T}_{\mathrm{B}}\enspace \hat{n}\left({t}_{1},{x}_{1}\right)\hat{n}\left({t}_{2},{x}_{2}\right)\hat{n}\left({t}_{3},{x}_{3}\right)\hat{n}\left({t}_{4},{x}_{4}\right)\vert \text{Fermi}\rangle ,\end{equation} \tag{ 7.29 }$$

we can express the connected piece as

$$\begin{equation}{\mathcal{Q}}_{\mathrm{c}}\left(\left\{{t}_{k},{x}_{k}\right\}\right)=-\frac{1}{4}\sum _{\sigma \in \mathrm{P}\mathrm{e}\mathrm{r}\mathrm{m}\left(4\right)}\enspace \prod _{k=1}^{4}{S}_{\mathrm{F}}\left({t}_{{\sigma }_{k}},{x}_{{\sigma }_{k}};{t}_{{\sigma }_{k+1}},{x}_{{\sigma }_{k+1}}\right).\end{equation} \tag{ 7.30 }$$

Fourier transforming, extracting the term without the ${\mathrm{e}}^{{\pm}\mathrm{i}{x}^{2}/4}$ oscillations, and using (7.16) we end up with

$$\begin{align}\hfill S\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\enspace \vert \enspace {\omega }_{4}\right)\enspace \equiv {\mathrm{e}}^{-\mathrm{i}{\omega }_{4}\mathrm{log}\left(-{\nu }_{\mathrm{F}}\right)}\enspace & \left[{\int }_{{\omega }_{2}+{\omega }_{3}}^{{\omega }_{1}+{\omega }_{2}+{\omega }_{3}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}}-{\int }_{{\omega }_{2}}^{{\omega }_{1}+{\omega }_{2}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}}\right.\hfill \\ \hfill & \left.-{\int }_{{\omega }_{3}}^{{\omega }_{1}+{\omega }_{3}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}}+{\int }_{0}^{{\omega }_{1}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}}\right].\hfill \end{align} \tag{ 7.31 }$$

Note that there is an overall minus sign arising from our original definition of the scattering amplitude (7.16) and (7.30). Expanding the above expression at large ν_F by making use of (E.1) we find

$$\begin{align}\hfill S\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\enspace \vert \enspace {\omega }_{4}\right)=& -\frac{1}{{\nu }_{\mathrm{F}}^{2}}\enspace {\omega }_{1}{\omega }_{2}{\omega }_{3}{\omega }_{4}\left({\omega }_{4}+i\right)+\frac{1}{24{\nu }_{\mathrm{F}}^{4}}\enspace {\omega }_{1}{\omega }_{2}{\omega }_{3}{\omega }_{4}\left({\omega }_{4}+i\right)\left({\omega }_{4}+2i\right)\left({\omega }_{4}+3i\right)\hfill \\ \hfill & {\times}\left({\omega }_{1}\left({\omega }_{1}-i\right)+{\omega }_{2}\left({\omega }_{2}-i\right)+{\omega }_{3}\left({\omega }_{3}-i\right)+1\right)+\mathcal{O}\left({\nu }_{\mathrm{F}}^{-6}\right).\hfill \end{align} \tag{ 7.32 }$$

One-to-(n − 1) scattering. From the above scattering amplitudes we observe

$$\begin{equation}S\left({\omega }_{1},\dots ,{\omega }_{n-1}\vert \enspace {\omega }_{n}\right)=\sum _{L=0}^{\infty }{\nu }_{\mathrm{F}}^{n-1+2L}{S}^{\left(L\right)}\left({\omega }_{1},\dots ,{\omega }_{n-1}\vert \enspace {\omega }_{n}\right),\end{equation} \tag{ 7.33 }$$

where for the case of one-to-(n − 1) scattering the left hand side is given by [106]

$$\begin{equation}S\left({\omega }_{1},\dots ,{\omega }_{n-1}\vert \enspace {\omega }_{n}\right)\equiv {\mathrm{e}}^{-\mathrm{i}{\omega }_{n}\mathrm{log}\left(-{\nu }_{\mathrm{F}}\right)}\enspace \sum _{{\Omega}\enspace \subseteq \left\{{\omega }_{1},{\ldotp\ldotp},\enspace {\omega }_{n-1}\right\}}\enspace {\left(-1\right)}^{\vert {\Omega}\vert +1}\enspace {\int }_{0}^{\omega \left({\Omega}\right)}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{n}}.\end{equation} \tag{ 7.34 }$$

where ω(Ω) is the sum of all elements in Ω. Furthermore, for the leading term in a large ν_F expansion of the one-to-(n − 1) scattering amplitude, we observe

$$\begin{equation}{S}^{\left(0\right)}\left({\omega }_{1},\dots ,{\omega }_{n-1}\vert \enspace {\omega }_{n}\right)=-{\mathrm{i}}^{n}\enspace {\omega }_{1}\cdots {\omega }_{n-1}\prod _{k=0}^{n-3}\enspace \left({\omega }_{n}+\mathrm{i}k\right),\quad n{\geqslant}3.\end{equation} \tag{ 7.35 }$$

For the the first sub-leading term and n ⩾ 3 we find

$$\begin{equation}{S}^{\left(1\right)}\left({\omega }_{1},\dots ,{\omega }_{n-1}\vert \enspace {\omega }_{n}\right)={\left(\mathrm{i}\right)}^{n}\enspace \frac{1}{24}\enspace {\omega }_{1}\cdots {\omega }_{n-1}\prod _{k=0}^{n-1}\enspace \left({\omega }_{n}+ki\right)\left(\sum _{\ell =1}^{n-1}{\omega }_{\ell }\left({\omega }_{\ell }-i\right)+1\right),\end{equation} \tag{ 7.36 }$$

where ${\omega }_{n}=-{\sum }_{k=1}^{n-1}{\omega }_{k}$. For both expressions we made use of the large ν_F expansion (E.1). In appendix F we obtain some higher-point scattering amplitudes that confirm the above expressions, at least for low enough particle number.

Two-to-two scattering. To obtain the two-to-two amplitude, instead of looking for terms of the form (7.16) we must extract the term proportional to

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}{\omega }_{1}{u}_{1}-\mathrm{i}{\omega }_{2}{u}_{2}+\mathrm{i}{\omega }_{3}{u}_{3}+\mathrm{i}{\omega }_{4}{u}_{4}}\enspace \delta \left({\omega }_{1}+{\omega }_{2}+{\omega }_{3}+{\omega }_{4}\right)\enspace S\left({\omega }_{1},{\omega }_{2}\enspace \vert \enspace {\omega }_{3},{\omega }_{4}\right)\enspace \end{equation} \tag{ 7.37 }$$

from the connected part of the four point correlation function (7.30). This defines for us the two-to-two amplitude S(ω₁, ω₂|ω₃, ω₄). In the above expression we assume that ω₁ and ω₂ are the positive frequencies of the outgoing waves, while ω₃ and ω₄ are the negative frequencies of the two incoming waves. Additionally we have to order the frequencies. We choose

$$\begin{equation}{\omega }_{2}{ >}{\omega }_{1},\qquad {\omega }_{4}{ >}{\omega }_{3},\qquad {\omega }_{2}{ >}-{\omega }_{3},\qquad {\omega }_{2}{ >}-{\omega }_{4},\end{equation} \tag{ 7.38 }$$

such that ω₂ = max{ω₁, ω₂, |ω₃|, |ω₄|}. Combining with (7.37) and (7.30) we obtain [106–108]

$$\begin{align}\hfill S\left({\omega }_{1},{\omega }_{2}\vert {\omega }_{3},{\omega }_{4}\right)\equiv & {\mathrm{e}}^{-\mathrm{i}\left({\omega }_{3}+{\omega }_{4}\right)\mathrm{log}\left(-{\nu }_{\mathrm{F}}\right)}\left[{\int }_{{\omega }_{3}}^{{\omega }_{1}+{\omega }_{3}}\mathrm{d}\upsilon {R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{3}}^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}}\right.\hfill \\ \hfill & +{\int }_{{\omega }_{2}}^{{\omega }_{1}+{\omega }_{2}}\mathrm{d}\upsilon {R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{1}+{\omega }_{2}}\hfill \\ \hfill & -{\int }_{{\omega }_{3}}^{{\omega }_{1}+{\omega }_{3}}\mathrm{d}\upsilon {R}_{{\nu }_{\mathrm{F}}-\upsilon }{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{3}}^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{1}+{\omega }_{3}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{2}}^{{\ast}}\hfill \\ \hfill & \left.-{\int }_{{\omega }_{4}}^{{\omega }_{1}+{\omega }_{4}}\mathrm{d}\upsilon {R}_{{\nu }_{\mathrm{F}}-\upsilon }{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{4}}^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{1}+{\omega }_{4}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{2}}^{{\ast}}\right]\hfill \\ \hfill =& -\frac{1}{{\nu }_{\mathrm{F}}^{2}}{\omega }_{1}{\omega }_{2}{\omega }_{3}{\omega }_{4}\left({\omega }_{2}-i\right)+\mathcal{O}\left({\nu }_{\mathrm{F}}^{-4}\right),\hfill \end{align} \tag{ 7.39 }$$

as the two-to-two scattering amplitude.

Two-to-three scattering. To obtain the two-to-three amplitude we consider the connected correlation function (F.1). We are then looking for terms in proportional to

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}{\omega }_{1}{u}_{1}-\mathrm{i}{\omega }_{2}{u}_{2}-\mathrm{i}{\omega }_{3}{u}_{3}+\mathrm{i}{\omega }_{4}{u}_{4}+\mathrm{i}{\omega }_{5}{u}_{5}}\delta \left({\omega }_{1}+{\omega }_{2}+{\omega }_{3}+{\omega }_{4}+{\omega }_{5}\right)\enspace S\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\enspace \vert \enspace {\omega }_{4},{\omega }_{5}\right)\end{equation} \tag{ 7.40 }$$

with ω₁, ω₂, ω₃ outgoing and positive and ω₄, ω₅ incoming negative. The above expression defines for us the scattering amplitude for two-to-three scattering. Additionally we assume

$$\begin{equation}{\omega }_{3}{ >}{\omega }_{2}{ >}{\omega }_{1},\qquad {\omega }_{5}{ >}{\omega }_{4},\qquad {\omega }_{3}{ >}-{\omega }_{4},\qquad {\omega }_{2}{< }-{\omega }_{5},\qquad {\omega }_{1}{< }-{\omega }_{5},\end{equation} \tag{ 7.41 }$$

and x₁ > x₂ > x₃ > x₄ > x₅. In particular (7.41) implies that ω₃ = max{ω₁, ω₂, ω₃, |ω₄|, |ω₅|}. We then find:

$$\begin{align}\hfill S\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\enspace \vert \enspace {\omega }_{4},{\omega }_{5}\right)\equiv & \enspace {\mathrm{e}}^{-\mathrm{i}\left({\omega }_{4}+{\omega }_{5}\right)\mathrm{log}\left(-{\nu }_{\mathrm{F}}\right)} \left[{\int }_{0}^{{\omega }_{1}} \enspace \enspace \enspace \mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{2}-{\omega }_{4}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{1}+{\omega }_{5}}^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{1}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{2}}^{{\ast}}\right.\hfill \\ \hfill & -{\int }_{{\omega }_{4}}^{{\omega }_{1}+{\omega }_{4}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{3}-{\omega }_{5}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{4}}^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon }{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{3}}^{{\ast}}\hfill \\ \hfill & -{\int }_{{\omega }_{2}+{\omega }_{3}}^{{\omega }_{1}+{\omega }_{2}+{\omega }_{3}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}-{\omega }_{5}}{R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}\hfill \\ \hfill & -{\int }_{{\omega }_{5}}^{{\omega }_{1}+{\omega }_{5}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{3}-{\omega }_{4}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{5}}^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon }{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{3}}^{{\ast}}\hfill \\ \hfill & +{\int }_{{\omega }_{3}}^{{\omega }_{1}+{\omega }_{3}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}-{\omega }_{5}}{R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}\hfill \\ \hfill & +{\int }_{0}^{{\omega }_{1}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{2}-{\omega }_{5}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{1}+{\omega }_{4}}^{{\ast}}{R}_{{\nu }_{\mathrm{F}}-\upsilon +{\omega }_{1}}{R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{2}}^{{\ast}}\hfill \\ \hfill & +{\int }_{0}^{{\omega }_{1}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}-{\omega }_{5}}{R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}-\left.{\int }_{{\omega }_{2}}^{{\omega }_{1}+{\omega }_{2}}\mathrm{d}\upsilon \enspace {R}_{{\nu }_{\mathrm{F}}-\upsilon -{\omega }_{4}-{\omega }_{5}}{R}_{{\nu }_{\mathrm{F}}-\upsilon }^{{\ast}}\right]\hfill \\ \hfill =& \frac{\mathrm{i}}{{\nu }_{\mathrm{F}}^{3}}{\omega }_{1}{\omega }_{2}{\omega }_{3}{\omega }_{4}{\omega }_{5}\left({\omega }_{3}-i\right)\left({\omega }_{1}+{\omega }_{2}+{\omega }_{3}-2i\right)+\mathcal{O}\left({\nu }_{\mathrm{F}}^{-5}\right).\hfill \end{align} \tag{ 7.42 }$$

We note that any two-to-n amplitude contains factors of either two or four reflection coefficients.

In such a way we can continue to compute higher amplitudes from our basic building blocks. We thus observe the emergence of a perturbative S-matrix admitting a small 1/ν_F expansion. So far, we obtained the S-matrix directly from the large N quantum mechanical matrix model. It is natural to ask whether the same S-matrix can be obtained directly from some weakly coupled system. Before turning to this we summarise the results of the (n − m)-to-m scattering amplitudes S(ω₁, ..., ω_m|ω_m+1, ..., ω_n) calculated in this section and appendix F (tables 1–3).

Table 1. Summary of the scattering amplitudes calculated in this section and appendix F. Each box contains the contribution to the scattering amplitude at the inverse power of ν_F in the leftmost box. An x indicates a nonzero result easily obtained by expanding the reflection coefficients to higher order in 1/ν_F. Within S(ω₁, ..., ω_n−1|ω_n), ω_n denotes the incoming (negative) frequency, whereas {ω₁, ..., ω_n−1} denote the n − 1 outgoing (positive) frequencies.

*	S(ω1|ω2)	S(ω1, ω2|ω3)	S(ω1, ω2, ω3|ω4)
${\nu }_{\mathrm{F}}^{0}$	ω1	0	0
${\nu }_{\mathrm{F}}^{-1}$	0	iω1ω2ω3	0
${\nu }_{\mathrm{F}}^{-2}$	$\frac{1}{24}\left(\mathrm{i}{\omega }_{1}^{2}-{\omega }_{1}^{4}\left({\omega }_{1}-2i\right)\right)$	0	−ω1ω2ω3ω4(ω4 + i)
${\nu }_{\mathrm{F}}^{-3}$	0	$-\frac{\mathrm{i}}{24}\enspace {\omega }_{1}{\omega }_{2}{\prod }_{k=0}^{2}\enspace \left({\omega }_{3}+ki\right)$	0
		${\times}\left({\sum }_{\ell =1}^{2}{\omega }_{\ell }\left({\omega }_{\ell }-i\right)+1\right)$
${\nu }_{\mathrm{F}}^{-4}$	x	0	$\frac{1}{24}\enspace {\omega }_{1}{\omega }_{2}{\omega }_{3}{\prod }_{k=0}^{3}\enspace \left({\omega }_{4}+ki\right)$
			${\times}\left({\sum }_{\ell =1}^{3}{\omega }_{\ell }\left({\omega }_{\ell }-i\right)+1\right)$
	0	x	0
	x	0	x

Table 2. Summary of the scattering amplitudes calculated in this section and appendix F. Each box contains the contribution to the scattering amplitude at the inverse power of ν_F in the leftmost box. An x indicates a nonzero result easily obtained by expanding the reflection coefficients to higher order in 1/ν_F. Within S(ω₁, ..., ω_n−1|ω_n), ω_n denotes the incoming (negative) frequency, whereas {ω₁, ..., ω_n−1} denote the n − 1 outgoing (positive) frequencies.

∗	S(ω1, ..., ω4|ω5)	S(ω1, ..., ω5|ω6)	S(ω1, ..., ω6|ω7)
	0	0	0
${\nu }_{\mathrm{F}}^{-3}$	−iω1ω2ω3ω4ω5	0	0
	× (ω5 + i)(ω5 + 2i)
${\nu }_{\mathrm{F}}^{-4}$	0	ω1⋯ω6	0
		× (ω6 + i)(ω6 + 2i)(ω6 + 3i)
${\nu }_{\mathrm{F}}^{-5}$	$\frac{\mathrm{i}}{24}\enspace {\omega }_{1}\dots {\omega }_{4}{\prod }_{k=0}^{4}\enspace \left({\omega }_{5}+ki\right)$	0	$\mathrm{i}{\omega }_{1}\cdots {\omega }_{7}$
	${\times}\left({\sum }_{\ell =1}^{4}{\omega }_{\ell }\left({\omega }_{\ell }-i\right)+1\right)$		${\times}{\prod }_{k=0}^{4}\left({\omega }_{7}+ki\right)$
${\nu }_{\mathrm{F}}^{-6}$	0	$-\frac{1}{24}\enspace {\omega }_{1}\cdots {\omega }_{5}{\prod }_{k=0}^{5}\enspace \left({\omega }_{6}+ki\right)$	0
		${\times}\left({\sum }_{\ell =1}^{5}{\omega }_{\ell }\left({\omega }_{\ell }-i\right)+1\right)$
${\nu }_{\mathrm{F}}^{-7}$	x	0	$-\frac{\mathrm{i}}{24}\enspace {\omega }_{1}\cdots {\omega }_{6}{\prod }_{k=0}^{6}\enspace \left({\omega }_{7}+ki\right)$
			${\times}\left({\sum }_{\ell =1}^{6}{\omega }_{\ell }\left({\omega }_{\ell }-i\right)+1\right)$
	0	x	0
	x	0	x

Table 3. Summary of the scattering amplitudes calculated in this section and appendix F. Each box contains the contribution to the scattering amplitude at the inverse power of ν_F in the leftmost box. An x indicates a nonzero result easily obtained by expanding the reflection coefficients to higher order in 1/ν_F. Within S(ω₁, ..., ω_m|ω_m+1, ..., ω_n), {ω_m+1, ..., ω_n} denote the n − m incoming (negative) frequencies and {ω₁, ..., ω_m} the m outgoing (positive) frequencies. Note in particular that we choose an ordering of the frequencies which implies that ω₂ is the maximum for the two-to-two and three-to-two amplitude, whereas ω₃ is the maximum for the two-to-three amplitude.

∗	S(ω1, ω2|ω3, ω4)	S(ω1, ω2|ω3, ω4, ω5)	S(ω1, ω2, ω3|ω4, ω5)
	0	0	0
${\nu }_{\mathrm{F}}^{-2}$	−ω1ω2ω3ω4(ω2 − i)	0	0
${\nu }_{\mathrm{F}}^{-3}$	0	−iω1ω2ω3ω4ω5	iω1ω2ω3ω4ω5
		× (ω2 − i)(ω2 − 2i)	× (ω3 − i)(ω1 + ω2 + ω3 − 2i)
${\nu }_{\mathrm{F}}^{-4}$	x	0	0
${\nu }_{\mathrm{F}}^{-5}$	0	x	x
	x	0	0
	0	x	x

We begin our discussion by considering a theory of pure geometry in two-dimensions. The question of interest is as follows. Let us endow a compact genus h Riemann surface, Σ_h, with a metric g_ij. The most general local, diffeomorphism invariant, two-derivative action is given by

$$\begin{equation}{S}_{\text{grav}}\left[{g}_{ij};{{\Sigma}}_{h}\right]=-\frac{\vartheta }{4\pi }{\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{g}\enspace R+{\Lambda}{\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{g},\end{equation} \tag{ 8.1 }$$

where ϑ and Λ are real parameters and R is the Ricci scalar. By the Gauss–Bonnet theorem, the first term is a topological invariant proportional to the Euler characteristic χ_h. The second term, which is the cosmological constant term, computes the area of the surface. Motivated by the continuum limit of the discrete picture discussed in section 3 and the resulting expression (3.80), we are prompted to study the path integral

$$\begin{equation}{\mathcal{Z}}_{{{\Sigma}}_{h}}\left[\vartheta ;{\Lambda}\right]={\mathrm{e}}^{\vartheta \enspace {\chi }_{h}}\enspace \int \frac{\mathcal{D}{g}_{ij}}{\mathrm{v}\mathrm{o}\mathrm{l}\enspace \mathrm{d}\mathrm{i}\mathrm{ff}}\enspace \mathrm{exp}\left(-{\Lambda}{\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{g}\right).\end{equation} \tag{ 8.2 }$$

At this stage, it is useful to recall that a two-dimensional metric is diffeomorphic, at least locally, to the following

$$\begin{equation}{g}_{ij}\left(x\right)={\mathrm{e}}^{\varphi \left(x\right)}{\tilde {g}}_{ij}\left(x\right),\end{equation} \tag{ 8.3 }$$

where ${\tilde {g}}_{ij}$ is a fixed reference metric, often referred to as the fiducial metric, and φ(x) is a local Weyl factor that remains unfixed. The parameterization (8.3) leaves a particular subgroup of diffeomorphisms unfixed. Indeed, if we consider the set of coordinate transformations that map ${\tilde {g}}_{ij}$ to itself times a local Weyl factor e^σ(x), we can return to the original expression by shifting φ(x) → φ(x) − σ(x). That the group of residual diffeomorphisms in the gauge (8.3) is the two-dimensional conformal group will play an important role in what comes. In fact, the transformation

$$\begin{equation}{\tilde {g}}_{ij}\left(x\right)\to {\mathrm{e}}^{\tilde {\sigma }\left(x\right)}{\tilde {g}}_{ij}\left(x\right),\qquad \varphi \left(x\right)\to \varphi \left(x\right)-\tilde {\sigma }\left(x\right)\enspace \end{equation} \tag{ 8.4 }$$

for general $\tilde {\sigma }\left(x\right)$ is a redundancy of the parameterization (8.3). As such, the resulting theory governing φ(x) and the accompanying Fadeev–Popov ghost fields should be invariant under the more general Weyl transformations (8.4). In the gauge (8.3), the ghost theory is described by a conformal field theory of bc-ghosts in a fixed background ${\tilde {g}}_{ij}$ whose central charge is c_g = −26. The resulting theory consisting of the bc-ghosts and the φ-field should be invariant under the redundancy (8.4). Thus, the theory governing the Weyl mode φ(x) must be a two-dimensional conformal field theory with central charge c_L = 26 such that the net conformal anomaly vanishes. If we further assume the φ-sector admits a local description, one can postulate an action of the following type [111, 112]

$$\begin{equation}{S}_{L}\left[\varphi ,{\tilde {g}}_{ij}\right]=\frac{1}{4\pi }{\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\left({\tilde {g}}^{ij}{\partial }_{i}\varphi {\partial }_{j}\varphi +Q\tilde {R}\varphi +4\pi {\Lambda}\enspace {\mathrm{e}}^{2b\varphi }\right).\end{equation} \tag{ 8.5 }$$

The above action is known as the Liouville action and has been studied extensively [33, 114, 121]. In order for the kinetic term to have a standard normalisation, we have allowed for a constant rescaling of φ by 2b, such that the resulting theory should be invariant under the following redundancy

$$\begin{equation}{\tilde {g}}_{ij}\left(x\right)\to {\mathrm{e}}^{\tilde {\sigma }\left(x\right)}{\tilde {g}}_{ij}\left(x\right),\qquad \varphi \left(x\right)\to \varphi \left(x\right)-\tilde {\sigma }\left(x\right)/2b.\end{equation} \tag{ 8.6 }$$

When Q = b + 1/b, the Liouville action describes a two-dimensional conformal field theory. It has some unusual features such as a continuous spectrum and the absence of a normalisable vacuum state. Nevertheless, it admits a consistent quantisation and many of its properties are known explicitly. For instance, its central charge is given by c_L = 1 + 6Q², and the theory contains spinless primary operators ${\mathcal{O}}^{\left(\alpha \right)}={\mathrm{e}}^{2\alpha \varphi }$ whose conformal dimensions are

$$\begin{equation}{{\Delta}}_{\alpha }=\alpha \left(Q-\alpha \right),\qquad {\overline{{\Delta}}}_{\alpha }=\alpha \left(Q-\alpha \right).\end{equation} \tag{ 8.7 }$$

These conformal dimensions can be used to fix b in (8.5), since conformal invariance requires ${\mathcal{O}}^{\left(b\right)}={\mathrm{e}}^{2b\varphi }$ to have conformal dimensions $\left({{\Delta}}_{b},{\overline{{\Delta}}}_{b}\right)=\left(1,1\right)$. One thus finds

$$\begin{equation}b=\frac{Q}{2}{\pm}\sqrt{\frac{{Q}^{2}}{4}-1}.\end{equation} \tag{ 8.8 }$$

Requiring c_L + c_g = 0 fixes $Q={\pm}5/\sqrt{6}$. Given the discrete symmetry Q → −Q, b → −b, and φ → −φ, we can pick Q > 0 without loss of generality. One then finds the two solutions b = b_± with ${b}_{+}=\sqrt{3/2}$ and ${b}_{-}=\sqrt{2/3}$. If one further requires the existence of a semi-classical limit, i.e. a limit where Q → ∞ gives rise to a sensible saddle-point approximation, one is forced to pick the negative root in (8.8). We will consider the negative root in what follows even in the absence of a semiclassical limit.

In this way we obtain a description of the path integral for a theory of pure two-dimensional geometry in terms of Liouville theory. A natural collection of observables is given by the set of diffeomorphism invariant functionals of g_ij. In the Weyl gauge, this is given by the set of conformal primaries built from φ with conformal dimensions $\left({\Delta},\overline{{\Delta}}\right)=\left(1,1\right)$. A particularly important one is the area operator

$$\begin{equation}{\mathcal{A}}_{h}={\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{g}={\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\enspace \enspace {\mathrm{e}}^{2b\varphi },\end{equation} \tag{ 8.9 }$$

where b is given by (8.8). Expectation values of ${\mathcal{A}}_{h}$ are calculated by taking derivatives of $-\mathrm{log}\enspace {\mathcal{Z}}_{{{\Sigma}}_{h}}\left[\vartheta ;{\Lambda}\right]$ with respect of Λ.

Generally speaking a local quantum field theory can be placed on an arbitrary curved background. If the theory has no diffeomorphism anomaly, we can further gauge the diffeomorphism group. Thus, it would seem we can couple arbitrary matter content to two-dimensional gravity. Following this reasoning, and going to the Weyl gauge (8.3), a consistent theory of matter coupled to two-dimensional gravity is described by the following action

$$\begin{equation}{S}_{\text{eff}}={S}_{L}\left[\varphi ,{\tilde {g}}_{ij}\right]+{S}_{\text{ghost}}\left[b,c,{\tilde {g}}_{ij}\right]+{S}_{\text{matter}}\left[X,\varphi ,{\tilde {g}}_{ij}\right].\end{equation} \tag{ 8.10 }$$

Invariance under the redundancy (8.4) and the residual diffeomorphisms implies that the above theory must be a two-dimensional conformal field theory with vanishing central charge. This statement must be true, curious as it may sound, even for a massive matter theory. Somehow, the Liouville mode must 'dress' the matter theory in such a way as to make the resulting system conformal. Of course, given that the original matter theory is not conformal, the resulting conformal theory obtained upon coupling to the Liouville mode must be strongly interacting even if the original matter theory is not [117].

In what follows we will consider the simpler situation, namely, a matter theory which is itself a two-dimensional conformal field theory with central charge c_m. In this case, the matter action S_matter becomes solely a function of the matter fields and ${\tilde {g}}_{ij}$. However, the path integration measures for the matter and ghost fields produce the Liouville action in the Weyl gauge [9].⁹ Thus, one can postulate that (8.5) continues to govern the φ sector, except now we must fix the parameters Q and b such that c_L + c_m − 26 = 0. Proceeding along similar lines to the previous sub-section, we end up with the expressions

$$\begin{equation}Q=\sqrt{\frac{25-{c}_{m}}{6}},\qquad b=\frac{\sqrt{25-{c}_{m}}{\pm}\sqrt{1-{c}_{m}}}{2\sqrt{6}}.\end{equation} \tag{ 8.11 }$$

For c_m = 0, we recover (8.8).

The expressions (8.11) demonstrate an important point. The resulting theory is highly sensitive to c_m. For c_m ⩽ 1, we have that b is real and positive. For c_m ∈ (1, 25) the parameter b becomes complex with non-vanishing real and imaginary parts, while Q remains real. Finally, for c_m ⩾ 25, we have that b and Q become pure imaginary. If we take c_m = 25 and tune Λ = 0, we observe that (8.5) reduces to the action of a free scalar, such that our resulting theory is nothing more than the critical bosonic string. All other cases describe non-critical string theories.

In the presence of matter, we can extend the space of admissible observables. For instance, a spinless weight (Δ, Δ) conformal primary ${\mathcal{O}}_{{\Delta}}$ of the matter theory can be combined with the Liouville exponential operator to produce a 'dressed' operator ${\mathcal{O}}_{{\Delta}}={\mathrm{e}}^{2\sigma \varphi }{\mathcal{O}}_{{\Delta}}$ with

$$\begin{equation}\sigma =\frac{\sqrt{25-{c}_{m}}{\pm}\sqrt{24{\Delta}+1-{c}_{m}}}{2\sqrt{6}}.\end{equation} \tag{ 8.12 }$$

In this way, we can obtain a significant collection of observables. It is worth noticing that for certain values of c_m and Δ, σ may become complex. For such situations, ${\mathcal{O}}^{\left(\sigma \right)}={\mathrm{e}}^{2\sigma \varphi }$ becomes a complex (rather than Hermitian) operator. To render it Hermitian we can consider linear combinations of ${\mathcal{O}}^{\left(\sigma \right)}$ with its complex conjugate.

Small and/or negative central charge. The most understood situation occurs when c_m < 1. For example, we could take the matter conformal field theory to be one of the minimal models [11, 12] with central charge

$$\begin{equation}{c}_{m}^{\left(p,q\right)}=1-\frac{6{\left(p-q\right)}^{2}}{pq},\qquad \left(p,q\right)\enspace \text{coprime},\end{equation} \tag{ 8.13 }$$

and operators of conformal dimension

$$\begin{equation}{{\Delta}}_{r,s}=\frac{{\left(rq-sp\right)}^{2}-{\left(p-q\right)}^{2}}{4pq},\quad r=1,\dots ,p-1,s=1,\dots ,q-1.\end{equation} \tag{ 8.14 }$$

The limit c_m → −∞ leads to a parameterically large Q and, for the negative branch, a parametrically small b. In this limit, the system is driven to a semiclassical regime where saddle point techniques may be employed. Though there are several known minimal models with large and negative c_m, they are all non-unitary.¹⁰

Large central charge. For c_m > 25, we could imagine a contour rotation φ → iφ that yields the action (8.5) real. This comes at the price (or perhaps reward) of changing the sign of the kinetic term of φ. In critical string theory, we would interpret a wrong sign scalar as a time-like direction in the target space. Perhaps the same is true for c_m > 25, in which case non-critical string theory might provide an interesting window into time-dependent backgrounds in string theory. In such circumstances caution must be exercised, as many of the techniques we understand require rotating time-like target space coordinates back to their spacelike values, which for c_m > 25 would take us back to a complex action.

Intermediate central charge. For c_m ∈ (1, 25) the situation is far less understood. In this case, the action is intrinsically complex, and new methods are required to deal with it. Curiously, and perhaps interestingly, this case includes a matter conformal field theory with c_m = 3 and c_m = 4.

8.2.1. Matter matters

At this stage we have equipped ourselves with the necessary tools to compute the following quantity

$$\begin{equation}{\tilde {\mathcal{Z}}}_{{{\Sigma}}_{h}}\left[\upsilon ;{\Lambda}\right]=\int \frac{\mathcal{D}{g}_{ij}}{\text{vol}\enspace \text{diff}}\enspace {\mathcal{D}}_{{g}_{ij}}X\enspace {\mathrm{e}}^{-{S}_{\text{grav}}\left[{g}_{ij}\right]-{S}_{\text{matter}}\left[X,{g}_{ij}\right]}\enspace \delta \enspace \left({\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{g}-\upsilon \right).\end{equation} \tag{ 8.15 }$$

The δ-function in the above expression fixes the area to the particular value υ. By fixing the area operator (8.9) to take a fixed value, the fixed area partition function (8.15) tames any potential divergences in the original path integral which can subsequently be analysed in a clearer fashion.

Once again, resorting to the Weyl gauge (8.3) and keeping in mind the rescaling of φ, we can make progress in evaluating (8.15). We obtain the path integral expression

$$\begin{equation}{\tilde {\mathcal{Z}}}_{{{\Sigma}}_{h}}\left[\upsilon ;{\Lambda}\right]={\mathrm{e}}^{\vartheta \enspace {\chi }_{h}}\enspace \int \frac{\mathcal{D}\varphi }{\text{vol}\enspace \mathcal{G}}\enspace {\mathrm{e}}^{-{S}_{L}\left[\varphi ,{\tilde {g}}_{ij}\right]}\enspace \delta \enspace \left({\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\enspace {\mathrm{e}}^{2b\varphi }-\upsilon \right).\end{equation} \tag{ 8.16 }$$

Any local ultraviolet divergences arising upon integrating out the matter and ghost fields are absorbed in the bare couplings of the gravitational theory. The group $\mathcal{G}$ is the residual gauge group upon fixing the Weyl gauge. For instance, when h = 0 we could fix ${\tilde {g}}_{ij}$ to be the round metric on the two-sphere in which case $\mathcal{G}=PSL\left(2,\mathbb{C}\right)$. In contrast to the critical string, dividing by the potentially divergent volume of $\mathcal{G}$ does not mean that the partition function vanishes. This is due to the fact that φ transforms non-trivially under the residual gauge transformations.

Given (8.16), the approach is to consider a constant shift φ → φ + log υ^1/2b that allows us to remove υ from the δ-function while only affecting the remaining Liouville theory in a simple way [111, 112]. Using δ(ζx) = δ(x)/|ζ| and (8.5) it is readily found that

$$\begin{equation}{\tilde {\mathcal{Z}}}_{{{\Sigma}}_{h}}\left[\upsilon ;{\Lambda}\right]={\mathcal{N}}_{h}\enspace {\mathrm{e}}^{\vartheta {\chi }_{h}}\enspace {\mathrm{e}}^{-{\Lambda}\upsilon }\enspace {\upsilon }^{-Q{\chi }_{h}/2b-1},\end{equation} \tag{ 8.17 }$$

where ${\mathcal{N}}_{h}$ is a normalisation constant that may depend on the genus h but is independent of υ, ϑ, and Λ. We recall that b is given by (8.11). The continuum partition function can be recovered by integrating over υ as follows

$$\begin{equation}{\mathcal{Z}}_{{{\Sigma}}_{h}}\left[\vartheta ;{\Lambda}\right]={\mathrm{e}}^{\vartheta {\chi }_{h}}{\int }_{0}^{\infty }\mathrm{d}\upsilon \enspace {\tilde {\mathcal{Z}}}_{{{\Sigma}}_{h}}\left[\upsilon ;{\Lambda}\right]={\mathcal{N}}_{h}\enspace {\mathrm{e}}^{\vartheta {\chi }_{h}}\enspace {{\Lambda}}^{\frac{Q{\chi }_{h}}{2b}}\enspace {\Gamma}\left(-Q{\chi }_{h}/2b\right),\quad h\ne 1,\end{equation} \tag{ 8.18 }$$

where we have regularised the integral near υ = 0 by analytic continuation. When h = 1 the resulting expression is logarithmic in Λ, regardless of the value of Q/2b.

It is tempting to view (8.18) from the perspective of critical systems near a phase transition, with Λ being a tuneable parameter driving the system towards criticality. To this end, it is customary in the literature to define a critical exponent Γ_str ≡ 2 − Q/b, often referred to as the string susceptibility. Combining Q and b in (8.11) and choosing the negative branch for b, we obtain the KPZ relation [110]

$$\begin{equation}{{\Gamma}}_{\text{str}}=\frac{1}{12}\left({c}_{m}-1-\sqrt{\left({c}_{m}-1\right)\left({c}_{m}-25\right)}\right).\end{equation} \tag{ 8.19 }$$

Some examples are given by [14–16, 82, 83, 87]

$$\begin{equation}\left({{\Gamma}}_{\text{str}},Q/2b\right)=\begin{cases}\left(-1/2,5/4\right)\quad \hfill & \mathrm{f}\mathrm{o}\mathrm{r}\enspace {c}_{m}=0,\hfill \\ \left(-1/3,7/6\right)\quad \hfill & \mathrm{f}\mathrm{o}\mathrm{r}\enspace {c}_{m}=1/2,\hfill \\ \left(0,1\right)\quad \hfill & \mathrm{f}\mathrm{o}\mathrm{r}\enspace {c}_{m}=1.\hfill \end{cases}\end{equation} \tag{ 8.20 }$$

At large negative c_m we have

$$\begin{equation}{\mathcal{Z}}_{{{\Sigma}}_{h}}\left[\vartheta ;{\Lambda}\right]={\tilde {\mathcal{N}}}_{h}\enspace {\mathrm{e}}^{\vartheta {\chi }_{h}}{\left(\frac{1}{\sqrt{{\Lambda}}}\right)}^{\frac{{c}_{m}{\chi }_{h}}{6}},\end{equation} \tag{ 8.21 }$$

where ${\tilde {\mathcal{N}}}_{h}$ is again independent of Λ. For genus zero (8.21) is reminiscent of the partition function of a CFT of central charge c_m on a round two-sphere with Ricci scalar $\tilde {R}=2{\Lambda}$ [116], and consequently an entanglement entropy [195, 196, 200]. For a genus h Riemann surface (8.21) is reminiscent of the partition function obtained from Euclidean AdS₃/CFT₂ considerations at large positive central charge (see for instance [115]).

8.3.1. A gravitational conformal weight

When coupling a conformal field theory to gravity, one is naturally led to the question of what replaces the conformal weights of a spinless conformal primary ${\mathcal{O}}_{{\Delta}}$. One way to quantify this is by computing the following quantity

$$\begin{align}\hfill {\langle {\mathcal{O}}_{{\Delta}}\rangle }_{\upsilon }& =\frac{1}{{\tilde {\mathcal{Z}}}_{{{\Sigma}}_{h}}\left[\upsilon ;{\Lambda}\right]}\hfill \\ \hfill & {\times}{\mathrm{e}}^{\vartheta \enspace {\chi }_{h}}\enspace \int \frac{\mathcal{D}\varphi }{\text{vol}\enspace \mathcal{G}}\mathcal{D}X\enspace {\mathrm{e}}^{-{S}_{L}\left[\varphi ,{\tilde {g}}_{ij}\right]-{S}_{\text{matter}}\left[X,{\tilde {g}}_{ij}\right]}\enspace \delta \enspace \left({\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\enspace {\mathrm{e}}^{2b\varphi }-\upsilon \right){\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\enspace {\mathrm{e}}^{2\sigma \varphi }{\mathcal{O}}_{{\Delta}},\hfill \end{align} \tag{ 8.22 }$$

where σ was given in (8.12) and b in (8.11). Once again, using the technique of [111, 112], we can shift φ → φ + log υ^1/(2b) to conclude

$$\begin{equation}{\langle {\mathcal{O}}_{{\Delta}}\rangle }_{\upsilon }={\mathcal{N}}_{h}\enspace {\upsilon }^{\sigma /b},\end{equation} \tag{ 8.23 }$$

with a υ-independent normalisation constant ${\mathcal{N}}_{h}$. Rescaling the area parameter υ → k²υ we find ${\langle {\mathcal{O}}_{{\Delta}}\rangle }_{\upsilon }\to {k}^{2\sigma /b}{\langle {\mathcal{O}}_{{\Delta}}\rangle }_{\upsilon }$. In an ordinary two-dimensional conformal field theory, an integrated spinless conformal operator with weights $\left(\tilde {{\Delta}},\tilde {{\Delta}}\right)$ would scale as ${k}^{2\left(1-\tilde {{\Delta}}\right)}$. We are thus prompted to define the quantity Δ_grav ≡ 1 − σ/b as the gravitational analogue of the conformal weight. Explicitly

$$\begin{equation}{{\Delta}}_{\text{grav}}=\frac{\sqrt{24{\Delta}+1-{c}_{m}}-\sqrt{1-{c}_{m}}}{\sqrt{25-{c}_{m}}-\sqrt{1-{c}_{m}}},\end{equation} \tag{ 8.24 }$$

which is our second KPZ relation. For the identity operator both Δ = 0 and Δ_grav = 0. We also have that Δ_grav = 1 when Δ = 1. For large and negative c_m we have ${{\Delta}}_{\text{grav}}={\Delta}+\mathcal{O}\left(1/{c}_{m}\right)$. More generally, the two differ.

We are now in a position to connect the continuum description developed in this section to the matrix models we examined throughout the previous ones.

Single matrix integrals. Let us consider first the integrals over single matrices considered in section 3. There, we argued for the existence of a continuum limit in which the expectation number ⟨n_h⟩ defined in (3.78) diverges as we approach a critical coupling. We presented evidence that near criticality and at large N

$$\begin{equation}{\mathcal{F}}_{N}\left(\alpha \right)\approx \sum _{h=0}^{\infty }{f}_{h}\enspace {\mathrm{e}}^{{\chi }_{h}\mathrm{log}N}{\left(\alpha -{\alpha }_{\mathrm{c}}\right)}^{5{\chi }_{h}/4}.\end{equation} \tag{ 8.25 }$$

How should we compare the above to our continuum expressions? We might expect that the sum (3.71) should be viewed as the discrete version of the pure geometry path integral (8.2), with ϑ = log N. In such a case, we should compare the critical limit (3.79) to (8.18). Using the c_m = 0 result for Q/2b, and recalling (3.79), we conclude that the cosmological constant of the continuum theory Λ is propotional to the deviation from criticality, i.e. Λ ∝ (α − α_c). Notice that we recover the logarithmic behaviour in (4.28) for h = 1.

That the continuum cosmological constant is defined as a deviation from the critical coupling α_c is part of a recurring theme. It is the same theme that led us to consider physics just above the filled Fermi sea.

Multicritical matrix model. It has been argued [44, 66, 73] that multicritical matrix models correspond to two-dimensional gravity coupled to a non-unitary minimal CFT (8.13) with (p, q)= (2, 2m − 1) and string susceptibility Γ_str = −1/m. For instance, the multicritical model discussed in (3.49) corresponds to m = 3: i.e. two-dimensional gravity coupled to the Lee–Yang minimal model with c_m = −22/5. For these models, one has to use a slightly modified version [73] of the KPZ relation due to the presence of negative weight operators in the spectrum of the matter CFT. It is remarkable that tuning the parameters of V(M) for a single matrix leads to a continuum description of gravity coupled to conformal field theories with varying central charges.

Double matrix integrals. Given (5.34) and using Λ ∝ (α − α_c) one is led to a theory with Q/(2b) = 7/6. Looking at (8.20) we see that this occurs for c_m = 1/2. It is natural to propose [82] that the two matrix model corresponds to c_m = 1/2 matter, and in particular two-dimensional gravity coupled to a free massless fermion. The free massless fermion in two dimensions is a minimal model which has three spinless conformal primaries, namely the identify operator with Δ₀ = 0, the σ-operator with Δ_σ = 1/16, and the energy operator with Δ = 1/2. These operators are sourced by the cosmological constant Λ, the magnetic field H and the temperature β. These sources make their appearance in the matrix model as simple functions of the couplings α, g, and c of the potential (5.40). Using (8.24) we infer the gravitational conformal weight of these primaries to be Δ_grav,0 = 0, Δ_grav,σ = 1/6, and Δ_grav, = 2/3. These were obtained from considerations of the generalised two-matrix model (5.40) in [73] where it was found that there are certain non-analyticities of the large N two-matrix integral with potential (5.40) upon tuning α, c, and g. The critical coefficients stemming from these non-analyticities are expressed in terms of the gravitational conformal weights Δ_grav,0, Δ_grav,, and Δ_grav,σ via the scaling relations of the two-dimensional Ising model. Concretely, the critical exponents associated to c and g, which we denote by α and β respectively, are related to the gravitational conformal weights as

$$\begin{equation}\boldsymbol{\alpha }=\frac{1-2{{\Delta}}_{\text{grav},{\epsilon}}}{1-{{\Delta}}_{\text{grav},{\epsilon}}}=-1,\qquad \boldsymbol{\beta }=\frac{{{\Delta}}_{\text{grav},\sigma }}{1-{{\Delta}}_{\text{grav},{\epsilon}}}=\frac{1}{2}.\end{equation} \tag{ 8.26 }$$

It might be worth recalling that a free massless fermion describes the critical behaviour of the two-dimensional Ising model. The Ising model has two states at every point on the lattice. This may remind us, somewhat, of the decorated Feynman diagrams discussed in section 5.3. This bears some truth, but caution should be exercised before decorating the Feynman diagrams too elaborately. For instance, the continuum limit of mulitcritical single matrix models corresponds to non-unitary minimal models with c_m < 0 coupled to gravity. It has been further argued [84, 85] that the whole remaining family of minimal models can be obtained from particularly selected two-matrix models.

Upon coupling a minimal model to two-dimensional quantum gravity one obtains a new set of dressed and integrated operators generalising to the area operator in (8.22). In addition to their gravitional conformal weights (8.24), one can also study the correlation functions of these operators and compare them to quantities obtained from matrix integral calculations. Work in this direction includes [125–129].

Non-perturbative features from the continuum. In our discussion of large N matrices we touched upon the possibility of non-perturbative corrections to the planar expansion. In the double scaling limit, these manifest themselves in terms of small WKB type corrections of the series expansion solutions to the various non-linear differential equations, such as the Painlevé I equation. We are naturally led to the question of whether such non-perturbative terms can be recovered in the continuum picture.

In [141] it was shown that Liouville theory can be studied on the disk topology, which is a manifold with a boundary. In this case, the Liouville field diverges near the boundary of the disk. Taking the metric ${\tilde {g}}_{ij}$ on the disk to be ds² = dρ² + ρ²dθ² with ρ ∈ (0, 1), in the semiclassical b → 0 limit the solution of φ is given by

$$\begin{equation}{\mathrm{e}}^{2b{\varphi }_{cl}}=\frac{1}{\pi {b}^{2}{\Lambda}}\frac{1}{{\left(1-{\rho }^{2}\right)}^{2}},\end{equation} \tag{ 8.27 }$$

such that the physical metric ${g}_{ij}={\mathrm{e}}^{2b\varphi }{\tilde {g}}_{ij}$ is the Poincaré disk. The authors of [141] showed the existence of a family of boundary states labelled by two integers (m, n). These configurations are known as ZZ-branes or instantons. Processes involving Liouville theory on the disk have been argued [101, 143, 144] to be related to non-perturbative features in the matrix models. For instance, [142] recovered the exponent (4.36) and (5.38) of the non-perturbative corrections in the double scaling limit by calculating the Liouville disk partition function. In [124], the ZZ-instantons were interpreted as single eigenvalues sitting on critical points of the potential away from the dense eigenvalue distribution.

There is another type of boundary condition that arises in the study of Liouville theory on a disk topology. Here, one adds a boundary interaction [139, 140]

$$\begin{equation}{S}_{\text{bdy}}={\int }_{{S}^{1}}\mathrm{d}u\sqrt{h}\left(\frac{Q}{2\pi }K\varphi +{{\Lambda}}_{\mathrm{B}}\enspace {\mathrm{e}}^{b\varphi }\right)\end{equation} \tag{ 8.28 }$$

to the Liouville action (8.5). The boundary cosmological constant Λ_B is a new continuous parameter labelling the states in this theory. The induced metric at the S¹ boundary of the disk is h and K is the extrinsic curvature at the boundary. Λ_B can be viewed as a chemical potential for the size of the boundary of the physical metric on a disk topology. For Λ_B > 0, the boundary action S_bdy suppresses configurations at large values of φ. These configurations are known as FZZT branes and can be viewed as extended across the weakly coupled region φ ≲ −(log Λ_B)/b. From the matrix model perspective, an interpretation of the FZZT brane should be an object that is parameterised by a continuous variable. Once again, one may ask whether FZZT branes are related to non-perturbative features of the matrix integral. In [124], the FZZT branes were argued to be described by the single trace $\mathrm{T}\mathrm{r}\enspace \mathrm{log}\enspace \left({{\Lambda}}_{\mathrm{B}}-M\right)$. Indeed, the relation of $\mathrm{T}\mathrm{r}\enspace \mathrm{log}\enspace \left({{\Lambda}}_{\mathrm{B}}-M\right)$ to triangulated Riemann surfaces with a boundary was explored in section 3.5.

Quantum mechanical matrices. In the case of quantum mechanical matrices, we uncovered a rich structure in the form of a two-dimensional S-matrix. We now turn to its realisation in the continuum picture.

We now discuss how to recover the scattering discussed in the previous section from the continuum theory, in particular the one-to-one and one-to-two amplitudes. In order to do so, we must first understand what is being scattered in the continuum picture. For instance the two-dimensional theory (8.10) lives on a Euclidean compact surface Σ_h which has no asymptotic regions to scatter to and from. If, on the other hand, we view the two-dimensional theory as the worldsheet theory of a string propagating some target space, there may be room for a scattering amplitude in the assymptotia of the target space. To this end we take the c_m = 1 theory to be a free scalar field denoted by X⁰. In particular for c_m = 1 we conclude using (8.11) that b = 1 and Q = 2. Our worldsheet theory is then given by

$$\begin{equation}{S}_{\text{w.s.}}=-\frac{1}{4\pi }{\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\enspace {\tilde {g}}^{ij}\enspace {\partial }_{i}{X}^{0}{\partial }_{j}{X}^{0}+\frac{1}{4\pi }{\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\left({\tilde {g}}^{ij}{\partial }_{i}\varphi {\partial }_{j}\varphi +2\tilde {R}\enspace \varphi +4\pi {\Lambda}\enspace \enspace {\mathrm{e}}^{2\varphi }\right).\end{equation} \tag{ 8.29 }$$

The free scalar X⁰ encodes the time coordinate of the target space, and as such comes with a kinetic term with the wrong sign. The Liouville mode φ is then associated to the spatial coordinate of the target space, such that the target space is two-dimensional. Taking the contribution of the ghost conformal field theory into account, the net theory has vanishing central charge.¹¹ A string propagating in two-dimensions has no transverse excitations. Consequently, its spectrum is given by a single scalar field encoding the target space centre of mass position of the string. The natural asymptotic region in the target space is given by the region φ → −∞, where the Liouville interaction is switched off.

We can now compute the S-matrix elements as is done for more familiar worldsheet theories. We calculate the expectation value of several vertex operator insertions. For a free scalar we can construct conformal operators by exponentiating X⁰. It follows from our discussion in section 8.2, that we can construct the dressed integrated vertex operator

$$\begin{equation}{\mathcal{V}}_{\omega }^{{\pm}}={g}_{s}{\int }_{{{\Sigma}}_{h}}{\mathrm{d}}^{2}x\sqrt{\tilde {g}}\enspace :{\mathrm{e}}^{{\pm}\mathrm{i}\omega {X}^{0}}:\mathcal{O}\left(\omega \right),\quad \omega {\geqslant}0.\end{equation} \tag{ 8.30 }$$

In the above, $\mathcal{O}\left(\omega \right)$ is taken to be a Hermitian Liouville operator of the form

$$\begin{equation}\mathcal{O}\left(\omega \right)\equiv S\left(\omega \right){\mathrm{e}}^{\left(2+\mathrm{i}\omega \right)\varphi }+S{\left(\omega \right)}^{{\ast}}{\mathrm{e}}^{\left(2-\mathrm{i}\omega \right)\varphi },\quad \omega {\geqslant}0.\end{equation} \tag{ 8.31 }$$

The phase factors have been chosen such that the two-point function of the Liouville operator $\mathcal{O}\left(\omega ;z,\overline{z}\right)$ has standard δ-function normalisation

$$\begin{equation}\langle \mathcal{O}\left({\omega }_{1};{z}_{2},{\overline{z}}_{2}\right)\mathcal{O}\left({\omega }_{2};{z}_{2},{\overline{z}}_{2}\right)\rangle =2\pi \delta \left({\omega }_{1}-{\omega }_{2}\right){\times}\frac{1}{\vert {z}_{1}-{z}_{2}{\vert }^{4+{\omega }_{1}^{2}}},\end{equation} \tag{ 8.32 }$$

where z and $\overline{z}$ label points on $\mathbb{C}$, and we adhere to the conventions of [29] with α' = 1. The phase factor is given by [114]

$$\begin{equation}S\left(\omega \right)=\sqrt{\frac{{\Gamma}\left(\mathrm{i}\omega \right)}{{\Gamma}\left(1-\mathrm{i}\omega \right)}\frac{{\Gamma}\left(1+\mathrm{i}\omega \right)}{{\Gamma}\left(-\mathrm{i}\omega \right)}}.\end{equation} \tag{ 8.33 }$$

To compute target space S-matrix elements we must calculate expectation values for the above vertex operators. The parameter ω corresponds to the energy of the asymptotic state, and the sign determines whether it is incoming (negative) or outgoing (positive).¹²

One-to-one from the continuum. We begin by considering the one-to-one amplitude in the continuum picture. For this, we require the expectation value for two vertex operator insertions on the two-sphere. Given the two-point function (8.32), and following the discussion in [135], one finds

$$\begin{equation}{A}_{{S}^{2}}\left({\omega }_{1}\vert {\omega }_{2}\right)=2\pi {\omega }_{1}\delta \left({\omega }_{1}-{\omega }_{2}\right){\times}\left(\frac{{g}_{s}^{2}{c}_{{S}^{2}}}{4\pi }\right),\end{equation} \tag{ 8.34 }$$

where ${c}_{{S}^{2}}$ is a coefficient which we will soon relate to the one-to-two amplitude. Comparing to (7.23) we fix ${c}_{{S}^{2}}=4\pi /{g}_{s}^{2}$.

One-to-two from the continuum. Next we consider three insertions on a genus zero surface

$$\begin{equation}{A}_{{S}^{2}}\left({\omega }_{1},{\omega }_{2}\enspace \vert \enspace {\omega }_{3}\right)\equiv {\langle c\tilde {c}{\mathcal{V}}_{{\omega }_{1}}^{+}\enspace c\tilde {c}{\mathcal{V}}_{{\omega }_{2}}^{+}\enspace c\tilde {c}{\mathcal{V}}_{{\omega }_{3}}^{-}\rangle }_{{S}^{2}}.\end{equation} \tag{ 8.35 }$$

We follow the presentation in [133]. To evaluate the above expression requires knowledge of the three-point function for the Liouville operator $\mathcal{O}\left(\omega ;z,\overline{z}\right)$

$$\begin{align}\hfill & \langle \mathcal{O}\left({\omega }_{1};{z}_{1},{\overline{z}}_{1}\right)\enspace \mathcal{O}\left({\omega }_{2};{z}_{2},{\overline{z}}_{2}\right)\enspace \mathcal{O}\left({\omega }_{3};{z}_{3},{\overline{z}}_{3}\right)\rangle \hfill \\ \hfill & \qquad =\frac{\mathcal{C}\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\right)}{\vert {z}_{12}{\vert }^{2+\left({\omega }_{1}^{2}+{\omega }_{2}^{2}-{\omega }_{3}^{2}\right)/2}\vert {z}_{23}{\vert }^{2+\left({\omega }_{2}^{2}+{\omega }_{3}^{2}-{\omega }_{1}^{2}\right)/2}\vert {z}_{31}{\vert }^{2+\left({\omega }_{3}^{2}+{\omega }_{1}^{2}-{\omega }_{2}^{2}\right)/2}}.\hfill \end{align} \tag{ 8.36 }$$

This quantity was studied by Dorn–Otto [120] and Zamolodchikov–Zamolodchikov [122] which gives it the name DOZZ formula. The structure constant $\mathcal{C}\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\right)$ is given by

$$\begin{equation}\mathcal{C}\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\right)\equiv {\left(S\left({\omega }_{1}\right)S\left({\omega }_{2}\right)S\left({\omega }_{3}\right)\right)}^{-1}C\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\right),\end{equation} \tag{ 8.37 }$$

with the DOZZ coefficient

$$\begin{align}\hfill C\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\right)=& \enspace \frac{{{\Upsilon}}_{1}^{\prime }\left(0\right)}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{3}+{\omega }_{2}+{\omega }_{1}\right)/2\right)}\frac{{{\Upsilon}}_{1}\left(2+\mathrm{i}{\omega }_{1}\right)}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{1}+{\omega }_{2}-{\omega }_{3}\right)/2\right)}\hfill \\ \hfill & {\times}\enspace \frac{{{\Upsilon}}_{1}\left(2+\mathrm{i}{\omega }_{2}\right)}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{2}+{\omega }_{3}-{\omega }_{1}\right)/2\right)}\frac{{{\Upsilon}}_{1}\left(2+\mathrm{i}{\omega }_{3}\right)}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{3}+{\omega }_{1}-{\omega }_{2}\right)/2\right)}.\hfill \end{align} \tag{ 8.38 }$$

In the literature the three-point function contains an additional phase has been absorbed in a shift of the Liouville field. The function ϒ_b(z) is a holomorphic function that admits an integral expression as

$$\begin{equation}\mathrm{log}\enspace {{\Upsilon}}_{b}\left(z\right)={\int }_{0}^{\infty }\frac{\mathrm{d}t}{t}\left[{\left(\frac{Q}{2}-z\right)}^{2}\enspace {\mathrm{e}}^{-t}-\frac{{\mathrm{sinh}}^{2}\left[\left(\frac{Q}{2}-z\right)\frac{t}{2}\right]}{\mathrm{sinh}\enspace \frac{tb}{2}\enspace \mathrm{sinh}\enspace \frac{t}{2b}}\right],\end{equation} \tag{ 8.39 }$$

where the real part of z is restricted to the interval z ∈ (0, Q) and b + b⁻¹ = Q. To evaluate ϒ_b(z) at complex values of z we must analytically continue the above expression. It is useful to note that we can express ϒ_b(z) in terms of the Barnes double Gamma-function Γ₂(z|b₁, b₂):

$$\begin{equation}{{\Upsilon}}_{b}\left(z\right)=\frac{1}{{{\Gamma}}_{b}\left(z\right){{\Gamma}}_{b}\left(Q-z\right)},\qquad {{\Gamma}}_{b}\left(z\right)\equiv \frac{{{\Gamma}}_{2}\left(z\enspace \vert \enspace b,{b}^{-1}\right)}{{{\Gamma}}_{2}\left(Q/2\enspace \vert \enspace b,{b}^{-1}\right)}.\end{equation} \tag{ 8.40 }$$

In particular we have Γ₂(z|1, 1) = (2π)^z/2/G(z) with G(z) being the Barnes G-function. ϒ_b(z) also satisfies

$$\begin{equation}{{\Upsilon}}_{b}\left(z+b\right)=\gamma \left(bz\right)\enspace {b}^{1-2bz}\enspace {{\Upsilon}}_{b}\left(z\right),\qquad {{\Upsilon}}_{b}\left(z\right)={{\Upsilon}}_{{b}^{-1}}\left(z\right),\qquad {{\Upsilon}}_{b}\left(Q-z\right)={{\Upsilon}}_{b}\left(z\right),\end{equation} \tag{ 8.41 }$$

with γ(z) ≡ Γ(z)/Γ(1 − z). The phase S(ω) in (8.33) can be expressed as S(ω) = γ(iω)^1/2γ(1 + iω)^1/2.

Additionally (8.30) contains a piece stemming from the operator ${\mathrm{e}}^{{\pm}\mathrm{i}\omega {X}^{0}}$. The structure constants for this piece (as may be familiar from the critical string) will not produce non-trivial functions of the ω_i. Nevertheless, due to the symmetry X⁰ → X⁰ + a, it yields a delta-function imposing the conservation of target space energy.

At this point we have asembled the tools required to evaluate (8.35). Using the relations (8.41) for (Q, b) = (2, 1) we can rewrite ϒ₁(2 + iω_k) as

$$\begin{equation}{{\Upsilon}}_{1}\left(2+\mathrm{i}{\omega }_{k}\right)={{\Upsilon}}_{1}{\left(2+\mathrm{i}{\omega }_{k}\right)}^{1/2}{{\Upsilon}}_{1}{\left(2+\mathrm{i}{\omega }_{k}\right)}^{1/2}={{\Upsilon}}_{1}{\left(-\mathrm{i}{\omega }_{k}\right)}^{1/2}{{\Upsilon}}_{1}{\left(\mathrm{i}{\omega }_{k}\right)}^{1/2}\gamma {\left(\mathrm{i}{\omega }_{k}\right)}^{1/2}\gamma {\left(-\mathrm{i}{\omega }_{k}\right)}^{-1/2},\end{equation} \tag{ 8.42 }$$

with γ(x) defined below (8.41) satisfying in particular $\gamma \left(1+\mathrm{i}{\omega }_{k}\right)=\gamma {\left(-\mathrm{i}{\omega }_{k}\right)}^{-1}$. This way we obtain for (8.37)

$$\begin{align}\hfill \mathcal{C}\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\right)& =\frac{1}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{1}+{\omega }_{2}+{\omega }_{3}\right)/2\right)}\left(\frac{{{\Upsilon}}_{1}{\left(\mathrm{i}{\omega }_{1}\right)}^{1/2}{{\Upsilon}}_{1}{\left(-\mathrm{i}{\omega }_{1}\right)}^{1/2}}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{1}+{\omega }_{2}-{\omega }_{3}\right)/2\right)}{\times}\mathrm{2}\enspace \mathrm{p}\mathrm{e}\mathrm{r}\mathrm{m}.\right)\hfill \\ \hfill & =\frac{1}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{1}+{\omega }_{2}+{\omega }_{3}\right)/2\right)}\left(\frac{{\omega }_{1}{{\Upsilon}}_{1}\left(1+\mathrm{i}{\omega }_{1}\right)}{{{\Upsilon}}_{1}\left(1+\mathrm{i}\left({\omega }_{1}+{\omega }_{2}-{\omega }_{3}\right)/2\right)}{\times}\mathrm{2}\enspace \mathrm{p}\mathrm{e}\mathrm{r}\mathrm{m}.\right),\hfill \end{align} \tag{ 8.43 }$$

where we made use of (8.40) to evaluate ϒ₁'(0) = 1 and in going to the second line we used

$$\begin{equation}{{\Upsilon}}_{1}\left(-\mathrm{i}{\omega }_{k}\right)={{\Upsilon}}_{1}\left(2+\mathrm{i}{\omega }_{k}\right)=\gamma \left(1+\mathrm{i}{\omega }_{k}\right){{\Upsilon}}_{1}\left(1+\mathrm{i}{\omega }_{k}\right).\end{equation} \tag{ 8.44 }$$

Putting everything together, and in particular applying the energy conservation obtained from the delta function as well as ϒ₁(1) = 1, we have that

$$\begin{equation}{A}_{{S}^{2}}\left({\omega }_{1},{\omega }_{2}\enspace \vert \enspace {\omega }_{3}\right)=\mathrm{i}\enspace {c}_{{S}^{2}}\enspace {g}_{s}^{3}\enspace \delta \left({\omega }_{1}+{\omega }_{2}-{\omega }_{3}\right)\mathcal{C}\left({\omega }_{1},{\omega }_{2},{\omega }_{3}\right)=\mathrm{i}\enspace {c}_{{S}^{2}}\enspace {g}_{s}^{3}\enspace \delta \left({\omega }_{1}+{\omega }_{2}-{\omega }_{3}\right)\enspace {\omega }_{1}{\omega }_{2}{\omega }_{3},\end{equation} \tag{ 8.45 }$$

where ${c}_{{S}^{2}}$ is a non-vanishing constant. Comparing to (7.28), we see that the amplitude matches the leading term in the large ν_F expansion [25–27], so long as we identify $1/\vert {\nu }_{\mathrm{F}}\vert ={c}_{{S}^{2}}{g}_{s}^{3}=4\pi {g}_{s}$.

The terms (7.28) which are subleading in 1/ν_F are predictions about higher genus contributions to the one-to-two scattering process in the continuum picture. We note that this identification is consistent with the powers of ν_F appearing in the various scattering amplitudes computed in section 7.2. A scattering amplitude on a Riemann surface of genus h and n vertex operator insertions scales as ${g}_{s}^{2h-2+n}$. A similar picture holds for other perturbative scattering amplitudes.

Further remarkable comparisons between the quantum mechanical matrix model S-matrix and Liouville theory have been obtained more recently in [133]. It is interesting to note that the whole structure of the perturbative string amplitudes is encoded in the reflection coefficient (7.5). Though manifest from the quantum mechanical matrix model, the explicit form of R_ν is far from obvious from the Liouville perspective.

8.5.1. Two-dimensional target space picture

To end the section, we would like to make some remarks about the target space picture of the non-critical string theory (8.29). For φ → −∞ the exponential interaction of the Liouville mode is switched off, and the theory is approximated by the combination of a free scalar and a linear dilaton CFT. From this perspective, at least in the limit of large and negative φ, we have a two-dimensional target space with a running dilaton Φ. The free scalar encodes the time direction of the target space, whereas φ encodes the spatial direction. The target space metric and dilaton field are given by

$$\begin{equation}\frac{\mathrm{d}{s}^{2}}{{\alpha }^{\prime }}=-{\left(\mathrm{d}{X}^{0}\right)}^{2}+\mathrm{d}{\varphi }^{2},\qquad {\Phi}\left({X}^{0},\varphi \right)=2\varphi ,\end{equation} \tag{ 8.46 }$$

such that the space-dependent string coupling g_s = e^Φ becomes parameterically weak as φ → −∞ and increasingly strong in the positive φ region. We also have reinstated the string tension α'.

Due to the space-dependence, the vacuum is not Poincaré invariant. As φ increases, so does the string coupling, and we moreover enter a regime where the Liouville interaction becomes significant. This regime is often referred to as the Liouville wall in the literature. As was mentioned, the spectrum of a string embedded in a two-dimensional target space consists of a single propagating degree of freedom, which is the 'tachyon field' τ. It is the target space field corresponding to the vertex operator (8.30). In two-dimensions, τ is in fact a massless field, so tachyon is somewhat of a misnomer. From the target space perspective, the S-matrix is given by sending an incoming tachyon wave from φ → −∞, which reflects off the Liouville wall, finally heading back to φ → −∞. Finally, we note that the background (8.46) is intrinsically stringy since the slope of Φ is the string length.

One of the main tools that has allowed for tremendous progress in many corners of modern theoretical physics, ranging from quantum field theory and mathematics all the way to string theory and black hole physics, is supersymmetry. Here we mention the supersymmetric generalisation of some quantum mechanical matrix models. A supersymmetric quantum mechanical theory has a Hamiltonian of the following form

$$\begin{equation}\left\{{\hat{\mathcal{Q}}}^{A},{\hat{\mathcal{Q}}}_{B}^{{\dagger}}\right\}=2{{\delta }^{A}}_{B}{\hat{H}}_{\text{SUSY}},\qquad \left\{{\hat{\mathcal{Q}}}^{A},{\hat{\mathcal{Q}}}_{B}\right\}=0,\end{equation} \tag{ 9.1 }$$

where the ${\hat{\mathcal{Q}}}^{A}$ are supercharges, with A and B being indices denoting additional symmetries or structures. In line with our previous considerations of matrix models, we further assume that the models have an SU(N) symmetry, which we may or may not gauge. We take the bosonic elements of the theory to be as in our previously studied example, namely a collection of Hermitian N × N matrices ${\hat{M}}_{IJ}^{a}$, with a = 1, 2, ..., n_B. In a supersymmetric theory, these are accompanied by adjoint valued fermionic operators ${\hat{{\Theta}}}_{IJ}^{\alpha }$, with α = 1, 2, ..., n_F, obeying the operator algebra

$$\begin{equation}\left\{{\hat{{\Theta}}}_{IJ}^{\alpha },{\hat{{\Theta}}}_{KL}^{\beta }\right\}={\delta }_{IK}{\delta }_{JL}{\delta }^{\alpha \beta }.\end{equation} \tag{ 9.2 }$$

As a reality condition we consider ${\left({\hat{{\Theta}}}_{IJ}^{\alpha }\right)}^{{\dagger}}={\hat{{\Theta}}}_{JI}^{\alpha }$. From the perspective of an action principle, the fermionic operators translate to Grassmann valued functions of time Θ_IJ(t) satisfying ${\left({{\Theta}}_{IJ}\left(t\right)\right)}^{{\ast}}={{\Theta}}_{JI}\left(t\right)$.

Marinari–Parisi models. The simplest class of models for a supersymmetric theory of quantum mechanical matrices was introduced by Marinari and Parisi [157]. They are built from the off-shell field content $\mathcal{M}=\left\{{M}_{IJ},{{\Psi}}_{IJ},{\overline{{\Psi}}}_{IJ},{F}_{IJ}\right\}$, where M_IJ and F_IJ are Hermitian matrices, while ${{\Psi}}_{IJ}={{\Theta}}_{IJ}^{1}+\mathrm{i}{{\Theta}}_{IJ}^{2}$ and ${\overline{{\Psi}}}_{IJ}={{\Theta}}_{IJ}^{1}-\mathrm{i}{{\Theta}}_{IJ}^{2}$ are Grassmann valued complex matrices. The action reads

$$\begin{equation}{S}_{\text{MP}}\left[M\left(t\right)\right]=N\enspace \mathrm{T}\mathrm{r}{\int }_{\mathbb{R}}\mathrm{d}t\left(\frac{1}{2}{\dot {M}}^{2}+\mathrm{i}\overline{{\Psi}}\dot {{\Psi}}-\vert W\left(M\right){\vert }^{2}-\overline{{\Psi}}\enspace {\partial }_{M}W\left(M\right)\enspace {\Psi}\right),\end{equation} \tag{ 9.3 }$$

where we have integrated out the auxiliary field F_IJ, and W(M) is the superpotential, which is itself a Hermitian matrix.

The models (9.3) are distinguished from those built from a single matrix M_IJ in an important way. Recall that the single matrix models admitted a description for the ground state sector that could be reduced to an N variable (rather than N²) problem. Moreover, had we considered gauging the SU(N) symmetry of the single matrix models we would end up again with N degrees of freedom. In the case at hand, Ψ_IJ adds additional degrees of freedom, and theories of the type (9.3) can no longer be mapped to a simple eigenvalue problem even upon gauging the SU(N). There are not enough symmetries. This conclusion remains true for essentially any model comprised of two or more interacting matrices and a single SU(N) symmetry, regardless of whether they are fermionic or bosonic. Thus, new techniques are required to understand such a situation.

To our knowledge, a complete description for the worldsheet theory describing the continuum limit of Marinari–Parisi models remains unknown. We now discuss a supersymmetric quantum mechanical matrix model for which a clearer picture of the emergent worldsheet is available.

D0-brane quantum mechanical matrices. An important supersymmetric quantum mechanical matrix model is the D0-brane quantum mechanical model. The model enjoys an SO(9) global symmetry and its field content is given by an SU(N) adjoint valued vector multiplet $\mathcal{V}=\left\{{M}_{IJ}^{a},{{\Theta}}_{IJ}^{\alpha },{A}_{IJ}\right\}$. The index a = 1, 2, ..., 9 is an SO(9) vector index, while α = 1, ..., 16 is a spinor index of the 16-dimensional representation of SO(9). The field A_IJ is a gauge field for the SU(N) symmetry. The role of A_IJ is solely to impose the gauge constraint, and carries no dynamics of its own. The system of interest has an interpretation stemming from considerations of type IIA superstring theory. It describes the low energy theory governing the open string sector of a collection of N D0-branes [159]. The theory is governed by the action

$$\begin{equation}{S}_{\text{D0}}=\frac{N}{{\lambda }^{2}}\enspace \mathrm{T}\mathrm{r}{\int }_{\mathbb{R}}\mathrm{d}t\left(\frac{1}{2}{\dot {M}}^{a}{\dot {M}}^{a}+\frac{\mathrm{i}}{2}{{\Theta}}^{T}\dot {{\Theta}}+\frac{1}{4}{\left[{M}^{a},{M}^{b}\right]}^{2}-\frac{1}{2}{{\Theta}}^{T}{{\Gamma}}^{a}\left[{\Theta},{M}^{a}\right]\right),\end{equation} \tag{ 9.4 }$$

where we have set A_IJ = 0, and must recall the gauge singlet constraint. The Γ^a are 16 × 16 dimensional matrices satisfying the Clifford algebra

$$\begin{equation}\left\{{{\Gamma}}^{a},{{\Gamma}}^{b}\right\}=2\enspace {\delta }^{ab}.\end{equation} \tag{ 9.5 }$$

In addition to the global SO(9) symmetry, the above model has 16 supercharges. The dimensionful parameter λ is the 't Hooft coupling of the theory. The theory is effectively weakly coupled at high energies, and strongly coupled at low energies.

There is a remarkable proposal originally developed by Banks, Fischler, Shenker, and Susskind (BFSS) relating the above theory in the large N limit to quantum gravity in the ten-dimensional geometry residing in the vicinity of a stack of N D0-branes [159]. The string frame metric and dilaton of this geometry are described by the SO(9) invariant solution

$$\begin{equation}\frac{\mathrm{d}{s}^{2}}{{\alpha }^{\prime }}=-\frac{{u}^{7/2}\mathrm{d}{t}^{2}}{\sqrt{a\lambda }}+\sqrt{a\lambda }\left(\frac{\mathrm{d}{u}^{2}}{{u}^{7/2}}+\frac{\mathrm{d}{{\Omega}}_{8}^{2}}{{u}^{3/2}}\right),\qquad {\mathrm{e}}^{{\Phi}}=\frac{{\left(2\pi \right)}^{2}{a}^{3/4}}{N}{\left(u\enspace {\lambda }^{-1/3}\right)}^{-21/4},\end{equation} \tag{ 9.6 }$$

with a = 240π⁵. In addition to the metric and dilaton, the solution contains N units Ramond–Ramond two-form flux sourced by the D0-branes. It is also worth noting there have been several variations of the D0-brane model, such as the BMN model [164] which describes a massive deformation of (9.4), whereby the flat direction of mutually commuting M^a's is lifted.

A particular incarnation of the BFSS proposal [162], which is in the same spirit as AdS/CFT, claims that quantum gravity in the background (9.6) is dual to the D0-brane quantum mechanics theory. If true, the proposal surpasses the c > 1 barrier in a remarkable way. Unravelling the details of the BFSS proposal, by whatever means need be [165–170], remains a fascinating challenge in string theory.

d-dimensional models? One might ask about simple supersymmetric generalisations of (9.4) which have an SO(d) symmetry. Let Γ^a be the s_d × s_d dimensional generators of the s_d dimensional, real representation of the Clifford algebra satisfying the Clifford relation

$$\begin{equation}\left\{{{\Gamma}}^{a},{{\Gamma}}^{b}\right\}=2\enspace {\delta }^{ab},\quad a,b=1,\dots ,d.\end{equation} \tag{ 9.7 }$$

The d-dimensional Clifford algebra contains irreducible representations with dimension

$$\begin{equation}{s}_{d}=\begin{cases}^{\lfloor \frac{d}{2}\rfloor },\quad d=0,1,2\enspace \mathrm{mod}\enspace 8\quad \hfill \\ {2}^{\lfloor \frac{d}{2}\rfloor +1},\quad \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\quad \hfill \end{cases}\end{equation} \tag{ 9.8 }$$

A simple degree of freedom counting informs us that any supersymmetric theory built from a vector multiplet $\mathcal{V}=\left\{{M}_{IJ}^{a},{{\Theta}}_{IJ}^{\alpha },{A}_{IJ}\right\}$ requires s_d = 2(d − 1). When combined with (9.8), this implies that only d = 2, 3, 5, 9, and correspondingly s_d = 2, 4, 8, 16 are valid dimensions [156, 163]. Perhaps the difficulty in constructing models with d > 9 is related to the absence of perturbative superstring theories in larger than ten spacetime dimensions.

Another natural class of matrix models consists of matrices living on a spacetime manifold $\mathcal{M}$, as opposed to being merely functions of time. An example of such a theory is Yang–Mills theory with SU(N) gauge group.

Perhaps the simplest example of such a theory is Yang–Mills theory in two spacetime dimensions, described by the Euclidean action

$$\begin{equation}{S}_{\text{YM}}=\frac{1}{4{g}^{2}}{\int }_{\mathcal{M}}{\mathrm{d}}^{2}x\sqrt{g}\enspace {F}_{ij}{F}^{ij},\end{equation} \tag{ 9.9 }$$

where F_ij is the field strength for the gauge field A_i and $\mathcal{M}$ is a two-dimensional manifold endowed with metric g_ij. The theory carries no local degrees of freedom, but it has a non-vanishing Hilbert space. The large N loop equations for two-dimensional Yang–Mills theory have been studied in [174].

If one quantises the theory on a spatial S¹, it can be shown that the theory is mapped to a quantum mechanical matrix model, but in this case the matrices are unitary rather than Hermitian [171–173, 178]. Such unitary matrices admit a treatment similar to that of Hermitian matrices, except the resulting fermion theory now lives on the circle rather than the half-line.

It has been proposed by Gross and Taylor [176] that two-dimensional Yang–Mills theory is described by a worldsheet theory in the large N limit. Although much remains to be understood about this worldsheet theory, it seems clear that it is not of the Liouville type.

Perhaps the second simplest quantum field theory built from adjoint matrices is the Chern–Simons theory in three-dimensions. The action of this theory is given by

$$\begin{equation}{S}_{\text{CS}}=\frac{k}{4\pi }\enspace \mathrm{T}\mathrm{r}{\int }_{\mathcal{M}}\left(A\wedge \mathrm{d}A+\frac{2}{3}A\wedge A\wedge A\right).\end{equation} \tag{ 9.10 }$$

The above theory has the remarkable property of being a topological theory, insensitive to the geometric features of $\mathcal{M}$. When quantised on a compact spatial two-manifold the above theory has a finite dimensional Hilbert space. For concreteness we consider an U(N) gauge group.

The standard perturbative expansion of Chern–Simons theory is given by taking the level $k\in \mathbb{Z}$ to be large while keeping N fixed. On the other hand, it is natural to explore the behaviour of Chern–Simons theory in the 't Hooft limit for which we take N large while keeping λ = N/(N + k) fixed. Indeed, as argued by Gopakumar and Vafa [180, 182], the 't Hooft limit of Chern–Simons theory on a three-sphere with U(N) gauge group is captured by the topological closed string on a resolved conifold background. The target space of the worldsheet theory has a Ramond–Ramond two-form B, a Kähler two-form J, and string coupling g_s. These are related to the Chern–Simons parameters by

$$\begin{equation}\lambda =-\frac{1}{2\pi }{\int }_{{S}^{2}}\left(B+\mathrm{i}J\right),\qquad {g}_{s}=\frac{2\pi \lambda }{N}.\end{equation} \tag{ 9.11 }$$

Notice that although λ was originally a real number, the above identifications require extending its domain to parts of the complex plane. It is also worth noting that the partition function of Chern–Simons theory can be mapped to an ordinary matrix intergral [184], rather than a quantum mechanical matrix model.

It is worth noting that the topological string is substantially different than the Liouville string discussed in the previous section. Moreover, the duality between Chern–Simons and the topological string does not require tuning the 't Hooft coupling λ to some critical value.

When discussing supersymmetric models we introduced a new set of fermionic operators ${\hat{{\Psi}}}_{IJ}^{\alpha }={\hat{{\Theta}}}_{IJ}^{\alpha }+\mathrm{i}\enspace {\hat{{\Upsilon}}}_{IJ}^{\alpha }$. This leads to yet another class of quantum mechanical matrix models built entirely from the fermions. As for Chern–Simons theory quantised on a compact Riemann surface, the Hilbert space of these modes is finite dimensional for any finite value of N—namely $\mathrm{dim}\enspace {\mathcal{H}}_{N}={2}^{{N}^{2}}$. The finiteness of the Hilbert space distinguishes these models from their bosonic and supersymmetric cousins. A simple class of models [188–191, 194] are described by a quartic Hamiltonian of the general type

$$\begin{equation}{\hat{H}}_{\text{F}}=\frac{N}{\lambda }\enspace \mathrm{T}\mathrm{r}\enspace {\hat{{\Psi}}}^{{\dagger}}{{\Gamma}}^{a}\hat{{\Psi}}\enspace {\hat{{\Psi}}}^{{\dagger}}{{\Gamma}}^{a}\hat{{\Psi}},\end{equation} \tag{ 9.12 }$$

where as before the Γ^a are the s_d × s_d dimensional Γ^a matrices.

Although the wordsheet description capturing the 't Hooft limit of the above models remains unknown, simple versions of the above model exhibit features that bear some similarity to those encountered in Chern–Simons theory, as well as two-dimensional Yang–Mills theory. The models (9.12) have been considered both for a gauged as well as ungauged SU(N) symmetry.

Recently, a novel type of large N limit has been exploited to analyse systems which have several features of interest from the perspective of the physics of black holes, as well as condensed matter systems. One unusual ingredient in these models, known as quenched disorder, is most often used in the context of glassy physics. Quenched disorder treats the couplings of the theory as random variables drawn from an ensemble. The simplest class of Hamiltonian's considered [213, 215, 216, 219] are again built purely from real quantum mechanical fermions ${\hat{{\Theta}}}_{I}$ satisfying the algebra $\left\{{\hat{{\Theta}}}_{I},{\hat{{\Theta}}}_{J}\right\}={\delta }_{IJ}$, with I = 1, ..., N. They are

$$\begin{equation}{\hat{H}}_{\text{SYK}}=\frac{{\mathrm{i}}^{q/2}}{q!}{\mathcal{J}}_{{I}_{1}\dots {I}_{q}}{\hat{{\Theta}}}_{{I}_{1}}{\hat{{\Theta}}}_{{I}_{2}}\dots {\hat{{\Theta}}}_{{I}_{q}},\end{equation} \tag{ 9.13 }$$

with totally anti-symmetric real couplings ${\mathcal{J}}_{{I}_{1}\dots {I}_{q}}$ drawn from a random Gaussian ensemble with variance

$$\begin{equation}\langle {\mathcal{J}}_{{I}_{1}\dots {I}_{q}}{\mathcal{J}}_{{I}_{1}^{\prime }\dots {I}_{q}^{\prime }}\rangle =\frac{{\mathcal{J}}^{2}\left(q-1\right)!}{{N}^{q-1}}{\times}{\delta }_{{I}_{1}\dots {I}_{q},{I}_{1}^{\prime }\dots {I}_{q}^{\prime }}\enspace .\end{equation} \tag{ 9.14 }$$

Notice that the theory is built from vector-like constituents, and the Hilbert space is 2^N dimensional. As such, the models do not exhibit a 't Hooft expansion, at least in the standard sense. Nevertheless, the perturbative structure exhibits an interesting structure at large N, often referred to as a melonic expansion. These theories have no continuous global symmetries, and as such bear some relation to matrix models with the SU(N) ungauged.

We can analyse the diagrammatic expansion most easily by expressing the theory in terms of an action. For any given realisation of the couplings, the action is given by

$$\begin{equation}{S}_{\text{SYK}}\left[{{\Theta}}_{I};{\mathcal{J}}_{{I}_{1}\dots {I}_{q}}\right]={\int }_{\mathbb{R}}\mathrm{d}t\left(\frac{\mathrm{i}}{2}{{\Theta}}_{I}{\dot {{\Theta}}}_{I}-\frac{{\mathrm{i}}^{q/2}}{q!}{\mathcal{J}}_{{I}_{1}\dots {I}_{q}}{{\Theta}}_{{I}_{1}}{{\Theta}}_{{I}_{2}}\dots {{\Theta}}_{{I}_{q}}\right).\end{equation} \tag{ 9.15 }$$

To build the diagrams for the theory with randomised couplings, one resorts to the replica trick whereby the theory is replicated to $n\in \mathbb{Z}$ copies and n is then allowed to take real values. This adds an additional 'replica' index to the degrees of freedom. Assuming that the symmetry permuting the replicas is unbroken at large N, one can simply treat the ${\mathcal{J}}_{{I}_{1}\dots {I}_{q}}$ as additional non-dynamical fields whose propagator is dictated by the Gaussian probability distribution (9.14). At large N, a particular class of diagrams dominate. This class is displayed in figure 16.

Zoom In Zoom Out Reset image size

The leading N contribution to the correlation functions can be recast in the form of a Schwinger–Dyson equation which, in turn, exhibits a one-dimensional conformal symmetry at sufficiently large time separations. The conformal symmetry is the space of maps from $\mathbb{R}$ (or S¹ when the theory is placed at finite temperature) to itself. It is suggestive of a two-dimensional anti-de Sitter spacetime black hole holographic dual. This interpretation is further supported by considerations of the symmetry breaking low energy sector which is governed by a Schwarzian mode [217]. It turns out that precisely the same diagrams are those that organise a class of models known as tensor models, which are built from objects with multiple indices [214, 218]. Tensor models have been previously introduced [208, 209] as an interesting generalisation to matrix models. The quantum mechanical tensor models do not have random coupling constants.

The SYK and tensor models are neither matrix models nor vector models. Rather, they are an entirely different class of large N theories. Unlike most quantum mechanical matrix models, they can be dealt with without imposing any gauge singlet constraint. This is perhaps related to their connection to black hole physics. Though not matrix models, it is interesting to note that both matrix integrals as well as two-dimensional gravity seem to stem from detailed considerations of these non-matrix models [221].

As a final note, we remark that (4.46) and (B.6) are part of a set of equations often referred to as the string equations [66]. For the mth multicritical polynomial we have

$$\begin{equation}z=\frac{{\left(-1\right)}^{m}{2}^{m+1}m!}{\left(2m-1\right){\Exclam}}\enspace {\mathcal{R}}_{m}\left[\delta r\left(z\right)\right],\end{equation} \tag{ B.10 }$$

where the ${\mathcal{R}}_{m}\left[\delta r\left(z\right)\right]$ are the Gelfand–Dikii polynomials [57]. These polynomials also appear in Witten's conjecture [19, 20] relating the intersection numbers of moduli spaces to an infinite set of partial differential equations (KdV hierarchy). They are obtained from the recursion formula:

$$\begin{equation}{\mathcal{R}}_{0}\left[\delta r\left(z\right)\right]=\frac{1}{2},\qquad {\mathcal{R}}_{m+1}^{\prime }\left[\delta r\left(z\right)\right]=\frac{1}{4}{\mathcal{R}}_{m}^{{\prime\prime\prime}}\left[\delta r\left(z\right)\right]-\delta r\left(z\right){\mathcal{R}}_{m}^{\prime }\left[\delta r\left(z\right)\right]-\frac{1}{2}\delta {r}^{\prime }\left(z\right){\mathcal{R}}_{m}\left[\delta r\left(z\right)\right].\end{equation} \tag{ B.11 }$$

We list them up to m = 4:

$$\begin{align}\hfill m=0:{\mathcal{R}}_{0}\left[\delta r\left(z\right)\right]& =\frac{1}{2},\hfill \\ \hfill m=1:{\mathcal{R}}_{1}\left[\delta r\left(z\right)\right]& =-\frac{1}{4}\delta r\left(z\right),\hfill \\ \hfill m=2:{\mathcal{R}}_{2}\left[\delta r\left(z\right)\right]& =\frac{1}{16}\left(3\delta r{\left(z\right)}^{2}-\delta {r}^{{\prime\prime}}\left(z\right)\right),\hfill \\ \hfill m=3:{\mathcal{R}}_{3}\left[\delta r\left(z\right)\right]& =-\frac{1}{64}\left(10\delta r{\left(z\right)}^{3}-10\delta {r}^{{\prime\prime}}\left(z\right)\delta r\left(z\right)-5\delta {r}^{\prime }{\left(z\right)}^{2}+\delta {r}^{\left(4\right)}\left(z\right)\right),\hfill \\ \hfill m=4:{\mathcal{R}}_{4}\left[\delta r\left(z\right)\right]& =\frac{1}{256}\left(35\delta r{\left(z\right)}^{4}-70\delta r{\left(z\right)}^{2}\delta {r}^{\prime }\left(z\right)+21\delta {r}^{{\prime\prime}}{\left(z\right)}^{2}+28\delta {r}^{\prime }\left(z\right)\delta {r}^{\left(3\right)}\left(z\right)\right.\hfill \\ \hfill & \left.\quad -70\delta r\left(z\right)\delta {r}^{\prime }{\left(z\right)}^{2}+14\delta {r}^{\left(4\right)}\left(z\right)-\delta {r}^{\left(6\right)}\left(z\right)\right).\hfill \end{align} \tag{ B.12 }$$

The cases m = 3 and m = 4 then clearly reproduce (4.46) and (B.6).

From (3.49) and (B.2) we observe that for multicritical models V(M) is a polynomial in M with alternating coeffficients. This holds true for general multicritical models $V\left(M\right)={\gamma }^{-1}{\sum }_{n=1}^{m}{v}_{n}{M}^{2n}$ with [45]

$$\begin{equation}{v}_{n}=\frac{{\left(-1\right)}^{n+1}}{4m}\left(\genfrac{}{}{0.0pt}{}{m}{n}\right)B\left(n,1/2\right),\quad n{\leqslant}m,\end{equation} \tag{ B.13 }$$

where B(x, y) is the Euler beta function. For m odd the multicritical V(M) has no real maximum whereas for m even we always find at least one real maximum. For the mth multicritical model with m odd the coefficients are Borel summable, while for m even the coefficients are not Borel summable [63]. The density of eigenvalues for the mth multicritical model with eigenvalues distributed in [−a, a] is given by [45]

$$\begin{equation}{\rho }_{\text{ext}}\left(\lambda \right)=\frac{{\left(1-{a}^{2}\right)}^{m}{\sqrt{{a}^{2}-{\lambda }^{2}}\enspace }_{2}{F}_{1}\left(1,m+\frac{1}{2},\frac{3}{2},\frac{{a}^{2}-{\lambda }^{2}}{1-{\lambda }^{2}}\right)}{2\pi \gamma \left(1-{\lambda }^{2}\right)},\end{equation} \tag{ B.14 }$$

where

$$\begin{equation}\gamma =\frac{1}{4m}-\frac{1}{4m}{\left(1-{a}^{2}\right)}^{m},\end{equation} \tag{ B.15 }$$

and γ_c = 1/4m (figure 18).

Zoom In Zoom Out Reset image size