text:
# Transient tunneling effects of resonance doublets in triple barrier systems

## I Introduction

In this work we address the issue of time-dependent tunneling phenomena in triple barrier resonant systems, aimed to study the transient behavior of the probability density near a resonance doublet. We shall refer to these structures as two-level open systems, in the sense that their finite-width barriers enable the system to interact with incident particles via a tunneling process. The dynamical properties of triple barrier structures have not drawn the attention they deserve. In this work we wish to emphasize that triple barriers involve novel dynamical aspects not present in double-barrier structures, where the tunneling dynamics near resonance energy is governed by a single resonance.hauge ; gcr97 ; apl1 ; apl2

The purpose of this paper is to demonstrate, based on an exact analytical approach, that the dynamics of the transient probability density is governed by three relevant frequencies that involve the resonance energies of the doublet and the incidence energy $`E`$. We find that in addition to the Bohr frequency, $`\omega _{21}=|_2_1|/\mathrm{}`$, which is an intrinsic property of the system, there are two additional frequencies, $`\omega _1=|E_1|/\mathrm{}`$ and $`\omega _2=|E_2|/\mathrm{}`$, where the resonance energies $`_1`$ and $`_2`$ are the real parts of the corresponding complex resonance energies $`E_n=_ni\mathrm{\Gamma }_n/2`$ (n=1,2) of the problem. This should be contrasted with the well known cohen  dynamical behavior of closed two-level systems, which is governed only by the Bohr frequency, $`\mathrm{\Omega }_{12}=(\mathrm{E}_2\mathrm{E}_1)/\mathrm{}`$ where $`\mathrm{E}_1`$ and $`\mathrm{E}_2`$ are the real energy eigenvalues of the system. As shown below, the above frequencies may be manipulated to produce a significant enhancement of the short-time transient behavior of the probability density.

The paper is organized as follows: Section II presents an overview of the formalism. In section III we discuss the transient behavior of the probability density through several numerical examples. Finally, section IV, provides some concluding remarks.

## II The Formalism

The model used in this work deals with an explicit solution to the time-dependent Schrödinger equation with cutoff initial conditions, gcr97  and consists of a generalization to tunneling problems of the free quantum shutter setup that predicts diffraction in time.mm ; kleber ; holland  The phenomenon of diffraction of matter in time has been recently experimentally verifieddalibard ; hils  and has also stimulated further studies.zeilinger  The setup used in this work may be visualized as a quantum shutterci  placed at $`x=0`$, just at the left edge of the resonant structure that extends over the interval $`0xL`$. Upon opening the shutteropening  at $`t=0`$, the incoming initial wave, represented by a cutoff plane wave,

$$\mathrm{\Psi }(x,k;t=0)=\{\begin{array}{cc}e^{ikx}e^{ikx},& x0\\ 0,& x>0,\end{array}$$
(1)
interacts with the internal region ($`0xL`$) of the potential. The wave solution to the time-dependent problem $`\mathrm{\Psi }(x,k;t)`$ for $`x>0`$ and $`t>0`$, is given by,gcr97

$$\mathrm{\Psi }=\mathrm{\Phi }_kM(y_k)\mathrm{\Phi }_kM(y_k)\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\rho _nM(y_{k_n}).$$
(2)
The quantities $`\mathrm{\Phi }_{\pm k}\mathrm{\Phi }(x,\pm k)`$ refer to the stationary wave solution, and the factors,

$$\rho _n(x,k)2iku_n(0)u_n(x)/(k^2k_n^2),$$
(3)
are given in terms of the resonant states $`\{u_n(x)\}`$ and the complex energy eigenvalues $`E_n=\mathrm{}^2k_n^2/2m`$ of the problem. The complex energy eigenvalues may be written in terms of the complex wave numbers $`k_n=a_nib_n`$, and correspond to the S-matrix poles of the problem. They are distributed in the third and fourth quadrants on the complex $`k`$-plane in a well known manner. The $`M`$-functions are defined as,gcr97

$$M(y_s)=\frac{1}{2}w(iy_s),$$
(4)
where the functions $`w(iy_s)`$ stand for the complex error function,wiz

$$w(iy_s)=\mathrm{e}^{y_s^2}\mathrm{erfc}(y_s)$$
(5)
and the argument $`y_s`$, reads,

$$y_s\mathrm{e}^{i3\pi /4}\left(\frac{\mathrm{}}{2m}\right)^{1/2}st^{1/2},$$
(6)
and $`s`$ stands for $`\pm k`$ or $`k_{\pm n}`$. As shown elsewhere,gcr97  the time-dependent solution given by Eq. (2) goes into the stationary solution $`\mathrm{\Phi }_k`$ at asymptotically long times.

For triple barrier systems, the resonance spectra typically corresponds to a succession of resonance doublets, formed by the coupling of the single resonances associated with each of the two wells of the system. We shall be interested in systems where the first doublet is isolated. The approximation to Eq. (2) then reads,

$`\mathrm{\Psi }`$ $``$ $`\mathrm{\Phi }_kM(y_k)\mathrm{\Phi }_k^{}M(y_k)`$ (7)
$`{\displaystyle \underset{n=1}{\overset{2}{}}}\left\{\rho _nM(y_{k_n})+\rho _nM(y_{k_n})\right\},`$
where we have used $`\mathrm{\Phi }_k=\mathrm{\Phi }_k^{}`$. For a resonance doublet the stationary function may also be written as the sum over the first two resonance terms,overlapping  namely,

$$\mathrm{\Phi }_k(x,k)\rho _1(x,k)+\rho _2(x,k),$$
(8)
and consequently,

$$|\mathrm{\Phi }_k|^2|\rho _1|^2+|\rho _2|^2+\rho _{12},$$
(9)
where $`\rho _{12}=2\mathrm{R}\mathrm{e}\{\rho _1\rho _2^{}\}`$. Although the time-dependence of Eq. (7) is contained in the $`M`$-functions, a considerable simplification of this two-level formula can be derived, in which the time dependence is explicitly given in terms of simple functions. Such a derivation is discussed in detail elsewhereyo  and we will recount it here briefly. The $`M`$-functions $`M(y_k)`$ and $`M(y_{k_n})`$ contained in Eq. (7), can be related to functions of the form $`M(y_k)`$ and $`M(y_{k_n})`$ by means of the symmetry relation,gcr97

$$M(y_s)=e^{y_s^2}M(y_s).$$
(10)
Using Eq. (10) we can rewrite Eq. (7) as,

$$\mathrm{\Psi }=\underset{n=1}{\overset{2}{}}\rho _n(x)\left[\mathrm{e}^{y_k^2}\mathrm{e}^{y_{k_n}^2}\right]+\mathrm{\Delta }(x,t),$$
(11)
where $`\mathrm{\Delta }(x,t)`$ accounts for all the terms containing $`M`$ functions of the form $`M(y_k)`$ and $`M(y_{k_n})`$, which behave as an inverse power of $`t`$, as follows from its series expansion,gcr97  $`M(y_s)1/2[1/(\pi ^{1/2}y_s)1/(\pi ^{1/2}y_s^3)+\mathrm{}]`$. Thus, except for extremely short or very long times compared with the lifetimes of the resonance levels of the doublet, the term $`\mathrm{\Delta }(x,t)`$ gives a negligible contribution to the solution and can be neglected. By doing this, we can obtain a simple expression for the probability density, valid for the internal region and an energy $`E`$ close to the doublet, namely,

$$\left|\mathrm{\Psi }(E,t)\right|^2=\varphi _1(E,t)+\varphi _2(E,t)+\varphi _{12}(E,t),$$
(12)
where $`\varphi _n(E,t)`$ and the interference terms $`\varphi _{mn}(E,t)`$ (n=1,2) are given respectively by,

$$\varphi _n(E,t)=\left|\rho _n\right|^2\chi _n(E,t),$$
(13)
and

$$\varphi _{mn}(E,t)=2Re\{\rho _m\rho _n^{}\xi _{mn}(E,t)\},$$
(14)
where $`E=\mathrm{}^2k^2/2m`$ is the incidence energy, and the functions $`\chi _n`$ and $`\xi _{mn}`$ have the following closed analytic expressions,

$$\chi _n(E,t)=12\mathrm{cos}(\widehat{\omega }_nt)e^{\mathrm{\Gamma }_nt/2\mathrm{}}+e^{\mathrm{\Gamma }_nt/\mathrm{}};$$
(15)
$`\xi _{mn}(E,t)`$ $`=`$ $`[1e^{i\widehat{\omega }_mt\mathrm{\Gamma }_mt/2\mathrm{}}e^{i\widehat{\omega }_nt\mathrm{\Gamma }_nt/2\mathrm{}}`$ (16)
$`+e^{i\widehat{\omega }_{mn}t(\mathrm{\Gamma }_m+\mathrm{\Gamma }_n)t/2\mathrm{}}].`$
In the above expressions, $`\widehat{\omega }_1`$, $`\widehat{\omega }_2`$ and $`\widehat{\omega }_{12}`$ are defined by $`\widehat{\omega }_n(E_n)/\mathrm{}`$ and $`\widehat{\omega }_{21}(_2_1)/\mathrm{}`$.

The formula given by Eq. (12) is an important analytical result since it explicitly reveals novel aspects of the quantum dynamics of tunneling structures with resonance doublets. According to Eq. (12), the time-dependent probability density is the superposition of the three oscillating contributions, $`\varphi _1(E,t)`$, $`\varphi _2(E,t)`$ and $`\varphi _{21}(E,t)`$, which have in general different amplitudes and frequencies. The three characteristic frequencies that govern the time evolution during the transient regime are: $`\omega _1`$, $`\omega _2`$ and $`\omega _{21}`$, which are respectively the absolute values of $`\widehat{\omega }_1`$, $`\widehat{\omega }_2`$ and $`\widehat{\omega }_{21}`$. Note that at asymptotically long times, it is easily seen from Eqs. (15) and (16), respectively, that $`\chi _n1`$ and $`\xi _{mn}1`$, and hence the probability density for the two-level formula, given by Eq. (12), goes into the stationary solution given by Eq. (9).

## III examples

We shall be interested in analyzing the transient tunneling effects of the probability density at the right-hand edge of the system, $`x=L`$, because there appear the largest transient effects along the transmitted region. We consider as a first example a periodic triple barrier system with parameters given as in Ref. strong, , namely: barrier heights $`V_0=0.12`$ eV, barrier widths $`b_0=3.0`$ nm, well widths $`w_0=16.0`$ nm; and effective mass of the electron $`m=0.067m_e`$. The corresponding resonance parameters of the first doublet are: energy positions, $`_1=11.512`$ meV and $`_2=14.387`$ meV; and resonance widths, $`\mathrm{\Gamma }_1=0.4089`$ $`meV`$ and $`\mathrm{\Gamma }_2=0.6365`$ meV.

Let us first illustrate the reliability of our approximate formula derived above for a single doublet. In Fig. 1 we compare the behavior of the probability density using both the formal solution, Eq. (2) (solid line), and Eq. (12) (dashed line), for an incidence energy near the first resonance, $`E=_1+2.0\mathrm{\Gamma }_1=12.33`$ meV. As can be appreciated, the two level approximation to Eq. (2) given by Eq. (12) gives an excellent description of the probability density. In this particular example, we have included in Eq. (2) the first four resonances of the systems, i.e. $`N=4`$, in order to illustrate that the contribution of far away resonances is negligible. The irregular behavior of $`|\mathrm{\Psi }|^2`$ observed in Fig. 1, arises from the interplay between $`\varphi _1`$, $`\varphi _2`$, and $`\varphi _{12}`$ of Eq. (12). This situation contrasts with the regular behavior observed in double barrier structures. As shown in a recent work,gcr97  in the case of a double barrier system the probability density grows exponentially for incidence at resonance, and exhibits regular oscillations with a single frequency if the incidence occurs near resonance.apl1  This is due to the fact that in the double barrier case, the one-level approximation stands, and hence only the term $`\varphi _1(E,t)`$ is important. An interesting feature of triple barrier systems, not present in double barrier structures, is that the frequencies can be tuned by a proper choice of the incidence energy $`E`$. This allows to manipulate the frequencies in such a way that the irregularities observed in Fig. 1 dissappear. This occurs at two special situations that depend on $`E`$. The first situation is when the incidence energy coincides with one of the two resonances, and the second one occurs when the incidence energy coincides with the middle point between the two resonances of the doublet. In the first case, only one of the three terms of Eq. (12) dominates over the remaining two, for example if $`E=_1`$, then $`\varphi _2(E,t)`$ and $`\varphi _{21}(E,t)`$ are negligible in comparison with $`\varphi _1(E,t)`$. Since in this case $`\widehat{\omega }_1=0`$, then we have that the probability density is governed by the following simple expression,

$$|\mathrm{\Psi }(E=_1)|^2T(_1)(1e^{t/\tau _1})^2,$$
(17)
where $`\tau _1=\mathrm{}/\mathrm{\Gamma }`$ is the lifetime of the resonance $`n=1`$ and $`T(_1)`$ is the peak value of the transmission coefficient, which is unity for this symmetrical system. The results of this resonant case are depicted in Fig. 2(a), where we show the calculations using Equations (7) (solid line) and (17) (dashed line). Both curves almost coincide, except by the very small oscillations of the exact two-level formula, i.e. Eq. (7) (solid line), which are due to the effect of the second resonance of the doublet. In the second case, when the incidence energy is chosen just at the middle of the two resonances i.e., $`E=\overline{E}(_1+_2)/2`$, we have $`\overline{\omega }\widehat{\omega }_1=\widehat{\omega }_2=\widehat{\omega }_{21}/2`$, that is, the dynamics is governed by a single frequency, $`\overline{\omega }`$, and the behavior of $`\left|\mathrm{\Psi }(L,k;t)\right|^2`$ vs $`t`$ is similar to a diffraction in time pattern,mm  see Fig. 2(b). Here the numerical value of $`\overline{E}=12.949`$ meV$`=_1+3.515\mathrm{\Gamma }_1`$.

Compared to the double barrier case, this constructive interference effect produces an important enhancement of the transient probability density. The comparison is shown in Fig. 3(a) in which we used the same triple barrier structure parameters of the previous figures, and the double barrier system with potential parameters: barrier heights $`V_0=0.23`$ eV, barrier widths $`b_0=5.0`$ nm, well width $`w_0=5.0`$ nm. The first resonant state of the system has energy position, $`_1=80.11`$ meV; and resonance width, $`\mathrm{\Gamma }_1=1.033`$ meV. The incidence energy was also chosen with the same deviation from resonance in units of the resonance width, that is, $`E=_1+3.515\mathrm{\Gamma }_1`$ whose numerical value is $`83.740`$ meV. Note that the scale in the time axis was normalized to the lifetime of the first resonance of each system, which for the triple barrier has the value $`\tau _1=1.61`$ ps, and for the double barrier, $`\tau _1=6.37`$ ps. Both curves tends to their correct asymptotic limits, the transmission coefficient, which for the double barrier system has the value T(E)=0.0229, and for the triple barrier system, T(E)=0.119. These values are indicated by the small arrows in Fig. 3(a).

As it is well known from time-independent studies in triple barrier systems,strong ; yama  for incidence energies at the center of an isolated doublet, the transmission coefficient increases to unity as we increase the width $`b_2`$ of the central barrier to about twice the value of the width of the external barriers. In view that the transmission coefficient is the asymptotic value of the time dependent probability density at $`x=L`$, it is expected that the latter can also be enhanced in the same fashion. In Fig. 3(b) we illustrate how we can manipulate this extra degree of freedom to enhance the amplitude of the oscillations of the transient probability density. Here we considered $`b_2=4.0`$ nm (dashed line), and $`b_2=5.0`$ nm (dotted line); the solid line corresponds to the same triple barrier system of Fig. 3(a) ($`b_2=3.0`$ nm), which is also included here for comparison. Also in this figure, the values of the transmission coefficient at the energies $`\overline{E}`$ for each system is indicated by arrows to illustrate how the time dependent probability density tends to the correct asymptotic behavior as $`t\mathrm{}`$.

## IV concluding remarks

In conclusion, the dynamics of the probability density for triple barrier resonant structures, which is a typical example of an open two-level system, has been explored. We have derived a simple analytic expression for the probability density that provides an accurate description for energies near the resonance doublet of the system. The two-level formula allows to identify three relevant frequencies that govern the transient behavior as a function of time. The derived formula goes into the stationary two-level solution at asymptotically long times, and thus establishes a link with the usual stationary approach. Our results suggest that the transient effects that we have discussed are of relevance at short times and distances from the interaction region. Hopefully our results may stimulate experimental work on this subject.

###### Acknowledgements.

The authors acknowledge financial support from Conacyt, México, through Contract No. 431100-5-32082E. One of the authors (G.G.C) acknowledges financial support of DGAPA-UNAM under Grant No. IN101301.


format: text

text:
# 1 Introduction

## 1 Introduction

Random matrix models  have a wide range of applications in mathematics and physics where they constitute a major field of activity. They are involved in condensed matter physics (quantum chaos , localization, crystal growths ,…etc), statistical physics  (on a 2d fluctuating surface, also called 2d euclidean quantum gravity, linked to conformal field theory), high energy physics (string theory , quantum gravity , QCD ,…), and they are very important in mathematics too: (they seem to be linked to the Riemann conjecture ), they are important in combinatorics, and provide a wide class of integrable systems .

In the 80’s, random matrix models were introduced as a toy model of 0-dimensional string theory and quantum gravity . By fine-tuning the parameters of the model, one can reach many multi-critical points, which are in relationship with the $`(p,q)`$ minimal conformal models . The 1-hermitian matrix models is limited to $`(p,2)`$ conformal minimal models, whereas the 2-matrix model can represent any $`(p,q)`$ conformal minimal model . The 2-matrix model is thus more general than the 1-matrix model, and it has been the source of considerable interest in the past few months .

The free energy of the 1-matrix model was conjectured  and also rigorously proven  to have a $`1/N^2`$ expansion called topological expansion ($`N`$ is the size of the matrix):

$$F=\underset{h=0}{\overset{\mathrm{}}{}}N^{2h}F^{(h)}$$
(1.1)
The authors of  invented an efficient method to compute recursively all the $`F^{(h)}`$’s, and they improved it in .

The 2-matrix model is conjectured to have a similar $`1/N^2`$ topological expansion. Indeed, this expansion is the main motivation for applications to 2-dimensional quantum gravity , because each $`F^{(h)}`$ is the partition function of a statistical physics model on a genus $`h`$ surface. There is at the present time no rigorous proof222The Riemann-Hilbert approach seems to be the best way to prove the existence of the $`1/N^2`$ topological extension as in . The Riemann-Hilbert problem for the 2-matrix model has been formulated , and seems to be on the verge of being solved . of the existence of such an expansion, but, assuming that it exists, the aim of the present work is to give a method to compute recursively the terms of the expansion, similar to that of .

The 2-matrix model was first introduced as a model for two-dimensional gravity, with matter, and in particular with an Ising field . The diagrammatic expansion of the 2-matrix model’s partition function is known to generate 2-dimensional statistical physics models on a random discrete surface :

$$N^2F=\mathrm{ln}Z=\underset{\mathrm{surfaces}}{}\underset{\mathrm{matter}}{}\mathrm{e}^{\mathrm{Action}}$$
(1.2)
where the Action is the matter action (like Ising’s nearest neighboor spin coupling) plus the gravity action (total curvature and cosmological constant) . The cosmological constant couples to the area of the surface, and $`N`$ (the size of the matrix) couples to the total curvature, i.e. the genus of the surface. The large $`N`$ expansion thus generates a genus expansion:

$$F=\underset{h=0}{\overset{\mathrm{}}{}}N^{2h}F^{(h)}$$
(1.3)
where $`F^{(h)}`$ is the partitrion function of the statistical physics model on a random surface of fixed genus $`h`$.

$$F^{(h)}=\underset{\mathrm{genus}h\mathrm{surfaces}}{}\underset{\mathrm{matter}}{}\mathrm{e}^{\mathrm{Action}}$$
(1.4)
The leading term $`F^{(0)}`$ is the planar contribution. It can be computed by many different methods, for instance the saddle point method (along the method invented by  and rigorously established by ) or the loop equation method which we will explain below. Our goal in this article is to compute $`F^{(1)}`$ and present an algorithmic method for computing $`F^{(h)}`$ for $`h1`$. We generalize the method of .

### 1.1 Reminder: the one-matrix model

For a given polynomial $`V(x)`$ of degree $`d+1`$, called the potential, and a given integer $`N`$, we define the partition function $`Z`$ and the free energy $`F`$ as:

$$Z=\mathrm{e}^{N^2F}=\text{d}M\mathrm{e}^{N\mathrm{tr}V(M)}$$
(1.5)
where the integral is over the set of hermitian matrices of size $`N`$, with the measure $`\text{d}M`$ equal to the product of Lebesgue measures of all real components of $`M`$.

The free energy has a $`1/N^2`$ power series expansion (under the 1-cut asumption, this will be discussed in more details in section 4.1):

$$F=\underset{h=0}{\overset{\mathrm{}}{}}N^{2h}F^{(h)}$$
(1.6)
In their pioneering work, the authors of  invented a method to compute recursivelly all the $`F^{(h)}`$’s. In particular they found that the genus one free energy for the 1-matrix model is:

$$F^{(1)}=\frac{1}{24}\mathrm{ln}\left((ba)^4M(a)M(b)\right)$$
(1.7)
where $`b,a`$ and the polynomial $`M(x)`$ are related to the large $`N`$ limit of the average density of eigenvalues of the matrix:

$$\rho (x)=\frac{1}{2\pi }M(x)\sqrt{(xa)(bx)},M(x)=\underset{x\mathrm{}}{Pol}V^{}(x)/\sqrt{(xa)(xb)}$$
(1.8)
This density has a compact support $`[a,b]`$, and $`M(x)`$ is a polynomial of degree $`d1`$ (the 1-matrix case is re-explained in section 7).

We are going to extend this kind of expression (eq. (1.7)) for the 2-matrix model.

### 1.2 Outline cof the article

* In section 2 we introduce the definitions and notations, in particular we define the 1-loop functions and 2-loop functions, and the loop-insertion operators.
* In section 3 we explain the loop equation method, and derive the “Master loop equation”.
* In section 4, we observe that, to leading order, the master loop equation is an algebraic equation of genus zero, and we study the geometry and the sheet-structure of the underlying algebraic curve.
* In section 5 we use the algebraic master loop equation to compute all the previously defined loop functions to leading order.
* In section 6 we include the previously neglected $`1/N^2`$ term in the loop equation, and we compute the 1-loop function $`Y(x)`$ to next to leading order. Then we derive the next to leading order free energy $`F^{(1)}`$ by integrating $`Y^{(1)}(x)`$. We also discuss how to compute higher order terms.
* In section 7 we compare with the results previously known for the one matrix case.
* section 8 is the conclusion.

## 2 The 2-matrix model

### 2.1 The partition function

For two given polynomials $`V_1(x)`$ of degree $`d_1+1`$ and $`V_2(y)`$ of degree $`d_2+1`$, called the potentials, and a given integer $`N`$, we define the partition function $`Z`$ and the free energy $`F`$ as:

$$Z=\mathrm{e}^{N^2F}=\text{d}M_1\text{d}M_2\mathrm{e}^{N\mathrm{tr}[V_1(M_1)+V_2(M_2)M_1M_2]}$$
(2.1)
where the integral is over the set of pairs of hermitian matrices $`M_1`$ and $`M_2`$ of size $`N`$, with the measure $`\text{d}M_1\text{d}M_2`$ equal to the product of Lebesgue measures of all real components of $`M_1`$ and $`M_2`$ (see  for definitions).

The potentials $`V_1`$ and $`V_2`$ are written:

$$V_1(x)=\underset{k=0}{\overset{d_1}{}}\frac{g_{k+1}}{k+1}x^{k+1},V_2(y)=\underset{k=0}{\overset{d_2}{}}\frac{g_{k+1}^{}}{k+1}y^{k+1}$$
(2.2)
Remark: eq. (2.1) is well defined only when the potentials are bounded from below, i.e. when:

$$\left\{d_1\mathrm{and}d_2\mathrm{are}\mathrm{odd}g_{d_1+1}>0,g_{d_2+1}^{}>0\right\}$$
(2.3)
However, it is known  how to give a meaning to 2.1 for potentials non bounded from below, in the large $`N`$ limit (somehow the tunnel effect allowing the matrices to “escape from the potential wells” is suppressed like $`O(\mathrm{e}^N)`$ in the large $`N`$ limit, it plays no role in the $`1/N^2`$ expansion). Therefore, we will consider arbitrary polynomial potentials which do not necessarily satisfy condition (2.3). For simplicity, we will assume that the potentials $`V_1`$ and $`V_2`$ are real333The derivation of the loop equations in section 3 should be slightly modified in the case of non-real potentials, because some change of variables we consider would be no longer hermitian. But most of the results derived in this article are still valid for complex potentials..

In addition, we will assume that the potentials $`V_1`$ and $`V_2`$ are generic, i.e. they are not critical (we will explain what is a critical potential in section 4.4.1, see ).

The partition function $`F`$ can be interpreted (through the Feynmann graphs expansion, see ) as the generating function of surfaces made of polygons with a spin: spin “up” (+) for the polygons generated by the vertices of $`M_1`$, “down” (-) for the polygons generated by the vertices of $`M_2`$.

$$F=\underset{S}{}N^{2h(S)}\frac{1}{\mathrm{\#}\mathrm{Aut}(S)}\frac{1}{g_2g_2^{}}\underset{k=1,k2}{\overset{d_1+1}{}}(g_k)^{n_k(S)}\underset{j=1,j2}{\overset{d_2+1}{}}(g_{}^{}{}_{j}{}^{})^{\stackrel{~}{n}_j(S)}$$
(2.4)
where $``$ is the ensemble of surfaces made of up-polygons with at most $`d_1+1`$ sides and down-polygons with at most $`d_2+1`$ sides. Aut$`(S)`$ is the group of automorphisms of the graph of edges of $`S`$, $`h(S)`$ is the genus of $`S`$, $`n_k(S)`$ ($`1kd_1+1`$) is the number of up $`k`$-gons of $`S`$ and $`\stackrel{~}{n}_j(S)`$ ($`1jd_2+1`$) is the number of down $`j`$-gons of $`S`$.

We write digrammatically:

$$F=\text{}+\frac{1}{N^2}\text{}+\mathrm{}$$

This is why the two-matrix model can represent the Ising model on a random discretized surface , and other 2d statistical physics models on a random surface carrying some type of matter . It is clear from eq. (2.4) that relevant applications to 2d statistical physics on a random surface involve potentials which violate condition 2.3 (all the $`g_k`$’s and $`\stackrel{~}{g}_l`$’s must be negative in order to have positive weights).

### 2.2 Definition of the loop functions

We introduce the following expectation values computed with the probability weight $`\frac{1}{Z}\mathrm{e}^{N\mathrm{tr}[V_1(M_1)+V_2(M_2)M_1M_2]}`$:

$$T_{k,l}=\frac{1}{N}\mathrm{tr}M_1^kM_2^l$$
(2.5)
We will be particularly interested in their large $`N`$ limits and large $`N`$ expansion.

$`T_{k,l}`$ can be interpreted as the partition function of a statistical physics model on surfaces with one boundary of lenght $`k+l`$ ($`k`$ contiguous + and $`j`$ contiguous -):

$$T_{k,l}=\text{}+\frac{1}{N^2}\text{}+\mathrm{}$$

It is more convenient to introduce generating functions which contain all the $`T_{kl}`$’s at once:

#### 2.2.1 One-loop functions

We define the following formal functions of the complex variables $`x`$ and $`y`$:

$$W(x)=\frac{1}{N}\mathrm{tr}\frac{1}{xM_1}=\underset{k=0}{\overset{\mathrm{}}{}}\frac{T_{k,0}}{x^{k+1}},Y(x)=V_1^{}(x)W(x)$$
(2.6)
$$\stackrel{~}{W}(y)=\frac{1}{N}\mathrm{tr}\frac{1}{yM_2}=\underset{k=0}{\overset{\mathrm{}}{}}\frac{T_{0,k}}{y^{k+1}},X(y)=V_2^{}(y)\stackrel{~}{W}(y)$$
(2.7)
These functions are formally defined only through their large $`x`$ and large $`y`$ expansions, they are only a convenient rewritting of the collection of $`T_{k,0}`$ and $`T_{0,k}`$. The sum is not necessarily convergent.

We will find below, that $`W^{(0)}(x)=\mathrm{lim}_N\mathrm{}W(x)`$ (resp. $`\stackrel{~}{W}(y)`$) satisfies an algebraic equation, and thus $`W^{(0)}(x)`$ (resp. $`\stackrel{~}{W}^{(0)}(y)`$) is an analytical function of $`x`$ (resp. $`y`$), in the complex plane, with a cut. The radius of convergence of eq. (2.6) is finite in the large $`N`$ limit. Their analytical structure will be dicussed in details in part 4.4.

The functions $`W(x)`$ and $`\stackrel{~}{W}(y)`$ are called the resolvents, they provide most of the information about the statistical properties of the spectra of $`M_1`$ and $`M_2`$. For instance the location of the cut of $`W^{(0)}(x)`$, is the support of the large $`N`$ average density of eigenvalues of $`M_1`$, and in the large $`N`$ limit, the average density of eigenvalues is proportional to the discontinuity of $`W^{(0)}(x)`$ along its cuts.

Their diagrammatic expansion is the generating function for statistical physics models coupled to gravity, on a surface with one boundary (a disc if genus zero), this is why they are called one-loop functions.

$$W(x)=\underset{l_1=0}{\overset{\mathrm{}}{}}x^{l_11}\underset{\genfrac{}{}{0pt}{}{\text{surfaces, }}{\text{boundary }=l_1\text{ up}}}{}\underset{\mathrm{matter}}{}\mathrm{e}^{\mathrm{Action}}=\text{}+\mathrm{}$$
(2.8)
where the sum carries on all discretized surfaces with one boundary made of $`l_1`$ spins +.

$$\stackrel{~}{W}(y)=\underset{l_2=0}{\overset{\mathrm{}}{}}x^{l_21}\underset{\genfrac{}{}{0pt}{}{\text{surfaces, }}{\text{boundary }=l_2\text{ down}}}{}\underset{\mathrm{matter}}{}\mathrm{e}^{\mathrm{Action}}=\text{}+\mathrm{}$$
(2.9)
where the sum carries on all discretized surfaces with one boundary made of $`l_2`$ spin -.

We also define:

$$W(x,y)=\frac{1}{N}\mathrm{tr}\frac{1}{xM_1}\frac{1}{yM_2}=\underset{k=0}{\overset{\mathrm{}}{}}\underset{l=0}{\overset{\mathrm{}}{}}\frac{T_{k,l}}{x^{k+1}y^{l+1}}$$
(2.10)
This function provides information on the correlation between the spectrum of $`M_1`$ and that of $`M_2`$. Its diagrammatic expansion generates surfaces with one “bi-colored” boundary, i.e. a boundary of lenght $`l_1+l_2`$ made of $`l_1`$ consecutive spins up, followed by $`l_2`$ down:

$$W(x,y)=\underset{l_1=0}{\overset{\mathrm{}}{}}\underset{l_2=0}{\overset{\mathrm{}}{}}x^{l_11}y^{l_21}\underset{\genfrac{}{}{0pt}{}{\text{surfaces, }}{\text{bnd }=l_1\text{ up + }l_2\text{ down}}}{}\underset{\mathrm{matter}}{}\mathrm{e}^{\mathrm{Action}}=\text{}$$
(2.11)
We also define:

$$U(x,y)=\frac{1}{N}\mathrm{tr}\frac{1}{xM_1}\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}$$
(2.12)
$$\stackrel{~}{U}(x,y)=\frac{1}{N}\mathrm{tr}\frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}\frac{1}{yM_2}$$
(2.13)
which are polynomials in one of the variable and analytical near $`\mathrm{}`$ for the other variable.

And we define

$$P(x,y)=\frac{1}{N}\mathrm{tr}\frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}$$
(2.14)
which is a polynomial in both variables, of degree $`(d_11,d_21)`$.

In appendix C we compute another one-loop function $`H(x,y,x^{},y^{})`$, which generates surfaces with one quadri-colored boundary, i.e. surfaces whose boundary is +-+-. Somehow it is the random surface analogus of a square with opposite sides of the same spin and adjacent sides of opposite spin. It maybe very usefull to compute that kind of function to study how spin percolates from one side to the other. It may have applications to compute boundary operators in 2d statistical physics on a random surface.

#### 2.2.2 Two-loop functions

We also need to define the following correlation functions, which contain the product of two traces. Their diagrammatic expansion is the generating function for statistical physics models coupled to gravity, on a surface with two boundaries (cylinder-like), this is why they are called two-loop functions.

$$\mathrm{\Omega }(x;x^{})=\frac{W(x)}{V_1(x^{})}=\mathrm{tr}\frac{1}{xM_1}\mathrm{tr}\frac{1}{x^{}M_1}N^2W(x)W(x^{}),$$
(2.15)
$$\stackrel{~}{\mathrm{\Omega }}(y;x^{})=\frac{\stackrel{~}{W}(y)}{V_1(x^{})}=\mathrm{tr}\frac{1}{x^{}M_1}\mathrm{tr}\frac{1}{yM_2}N^2W(x^{})\stackrel{~}{W}(y),$$
(2.16)
$$U(x,y;x^{})=\frac{U(x,y)}{V_1(x^{})}=\mathrm{tr}\frac{1}{xM_1}\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}\mathrm{tr}\frac{1}{x^{}M_1}N^2U(x,y)W(x^{}),$$
(2.17)
$$\stackrel{~}{U}(x,y;x^{})=\frac{\stackrel{~}{U}(x,y)}{V_1(x^{})}=\mathrm{tr}\frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}\frac{1}{yM_2}\mathrm{tr}\frac{1}{x^{}M_1}N^2\stackrel{~}{U}(x,y)W(x^{}),$$
(2.18)
$`P(x,y;x^{})={\displaystyle \frac{P(x,y)}{V_1(x^{})}}`$ $`=`$ $`\mathrm{tr}{\displaystyle \frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}}{\displaystyle \frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}}\mathrm{tr}{\displaystyle \frac{1}{x^{}M_1}}`$ (2.20)
$`N^2P(x,y)W(x^{}),`$
where the “loop-insertion” operators (see ) are formally444By formally, we mean they make sense only through their large $`x`$ (resp. $`y`$) expansion. If $`k`$ is larger than $`d_1`$, $`/g_k`$ means: take a potential of degree $`k`$, and then take the derivative w.r.t. $`g_k`$ at $`g_j=0`$ for $`j>d_1`$. defined as:

$$\frac{}{V_1(x^{})}=\underset{k=1}{\overset{\mathrm{}}{}}\frac{k}{x^k}\frac{}{g_k},\frac{}{V_2(y^{})}=\underset{k=1}{\overset{\mathrm{}}{}}\frac{k}{y^k}\frac{}{g_k^{}}$$
(2.21)
their action on the partition function or on expectation values inserts a trace, i.e. a boundary (a “loop”) for the ensemble of surfaces. For instance we have:

$$W(x)=\frac{}{V_1(x)}F,\mathrm{\Omega }(x;x^{})=\frac{}{V_1(x^{})}W(x)=\frac{}{V_1(x)}\frac{}{V_1(x^{})}F=\mathrm{\Omega }(x^{};x)$$
(2.22)
The diagrammatic expansion of $`F`$ gives surfaces with no boundary, $`W(x)`$ gives surfaces with one boundary, $`\mathrm{\Omega }(x,x^{})`$ gives surfaces with two (spin up) boundaries, and $`\stackrel{~}{\mathrm{\Omega }}(x,y)`$ gives surfaces with one spin up boundary and one spin down boundary. And so on: $`\frac{}{V_1(x^{})}`$ inserts a + loop, $`\frac{}{V_2(y^{})}`$ inserts a - loop.

$$\mathrm{\Omega }(x;x^{})=\text{},\stackrel{~}{\mathrm{\Omega }}(y;x^{})=\text{}.$$

## 3 The loop equations

### 3.1 Generalities

The loop equations are a consequence of the reparametrization invariance of the matrix integral. Basically, one writes that under an infinitesimal change of variable of the type:

$$M_1\stackrel{~}{M}_1=M_1+ϵf(M_1,M_2)$$
(3.1)
the partition function is left unchanged to order 1 in $`ϵ`$. We have to compute the variation of the measure and the variation of the potential:

$`Z`$ $`=`$ $`{\displaystyle \text{d}\stackrel{~}{M}_1\text{d}M_2\mathrm{e}^{N\mathrm{tr}[V_1(\stackrel{~}{M}_1)+V_2(M_2)\stackrel{~}{M}_1M_2]}}`$ (3.2)
$`=`$ $`{\displaystyle \text{d}M_1\text{d}M_2(1+ϵJ(M_1,M_2))\mathrm{e}^{N\mathrm{tr}[V_1(M_1)+V_2(M_2)M_1M_2]}(1ϵK(M_1,M_2))}`$ (3.3)
where $`J(M_1,M_2)`$ is the jacobian of the change of variable computed to order $`1`$ in $`ϵ`$:

$$det\frac{\stackrel{~}{M_1}}{M_1}=1+ϵJ(M_1,M_2)+O(ϵ^2)$$
(3.5)
and $`K(M_1,M_2)`$ is the variation of the action to order $`1`$ in $`ϵ`$. We must have:

$$J(M_1,M_2)=K(M_1,M_2)$$
(3.6)
We are going to give the receipe555Easy to prove directly order by order in $`x`$. We leave the proof to the reader. A similar method is used in . to compute $`J(M_1,M_2)`$ for various change of variables $`f(M_1,M_2)`$.

Remarks:

\- It is important to take into account the non-commutativity of matrices.

\- It is also important to choose a change of variable which is a hermitian matrix.

\- We are presenting the method for variations of $`M_1`$ but of course the same applies to $`M_2`$ and we will use both changes of variables in the body of this article.

### 3.2 split and merge rules

We will be interested in change of variables of the following form666As before, $`\frac{1}{xM_1}=_kx^{k1}M_1^k`$ is a “formal” notation for the collection of $`M_1^k`$, $`k=0\mathrm{}\mathrm{}`$. :

* $`f(M_1,M_2)=A\frac{1}{xM_1}B`$:

where the matrices $`A`$ and $`B`$ are functions of $`M_1`$ and $`M_2`$. The variation of the measure is:

$$J(M_1,M_2)=\mathrm{tr}\left(A\frac{1}{xM_1}\right)\mathrm{tr}\left(\frac{1}{xM_1}B\right)+\mathrm{contributions}\mathrm{from}A,\mathrm{and}B.$$
(3.7)
somehow, you split the trace into two traces whenever you meet a $`\frac{1}{xM_1}`$ term outside a trace, and each side receives a $`\frac{1}{xM_1}`$.
* $`f(M_1,M_2)=A\mathrm{tr}\left(B\frac{1}{xM_1}\right)`$:

where $`A`$, $`B`$ are functions of $`M_1`$ and $`M_2`$. The variation of the measure is:

$$J(M_1,M_2)=\mathrm{tr}\left(A\frac{1}{xM_1}B\frac{1}{xM_1}\right)+\mathrm{contributions}\mathrm{from}A,\mathrm{and}B.$$
(3.8)
somehow, you merge the traces whenever you meet a $`\frac{1}{xM_1}`$ term inside a trace, and duplicate the $`\frac{1}{xM_1}`$ term (one before $`B`$, one after $`B`$).
* The variations of $`A`$ and $`B`$ are computed recursively by the chain rule.

The variation of the action is simply:

$$\frac{1}{N}K(M_1,M_2)=\mathrm{tr}\left(V_1^{}(M_1)f(M_1,M_2)M_2f(M_1,M_2)\right)$$
(3.9)
and we will use the following trick very often:

$$\mathrm{tr}\left(V_1^{}(M_1)\frac{1}{xM_1}A\right)=V_1^{}(x)\mathrm{tr}\left(\frac{1}{xM_1}A\right)\mathrm{tr}\left(\frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}A\right)$$
(3.10)
where it is important to notice that the second term is polynomial in $`x`$.

### 3.3 Factorization “theorem”

Then, we will make the asumption that in the large $`N`$ limit, the average of product of traces decouple:

$$\frac{1}{N^2}\mathrm{tr}A\mathrm{tr}B\frac{1}{N}\mathrm{tr}A\frac{1}{N}\mathrm{tr}B+O(1/N^2)$$
(3.11)
more precisely that the connected correlation function

$$\mathrm{tr}A\mathrm{tr}B_\mathrm{c}=\mathrm{tr}A\mathrm{tr}B\mathrm{tr}A\mathrm{tr}B$$
(3.12)
has a large $`N`$ limit.

The factorization “theorem” is not a theorem, it was never proven for the 2-matrix-model777It was proven rigorously for the one-matrix model in ., but there are plenty of heuristic reasons to believe that it holds, and it was used by nearly all authors in this field. As we said in the introduction, we are not going to prove the existence of the $`1/N^2`$ expansion, but assuming it, we are going to give an algorithmic method to compute it.

### 3.4 The master loop equation

Consider the change of variable:

$$M_2M_2+ϵ\frac{1}{xM_1}$$

which is indeed a hermitian matrix. It leads to (there is no variation of the measure since $`\delta M_2`$ is independent of $`M_2`$, the RHS is the variation of the potential):

$$0=\mathrm{tr}\frac{1}{xM_1}V_2^{}(M_2)\frac{M_1}{xM_1}$$
(3.13)
i.e.

$$\frac{1}{N}\mathrm{tr}\frac{1}{xM_1}V_2^{}(M_2)=xW(x)1$$
(3.14)
Now, consider the change of variables (which is indeed a hermitian matrix):

$$M_1M_1+ϵ\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}\frac{1}{xM_1}+ϵ\frac{1}{xM_1}\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}$$

it leads to (the LHS is the variation of the measure, the RHS of the potential, and we use the cyclicity of the trace):

$`\mathrm{tr}{\displaystyle \frac{1}{xM_1}}\mathrm{tr}{\displaystyle \frac{1}{xM_1}}{\displaystyle \frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}}=`$ (3.15)
$`N\mathrm{tr}{\displaystyle \frac{V_1^{}(M_1)}{xM_1}}{\displaystyle \frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}}\mathrm{tr}{\displaystyle \frac{1}{xM_1}}M_2{\displaystyle \frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}}`$ (3.16)
i.e. (using eq. (3.14)):

$$(yY(x))U(x,y)=V_2^{}(y)W(x)P(x,y)xW(x)+1\frac{1}{N^2}U(x,y;x)$$
(3.17)
Notice that $`U(x,y)`$ is a polynomial in $`y`$, therefore it is finite for $`y=Y(x)`$. If we choose $`y=Y(x)`$ we get (using $`W(x)=V_1^{}(x)Y(x)`$):

$$\{\begin{array}{c}y=Y(x)\hfill \\ (V_2^{}(y)x)(V_1^{}(x)y)P(x,y)+1=\frac{1}{N^2}U(x,y;x)\hfill \end{array}$$
(3.18)
which we write:

$`E(x,Y(x))={\displaystyle \frac{1}{N^2}}U(x,Y(x);x)`$      (3.19)
where $`E(x,y)`$ is a polynomial in $`x`$ (degree $`d_1+1`$) and $`y`$ (degree $`d_2+1`$):

$`E(x,y)=(V_1^{}(x)y)(V_2^{}(y)x)P(x,y)+1`$      (3.20)
eq. (3.19) is called the “master loop equation” , we will see below that it allows to determine the function $`Y(x)`$ and all other functions in the problem.

### 3.5 $`1/N^2`$ expansion

We assume that the function $`Y(x)`$ can be written as a power series in $`1/N^2`$:

$$Y(x)=Y^{(0)}(x)+\frac{1}{N^2}Y^{(1)}(x)+\mathrm{}$$
(3.21)
Such an expansion is believed to hold in the “1-cut case”, which will be explained below (section 4.1).

* To leading order, the loop equation is an algebraic equation:

$$E^{(0)}(x,Y^{(0)}(x))=0$$
(3.22)
where

$$E^{(0)}(x,y)=(V_1^{}(x)y)(V_2^{}(y)x)P^{(0)}(x,y)+1$$
(3.23)
The polynomial $`P^{(0)}(x,y)`$ has degrees $`(d_11,d_21)`$, it is determined by the large $`x`$ and large $`y`$ behaviors of $`Y(x)`$ and $`X(y)`$, as well as some analytical requirement, namely the “one-cut asumption”. As we will see below, in the 2-matrix model, this asumption is replaced by the asumption that the algebraic curve $`E(x,y)=0`$ has genus zero. In the next section we will study this algebraic curve, and we will explain how to compute the polynomial $`P^{(0)}(x,y)`$.

Let us note that this algebraic equation was derived by amny authors and various methods in the literature: withe the loop equation method, the existence of an algebraic equation is proven (and computed for small degree potentials) in  , and was presented for arbitrary polynomial potentials in the appendix of  and . The authors of  and  derived it by the “saddle point method” inspired from , and from the 2-toda hierarchy dispersionless limit, and that algebraic equation was used by many other authors in particular cases. And notice that it is the same equation as the spectral curve of any of the 4 dual differential systems in .
* To order $`1/N^2`$ we expand:

$$E(x,y)=E^{(0)}(x,y)\frac{1}{N^2}P^{(1)}(x,y)+\mathrm{}$$
(3.24)
where $`P^{(1)}(x,y)`$ has degrees $`(d_11,d_21)`$. We will see how to compute it in section 6.

We easily get:

$$Y^{(1)}(x)=\frac{U^{(0)}(x,Y^{(0)}(x);x)+P^{(1)}(x,Y^{(0)}(x))}{_yE^{(0)}(x,Y^{(0)}(x))}$$
(3.25)
In this expression, all terms in the RHS are computed to leading order, except the polynomial $`P^{(1)}(x,y)`$. But $`P^{(1)}(x,y)`$ is completely determined (in terms of leading order functions only) by the condition that all the unwanted zeroes of the denominator (those which don’t correspond branch-points) cancel. In other words, once we know the functions $`Y(x)`$ and $`U(x,y;x^{})`$ to leading order, we can compute $`Y(x)`$ to order $`1`$. Similarly, we can compute every loop function to order $`1`$.
* We can repeat the procedure to find higher orders. Once we know all the loop functions up to order $`h`$, we can compute them to order $`h+1`$. Of course the expressions become more and more complicated with increasing orders, but the method is the same. We will illustrate our method by computing the order $`1`$ only.
* once we know $`Y^{(h)}(x)`$, we can find the order $`h`$ free energy $`F^{(h)}`$, by integration: $`Y^{(h)}(x)=F^{(h)}/V_1(x)`$.

## 4 Leading order: algebraic equation

While we study the leading order, we will drop the subscripts $`(0)`$ for lisibility.

### 4.1 Reminder: one-cut asumption for the 1-matrix model

In the 1-matrix model, the existence of the $`1/N^2`$ topological expansion is related to the 1-cut asumption. We refer the reader to  for more detailed explanations, and we just summarize the main idea.

When there is more than one cut, i.e. when the density of eigenvalues has a disconnected support, there are tunnel effects between the different cuts which generate an oscillating function of $`N`$ (a theta-function). Such a periodic function cannot have a $`1/N^2`$ expansion. The frequencies of oscillations are the filling fraction of eigenvalues in the cuts, in other words, the frequencies are equal to the integrals of $`W(x)`$ along the non-trivial cycles of the hyper-elliptical surface generated by $`(x,W(x))`$ (In the 1-matrix case the surface is hyperelliptical, i.e. $`W(x)`$ obeys a quadratic equation). It is only in the one-cut case that there is no oscillations, and the $`1/N^2`$ expansion can exist.

We expect the same idea to hold for the 2-matrix model. The method of  can be adapted to the 2-matrix case: the frequencies are the contour integrals of $`W(x)dx`$, which are the same as the contour integrals of $`Y(x)dx`$ (which are also the same as the integrals of $`X(y)dy`$ by integration by parts), along the non-trivial cycles of the algebraic curve $`(x,Y(x))`$. The $`1/N^2`$ expansion cannot exist if there are oscillations. Therefore we will require that there is no oscillation, i.e. that there is no non-trivial cycle, i.e. that the algebraic curve has genus zero.

Notice a difference between the 1-matrix case and the 2-matrix case:

The number of connected parts of the support of the density is actually the number of cuts in the first sheet only. There may be cuts in other sheets.

\- In the 1-matrix case, the curve is hyper-elliptical, it has only two sheets, thus the first sheet contains all the cuts. The number of non-trivial cycles is equal to the number of cuts minus one. Therefore in the 1-matrix model, the genus zero asumption is equivalent to the 1-cut asumption and is equivalent to a density with connected support.

\- In the 2-matrix case, the curve is only algebraic. It can have more than two sheets, so that the non-trivial cycles are not necessarily in the first sheet. One cut in the first sheet is not equivalent to genus zero. In fact, there are $`d_2+1`$ sheets and at least $`d_2`$ cuts, $`d_21`$ of them are NOT in the first sheet. In the 2-matrix case, the connectedness of the support of the density is a necessary condition, but is not sufficient. The necessary and sufficient condition is the genus zero condition.

Therefore, in the 2-matrix model, we claim that the existence of the $`1/N^2`$ expansion is related to the asumption that the algebraic curve $`E^{(0)}(x,y)=0`$ has genus zero.

From now on, we will consider only the genus zero case. We would like to emphasize that many of the results derived below for the leading order loop functions, can be derived for arbitrary genus with no difficulty (for instance the functions $`\mathrm{\Omega }`$ are always the Bergmann kernels), but are out of the scope of this article. Many of them are presented in .

The previous discussion applies to the case where the partition function eq. (2.1) is well defined, i.e. the potentials satisfy condition eq. (2.3). However, there is another case of “existence” of the $`1/N^2`$ expansion: when the genus $`h`$ free energy $`F^{(h)}`$ is formally defined by its diagrammatic expansion eq. (1.4) and eq. (2.4). The formal series for $`F^{(h)}`$ is obtained by expanding in the vicinity of a gaussian integral , i.e. by considering that the potentials are “small deformations” of quadratic potentials. When both potentials are quadratic, the algebraic curve $`E(x,y)=0`$ clearly has genus zero. The formal expansion in power series of the $`g_k`$’s and $`g_l^{}`$’s can not change the genus. Therefore, the formal diagrammatic expansion of $`F^{(h)}`$ is also obtained by a genus zero asumption.

### 4.2 Genus zero asumption

Since we require that $`E(x,y)=0`$ be a genus zero algebraic curve, there must exist a rational uniformization. Such a rational uniformization appears naturally in the framework of bi-orthogonal polynomials as shown in . That leads us to choose it of the following form:

$$x=𝒳(s)=\gamma s+\underset{k=0}{\overset{d_2}{}}\frac{\alpha _k}{s^k},y=𝒴(s)=\frac{\gamma }{s}+\underset{j=0}{\overset{d_1}{}}\beta _js^j$$
(4.1)
This means that for every $`(x,y)`$ which satisfy $`E(x,y)=0`$, there exists at least one $`s`$ such that $`x=𝒳(s)`$ and $`y=𝒴(s)`$.

We note the inverse functions:

$$x=𝒳(s)s=\sigma (x),y=𝒴(s)s=\stackrel{~}{\sigma }(y)$$
(4.2)
The functions $`\sigma (x)`$ and $`\stackrel{~}{\sigma }(y)`$ are multivalued, we will discuss their sheet structure below. The functions $`Y(x)`$ and $`X(y)`$ are:

$$Y(x)=𝒴(\sigma (x)),X(y)=𝒳(\stackrel{~}{\sigma }(y))$$
(4.3)
They are multivalued too, and their sheet structure will be discussed below.

### 4.3 The parameters $`\alpha _k`$ and $`\beta _j`$

The $`\alpha `$’s and $`\beta `$’s are merely a reparametrization of the coefficients of the potentials.

Indeed, since for large $`x`$ we have: $`V_1^{}(x)Y(x)=W(x)\frac{1}{x}+O(1/x^2)`$ (this corresponds to $`s\mathrm{}`$), and for large $`y`$ we have: $`V_2^{}(y)X(y)=\stackrel{~}{W}(y)\frac{1}{y}+O(1/y^2)`$ ($`s0`$), we must have:

$`𝒴(s)V_1^{}(𝒳(s))`$ $`\underset{s\mathrm{}}{}`$ $`{\displaystyle \frac{1}{\gamma s}}+O(1/s^2)`$ (4.4)
$`𝒳(s)V_2^{}(𝒴(s))`$ $`\underset{s0}{}`$ $`{\displaystyle \frac{s}{\gamma }}+O(s^2)`$ (4.5)
Both these two equations imply that (compute $`\mathrm{Res}𝒴(s)𝒳^{}(s)\text{d}s=\mathrm{Res}𝒳(s)𝒴^{}(s)\text{d}s`$):

$`\gamma ^2=1+{\displaystyle \underset{k=1}{\overset{\mathrm{min}(d_1,d_2)}{}}}k\alpha _k\beta _k`$      (4.6)
And eqs.(4.4) allow to write the coupling constants $`g_l`$ and $`g_l^{}`$ as functions of the $`\alpha _k`$’s and $`\beta _k`$’s (compute $`\mathrm{Res}V_1^{}(x)x^{l1}\text{d}x`$ and $`\mathrm{Res}V_2^{}(y)y^{l1}\text{d}y`$):

$$g_{l+1}=\frac{1}{\gamma ^{l+1}}\underset{k=1}{\overset{d_1l1}{}}\left(\underset{p=0}{\overset{d_2}{}}\underset{i_p=0}{\overset{[d_1/(p+1)]}{}}\right)k\frac{(l+\underset{p}{}i_p)!}{l!_pi_p!}\alpha _k\beta _{(l+k+1+_p(p+1)i_p)}\underset{p=0}{\overset{d_2}{}}\left(\frac{\alpha _p}{\gamma }\right)^{i_p}$$
(4.7)
$$g_{l+1}^{}=\frac{1}{\gamma ^{l+1}}\underset{k=1}{\overset{d_2l1}{}}\left(\underset{p=0}{\overset{d_1}{}}\underset{i_p=0}{\overset{[d_2/(p+1)]}{}}\right)k\frac{(l+\underset{p}{}i_p)!}{l!_pi_p!}\beta _k\alpha _{(l+k+1+_p(p+1)i_p)}\underset{p=0}{\overset{d_1}{}}\left(\frac{\beta _p}{\gamma }\right)^{i_p}$$
(4.8)
where we take the convention that $`\beta _j=0`$ if $`j>d_1`$ and $`\alpha _j=0`$ if $`j>d_2`$, and $`\alpha _1=\beta _1=\gamma `$.

As an example, for $`l=d_1`$ and $`l=d_11`$, the formula reduces to:

$$g_{d_1+1}=\frac{\beta _{d_1}}{\gamma ^{d_1}},g_{d_1}=\frac{\beta _{d_11}}{\gamma ^{d_11}}d_1\alpha _0\frac{\beta _{d_1}}{\gamma ^{d_1}},\mathrm{}$$
(4.9)
In principle, it should be possible to revert these formula, and compute the $`\alpha `$’s and $`\beta `$’s as functions of the coupling constants. This can be done at least numerically.

Remark:

Notice that there might exist more than one solution. Only the solution which leads to an absolute minimum of the leading order free energy $`F^{(0)}`$ should be considered. If the potentials $`V_1`$ and $`V_2`$ are real (which is the case for most physical models), the density of eigenvalues of the first (resp. second) matrix is the imaginary part of $`Y(x)`$ (resp. $`X(y)`$) along the cut, and the density must be positive. Any solution for the $`\alpha `$’s and $`\beta `$’s which does not satisfy that is unacceptable. If there is no acceptable solution, this means that our choice of $`V_1`$ and $`V_2`$ does not correspond to a genus zero case. In general, if $`V_1`$ and $`V_2`$ have only one well, or if they are close to quadratic potentials, we have a genus zero solution.

### 4.4 Sheet structure of the function $`Y(x)`$

The function $`\sigma (x)`$ is multivalued, indeed the equation

$$𝒳(s)=x$$
(4.10)
has $`d_2+1`$ solutions, which we note

$$\{\sigma _0(x),\sigma _1(x),\mathrm{},\sigma _{d_2}(x)\}$$
(4.11)
and therefore, the function $`Y(x)=𝒴(\sigma (x))`$ is multivalued with $`d_2+1`$ values which we note:

$$\{Y_0(x),Y_1(x),\mathrm{},Y_{d_2}(x)\}\mathrm{where}Y_k(x)=𝒴(\sigma _k(x))$$
(4.12)
We define $`\sigma (x)=\sigma _0(x)`$ to be the ”physical sheet” solution, i.e. the one such that $`\sigma (x)\mathrm{}`$ when $`x\mathrm{}`$, i.e. the one for which $`W(x)=V_1^{}(x)Y(x)O(1/x)`$. The other solutions $`\sigma _k(x)`$ with $`k=1,\mathrm{},d_2`$ all tend to zero when $`x\mathrm{}`$.

#### 4.4.1 Endpoints and cuts

The functions $`Y(x)`$ and $`\sigma (x)`$ have cuts, whith a generic square root behaviour near the endpoints, which means that the derivative with respect to $`x`$ diverges at the endpoints. Since $`Y^{}(x)=\frac{𝒴^{}(s)}{𝒳^{}(s)}`$, the endpoints are solutions of $`s=0`$ or solutions of the equation:

$$𝒳^{}(s)=0$$
(4.13)
which has $`d_2+1`$ solutions. We note them:

$$\{e_0,e_1,\mathrm{},e_{d_2}\}$$
(4.14)
the endpoints are thus located at the $`𝒳(e_k)`$’s, as well as at infinity. In the physical sheet, the function $`Y(x)`$ is analytical at infinity, and the only cut in the physical sheet is a line joining two endpoints of the type $`𝒳(e_k)`$. In the other sheets, the function $`Y(x)`$ behaves as $`x^{1/d_2}`$ near infinity, the cuts are lines joining some $`𝒳(e_k)`$ to $`\mathrm{}`$. There are exactly $`d_2`$ cuts, in agreement with the fact that the genus is zero.

Similarly, the function $`X(y)`$ has exactly $`d_1`$ cuts, whose endpoints are $`𝒴(\stackrel{~}{e}_j)`$’s and $`\mathrm{}`$, where

$$\{\stackrel{~}{e}_0,\stackrel{~}{e}_1,\mathrm{},\stackrel{~}{e}_{d_1}\}$$
(4.15)
are the $`d_1+1`$ roots of

$$𝒴^{}(s)=0.$$
(4.16)
The cut in the physical sheet is a line joining two endpoints of the type $`𝒴(\stackrel{~}{e}_k)`$, and the cuts in the other sheets are lines joining some $`𝒴(\stackrel{~}{e}_k)`$ to $`\mathrm{}`$.

As said in the introduction, we assume that the potentials are not critical, which means that the algebraic curve has no singular points (other than $`s=0`$ and $`s=\mathrm{}`$). That means that all the $`e_k`$ and $`\stackrel{~}{e}_l`$ must be distinct (examples of critical potentials are given in : those which produce the rational minimal conformal models $`(p,q)`$). The non-criticality asumption, is crucial for the derivations we are going to present below: we will often have $`𝒳^{\prime \prime }(e_k)`$ or $`𝒴^{}(e_k)`$ in the denominators.

#### 4.4.2 Sheets in the $`s`$-plane

We define the $`x`$-sheets (resp. $`y`$-sheets) in the $`s`$-plane, as domains which contain only one $`\sigma _k(x)`$ (resp. $`\stackrel{~}{\sigma }_k(y)`$) when $`x`$ (resp. $`y`$) sweeps the complex plane. In other words, the function $`𝒳(s)`$ (resp. $`𝒴(s)`$) restricted to one domain is one to one.

The domains are separated by closed curves in the $`s`$-plane, which are the projections of the cuts in the $`xplane`$ (resp. $`y`$-plane). There is some arbitrariness in the choice of domains and cuts (one could take the cuts as straight lines or any other curves which link the endpoints). The only constraint is that all the cuts go through the endpoints.

The $`x`$-sheets (resp. $`y`$-sheets) have the following properties:
\- The physical sheet is the sheet which contains $`\mathrm{}`$ (resp. $`0`$).
\- All the $`d_2`$ (resp. $`d_1`$) other sheets contain $`0`$ (resp. $`\mathrm{}`$).

And the cuts:
\- The cut in the physical sheet corresponds to a contour in the $`s`$-plane which encircles $`0`$. If the potentials are real, we expect the cut in the $`x`$-plane (reps. $`y`$-plane) to be on the real axis, and its $`s`$-plane image to be a contour which crosses the real axis, and symmetric with respect to the real axis. The two endpoints are solutions of $`𝒳^{}(s)=0`$ (resp. $`𝒴^{}(s)=0`$).
\- The other $`d_21`$ (reps. $`d_11`$) cuts are contours in the other sheets. They all meet in $`0`$ (resp. $`\mathrm{}`$). One of their endpoint is a solution of $`𝒳^{}(s)=0`$ (resp. $`𝒴(s)=0`$) and the other endpoint is at $`s=0`$ (reps. $`s=\mathrm{}`$).

notational remark: when we write $`\sigma (x)`$ this means that $`\sigma (x)`$ is considered as a multivalued function. For instance the multivalued function $`Y(x)`$ is equal to

$$Y(x)=𝒴(\sigma (x))$$
(4.17)
We write $`\sigma _k(x)`$ or $`Y_k(x)`$ only when we want to specify in which sheet is $`x`$. By default, we will consider that $`x`$ is in the physical sheet, and therefore $`\sigma _0(x)`$ is in the domain which contains $`s=\mathrm{}`$, and $`\sigma _k(x)`$ with $`k1`$ is in a domain which contains $`s=0`$. Note that any symmetric function of $`(\sigma _0(x),\mathrm{},\sigma _{d_2}(x))`$ is a regular function of $`x`$, i.e. it has no cut.

### 4.5 The polynomial $`P(x,y)`$

The polynomial $`P(x,y)`$ is found from $`E(𝒳(s),𝒴(s))=0`$, i.e. :

$$P(𝒳(s),𝒴(s))=(V_1^{}(𝒳(s))𝒴(s))(V_2^{}(𝒴(s))𝒳(s))+1$$
(4.18)
both sides are Laurent polynomials of $`s`$ of the same degree. By indentifying the coefficients of powers of $`s`$ on both sides, one can compute explicitely all the coefficients of the polynomial $`P(x,y)`$. The expression of $`P(x,y)`$ in terms of $`\alpha _k`$ and $`\beta _j`$ can be written explicitely with the determinant formula:

$$E(x,y)=\frac{1}{\gamma ^{d_1+d_2}}det\left(\begin{array}{cccccccc}\gamma & \alpha _0x& \alpha _1& \mathrm{}& \alpha _{d_2}& 0& \mathrm{}& 0\\ 0& \mathrm{}& & & & \mathrm{}& & 0\\ \mathrm{}& & \mathrm{}& & & & \mathrm{}& \mathrm{}\\ 0& \mathrm{}& 0& \gamma & \alpha _0x& \alpha _1& \mathrm{}& \alpha _{d_2}\\ \beta _{d_1}& \mathrm{}& \beta _1& \beta _0y& \gamma & 0& \mathrm{}& 0\\ 0& \mathrm{}& & & & \mathrm{}& & \mathrm{}\\ \mathrm{}& & \mathrm{}& & & & \mathrm{}& 0\\ 0& \mathrm{}& 0& \beta _{d_1}& \mathrm{}& \beta _1& \beta _0y& \gamma \end{array}\right)$$
(4.19)
This determinant of size $`d_1+d_2+2`$ is the discriminant which vanishes if and only if the Laurent polynomials in $`s`$: $`𝒳(s)x`$ and $`𝒴(s)y`$ have a common root.

There is also the formula:

$$E(x,y)=(g_{d_1+1})^{d_2+1}\frac{1}{\gamma ^{2d_2}}\underset{i=0}{\overset{d_2}{}}\underset{j=0}{\overset{d_1}{}}(\gamma \sigma _i(x)\gamma \stackrel{~}{\sigma }_j(y))$$
(4.20)
indeed, it is a symmetric function of the $`\sigma _i(x)`$ and $`\stackrel{~}{\sigma }_j(y)`$, with no poles, therefore it is a polynomial in $`x`$ and $`y`$, and it vanishes when $`𝒳(s)x`$ and $`𝒴(s)y`$ have a common root.

Remark: Had we not made the genus zero asumption, we could determine $`P(x,y)`$ by a determinant formula888Thanks to J. Hurtubise for that remark. very similar to eq. (4.19). We can also determine $`P(x,y)`$ by requiring that $`E(x,y)=0`$ is a genus $`g`$ curve, and that all the b-cycle integrals of $`Y(x)\text{d}x`$ vanish. There is one (or more) solution for all $`g`$. For any $`g`$ between $`0`$ and $`\left[\frac{d_1d_21}{2}\right]`$ we can compute all the loop functions (to leading order), and therefore we can compute the free energy $`F^{(0)}(g)`$ as a function of the genus $`g`$. The actual value of $`g`$ is the one which minimizes $`F^{(0)}`$. The genus zero asumption is thus that the minimum is obtained for $`g=0`$.

## 5 Leading order loop functions

### 5.1 Some leading order Loop equations

The change of variables (we drop the $`O(1/N^2)`$ terms):

* $`M_1M_1+ϵ\frac{1}{xM_1}\frac{1}{yM_2}+ϵ\frac{1}{yM_2}\frac{1}{xM_1}`$ implies:

$$(yY(x))W(x,y)=W(x)\stackrel{~}{U}(x,y)=V_1^{}(x)Y(x)\stackrel{~}{U}(x,y)$$
(5.1)
* $`M_2M_2+ϵ\frac{1}{xM_1}\frac{1}{yM_2}+ϵ\frac{1}{yM_2}\frac{1}{xM_1}`$ implies:

$$(xX(y))W(x,y)=\stackrel{~}{W}(y)U(x,y)=V_2^{}(y)X(y)U(x,y)$$
(5.2)
* $`M_1M_1+ϵ\frac{1}{xM_1}\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}+ϵ\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}\frac{1}{xM_1}`$ implies:

$$(yY(x))U(x,y)=P(x,y)+V_2^{}(y)W(x)\frac{1}{N}<\mathrm{tr}\frac{1}{xM_1}V_2^{}(M_2)>$$
(5.3)
* $`M_2M_2+ϵ\frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}\frac{1}{yM_2}+ϵ\frac{1}{yM_2}\frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}`$ implies:

$$(xX(y))\stackrel{~}{U}(x,y)=P(x,y)+V_1^{}(x)\stackrel{~}{W}(y)\frac{1}{N}<\mathrm{tr}\frac{1}{yM_2}V_1^{}(M_1)>$$
(5.4)
* $`M_2M_2+ϵ\frac{1}{xM_1}`$ implies:

$$\frac{1}{N}<\mathrm{tr}\frac{1}{xM_1}V_2^{}(M_2)>=xW(x)1$$
(5.5)
* $`M_1M_1+ϵ\frac{1}{yM_2}`$ implies:

$$\frac{1}{N}<\mathrm{tr}\frac{1}{yM_2}V_1^{}(M_1)>=y\stackrel{~}{W}(y)1$$
(5.6)

### 5.2 The function $`U(x,y)`$

Using eq. (5.3) and eq. (5.5), we find that to leading order, the function $`U(x,y)`$ is simply:

$`U(x,y)=V_2^{}(y)x+{\displaystyle \frac{E(x,y)}{yY(x)}}`$      (5.7)
Notice that this function has no pole, and it is indeed a polynomial in $`y`$.

### 5.3 $`W(x,y)`$

Using eq. (5.2) and eq. (5.7), we find that, to leading order, the function $`W(x,y)`$ is:

$`W(x,y)=1{\displaystyle \frac{E(x,y)}{(xX(y))(yY(x))}}`$      (5.8)
It can also be written:

$$W(x,y)=1\left(\frac{\beta _{d_1}}{\gamma }\right)^{d_2}\frac{\sigma _0(x)\stackrel{~}{\sigma }_0(y)^{d_2}}{\sigma _0(x)\stackrel{~}{\sigma }_0(y)}\underset{i=1}{\overset{d_2}{}}\underset{j=1}{\overset{d_1}{}}(\sigma _i(x)\stackrel{~}{\sigma }_j(y))$$
(5.9)
it has a pole when $`\sigma _0(x)=\stackrel{~}{\sigma }_0(y)`$.

We are now going to compute the two-loop functions.

### 5.4 The functions $`\mathrm{\Omega }(x;x^{})`$ and $`\stackrel{~}{\mathrm{\Omega }}(y;x^{})`$

These functions were first computed by  using the orthogonal polynomial’s technics, they can also be found by another method, namely: they must be rational functions of $`\sigma (x)`$ (resp. $`\stackrel{~}{\sigma }(y)`$) and $`\sigma (x^{})`$, and they are determined completely by their poles and behaviours near $`\mathrm{}`$. The result is the following:

$`\mathrm{\Omega }(x;x^{})`$ $`=`$ $`{\displaystyle \frac{}{V_1(x^{})}}W(x)=_x_x^{}\mathrm{ln}{\displaystyle \frac{\sigma (x)\sigma (x^{})}{xx^{}}}={\displaystyle \frac{1}{(xx^{})^2}}+{\displaystyle \frac{\sigma ^{}(x)\sigma ^{}(x^{})}{(\sigma (x)\sigma (x^{}))^2}}`$ (5.10)
$`=`$ $`{\displaystyle \frac{1}{(xx^{})^2}}+{\displaystyle \frac{1}{𝒳^{}(\sigma (x))𝒳^{}(\sigma (x^{}))(\sigma (x)\sigma (x^{}))^2}}`$ (5.11)
and since $`Y(x)=V_1^{}(x)W(x)`$ we have:

$$\frac{}{V_1(x^{})}Y(x)=\frac{\sigma ^{}(x)\sigma ^{}(x^{})}{(\sigma (x)\sigma (x^{}))^2}=_x_x^{}\mathrm{ln}(\sigma (x)\sigma (x^{}))$$
(5.12)
i.e. $`\frac{}{V_1(x^{})}Y(x)`$ is the Bergmann kernel999 Without the genus zero asumption, one has $`\frac{}{V_1(x^{})}Y^{(0)}(x)=`$ Bergman kernel, but $`\frac{}{V_1(x^{})}Y(x)=`$ square of the Szego kernel $`+O(1/N)`$..

At $`x=x^{}`$ we get the Schwarzian derivative of $`\sigma (x)`$:

$$\mathrm{\Omega }(x;x)=\frac{1}{6}\frac{\sigma ^{\prime \prime \prime }(x)}{\sigma ^{}(x)}\frac{1}{4}\frac{\sigma ^{\prime \prime }(x)^2}{\sigma ^{}(x)^2}$$
(5.13)
Similarly , the function $`\stackrel{~}{\mathrm{\Omega }}(y;x^{})`$ is:

$$\stackrel{~}{\mathrm{\Omega }}(y;x^{})=\frac{}{V_1(x^{})}\stackrel{~}{W}(y)=x^{}_y\mathrm{ln}(\sigma (x^{})\stackrel{~}{\sigma }(y))=\frac{\sigma ^{}(x^{})\stackrel{~}{\sigma }^{}(y)}{(\sigma (x^{})\stackrel{~}{\sigma }(y))^2}$$
(5.14)
i.e.:

$$\frac{}{V_1(x^{})}X(y)=\frac{\stackrel{~}{\sigma }^{}(y)\sigma ^{}(x^{})}{(\stackrel{~}{\sigma }(y)\sigma (x^{}))^2}=_y_x^{}\mathrm{ln}(\stackrel{~}{\sigma }(y)\sigma (x^{}))$$
(5.15)
### 5.5 More loop equations

* The change of variables $`\delta M_2=\frac{1}{xM_1}\mathrm{tr}\frac{1}{x^{}M_1}`$ implies:

$$\mathrm{tr}\frac{1}{xM_1}V_2^{}(M_2)\mathrm{tr}\frac{1}{x^{}M_1}_\mathrm{c}=x\mathrm{\Omega }(x;x^{})$$
(5.16)
* Using this result, the change of variables $`\delta M_1=\frac{1}{xM_1}\frac{V_2^{}(y)V_2^{}(M_2)}{yM_2}\mathrm{tr}\frac{1}{x^{}M_1}+\mathrm{h}.\mathrm{c}.`$ implies:

$`(yY(x))U(x,y;x^{})`$ $`=`$ $`(V_2^{}(y)xU(x,y))\mathrm{\Omega }(x;x^{})P(x,y;x^{})`$ (5.18)
$`_x^{}{\displaystyle \frac{U(x,y)U(x^{},y)}{xx^{}}}`$
i.e., using eq. (5.7):

$$(yY(x))U(x,y;x^{})=\mathrm{\Omega }(x;x^{})\frac{E(x,y)}{yY(x)}P(x,y;x^{})_x^{}\frac{U(x,y)U(x^{},y)}{xx^{}}$$
(5.19)
which can also be written, using eq. (5.7) again:

$`(yY(x))U(x,y;x^{})`$ $`=`$ $`(\mathrm{\Omega }(x;x^{})+{\displaystyle \frac{1}{(xx^{})^2}}){\displaystyle \frac{E(x,y)}{yY(x)}}P(x,y;x^{})`$ (5.21)
$`+_x^{}{\displaystyle \frac{E(x^{},y)}{(xx^{})(yY(x^{}))}}`$
* The change of variables $`\delta M_1=\frac{1}{yM_2}\mathrm{tr}\frac{1}{x^{}M_1}`$ implies:

$$\mathrm{tr}V_1^{}(M_1)\frac{1}{yM_2}\mathrm{tr}\frac{1}{x^{}M_1}_\mathrm{c}=y\stackrel{~}{\mathrm{\Omega }}(y;x^{})_x^{}W(x^{},y)$$
(5.22)
* Using this result, the change of variables $`\delta M_2=\frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}\frac{1}{yM_2}\mathrm{tr}\frac{1}{x^{}M_1}+\mathrm{h}.\mathrm{c}.`$ implies:

$$(xX(y))\stackrel{~}{U}(x,y;x^{})=\stackrel{~}{\mathrm{\Omega }}(y;x^{})(V_1^{}(x)\stackrel{~}{U}(x,y)y)P(x,y;x^{})+_x^{}W(x^{},y)$$
(5.23)
which can also be written:

$$(xX(y))\stackrel{~}{U}(x,y;x^{})=P(x,y;x^{})\stackrel{~}{\mathrm{\Omega }}(y;x^{})\frac{E(x,y)}{xX(y)}+_x^{}W(x^{},y)$$
(5.24)

### 5.6 Determination of the Polynomial $`P(x,y;x^{})`$

In particular, if for any $`k=0,\mathrm{},d_2`$, we choose $`y=Y_k(x)`$ in eq. (5.24) , we have $`X(y)=x`$, and using eq. (5.8):

$$P(x,Y_k(x);x^{})=_x^{}\left[\frac{E(x^{},Y_k(x))}{(xx^{})(Y_k(x)Y(x^{}))}\right]\stackrel{~}{\mathrm{\Omega }}(Y_k(x);x^{})E_x(x,Y_k(x))$$
(5.25)
i.e.

$$P(x,Y_k(x);x^{})=_x^{}\left[\frac{E(x^{},Y_k(x))}{(xx^{})(Y_k(x)Y(x^{}))}+\frac{E_x(x,Y_k(x))}{𝒴^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))}\right]$$
(5.26)
We know that $`P`$ is a polynomial in $`y`$ of degree $`d_21`$, and we know its value in $`d_2+1`$ points, therefore we can determine it with the interpolation formula (see appendix A), which, in this case, reduces to:

$$P(x,y;x^{})=\underset{k=0}{\overset{d_2}{}}\frac{P(x,Y_k(x);x^{})}{(yY_k(x))}\frac{E(x,y)}{E_y(x,Y_k(x))}$$
(5.27)
The latter expression seems to be a polynomial of degree $`d_2`$ in $`y`$, but one can check that the leading term vanishes (see appendix B). Thus we may also write:

$$P(x,y;x^{})=\underset{k=1}{\overset{d_2}{}}\frac{P(x,Y_k(x);x^{})(Y_k(x)Y_0(x))}{(yY_k(x))(yY_0(x))}\frac{E(x,y)}{E_y(x,Y_k(x))}$$
(5.28)
Both expressions give either:

$`P(x,y;x^{})=`$ $`_x^{}{\displaystyle \underset{k=1}{\overset{d_2}{}}}{\displaystyle \frac{E_x(x,Y_k(x))E(x,y)(Y_k(x)Y_0(x))}{𝒴^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))E_y(x,Y_k(x))(yY_k(x))(yY_0(x))}}`$ (5.30)
$`+{\displaystyle \frac{E(x^{},Y_k(x))}{(xx^{})(Y_k(x)Y(x^{}))}}{\displaystyle \frac{E(x,y)}{E_y(x,Y_k(x))}}{\displaystyle \frac{Y_k(x)Y_0(x)}{(yY_k(x))(yY_0(x))}}`$
Using $`E(𝒳(s),𝒴(s))=0`$ we have $`𝒴^{}(\sigma _k(x))E_y(x,Y_k(x))=𝒳^{}(\sigma _k(x))E_x(x,Y_k(x))`$, and we can write:

$`P(x,y;x^{})=`$ $`_x^{}{\displaystyle \underset{k=1}{\overset{d_2}{}}}{\displaystyle \frac{E(x,y)}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))}}{\displaystyle \frac{Y_0(x)Y_k(x)}{(yY_k(x))(yY_0(x))}}`$ (5.33)
$`+{\displaystyle \frac{E(x^{},Y_k(x))}{(xx^{})(Y_k(x)Y(x^{}))}}{\displaystyle \frac{E(x,y)}{E_y(x,Y_k(x))}}{\displaystyle \frac{Y_k(x)Y_0(x)}{(yY_k(x))(yY_0(x))}}`$
or, using eq. (5.27):

$`P(x,y;x^{})`$ $`=_x^{}{\displaystyle \underset{k=0}{\overset{d_2}{}}}`$ $`{\displaystyle \frac{E(x,y)}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))}}{\displaystyle \frac{1}{(yY_k(x))}}`$ (5.36)
$`{\displaystyle \frac{E(x^{},Y_k(x))}{(xx^{})(Y_k(x)Y(x^{}))}}{\displaystyle \frac{E(x,y)}{E_y(x,Y_k(x))}}{\displaystyle \frac{1}{(yY_k(x))}}`$
Eq.(5.36) is clearly a polynomial in $`y`$, let us explain how one can check that it is also a polynomial in $`x`$: it is symmetric in the $`\sigma _k(x)`$’s, therefore it has no cut (it takes the same value in all sheets), it must be a rational function of $`x`$. The only possible poles could be at $`x=x^{}`$, $`𝒳^{}(\sigma _k(x))=0`$ or $`E_y(x,Y_k(x))=0`$.

\- Take $`x=x^{}`$, more precisely, for some $`k[0,d_2]`$, $`\sigma _k(x)=\sigma (x^{})`$, because of the symmetry one may assume that $`k=0`$, and one sees that expression eq. (5.33) has no pole at $`\sigma _0(x)=\sigma (x^{})`$.

\- Now choose some $`k`$ and $`x`$ such that $`𝒳^{}(\sigma _k(x))=ϵ`$ is very small. There must exist some (unique because our potentials are non critical) $`jk`$ such that $`\sigma _j(x)\sigma _k(x)=O(ϵ)`$, and one has $`𝒳^{}(\sigma _j(x))=𝒳^{}(\sigma _k(x))+O(ϵ^2)`$. The poles for $`j`$ and $`k`$ thus cancel each other.

\- Now choose some $`k`$ and $`x`$ such that $`E_y(x,Y_k(x))=ϵ`$ is very small. This means there must exist some (unique because our potentials are non critical) $`jk`$ such that $`Y_j(x)Y_k(x)=O(ϵ)`$, and one has $`E_y(x,Y_j(x))=E_y(x,Y_k(x))+O(ϵ^2)`$. The poles for $`j`$ and $`k`$ thus cancel each other in eq. (5.36).

### 5.7 Determination of $`U(x,y;x^{})`$

We can now compute $`U(x,y;x^{})`$ using eq. (5.21):

$`(yY_0(x))U(x,y;x^{})`$ $`=_x^{}`$ $`{\displaystyle \underset{k=1}{\overset{d_2}{}}}{\displaystyle \frac{E(x,y)}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))(yY_k(x))}}`$ (5.40)
$`{\displaystyle \frac{E(x,y)}{(xx^{})}}{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{E(x^{},Y_k(x))}{(Y_k(x)Y(x^{}))E_y(x,Y_k(x))(yY_k(x))}}`$
$`+{\displaystyle \frac{E(x^{},y)}{(xx^{})(yY(x^{}))}}`$
which simplifies to:

$`(yY_0(x))U(x,y;x^{})`$ $`=_x^{}`$ $`{\displaystyle \underset{k=1}{\overset{d_2}{}}}{\displaystyle \frac{E(x,y)}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))(yY_k(x))}}`$ (5.43)
$`+{\displaystyle \frac{E(x,y)}{(xx^{})(yY(x^{}))}}{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{E(x^{},Y_k(x))}{E_y(x,Y_k(x))(Y_k(x)y)}}`$
$`+{\displaystyle \frac{E(x^{},y)E(x,y)}{(xx^{})(yY(x^{}))}}`$
Note the following identity (see appendix B, Lemma B.5):

$$\underset{k=0}{\overset{d_2}{}}\frac{\underset{j=0}{\overset{d_2}{}}(Y_k(x)Y_j(x^{}))}{(Y_k(x)y)_{jk}(Y_k(x)Y_j(x))}=1\frac{\underset{j}{}(yY_j(x^{}))}{_j(yY_j(x))}$$
(5.44)
which implies:

$$\underset{k=0}{\overset{d_2}{}}\frac{E(x^{},Y_k(x))}{(Y_k(x)y)E_y(x,Y_k(x))}=1\frac{E(x^{},y)}{E(x,y)}$$
(5.45)
and thus:

$`U(x,y;x^{})=_x^{}{\displaystyle \underset{k=1}{\overset{d_2}{}}}{\displaystyle \frac{E(x,y)}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))(yY_k(x))(yY_0(x))}}`$      (5.46)
It is easy to check that this expression is a polynomial in $`y`$ and goes to $`0`$ when $`x\mathrm{}`$ in the physical sheet ($`s_0\mathrm{}`$) or when $`x^{}\mathrm{}`$ in the physical sheet. It has inverse square root singularities near the endpoints where $`𝒳^{}(\sigma _k(x))=0`$.

Computation at $`y=Y_0(x)`$

$$U(x,Y_0(x);x^{})=_x^{}\underset{k=1}{\overset{d_2}{}}\frac{E_y(x,Y_0(x))}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))(Y_0(x)Y_k(x))}$$
(5.47)
Computation at $`x^{}=x`$

$$U(x,Y_0(x);x)=\underset{k=1}{\overset{d_2}{}}\frac{E_y(x,Y_0(x))}{𝒳^{}(\sigma _k(x))𝒳^{}(\sigma _0(x))(\sigma _k(x)\sigma _0(x))^2(Y_0(x)Y_k(x))}$$
(5.48)
We are now equipped to compute the next to leading order functions…

## 6 Next to leading order

### 6.1 The one-loop function: computation of $`Y^{(1)}`$

We have (see eq. (3.25)):

$$Y^{(1)}(x)=\frac{P^{(1)}(x,Y_0(x))}{E_y(x,Y_0(x))}+\underset{k=1}{\overset{d_2}{}}\frac{1}{𝒳^{}(\sigma _k(x))𝒳^{}(\sigma _0(x))(\sigma _k(x)\sigma _0(x))^2(Y_0(x)Y_k(x))}$$
(6.1)
where $`P^{(1)}`$ has degrees $`(d_11,d_21)`$, and the coefficient of $`x^{d_11}y^{d_21}`$ vanishes. We thus have $`d_1d_21`$ unknown coefficients to determine. On the other hand, there are $`2d_1d_22`$ values of $`s`$ for which $`E_y(𝒳(s),𝒴(s))=0`$ which are not endpoints (see appendix D, eq. (D.10)). If we write that the poles at all these points cancel in equation 6.1, we can determine $`P^{(1)}(x,y)`$. However, we don’t need to determine $`P^{(1)}(x,y)`$, we will follow another method which allows to determine $`Y^{(1)}(x)`$ directly.

Eq.(6.1) is a rational function whose only poles are the endpoints, with degree up to 5. Since eq. (6.1) behaves as $`O(s^2)`$ when $`s\mathrm{}`$, $`Y^{(1)}`$ has no pole at $`\mathrm{}`$, and since eq. (6.1) behaves as $`O(s^{1+d_2})`$ when $`s0`$, $`𝒳^{}(s)Y^{(1)}`$ has no pole at $`0`$.

Therefore we may write:

$$Y^{(1)}(𝒳(s))=\frac{1}{𝒳^{}(s)}\underset{k=0}{\overset{d_2}{}}\frac{A_k}{(se_k)^4}+\frac{B_k}{(se_k)^3}+\frac{C_k}{(se_k)^2}+\frac{D_k}{(se_k)}$$
(6.2)
The coefficients $`A_k`$, $`B_k`$, $`C_k`$, $`D_k`$ are determined by matching the poles in eq. (6.1) ($`P^{(1)}`$ doesnot contribute to them). We will see below that $`D_k=0`$. This is not surprising since we expect $`Y^{(1)}(x)`$ to be a derivative (indeed $`Y^{(1)}(x)=/V_1(x)F^{(1)}=d/dx/V_1^{}(x)F^{(1)}`$).

Let $`k[0,d_2]`$, and choose $`s`$ close to $`e_k`$:

$$s=e_k+ϵ$$
(6.3)
there must exist $`\stackrel{~}{s}`$ (unique because the potentials are non-critical) such that $`𝒳(\stackrel{~}{s})=𝒳(s)`$ and $`\stackrel{~}{s}`$ is close to $`e_k`$ ($`\stackrel{~}{s}`$ is the $`\sigma _k(x)`$ in eq. (6.1)):

$$\stackrel{~}{s}=e_k\eta ,\eta =O(ϵ)$$
(6.4)
By solving $`𝒳(s)=𝒳(\stackrel{~}{s})`$ order by order in $`ϵ`$ we get:

$$\eta =\lambda ϵ,\lambda =1+r_kϵ+r_k^2ϵ^2+(2r_k^3+t_k)ϵ^3+O(ϵ^4)$$
(6.5)
where

$$r_k=\frac{1}{3}\frac{𝒳^{\prime \prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)},s_k=\frac{1}{6}\frac{𝒳^{IV}(e_k)}{𝒳^{\prime \prime }(e_k)},t_k=\frac{1}{60}\frac{𝒳^V(e_k)}{𝒳^{\prime \prime }(e_k)}r_ks_k$$
(6.6)
From eq. (6.1) we must have:

$$A_k+B_kϵ+C_kϵ^2+D_kϵ^3=\frac{ϵ^4}{𝒳^{}(e_k\eta )(ϵ+\eta )^2(𝒴(e_k+ϵ)𝒴(e_k\eta ))}+O(ϵ^4)$$
(6.7)
We note the 3rd degree polynomial:

$$P_k(ϵ)=A_k+B_kϵ+C_kϵ^2+D_kϵ^3$$
(6.8)
i.e.

$`P_k(ϵ)`$ $`=`$ $`(1+\lambda )^2`$ (6.12)
$`(\lambda 𝒳^{\prime \prime }(e_k)+ϵ{\displaystyle \frac{\lambda ^2}{2}}𝒳^{\prime \prime \prime }(e_k)ϵ^2{\displaystyle \frac{\lambda ^3}{6}}𝒳^{IV}(e_k)+ϵ^3{\displaystyle \frac{\lambda ^4}{24}}𝒳^V(e_k))^1`$
$`((1+\lambda )𝒴^{}(e_k)+{\displaystyle \frac{ϵ}{2}}(1\lambda ^2)𝒴^{\prime \prime }(e_k)+{\displaystyle \frac{ϵ^2}{6}}(1+\lambda ^3)𝒴^{\prime \prime \prime }(e_k)`$
$`+{\displaystyle \frac{ϵ^3}{24}}(1\lambda ^4)𝒴^{IV}(e_k))^1+O(ϵ^4)`$
$`=`$ $`(𝒳^{\prime \prime }(e_k)𝒴^{}(e_k))^1(1+\lambda )^3\lambda ^1`$ (6.16)
$`(1{\displaystyle \frac{3}{2}}ϵ\lambda r_k+ϵ^2\lambda ^2s_k{\displaystyle \frac{5}{2}}ϵ^3(t_k+r_ks_k))^1`$
$`(1+{\displaystyle \frac{ϵ}{2}}(1\lambda ){\displaystyle \frac{𝒴^{\prime \prime }(e_k)}{𝒴^{}(e_k)}}+{\displaystyle \frac{ϵ^2}{6}}(1\lambda +\lambda ^2){\displaystyle \frac{𝒴^{\prime \prime \prime }(e_k)}{𝒴^{}(e_k)}})^1`$
$`+O(ϵ^4)`$
It is easy to see that $`t_k`$ as well as $`𝒴^{IV}(e_k)`$ disappear, and we find that $`D_k=0`$.

After a straightforward calculation, one finds:

$`8𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(A_k+B_kϵ+C_kϵ^2+D_kϵ^3)`$ (6.17)
$`=1{\displaystyle \frac{1}{3}}ϵ{\displaystyle \frac{𝒳^{\prime \prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)}}+{\displaystyle \frac{1}{6}}ϵ^2\left({\displaystyle \frac{𝒳^{\prime \prime \prime }(e_k)^2}{𝒳^{\prime \prime }(e_k)^2}}{\displaystyle \frac{𝒳^{IV}(e_k)}{𝒳^{\prime \prime }(e_k)}}+{\displaystyle \frac{𝒳^{\prime \prime \prime }(e_k)𝒴^{\prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)}}{\displaystyle \frac{𝒴^{\prime \prime \prime }(e_k)}{𝒴^{}(e_k)}}\right)`$ (6.18)
After substitution into eq. (6.2), we find the genus one correction to the resolvent:

$`Y^{(1)}(𝒳(s))`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{1}{8𝒳^{}(s)𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^4}}`$ (6.22)
$`{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{𝒳^{\prime \prime \prime }(e_k)}{24𝒳^{}(s)𝒳^{\prime \prime 2}(e_k)𝒴^{}(e_k)(se_k)^3}}`$
$`{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{\frac{𝒳^{\prime \prime \prime }(e_k)^2}{𝒳^{\prime \prime }(e_k)^2}\frac{𝒳^{IV}(e_k)}{𝒳^{\prime \prime }(e_k)}+\frac{𝒳^{\prime \prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)}\frac{𝒴^{\prime \prime }(e_k)}{𝒴^{}(e_k)}\frac{𝒴^{\prime \prime \prime }(e_k)}{𝒴^{}(e_k)}}{48𝒳^{}(s)𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^2}}`$
That function represents the partition function of a statistical physics model on a genus one surface with one + boundary.

### 6.2 The free energy

We are now going to find the free energy $`F^{(1)}`$ such that:

$$Y^{(1)}(x)=\frac{F^{(1)}}{V_1(x)}$$
(6.24)
For that purpose, we need to compute the derivatives of various quantities with respect to $`V_1(x)`$. From our knowledge of the one matrix case , we guess that $`F^{(1)}`$ will be related to the “moments” $`𝒴^{}(e_k)`$. Therefore, we shall compute the action of $`/V_1(x)`$ on such quantities.

#### 6.2.1 The derivatives of $`𝒴(u)`$ and $`𝒳(u)`$

In this section, we note $`x=𝒳(s)`$ and:

$$\dot{𝒳}(u)=𝒳^{}(s)\frac{𝒳(u)}{V_1(x)},\dot{𝒴}(u)=𝒳^{}(s)\frac{𝒴(u)}{V_1(x)}$$
(6.25)
Notice that (using eq. (5.12) and $`\frac{𝒴(u)}{V_1(x)}=\frac{Y(𝒳(u))}{V_1(x)}=\frac{𝒳(u)}{V_1(x)}Y^{}(𝒳(u))+\frac{Y(x^{})}{V_1(x)}|_{x^{}=𝒳(u)}`$):

$$𝒴^{}(u)\dot{𝒳}(u)𝒳^{}(u)\dot{𝒴}(u)=\frac{1}{(su)^2}$$
(6.26)
In particular, if $`e_k`$, $`k=0,\mathrm{},d_2`$ is an endpoint $`𝒳^{}(e_k)=0`$ we have:

$$\dot{𝒳}(e_k)=\frac{1}{𝒴^{}(e_k)(se_k)^2}$$
(6.27)
Now, notice that $`\dot{𝒳}(u)`$ is a Laurent Polynomial101010Indeed $`d_2`$ is independent of $`V_1(x)`$. $`\dot{𝒴}(u)`$ is not a Laurent polynomial, because the degree of $`𝒴(u)`$ depends on $`V_1`$. in $`u`$:

$$\dot{𝒳}(u)=\dot{\gamma }u+\underset{k=0}{\overset{d_2}{}}\dot{\alpha }_ku^k$$
(6.28)
where

$$\dot{\gamma }=𝒳^{}(s)\frac{\gamma }{V_1(x)},\dot{\alpha }_k=𝒳^{}(s)\frac{\alpha _k}{V_1(x)}k=0,\mathrm{},d_2$$
(6.29)
We know its value at all $`e_k`$’s, $`k=0\mathrm{}d_2`$, therefore (see appendix A):

$$\dot{𝒳}(u)=u𝒳^{}(u)\left[\frac{\dot{\gamma }}{\gamma }+\underset{k=0}{\overset{d_2}{}}\frac{1}{e_k(se_k)^2(ue_k)𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)}\right]$$
(6.30)
Since $`\alpha _{d_2}=g_{d_2+1}^{}\gamma ^{d_2}`$, by comparing the coefficient of $`u^1`$ and $`u^{d_2}`$, we find:

$$\frac{\dot{\gamma }}{\gamma }=\frac{1}{d_2}\frac{\dot{\alpha }_{d_2}}{\alpha _{d_2}}=\frac{1}{2}\underset{k=0}{\overset{d_2}{}}\frac{1}{(se_k)^2e_k^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)}$$
(6.31)
and thus:

$$𝒳^{}(s)\frac{𝒳(u)}{V_1(x)}=\dot{𝒳}(u)=\frac{u𝒳^{}(u)}{2}\underset{k=0}{\overset{d_2}{}}\frac{u+e_k}{ue_k}\frac{1}{(se_k)^2e_k^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)}$$
(6.32)
Notice that, as a function of $`s`$, this expression has double poles at the endpoints $`e_k`$.

For any $`e_i`$, $`i=0,\mathrm{},d_2`$, $`𝒳^{}(e_i)=0`$ implies that $`𝒳^{}(s)\frac{e_i}{V_1(x)}=\frac{\dot{𝒳}^{}(e_i)}{𝒳^{\prime \prime }(e_i)}`$ and thus:

$`\dot{e_i}=𝒳^{}(s){\displaystyle \frac{e_i}{V_1(x)}}={\displaystyle \frac{\dot{𝒳}^{}(e_i)}{𝒳^{\prime \prime }(e_i)}}`$ $`=`$ $`{\displaystyle \frac{e_i}{2}}{\displaystyle \underset{ki}{}}{\displaystyle \frac{e_i+e_k}{e_ie_k}}{\displaystyle \frac{1}{(se_k)^2e_k^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)}}`$ (6.34)
$`{\displaystyle \frac{3+e_i\frac{𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)}}{2e_i(se_i)^2𝒳^{\prime \prime }(e_i)𝒴^{}(e_i)}}`$
From eq. (6.26) and eq. (6.32) we have:

$$\dot{𝒴}(u)=\frac{1}{𝒳^{}(u)}\left(𝒴^{}(u)\dot{𝒳}(u)\frac{1}{(su)^2}\right)$$
(6.35)
which implies:

$`𝒳^{}(s){\displaystyle \frac{𝒴^{}(e_i)}{V_1(x)}}=`$ $`\dot{𝒴}^{}(e_i)+`$ $`\dot{e_i}𝒴^{\prime \prime }(e_i)`$ (6.36)
$`=`$ $`{\displaystyle \frac{1}{2𝒳^{\prime \prime }(e_i)}}(`$ $`\dot{𝒳}(e_i)𝒴^{\prime \prime \prime }(e_i)+\dot{𝒳}^{\prime \prime }(e_i)𝒴^{}(e_i){\displaystyle \frac{6}{(se_i)^4}}`$ (6.39)
$`\dot{𝒳}(e_i)𝒴^{\prime \prime }(e_i){\displaystyle \frac{𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)}}\dot{𝒳}^{}(e_i)𝒴^{}(e_i){\displaystyle \frac{𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)}}`$
$`+{\displaystyle \frac{2}{(se_i)^3}}{\displaystyle \frac{𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)}})`$
This expression has poles of degree $`2`$, $`3`$ and $`4`$, at the $`e_k`$’s. We write:

$$𝒳^{}(s)\frac{𝒴^{}(e_i)}{V_1(x)}=\underset{k=0}{\overset{d_2}{}}\frac{A_{k,i}}{(se_k)^4}+\frac{B_{k,i}}{(se_k)^3}+\frac{C_{k,i}}{(se_k)^2}$$
(6.41)
with:

$$A_{k,i}=3\frac{\delta _{k,i}}{𝒳^{\prime \prime }(e_i)},B_{k,i}=\frac{\delta _{k,i}𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)^2}$$
(6.42)
and:

$$C_{k,i}=\frac{\left(g_k(e_i)(𝒴^{\prime \prime \prime }(e_i)𝒴^{\prime \prime }(e_i)\frac{𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)})g_k^{}(e_i)𝒴^{}(e_i)\frac{𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)}+g_k^{\prime \prime }(e_i)𝒴^{}(e_i)\right)}{2𝒳^{\prime \prime }(e_i)e_k^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)}$$
(6.43)
where

$$g_k(u)=\frac{u𝒳^{}(u)}{2}\frac{u+e_k}{ue_k}$$
(6.44)
in particular:

$$g_i(e_i)=e_i^2𝒳^{\prime \prime }(e_i),g_i^{}(e_i)=\frac{3}{2}e_i𝒳^{\prime \prime }(e_i)+\frac{1}{2}e_i^2𝒳^{\prime \prime \prime }(e_i)$$
(6.45)
$$g_i^{\prime \prime }(e_i)=𝒳^{\prime \prime }(e_i)+\frac{3}{2}e_i𝒳^{\prime \prime \prime }(e_i)+\frac{1}{3}e_i^2𝒳^{IV}(e_i)$$
(6.46)
and if $`ki`$:

$$g_k(e_i)=0,g_k^{}(e_i)=\frac{e_i}{2}\frac{e_i+e_k}{e_ie_k}𝒳^{\prime \prime }(e_i)$$
(6.47)
$$g_k^{\prime \prime }(e_i)=\frac{e_i+e_k}{e_ie_k}𝒳^{\prime \prime }(e_i)2\frac{e_ie_k}{(e_ie_k)^2}𝒳^{\prime \prime }(e_i)+\frac{e_i}{2}\frac{e_i+e_k}{e_ie_k}𝒳^{\prime \prime \prime }(e_i)$$
(6.48)
We thus have:

$$C_{i,i}=\frac{𝒴^{\prime \prime \prime }(e_i)𝒴^{\prime \prime }(e_i)\frac{𝒳^{\prime \prime \prime }(e_i)}{𝒳^{\prime \prime }(e_i)}}{2𝒳^{\prime \prime }(e_i)𝒴^{}(e_i)}\frac{𝒳^{\prime \prime \prime }(e_i)^2}{4𝒳^{\prime \prime }(e_i)^3}+\frac{1}{2e_i^2𝒳^{\prime \prime }(e_i)}+\frac{𝒳^{IV}(e_i)}{6𝒳^{\prime \prime }(e_i)^2}$$
(6.49)
and if $`ki`$:

$$C_{k,i}=\frac{(e_i^2e_k^22e_ie_k)𝒴^{}(e_i)}{2e_k^2(e_ie_k)^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)}$$
(6.50)
Now, let us compute the following quantity (our guess is that it is related to $`Y^{(1)}`$):

$`𝒳^{}(s){\displaystyle \underset{i=0}{\overset{d_2}{}}}{\displaystyle \frac{\mathrm{ln}𝒴^{}(e_i)}{V_1(x)}}`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{A_{k,k}}{𝒴^{}(e_k)(se_k)^4}}+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{B_{k,k}}{𝒴^{}(e_k)(se_k)^3}}`$ (6.52)
$`+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{C_{k,k}}{𝒴^{}(e_k)(se_k)^2}}+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \underset{ik}{}}{\displaystyle \frac{C_{k,i}}{𝒴^{}(e_i)(se_k)^2}}`$
substituting eqs.(6.42, 6.49, 6.50) into eq. (6.52), we get:

$`𝒳^{}(s){\displaystyle \underset{i=0}{\overset{d_2}{}}}{\displaystyle \frac{\mathrm{ln}𝒴^{}(e_i)}{V_1(x)}}`$ $`=`$ $`3{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{1}{𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^4}}`$ (6.58)
$`+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{𝒳^{\prime \prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)^2𝒴^{}(e_k)(se_k)^3}}`$
$`+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{\frac{𝒴^{\prime \prime \prime }(e_k)}{𝒴^{}(e_k)}\frac{𝒴^{\prime \prime }(e_k)}{𝒴^{}(e_k)}\frac{𝒳^{\prime \prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)}\frac{1}{2}\frac{𝒳^{\prime \prime \prime }(e_k)^2}{𝒳^{\prime \prime }(e_k)^2}+\frac{1}{3}\frac{𝒳^{IV}(e_k)}{𝒳^{\prime \prime }(e_k)}}{2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^2}}`$
$`+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{1}{2e_k^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^2}}`$
$`+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \underset{ik}{}}{\displaystyle \frac{(e_i^2e_k^22e_ie_k)}{2e_k^2(e_ie_k)^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^2}}`$
Notice the following identity111111One proves it by taking the log derivative of $`u^{d_2+1}𝒳^{}(u)=_k(ue_k)`$:

$`1+{\displaystyle \underset{ik}{}}{\displaystyle \frac{(e_i^2e_k^22e_ie_k)}{(e_ie_k)^2}}`$ $`=`$ $`d_2+12e_k^2{\displaystyle \underset{ik}{}}{\displaystyle \frac{1}{(e_ie_k)^2}}`$ (6.60)
$`=`$ $`(d_2+1)+{\displaystyle \frac{2}{3}}e_k^2{\displaystyle \frac{𝒳^{IV}(e_k)}{𝒳^{\prime \prime }(e_k)}}{\displaystyle \frac{e_k^2}{2}}{\displaystyle \frac{𝒳^{\prime \prime \prime }(e_k)^2}{𝒳^{\prime \prime }(e_k)^2}}`$ (6.61)
therefore:

$`𝒳^{}(s){\displaystyle \underset{i=0}{\overset{d_2}{}}}{\displaystyle \frac{\mathrm{ln}𝒴^{}(e_i)}{V_1(x)}}`$ $`=`$ $`3{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{1}{𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^4}}`$ (6.66)
$`+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{𝒳^{\prime \prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)^2𝒴^{}(e_k)(se_k)^3}}`$
$`+{\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{\frac{𝒴^{\prime \prime \prime }(e_k)}{𝒴^{}(e_k)}\frac{𝒴^{\prime \prime }(e_k)}{𝒴^{}(e_k)}\frac{𝒳^{\prime \prime \prime }(e_k)}{𝒳^{\prime \prime }(e_k)}\frac{𝒳^{\prime \prime \prime }(e_k)^2}{𝒳^{\prime \prime }(e_k)^2}+\frac{𝒳^{IV}(e_k)}{𝒳^{\prime \prime }(e_k)}}{2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^2}}`$
$`(d_2+1){\displaystyle \underset{k=0}{\overset{d_2}{}}}{\displaystyle \frac{1}{2e_k^2𝒳^{\prime \prime }(e_k)𝒴^{}(e_k)(se_k)^2}}`$
We recognize:

$$\underset{i=0}{\overset{d_2}{}}\frac{\mathrm{ln}𝒴^{}(e_i)}{V_1(x)}=24Y^{(1)}(x)(d_2+1)\frac{\mathrm{ln}\gamma }{V_1(x)}$$
(6.68)
#### 6.2.2 The free energy

The free energy is such that:

$$W(x)=\frac{F}{V_1(x)}=V_1^{}(x)Y(x)$$
(6.69)
therefore to order $`1/N^2`$ we have:

$$\frac{F^{(1)}}{V_1(x)}=Y^{(1)}(x)$$
(6.70)
which implies (up to a constant $`K`$ independent of $`V_1`$):

$$F^{(1)}=\frac{1}{24}\mathrm{ln}\left(\gamma ^4D\right)+K$$
(6.71)
where:

$$D=\frac{1}{\gamma ^{d_1+d_2+2}}det\left(\begin{array}{cccccccc}\gamma & 0& \alpha _1& \mathrm{}& d_2\alpha _{d_2}& 0& \mathrm{}& 0\\ 0& \mathrm{}& & & & \mathrm{}& & 0\\ \mathrm{}& & \mathrm{}& & & & \mathrm{}& \mathrm{}\\ 0& \mathrm{}& 0& \gamma & 0& \alpha _1& \mathrm{}& d_2\alpha _{d_2}\\ d_1\beta _{d_1}& \mathrm{}& \beta _1& 0& \gamma & 0& \mathrm{}& 0\\ 0& \mathrm{}& & & & \mathrm{}& & \mathrm{}\\ \mathrm{}& & \mathrm{}& & & & \mathrm{}& 0\\ 0& \mathrm{}& 0& d_1\beta _{d_1}& \mathrm{}& \beta _1& 0& \gamma \end{array}\right)$$
(6.72)
$`D`$ is such that:

$$D=(1)^{d_1+1}\left(\frac{d_2\alpha _{d_2}}{\gamma }\right)^2\underset{𝒳^{}(e)=0}{}\frac{𝒴^{}(e)}{\gamma }=(1)^{d_1+d_2}\gamma ^2\underset{𝒴^{}(\stackrel{~}{e})=0}{}\left(\frac{𝒳^{}(\stackrel{~}{e})}{\gamma }\right)$$
(6.73)
The Constant $`K`$ is independent of $`V_1`$, and if we repeat the calculation by exchanging the roles of $`M_1`$ and $`M_2`$, we find the same expression for $`F^{(1)}`$ with $`K`$ independent of $`V_2`$ as well. By comparing with the gaussian case we get $`K=0`$.

We have found the genus one free energy:

$`F^{(1)}={\displaystyle \frac{1}{24}}\mathrm{ln}\left(\gamma ^4D\right)`$      (6.74)
Remark: $`D`$ is a determinant which vanishes everytime the surface is singular, i.e. when the potentials are critical. We see that near a critical point, the genus one free energy has a logarithmic divergence. This is in agreement with known double-scaling limit exponents.

### 6.3 Higher genus expansion

The procedure can be repeated to higher orders. We have the three equations:

$$E(x,Y(x))=(V_1^{}(x)Y(x))(V_2^{}(Y(x))x)P(x,Y(x))+1=\frac{1}{N^2}U(x,Y(x);x)$$
(6.75)
$$(yY(x))U(x,y)=(V_1^{}(x)Y(x))(V_2^{}(y)x)P(x,y)+1\frac{1}{N^2}U(x,y;x)$$
(6.76)
$$U(x,y;x^{})=\frac{U(x,y)}{V_1(x^{})}$$
(6.77)
Imagine that we already know the functions $`Y(x)`$, $`E(x,y)`$ (i.e. $`P(x,y)`$), $`U(x,y)`$ and $`U(x,y;x^{})`$ up to order $`h1`$. Let us write eq. (6.75) to order $`h`$. we know all terms but $`Y^{(h)}(x)`$ and $`P^{(h)}(x,y)`$. The factor of the $`Y^{(h)}(x)`$ term is $`E_y^{(0)}(x,Y^{(0)}(x))`$, it vanishes for all $`(x,y)`$ such that $`y=Y^{(0)}(x)`$ and $`E_y^{(0)}(x,Y^{(0)}(x))`$, i.e. at $`2d_1d_22`$ points, therefore we know $`P^{(h)}(x,y)`$ in $`2d_1d_22`$ points. Since $`P^{(h)}(x,y)`$ is a polynomial of degree $`(d_11,d_21)`$, with $`d_1d_21`$ unknown coefficients, $`P^{(h)}(x,y)`$ can be determined completely. Then, from eq. (6.75), $`Y^{(h)}(x)`$ can be determined. Using eq. (6.76), we can determine $`U^{(h)}(x,y)`$, and using eq. (6.77) we can determine $`U^{(h)}(x,y;x^{})`$ to order $`h`$.

We clearly have an algorithmic recursive procedure to find $`Y^{(h)}(x)`$ to any order $`h`$.

Finding the free energy $`F^{(h)}`$ is more difficult, as one needs to integrate with respect to the potentials. The easiest is to have an ansatz for $`F^{(h)}`$ and check that its derivative with respect to $`V_1(x)`$ is indeed $`Y^{(h)}(x)`$, as in .

Define:

$$M_{0,i}=𝒴^{}(e_i),M_{k,i}=\frac{𝒴^{(k+1)}(e_i)}{𝒴^{}(e_i)}k1$$
(6.78)
$$J_{0,i}=𝒳^{\prime \prime }(e_i),J_{k,i}=\frac{𝒳^{(k+2)}(e_i)}{𝒳^{\prime \prime }(e_i)}k1$$
(6.79)
We expect the genus $`h>1`$ free energy to be a rational function of these quantities, of the form:

$$F^{(h)}=\underset{r,s,p,q,m0}{}\underset{k_1,\mathrm{},k_s}{}\underset{l_1,\mathrm{},l_s}{}A_{k_1,\mathrm{},k_s,l_1,\mathrm{},l_s,p,q,m}\frac{1}{\gamma ^m}\underset{i=0}{\overset{d_2}{}}\frac{\underset{j=1}{\overset{r}{}}\underset{u=1}{\overset{s}{}}M_{k_j,i}J_{l_u,i}}{M_{0,i}^pJ_{0,i}^q}$$
(6.80)
with the restrictions that:

$$p+q=\underset{k}{}k_j+\underset{u}{}l_u$$
(6.81)
The selection rules (of the form $`m+p+q4h4`$ as in ) are still to be understood.

## 7 The one matrix case: $`d_2=1`$

In this section we check that for $`d_2=1`$, our results reduce to those of .

Consider $`d_2=1`$ and $`V_2(y)=\frac{1}{2}g_2^{}y^2+g_1^{}y`$. $`V_2`$ is thus quadratic and $`M_2`$ can easily be integrated out. The 2-matrix model in this case reduces to a one matrix model with the potential:

$$V(x)=V_1(x)\frac{1}{2g_2^{}}x^2+\frac{g_1^{}}{g_2^{}}x$$
(7.1)
The function $`𝒳(s)`$ and $`𝒴(s)`$ are:

$$𝒳(s)=\gamma s+\alpha _0+\frac{\alpha _1}{s},𝒴(s)=\frac{\gamma }{s}+\underset{k=0}{\overset{d_1}{}}\beta _ks^k$$
(7.2)
From eq. (4.4) we find that:

$$\alpha _1=g_2^{}\gamma ,\alpha _0=g_2^{}\beta _0+g_1^{}$$
(7.3)
and:

$$V_1^{}(𝒳(s))=𝒴(s)+𝒴(g_2^{}/s)\beta _0\frac{\gamma }{s}\frac{\gamma }{g_2^{}}s$$
(7.4)
indeed this expression is a polynomial in $`x`$ and it satisfies eq. (4.4). It follows from eq. (7.1):

$$V^{}(𝒳(s))=𝒴(s)+𝒴(g_2^{}/s)2(\beta _0+\frac{\gamma }{s}+\frac{\gamma }{g_2^{}}s)$$
(7.5)
The endpoints $`a`$ and $`b`$ are the roots of $`𝒳^{}(s)`$ i.e. they correspond to:

$$e=\sqrt{g_2^{}},a=𝒳(e),b=𝒳(+e)$$
(7.6)
The resolvent (defined in eq. (2.6)) can be written :

$$W(x)=V_1^{}(x)Y(x)=\frac{1}{2}\left(V^{}(x)M(x)\sqrt{(xa)(xb)}\right)$$
(7.7)
where $`M(x)`$ is a polynomial in $`x`$ (the same as in ). Using eq. (7.5) and eq. (7.2) we find that:

$$M(𝒳(s))=\frac{𝒴(s)𝒴(g_2^{}/s)}{s𝒳^{}(s)},\sqrt{(𝒳(s)a)(𝒳(s)b)}=s𝒳^{}(s)$$
(7.8)
In particular at the endpoints we have:

$$M(a)=\frac{𝒴^{}(e)}{\gamma },M(b)=\frac{𝒴^{}(+e)}{\gamma }$$
(7.9)
It is clear now, that our result eq. (6.74) coincides with the genus one free energy found by the authors of :

$$F^{(1)}=\frac{1}{24}\mathrm{ln}\left(\left(\frac{ba}{4}\right)^4M(a)M(b)\right)$$
(7.10)
## 8 Conclusion and prospects

The recursive method presented in this article, allows one to compute the resolvents as well as other loop-functions, to any order in the topological expansion of the 2-matrix model. We have found the expression of the genus one free energy, but an algorithmic method to compute the free energy to higher orders is still incomplete. For that purpose, it would be usefull to define some “moments” and basis functions, as in , onto which one could decompose the loop-functions, and selection rules to know what set of moments should appear in the free energy at a given order.

It would be interesting also to understand what terms survive in the double scaling limit when the potentials are near a critical point, and to develop a recursive method to compute directly the free energy and resolvents in the double scaling limit.

As we said in the introduction, we have only developped a method to compute the terms in the $`1/N^2`$ expansion, but we have not proven the existence of that expansion. A rigorous proof of that existence can probably be obtained by Riemann-Hilbert technics .

It would also be usefull to understand what the $`1/N^2`$ expansion becomes in the non genus-zero case. It is very likely that it would involve the $`\theta `$-functions associated to the algebraic curve, as in .

Aknowledgements: The author would like to thank the CRM, Prof. J. Harnad and Prof. J. Hurtubise for their support when this work was completed.

## Appendix AThe interpolation formula

Let $`d`$ be an integer, and $`x_0,\mathrm{},x_d`$ be $`d+1`$ distinct complex numbers. We note:

$$\pi (x)=\underset{k=0}{\overset{d}{}}(xx_k)$$
(A.1)
Let $`P_0,\mathrm{},P_d`$ be $`d+1`$ complex numbers.

There is a unique polynomial $`P`$ of degree $`d`$ such that:

$$P(x_k)=P_kk=0,\mathrm{},d$$
(A.2)
That polynomial is:

$$P(x)=\underset{k=0}{\overset{d}{}}P_k\frac{\pi (x)}{(xx_k)\pi ^{}(x_k)}$$
(A.3)
Any polynomial $`P(x)`$ of degree $`>d`$ such that:

$$P(x_k)=P_kk=0,\mathrm{},d$$
(A.4)
can be written:

$$P(x)=Q(x)\pi (x)+\underset{k=0}{\overset{d}{}}P_k\frac{\pi (x)}{(xx_k)\pi ^{}(x_k)}$$
(A.5)
where $`Q(x)`$ is a polynomial of degree $`\mathrm{deg}Pd1`$.

## Appendix BProof that eq. (5.27) is a polynomial of degree $`d_21`$

The leading term of eq. (5.27) is (we use eq. (5.36)):

$$x^{}\frac{1}{xx^{}}\underset{k=0}{\overset{d_2}{}}\left[\frac{E(x^{},Y_k(x))}{(Y_k(x)Y(x^{}))}\frac{1}{E_y(x,Y_k(x))}\frac{xx^{}}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))}\right]$$
(B.1)
We have (see eq. (D.5)):

$$\underset{k=0}{\overset{d_2}{}}\frac{1}{𝒳^{}(\sigma _k(x))(\sigma _k(x)\sigma (x^{}))}$$
(B.2)
Therefore we have to prove that:

$$\underset{k=0}{\overset{d_2}{}}\frac{E(x^{},Y_k(x))}{(Y_k(x)Y(x^{}))}\frac{1}{E_y(x,Y_k(x))}=1$$
(B.3)
Note that:

$$\underset{k=0}{\overset{d_2}{}}\frac{E(x^{},Y_k(x))}{(Y_k(x)Y(x^{}))}\frac{1}{E_y(x,Y_k(x))}=\underset{k=0}{\overset{d_2}{}}\frac{\underset{j0}{}(Y_k(x)Y_j(x^{}))}{_{jk}(Y_k(x)Y_j(x))}$$
(B.4)
Lemma:

For any $`Y_0,\mathrm{},Y_{d_2}`$ and for any $`Y_1^{},\mathrm{},Y_{d_2}^{}`$ we have:

$$\underset{k=0}{\overset{d_2}{}}\frac{\underset{j0}{}(Y_kY_j^{})}{_{jk}(Y_kY_j)}=1$$
(B.5)
Proof:

$$\underset{k=0}{\overset{d_2}{}}\frac{\underset{j0}{}(Y_kY_j^{})}{_{jk}(Y_kY_j)}=\underset{k=0}{\overset{d_2}{}}(1)^kf(Y_k)\frac{\mathrm{\Delta }_{d_2}(Y_0,\mathrm{},Y_{k1},Y_{k+1},\mathrm{},Y_{d_2})}{\mathrm{\Delta }_{d_2+1}(Y_0,\mathrm{},Y_{d_2})}$$
(B.6)
where $`f(Y)=_{j0}(YY_j^{})`$ is a monic polynomial of degree $`d_2`$ in $`Y`$, and $`\mathrm{\Delta }_n`$ is the Vandermonde determinant of $`n`$ variables.

We recognize a determinant (which is $`1`$ by linearly combining columns):

$$\underset{k=0}{\overset{d_2}{}}\frac{\underset{j0}{}(Y_kY_j^{})}{_{jk}(Y_kY_j^{})}=\frac{det\left(\begin{array}{cccccc}1& Y_0& Y_0^2& \mathrm{}& Y_0^{d_21}& f(Y_0)\\ 1& Y_1& Y_1^2& \mathrm{}& Y_1^{d_21}& f(Y_1)\\ \mathrm{}& & & & & \mathrm{}\\ 1& Y_0& Y_{d_2}^2& \mathrm{}& Y_{d_2}^{d_21}& f(Y_{d_2})\end{array}\right)}{\mathrm{\Delta }_{d_2+1}(Y_0,\mathrm{},y_{d_2})}=1$$
(B.7)
## Appendix COther loop-functions

Define:

$`H(x,y,x^{},y^{})`$ $`=`$ $`{\displaystyle \frac{1}{2N}}\mathrm{tr}{\displaystyle \frac{1}{xM_1}}{\displaystyle \frac{1}{yM_2}}{\displaystyle \frac{1}{x^{}M_1}}{\displaystyle \frac{1}{y^{}M_2}}`$ (C.3)
$`+{\displaystyle \frac{1}{2N}}\mathrm{tr}{\displaystyle \frac{1}{y^{}M_2}}{\displaystyle \frac{1}{x^{}M_1}}{\displaystyle \frac{1}{yM_2}}{\displaystyle \frac{1}{xM_1}}`$
$`=\text{}`$
Notice that we have:

$$H(x,y,x^{},y^{})=H(x^{},y^{},x,y)=H(x,y^{},x^{},y)$$
(C.4)
We also define:

$`F(x,y,x^{},y^{})`$ $`=`$ $`{\displaystyle \frac{1}{2N}}\mathrm{tr}{\displaystyle \frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}}{\displaystyle \frac{1}{yM_2}}{\displaystyle \frac{1}{x^{}M_1}}{\displaystyle \frac{1}{y^{}M_2}}`$ (C.6)
$`+{\displaystyle \frac{1}{2N}}\mathrm{tr}{\displaystyle \frac{1}{y^{}M_2}}{\displaystyle \frac{1}{x^{}M_1}}{\displaystyle \frac{1}{yM_2}}{\displaystyle \frac{V_1^{}(x)V_1^{}(M_1)}{xM_1}}`$
we have:

$$F(x,y,x^{},y^{})=F(x,y^{},x^{},y)$$
(C.7)
Loop equations:

The invariance under the change of variable
$`\delta M_1=\frac{1}{xM_1}\frac{1}{yM_2}\frac{1}{x^{}M_1}\frac{1}{y^{}M_2}+\mathrm{h}.\mathrm{c}.`$ implies:

$$(y^{}Y(x))H(x,y,x^{},y^{})=\frac{W(x,y)W(x^{},y)}{xx^{}}(1W(x^{},y^{}))F(x,y,x^{},y^{})$$
(C.8)
The same equation with $`yy^{}`$ gives:

$`H(x,y,x^{},y^{})={\displaystyle \frac{1}{(xx^{})(yy^{})}}[{\displaystyle \frac{E(x^{},y)E(x,y^{})}{(x^{}X)(yY^{})(xX^{})(y^{}Y)}}`$ (C.9)
$`{\displaystyle \frac{E(x,y)E(x^{},y^{})}{(xX)(yY)(x^{}X^{})(y^{}Y^{})}}]`$ (C.10)
it follows:

$`F(x,y,x^{},y^{})`$ $`=`$ $`{\displaystyle \frac{E(x,y)E(x^{},y^{})}{(xx^{})(yy^{})(xX)(x^{}X^{})(y^{}Y^{})}}`$ (C.13)
$`{\displaystyle \frac{E(x^{},y)E(x,y^{})}{(xx^{})(yy^{})(xX^{})(x^{}X)(yY^{})}}`$
$`{\displaystyle \frac{E(x^{},y)E(x^{},y^{})}{(xx^{})(x^{}X)(x^{}X^{})(yY^{})(y^{}Y^{})}}`$
## Appendix DSome relations

The following formula are usefull for many calculations in this article. They are easy to prove.

$$E(x,y)=A\underset{i,j}{}(\sigma _i(x)\stackrel{~}{\sigma }_j(y))$$
(D.1)
with

$$A=(1)^{d_2+1}\gamma ^{1+d_1d_2}\beta _{d_1}^{d_2}g_{d_1+1}=(1)^{d_2+1}\frac{\beta _{d_1}^{d_2+1}}{\gamma ^{d_21}}$$
(D.2)
$$𝒳(s)x=\frac{\gamma }{s^{d_2}}\underset{i=0}{\overset{d_2}{}}(s\sigma _i(x)),𝒴(s)y=\frac{\beta _{d_1}}{s}\underset{j=0}{\overset{d_1}{}}(s\stackrel{~}{\sigma }_j(y))$$
(D.3)
$$𝒳^{}(\sigma _k(x))=\frac{\gamma }{\sigma _k^{d_2}(x)}\underset{ik}{}(\sigma _k(x)\sigma _i(x)),𝒴^{}(\stackrel{~}{\sigma }_k(y))=\frac{\beta _{d_1}}{\stackrel{~}{\sigma }_k(y)}\underset{jk}{}(\stackrel{~}{\sigma }_k(y)\stackrel{~}{\sigma }_j(y))$$
(D.4)
$$\underset{i=0}{\overset{d_2}{}}\frac{1}{𝒳^{}(\sigma _i(x))(s\sigma _i(x))}=\frac{1}{𝒳(s)x},\underset{j=0}{\overset{d_1}{}}\frac{1}{𝒴^{}(\stackrel{~}{\sigma }_j(y))(s\stackrel{~}{\sigma }_j(y))}=\frac{1}{𝒴(s)y}$$
(D.5)
$`E_y(x,Y_0(x))`$ $`=`$ $`A{\displaystyle \frac{1}{𝒴^{}(\sigma _0(x))}}{\displaystyle \underset{(i,j)(0,0)}{}}(\sigma _i(x)\stackrel{~}{\sigma }_j(Y_0(x)))`$ (D.6)
$`=`$ $`{\displaystyle \frac{\beta _{d_1}^{d_2}}{\gamma ^{d_2}}}\sigma _0^{d_2+1}(x)𝒳^{}(\sigma _0(x)){\displaystyle \underset{i0}{}}{\displaystyle \underset{j0}{}}(\sigma _i(x)\stackrel{~}{\sigma }_j(Y_0(x)))`$ (D.7)
$$Y_k(x)Y_0(x)=\frac{\beta _{d_1}}{\sigma _k(x)}\underset{j}{}(\sigma _k(x)\stackrel{~}{\sigma }_j(Y_0(x)))$$
(D.8)
Limits:

$$E_y(𝒳(s),𝒴(s))\underset{s\mathrm{}}{}\frac{\alpha _{d_2}\beta _{d_1}^{d_2}}{\gamma ^{d_2}}s^{d_1d_2},E_y(𝒳(s),𝒴(s))\underset{s0}{}d_2\frac{\beta _{d_1}\alpha _{d_2}^{d_1+1}}{\gamma ^{d_1+1}}s^{1d_1d_2d_2}$$
(D.9)
$`E_y(𝒳(s),𝒴(s))`$ is a Laurent polynomial:

$$E_y(𝒳(s),𝒴(s))=\frac{𝒳^{}(s)}{s^{2d_1d_2}}P_{2d_1d_22}(s)$$
(D.10)
where $`P_{2d_1d_22}(s)`$ is a polynomial of degree $`2d_1d_22`$. Therefore, $`E_y(𝒳(s),𝒴(s))=0`$ has $`2d_1d_22`$ zeroes.