Schwinger-Dyson and loop equations for a product of square Ginibre random matrices

TL;DR

This paper derives loop equations for the resolvents of a product of two Ginibre matrices using Schwinger-Dyson techniques, providing insights into their spectral properties and corrections beyond the large N limit.

Contribution

It introduces a novel approach to derive loop equations for matrix products without reformulating in terms of singular values, and analyzes the spectral curve geometry.

Findings

Large N limit of the 2-point correlation function's Stieltjes transform

First correction to the density's Stieltjes transform

Explicit results for various resolvents and their corrections

Abstract

In this paper, we study the product of two complex Ginibre matrices and the loop equations satisfied by their resolvents (i.e. the Stieltjes transform of the correlation functions). We obtain using Schwinger-Dyson equation (SDE) techniques the general loop equations satisfied by the resolvents. In order to deal with the product structure of the random matrix of interest, we consider SDEs involving the integral of higher derivatives. One of the advantage of this technique is that it bypasses the reformulation of the problem in terms of singular values. As a byproduct of this study we obtain the large $N$ limit of the Stieltjes transform of the $2$-point correlation function, as well as the first correction to the Stieltjes transform of the density, giving us access to corrections to the smoothed density. In order to pave the way for the establishment of a topological recursion formula we also study the geometry of the corresponding spectral curve. This paper also contains explicit results for different resolvents and their corrections.

Tables

Table 1. Table 1: Table of the first few cumulants c i , j [ 0 ] = lim N → ∞ 𝔼 ​ ( Tr ​ ( S 2 i ) ​ Tr ​ ( S 2 j ) ) − 1 N 2 ​ 𝔼 ​ ( Tr ​ ( S 2 i ) ) ​ 𝔼 ​ ( Tr ​ ( S 2 j ) ) subscript superscript 𝑐 delimited-[] 0 𝑖 𝑗 subscript → 𝑁 𝔼 Tr superscript subscript 𝑆 2 𝑖 Tr superscript subscript 𝑆 2 𝑗 1 superscript 𝑁 2 𝔼 Tr superscript subscript 𝑆 2 𝑖 𝔼 Tr superscript subscript 𝑆 2 𝑗 c^{[0]}_{i,j}=\lim_{N\rightarrow\infty}\mathbb{E}\left(\mathrm{Tr}(S_{2}^{i})\mathrm{Tr}(S_{2}^{j})\right)-\frac{1}{N^{2}}\mathbb{E}\bigl{(}\mathrm{Tr}(S_{2}^{i})\bigr{)}\mathbb{E}\bigl{(}\mathrm{Tr}(S_{2}^{j})\bigr{)} .

	1	2	3	4	5	6	7
1	3	20	126	792	5005	31824	203490
2	**	150	1008	6600	42900	278460	1808800
3	**	**	7056	47520	315315	2079168	13674528
4	**	**	**	326700	2202200	14702688	97675200
5	**	**	**	**	15030015	101359440	678978300
6	**	**	**	**	**	689244192	4649339520
7	**	**	**	**	**	**	31549089600

Equations

$$\prod _{l=1}^N{e}^{-V({\lambda }_l)}\prod _{1\le i<j\le N}({\lambda }_i-{\lambda }_j)({\lambda }_i^{\theta }-{\lambda }_j^{\theta }),$$

$$\frac{N}{2i\pi }{e}^{-N|X_{i,j}{|}^2}\mathrm{d}{\bar{X}}_{i,j}\mathrm{d}{X}_{i,j}.$$

$$\mathrm{d}{X}^{†}\mathrm{d}X=\prod _{i,j}\mathrm{d}{\bar{X}}_{i,j}\mathrm{d}{X}_{i,j},$$

$$\mathrm{d}\mu (X)=\frac{N^{N^2}}{(2i\pi {)}^{N^2}}{e}^{-N\mathrm{Tr}(X{X}^{†})}\mathrm{d}{X}^{†}\mathrm{d}X.$$

$$m_k=\mathbb{E}\left(\mathrm{Tr}(S_1^k)\right).$$

$$m_{k_1,\dots ,k_n}=\mathbb{E}\left(\prod _{i=1}^n\mathrm{Tr}(S_1^{k_i})\right).$$

$$m_k=\sum _{g\ge 0}{N}^{1-2g}{m}_k^{[g]},$$

$$m_k=\sum _{\mathcal{M}\in {\mathbb{M}}_k}{N}^{V_{\circ }(\mathcal{M})-k+F(\mathcal{M})},$$

$${\mathbb{M}}_{k_1,\dots ,k_n}=\{\mathcal{M}=(E,{\sigma }_{\bullet },{\sigma }_{\circ })\mid {\sigma }_{\bullet }={\gamma }_{k_1,\dots ,k_n}\}.$$

$${\mathbb{M}}_{k_1,\dots ,k_n}^c=\{\mathcal{M}=(E,{\sigma }_{\bullet },{\sigma }_{\circ })\mid {\sigma }_{\bullet }={\gamma }_{k_1,\dots ,k_n},⟨{\sigma }_{\bullet },{\sigma }_{\circ }⟩\text{ acts transitively on }E\}.$$

$$m_{k_1,\dots ,k_n}=\sum _{\mathcal{M}\in {\mathbb{M}}_{k_1,\dots ,k_n}}{N}^{V_{\circ }(\mathcal{M})-p+F(\mathcal{M})},$$

$$m_{k_1,\dots ,k_n}=\sum _{K⊢\{k_1,\dots ,k_n\}}\prod _{{\kappa }_i\in K}{c}_{{\kappa }_i}.$$

$$c_{k_1,\dots ,k_n}=\sum _{\mathcal{M}\in {\mathbb{M}}_{k_1,\dots ,k_n}^c}{N}^{V_{\circ }(\mathcal{M})-p+F(\mathcal{M})}.$$

$$c_{k_1,\dots ,k_n}=\sum _{g\ge 0}{N}^{2-n-2g}{c}_{k_1,\dots ,k_n}^{[g]}.$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}{m}_k=m_k^{[0]}.$$

$$W_{0,1}(x):=\underset{N\to \infty }{\lim }\frac{1}{N}\mathbb{E}\left(\mathrm{Tr}\left((x-W{)}^{-1}\right)\right)=\sum _{p\ge 0}\frac{m_p^{[0]}}{x^{p+1}}=\frac{x-\sqrt{x^2-4x}}{2x}.$$

$$\sum _{a,b=1}^N\int \frac{N^{N^2}}{(2i\pi {)}^{N^2}}\mathrm{d}{X}^{†}\mathrm{d}X{\partial }_{X_{ab}^{†}}\left((X^{†}{S}_1^k{)}_{ab}{e}^{-N\mathrm{Tr}(X{X}^{†})}\right)=0,$$

$$\sum _{\begin{array}{c}{p}_1,p_2\ge 0\\ p_1+p_2=k\end{array}}{m}_{p_1,p_2}-N{m}_{k+1}=0.$$

$${\overline{W}}_n(x_1,\dots ,x_n)$$

$$W_n(x_1,\dots ,x_n)$$

$${\overline{W}}_n(x_1,\dots ,x_n)=\sum _{K⊢\{1,\dots ,n\}}\prod _{K_i\in K}{W}_{\mid K_i\mid }(x_{K_i}),$$

$$\sum _{k\ge 0}\frac{1}{x^{k+1}}\left(\sum _{\begin{array}{c}{p}_1,p_2\ge 0\\ p_1+p_2=k\end{array}}{m}_{p_1,p_2}-N{m}_{k+1}\right)=0,$$

$${\overline{W}}_2(x,x)-N{W}_1(x)+N^2/x=0,$$

$$W_1(x{)}^2+W_2(x,x)-N{W}_1(x)+N^2/x=0.$$

$$W_n(x_1,x_2,\dots ,x_n)=\sum _{g\ge 0}{N}^{2-2g-n}{W}_{g,n}(x_1,x_2,\dots ,x_n)$$

$$W_1(x)=\sum _{g\ge 0}{N}^{1-2g}{W}_{g,1}(x),W_2(x,x)=\sum _{g\ge 0}{N}^{-2g}{W}_{g,2}(x,x).$$

$$x{W}_{0,1}(x{)}^2-x{W}_{0,1}(x)+1=0.$$

$$\delta f(x)=\underset{ϵ\to 0^{+}}{\lim }f(x+iϵ)-f(x-iϵ)$$

$$sf(x)=\underset{ϵ\to 0^{+}}{\lim }f(x+iϵ)+f(x-iϵ).$$

$$\delta (f_1{f}_2)(x)=\frac{1}{2}(\delta f_1(x)s{f}_2(x)+s{f}_1(x)\delta f_2(x))$$

$$s(f_1{f}_2)(x)=\frac{1}{2}(\delta f_1(x)\delta f_2(x)+s{f}_1(x)s{f}_2(x))$$

$$\delta (x{W}_{0,1}(x{)}^2-x{W}_{0,1}(x)+1)=0$$

Full text

Keywords: random matrices, product of Ginibre matrices, loop equations, Schwinger-Dyson equations, spectral curve

1 Introduction

The study of random matrices in mathematics can be traced back to the work of Hurwitz on the invariant measure for the matrix groups $U(N)$ and $SO(N)$ [Hur97, DF17]. In multivariate statistics another stream of random matrix theory was initiated with the work of Wishart [Wis28] on estimating the covariance matrices of multivariate statistics when the number of variables is large. In theoretical physics Wigner [Wig55] used random matrices to model energy spectrum of Hamiltonians of highly excited states of heavy nuclei. The works of physicists [tH74] on the large $N$ limit of $U(N)$ gauge theory provided yet another application to random matrices (and their generalized version often referred to as matrix models). Since then random matrix theory and matrix models have been found useful in an overwhelming number of contemporary fields, for example communication engineering [TV04], the analysis of algorithms [Tro15], and deep learning [PW17]. Many tools have been developed to understand the properties of different models and ensembles. One of these tools is called loop equations, and has led to the now well-known Chekhov-Eynard-Orantin topological recursion formula [Eyn04, CEO06, CE06]. In the realm of random matrix theory this formula allows for the systematic computation of correlation functions of random matrices, as series in $1/N$.

However some random matrix ensembles are, in the existing literature, still out of the scope of these loop equations. These are product ensembles, that is they are random matrices constructed out of a product of several random matrices. In this paper we describe the loop equations for such a product ensemble, specifically considering the case of a random matrix constructed out of the product of two complex Ginibre matrices. Such an ensemble was for instance considered in [BLMP07], with applications to the study of financial data, while a closely related product ensemble with applications to low energy QCD, was studied in [Osb04] (see also the text book treatment [For10, §15.11]), allowing for insight into the poorly understood regime of non-zero baryon chemical potential.

More generally the product ensembles are found to have many applications. Some of these applications are described in the thesis [Ips15]. Among those, one finds applications to telecommunication problems where product ensembles provide a model of communication channels where the signal has to pass through different media [Mul02]. One also finds applications to the study of spin chains with disorder [CPV93], quantum transport [Bee97], quantum information and random graph states [CNŻ10, CNŻ13]. The product ensembles also relate to the study of neural networks. Indeed information about the asymptotic behavior of such ensembles allows one to draw results about stability of gradient in a deep neural network with randomly initialized layers [HN18]. These product ensembles are also of interest for the study of the stability of large dynamical systems [Ben84, IF18]. As a consequence, finding mathematical and technical tools for investigating the properties of these ensembles can enable progress in these fields of study.

Yet another problem of importance is the one of Muttalib-Borodin ensembles. These ensembles were first defined as invariant ensembles, via their eigenvalue probability density function (PDF) [Mut95], and latter realized in terms of ensembles of random matrices with independent entries [Che18, FW17]. Their joint PDF is proportional to,

[TABLE]

where $\theta>0$ is a parameter and $V(\lambda_{l})$ can be interpreted as a confining potential. For general potential $V$ and $\theta=2$, this model relates to the $\mathcal{O}(\mathfrak{n})$ matrix model with $\mathfrak{n}=-2$, see [BE11], and it also relates to a particular model of disordered bosons [LSZ06]. A key structural interest in the Muttalib-Borodin ensembles is that they are biorthogonal ensembles. That is they admit a family of biorthogonal polynomials and their correlation functions can be expressed in determinantal form, with a kernel that can be expressed in terms of the biorthogonal polynomials; see [Bor98]. Although it is not immediately obvious, the singular values for the product of $M$ complex Ginibre matrices also give rise to biorthogonal ensembles [AIK13, KZ14]. Moreover, in the asymptotic regime of large separation, the PDF for the squared singular values reduces to (1) with $\theta=1/M$, and $V$ having the leading form $V(x)=-Mx^{1/M}$ [FLZJ15].

One attractive feature of both the Muttalib-Borodin ensemble, and the squared singular values of products of complex Ginibre matrices, is that in the global density limit the moments of spectral density are given by the Fuss-Catalan family of combinatorial numbers; see [PŻ11, FW17]. Another is the special role played by particular special functions of the Meijer-G and Wright Bessel function class. Underlying these special functions is a linear differential equation of degree $M+1$. Less well understood is the nonlinear differential system implied by the correlation kernel based on these special functions. These are relevant to the study of gap probabilities; see [WF17, MF18].

Other questions about products of random matrices have been investigated for instance in [DP19]. In this work, the authors are concerned about the behavior of traces of general words of Ginibre matrices. In particular they show that the limiting square singular values distribution is a Fuss-Catalan distribution for any words. In the work [DLN18], the authors study the traces of the general words in an alphabet of random matrices constructed out of the marginals of a random tensor. Using combinatorial techniques it is possible to show freeness of some marginals or to describe entirely the free cumulants when there is no freeness of the different marginals in the limit. One interesting aspect is that using these products of marginals it is possible to find distribution interpolating between the square of a Marc̆enko-Pastur law and the free multiplicative square of a Marc̆enko-Pastur law.

However there is in general little technical tools to describe the lower order in $N$ observables of product ensembles. Indeed free probability provides us with some useful techniques (free additive and multiplicative convolution), but those are restricted to the large $N$ limit, and comes in handy only for the study of the large $N$ density or the behavior of the large $N$ limit of the moments (with some extension to the fluctuations of the linear statistics via [CMSS07]).

In this paper we focus on describing the loop equations for the random matrix $S_{2}=X_{1}X_{1}^{\dagger}X_{2}^{\dagger}X_{2}$, where $X_{1},X_{2}$ are square complex Ginibre matrices. In order to obtain these loop equations we start with Schwinger-Dyson identities and use them to obtain relations between moments, later translated in terms of equations on the resolvents of $S_{2}$. These equations on the resolvents are the loop equations. One of the new features of the method presented here is that the starting point Schwinger-Dyson identities involve higher order derivatives. This allows us to obtain relations between moments of the matrix $S_{2}$ only without having to deal with mixed quantities. Thanks to the combinatorial interpretation of the moments of the matrix $S_{2}$ (that we also shortly describe), we show that the (connected) resolvents possess a $1/N$ expansion, which is the unique additional ingredient we need to be able to solve the loop equations recursively.

Using this data we illustrate the use of the obtained loop equations by computing the large $N$ limit of the resolvent $W_{0,1}(x)$, thus recovering known results relating to the generating function of the moments. We also compute $W_{0,2}(x_{1},x_{2})$ (that is the Stieltjes transform of the $2$-point correlation function) and show that it takes the expected universal form once expressed in the correct variables, thus relating to the Bergmann kernel on the sphere. We give explicit results for $W_{1,1}(x),W_{2,1}(x)$ (first and second correction to the large $N$ limit of the resolvent), $W_{1,2}(x_{1},x_{2})$ (first correction to $W_{0,2}(x_{1},x_{2})$), as well as $W_{0,3}(x_{1},x_{2},x_{3})$. One interesting aspect of the obtained loop equations are their structural properties, that seem to generalize in a very natural way the usual bilinear loop equations for random matrices or matrix models. In particular, the family of loop equations we obtain for this product of matrices are trilinear in the resolvents $W_{g,n}$. This is at the root of the appearance of the double ramification point of $W_{0,1}(x)$ and we expect that a topological recursion formula similar to the one obtained in [BE13] applies. Moreover they contain generalizations of the derivative difference term usually appearing in the bilinear setting, as well as derivatives of first and second order. Motivated by these interesting structural properties, we use the explicit computations to explore the analytical properties of the $W_{g,n}$ (or rather their analytic continuation on the associated spectral curve). These explorations give further hint that there is a topological recursion formula to compute them systematically. We expect that a similar technique allows to describe the loop equations for the product of $p\geq 2$ rectangular Ginibre matrices $S_{p}=X_{1}X_{2}\ldots X_{p}(X_{1}X_{2}\ldots X_{p})^{\dagger}$; we leave this study, as well as the one of a topological recursion formula, to further works. Note that, as a byproduct, we also expect that this technique applies to the interesting matrix models introduced in [AC14, AC18] to generate hypergeometric Hurwitz numbers.

Organisation of the paper.

The paper is organized as follows. In section 2, we use the Wishart case (that is the case of one Ginibre matrix) as a pedagogical example. It is used to sketch the combinatorial arguments allowing to show the existence of the $1/N$ expansion and to illustrate the Schwinger-Dyson equation technique in a simpler context. The reader already accustomed to Schwinger-Dyson equations obtained using the matrix elements variables and knowledgeable on the associated combinatorics may consider skipping this section.

In section 3, we describe the heart of this paper, that is the derivation of the Schwinger-Dyson equations and loop equations for a product matrix of the form $S_{2}=X_{1}X_{1}^{\dagger}X_{2}^{\dagger}X_{2}$. The loop equations take the form of a family of equations on the resolvents, that is the Stieltjes transforms (denoted $W_{n}(x_{1},\ldots,x_{n})$) of the $n$-point correlation functions. We present the results step by step to make the method transparent to the reader and the first few special cases that are the loop equations for $W_{1}(x)$, $W_{2}(x_{1},x_{2})$ and $W_{3}(x_{1},x_{2},x_{3})$ are presented in details. This section ends with the main result, that is the loop equations satisfied by any $W_{g,n}(x_{1},\ldots,x_{n})$ as shown on equation (96), where $W_{g,n}(x_{1},\ldots,x_{n})$ is the coefficient of order $g$ of the $1/N$ expansion of $W_{n}(x_{1},\ldots,x_{n})$.

In section 4, we take on a geometrical point of view in order to compute the $W_{g,n}$ more effectively from the loop equations. We describe in details the spectral curve geometry associated to the problem. We compute after a change of variables, $W_{0,2}(x_{1},x_{2})$, $W_{1,1}(x)$, $W_{2,1}(x)$, $W_{1,2}(x_{1},x_{2})$ and $W_{0,3}(x_{1},x_{2},x_{3})$ (see equations (141), (139), (140), (142), (144)). We use these explicit computations to explore the analytic properties of the loop equations. These properties are expected to be of importance to establish a topological recursion formula allowing to systematically compute every $W_{g,n}$.

Acknowledgments

Stephane Dartois would like to thank Valentin Bonzom, Alexandr Garbali, Jesper Ipsen and Paul Zinn-Justin for useful discussions and technical help related to this work as well as for references. This work was supported by the Australian Research Council grant DP170102028.

2 One matrix case, Wishart ensemble

In this section, we illustrate the problem that is our interest in this paper on a simpler case, that is the (trivial) product of one matrix. This is the case of a Wishart matrix. We first recall the combinatorial representation of moments of a Wishart ensemble matrix. We then show how we can compute the average resolvent of a Wishart matrix using the Schwinger-Dyson equation method. It is only in the next section that we consider the case of the product of two Ginibre matrices. Thus the technically knowledgeable reader can skip this section and start reading section 3.

2.1 Random Wishart matrices

In this paper we always consider square matrices. In the Wishart matrices case it corresponds to setting the asymptotic size ratio parameter $c$ to $1$. Let $X\in\mathcal{M}_{N\times N}(\mathbb{C})$ be a Ginibre random matrix. More concretely, $X$ is a random matrix whose entries are i.i.d. complex Gaussian with zero mean, or more formally, the entries $X_{i,j}$ are distributed according to the density

[TABLE]

In particular we denote,

[TABLE]

so that $X$ has the distribution

[TABLE]

A (complex) Wishart random matrix is the random variable defined as the product $S_{1}=XX^{\dagger}$.

Combinatorics of moments. The moments $m_{k}$ of order $k$ of a Wishart random matrix are defined as

[TABLE]

Further, for any sequence of positive integers $k_{1},\ldots,k_{n}$ we can define moments $m_{k_{1},\ldots,k_{n}}$ of order $k_{1},\ldots,k_{n}$. Similarly to the moments of order $k$ they are defined as the expectation of products of traces of powers of $S_{1}$

[TABLE]

As is for instance explained in [DLN18], the moments of order $k$ can be computed as a sum over labeled bicolored combinatorial maps $\mathcal{M}$ with one black vertex. This combinatorial representation of moments implies that the moments have a $1/N$ expansion. That is

[TABLE]

where $m_{k}^{[g]}$ are the coefficients of this expansion. This is a crucial point that allows one to solve the loop equations recursively. Note also that this expansion is finite, that is here $g<k/2$. Let us be a bit more explicit on this point.

We recall the definition of labeled bicolored combinatorial maps with possibly more than one black vertex.

Definition 1**.**

A labeled bicolored combinatorial map is a triplet $\mathcal{M}=(E,\sigma_{\bullet},\sigma_{\circ})$ where,

• •

$E$* is the set of edges of $\mathcal{M}$*

• •

$\sigma_{\bullet},\sigma_{\circ}$* are permutations on $E$*

• •

$\mathcal{M}$* is said to be connected if and only if the group $\langle\sigma_{\bullet},\sigma_{\circ}\rangle$ acts transitively on $E$.*

The cycles of $\sigma_{\circ}$ are called white vertices, the cycles of $\sigma_{\bullet}$ are called black vertices, and the cycles of $\sigma_{\bullet}\sigma_{\circ}$ are called faces. Combinatorial maps can be represented graphically [DLN18, Eyn16] as they encode embeddings of graphs on surfaces. We give a few examples in Fig. 1.

We define the set of combinatorial maps $\mathbb{M}_{p}=\{\mathcal{M}=(E,\sigma_{\bullet},\sigma_{\circ})\mid E=\{1,\ldots,p\},\sigma_{\bullet}=\gamma=(123\ldots p)\}$. One shows, using Wick-Isserlis theorem [Wic50, Iss18], that the moments of order $k$ can be written as a sum over combinatorial maps $\mathcal{M}\in\mathbb{M}_{p}$ (see [DLN18] for details)

[TABLE]

where $V_{\circ}(\mathcal{M})$ is the number of white vertices of $\mathcal{M}$ and $F(\mathcal{M})$ is the number of faces of $\mathcal{M}$. Using the fact that $V_{\bullet}+V_{\circ}(\mathcal{M})-k+F(\mathcal{M})=2-2g(\mathcal{M})$, where $g(\mathcal{M})$ is the genus of the combinatorial map (that is the genus of the surface in which the corresponding graph embedds), one can show equation (7).

Remark 1**.**

Note that elements of $\mathbb{M}_{p}$ are necessarily connected as $\gamma$ acts transitively of $\{1,\ldots,p\}$.

We now define the relevant set of maps for studying the moments of order $k_{1},\ldots,k_{n}$. In this case we denote $p=\sum_{i=1}^{n}k_{i}$, $E=\{1,\ldots,p\}$ and $\gamma_{k_{1},\ldots,k_{n}}=(12\ldots k_{1})(k_{1}+1\ldots k_{2})\ldots(k_{n-1}+1\ldots k_{n})$

[TABLE]

The maps in $\mathbb{M}_{k_{1},\ldots,k_{n}}$ are possibly non-connected as $\gamma_{k_{1},\ldots,k_{n}}$ does not act transitively on the set of edges. Consequently we define the corresponding set of connected maps

[TABLE]

We state without proof111The proof is very similar to the one black vertex case, already appearing in [DLN18]. that

[TABLE]

where $p=\sum_{i}k_{i}$. We can define the associated cumulants $c_{k_{1},\ldots,k_{n}}$ of the moments, through their relation to moments

[TABLE]

This relation is just the moment-cumulant relation for the family of random variables \bigl{\{}R_{k_{i}}:=\mathrm{Tr}(S_{1}^{k_{i}})\bigr{\}}. These cumulants can be expressed as sums over connected combinatorial maps

[TABLE]

Thanks to the connected condition, this sum is a polynomial in $1/N$ as long as $n>1$. That is to say we have

[TABLE]

This last equation is shown starting from (13) and again using $V_{\bullet}+V_{\circ}(\mathcal{M})-k+F(\mathcal{M})=2-2g(\mathcal{M})$ with $V_{\bullet}=n$.

Large $N$ limit of moments of a Wishart matrix. Using (7), one can study the large $N$ limit of the moments of order $k$ of a Wishart matrix, that is one can compute the limit

[TABLE]

This limit is given by the number of planar, labeled, bicolored combinatorial maps with one black vertex and $k$ edges. The number of such maps is given by the Catalan number222Note that one obtains Catalan numbers when the ratio parameter is set to $c=1$, however for general values of $c$ one obtains the Narayana statistics on trees, that is polynomials in $c$ whose coefficients are Narayana numbers [DR03]. $C_{k}$ so that $m^{[0]}_{k}=C_{k}=\frac{1}{k+1}\binom{2k}{k}$. This allows to compute the large $N$ limit $W_{0,1}(x)$ of the moment generating function of the Wishart matrix

[TABLE]

This last quantity is the Stieltjes transform of the limiting eigenvalues density of the Wishart matrix. The knowledge of $W_{0,1}(x)$ allows in principle333In this specific case one can recover explicitly the limiting eigenvalue density via the inverse transformation. However in general it can be more tedious to compute the inverse transform. In the cases where the equation determining $W_{0,1}$ is an algebraic equation, one can deduce a system of polynomial equations on two quantities $u(x),v(x)$, one of them being (proportional to) the large $N$ limit of the eigenvalue density $\rho_{0,1}(x)$. We illustrate this fact in the later Remarks 3, 4. to recover the limiting eigenvalues density via the inverse transformation.

Schwinger-Dyson equation method. In this part we use an alternative method to compute $W_{0,1}(x)$. We use the Wishart case as a pedagogical example. The Schwinger-Dyson equation method relies on the use of the simple identity

[TABLE]

after computing the derivatives explicitly we obtain the following set of relations between moments

[TABLE]

In order to continue this computation we define the $n$-points resolvents $\overline{W}_{n}(x_{1},\ldots,x_{n})$ and their connected counterpart $W_{n}(x_{1},\ldots,x_{n})$

[TABLE]

Note that we will often name the $n$-points resolvents and their connected counterpart simply resolvents, unless the context makes it unclear which object we are discussing. $W_{0,1}(x)$ is (up to normalization) the large $N$ limit of $W_{1}(x)$. We have the relation

[TABLE]

where we used the notation $x_{K_{i}}=\{x_{j}\}_{j\in K_{i}}$. The above relation is inherited from the moment-cumulant relation of equation (12).

Remark 2**.**

Note that $\overline{W}_{1}(x)=W_{1}(x)$.

With these definitions in mind, one considers the equality

[TABLE]

leading after some rewriting to

[TABLE]

or only in terms of the connected resolvents

[TABLE]

The (connected) resolvents inherit a $1/N$ expansion from the expansion of the cumulants,

[TABLE]

and thus we have

[TABLE]

In the large $N$ limit equation (24) reduces to an equation on $W_{0,1}(x)$,

[TABLE]

From which we select the solution which is analytic at infinity thus recovering expression (16).

Remark 3**.**

From this last equation we can obtain a polynomial equation on $\rho_{0,1}(x)$, that is the corresponding limiting eigenvalue density. To this aim, one introduces the two following operators acting on functions,

[TABLE]

We have the following polarization property, that is for two functions $f_{1},f_{2}$, we have

[TABLE]

Starting from equation (27) one deduces the two equalities

[TABLE]

After using the polarization formula, these equations boil down to the system on $u(x):=sW_{0,1}(x)$ and $v(x):=\delta W_{0,1}(x)$

[TABLE]

This in turn leads to $\rho_{0,1}(x)=\frac{1}{2i\pi}v(x)=\frac{1}{2\pi}\sqrt{\frac{x-4}{x}}$, where we choose the solution $v(x)$ that leads to a positive and normalized density.

3 Loop equations for the product of two Ginibre matrices

In this section we consider the problem of computing $W_{0,1}(x)$, $W_{0,2}(x_{1},x_{2})$ and $W_{1,1}(x)$ for a matrix $S_{2}=X_{1}X_{1}^{\dagger}X_{2}^{\dagger}X_{2}$ with $X_{1},X_{2}$ two random $N\times N$ complex matrices with normal entries of mean zero. We compute these quantities by exclusive use of Schwinger-Dyson equation techniques. More generally, we obtain the general equations satisfied by any $W_{g,n}$ for $(g,n)\geq(0,1)$.

In the first subsection, we briefly explain the combinatorics underlying the computation of the moments of the matrix $S_{2}$ that justifies the existence of a $1/N$ expansion for the $W_{g,n}$. In the second subsection we study in details the corresponding Schwinger-Dyson equations and obtain the loop equations satisfied by $W_{0,1}(x)$, $W_{0,2}(x_{1},x_{2})$ and $W_{1,1}(x)$ in this context. We show in particular that the loop equation satisfied by $W_{0,1}(x)$ is an algebraic equation of degree $3$ in $W_{0,1}$. Finally we describe the loop equations satisfied by any $W_{g,n}$.

3.1 Combinatorics of the moments of $S_{2}$ and existence of $1/N$ expansion

We describe here the combinatorics of the moments of the matrix $S_{2}$. This is a crucial point as this underlying combinatorics allows us to show that the cumulants of the random variables $\left\{\mathrm{Tr}(S_{2}^{i})\right\}_{i=0}^{\infty}$ have a $1/N$ expansion. In the subsequent developments, we keep the same notation for the moments $m_{k}$, $m_{k_{1},\ldots,k_{n}}$ but it should be clear that in this section and the following, the moments we consider are the moments of the matrix $S_{2}$, and that is so, in both the one trace case, and the multiple traces case. We have

[TABLE]

where the expectation is taken with respect to the density

[TABLE]

By using the Wick-Isserlis theorem, it is possible to give a combinatorial interpretation to the moments of $S_{2}$ (see for instance [DLN18]). The moments $m_{k}$ of $S_{2}$ write as a sum over combinatorial maps with one black vertex, $2k$ edges of two different types, type I and type II, such that there are $k$ edges of type I and $k$ edges of type II. Moreover the type of the edge alternates when going around the black vertex. Finally the white vertices can only be incident to edges of one given type. See Fig. 2 for examples.

We denote the set made of these maps by $\mathbb{M}_{2k}(2)$. In terms of permutations, these maps are such that $\sigma_{\bullet}=(12\ldots 2k)$ and the action of $\sigma_{\circ}$ on the set of edges $E=\{1,2,3,4,\ldots,2k\}$ factorizes over the odd and even subsets $E_{o}=\{1,3,5,\ldots,2k-1\},E_{e}=\{2,4,6,\ldots,2k\}$. More formally we have the decomposition

[TABLE]

Similarly, for moments of order $k_{1},\ldots,k_{n}$, we have the set of maps $\mathbb{M}_{2k_{1},2k_{2},\ldots,2k_{n}}(2)$, such that there are $n$ black vertices with degree distribution $2k_{1},2k_{2},\ldots,2k_{n}$ and a total of $p=2\sum_{i}k_{i}$ edges. Types of edge alternate around each black vertex, and white vertices can only be incident to edges of the same type see Fig. 2 for examples. We then have the decomposition

[TABLE]

Similarly we can express the cumulants $c_{k_{1},\ldots,k_{n}}$ for the family of random variables $\left\{\mathrm{Tr}(S_{2}^{i})\right\}_{i=0}^{\infty}$ as a sum over the set of connected maps $\mathbb{M}^{c}_{2k_{1},2k_{2},\ldots,2k_{n}}(2)$

[TABLE]

The connected condition ensures that the $c_{k_{1},\ldots,k_{n}}$ have a $1/N$ expansion for $n\geq 1$. This $1/N$ expansion as well as the definition of $c_{k_{1},\ldots,k_{n}}$ as the cumulants of the family $\left\{\mathrm{Tr}(S_{2}^{i})\right\}_{i=0}^{\infty}$ ensure that the resolvents for the matrix $S_{2}$ have the same structural properties than the resolvents of the Wishart matrix in equations (21), (25), that is we also have for the matrix $S_{2}$

[TABLE]

3.2 Equation on $W_{1}$ and $W_{0,1}$

We now want to write Schwinger-Dyson equations for the moments of the matrix $S_{2}$ in order to obtain the loop equations for the resolvents. We start with the set of identities

[TABLE]

After evaluating explicitly the action of the derivatives, we obtain relations,

[TABLE]

where for both equation, the first term comes from the evaluation of the derivative on the monomial, while the second term comes from the evaluation of the derivative on the exponential factor. Note however that these equations contain mixed terms of the form $\mathbb{E}\left(\mathrm{Tr}(S_{2}^{p_{1}})\mathrm{Tr}(S_{2}^{p_{2}}X_{2}^{*}X_{2})\right)$ and $\mathbb{E}\left(\mathrm{Tr}(S_{2}^{p_{1}}X_{1}X_{1}^{*})\mathrm{Tr}(S_{2}^{p_{2}})\right)$ that cannot be expressed in terms of the moments of $S_{2}$. Thus these two equations do not close on the set of moments of $S_{2}$. In order to obtain a set of relations that closes over the set of moments of $S_{2}$, we consider another identity involving higher derivatives. This is,

[TABLE]

where we sum over repeated indices. After some additional algebra to evaluate the action of both derivative operators, one gets relations between moments and additional mixed quantities

[TABLE]

where the first and second terms are obtained from the action of both derivatives operators on the monomial \bigl{[}X_{1}^{\dagger}X_{2}^{\dagger}X_{2}S_{2}^{k}X_{1}X_{1}^{\dagger}X_{2}^{\dagger}\bigr{]}_{ac}. The third term that involves mixed quantities is obtained by acting with one derivative operator on the monomial, while acting with the other derivative operator on the exponential factor. The last term is obtained from the action of both derivative operator on the exponential factor. These equations contain the mixed quantities already present in (43). Thus we can use (43) to get rid of these terms in (48). This leads to the equations on moments

[TABLE]

which is trilinear in the traces of $S_{2}$. Notice that this family of equations extends to the value “$k=-1$” by replacing the monomial \bigl{[}X_{1}^{\dagger}X_{2}^{\dagger}X_{2}S_{2}^{k}X_{1}X_{1}^{\dagger}X_{2}^{\dagger}\bigr{]}_{ac} by \bigl{[}X_{1}^{\dagger}X_{2}^{\dagger}\bigr{]}_{ac}. Therefore we allow ourselves to set $k=k-1$ and to use our moments notation to get

[TABLE]

We then multiply the above equation by $\frac{1}{x^{k+1}}$ and sum over $k\geq 0$ in order to get an equation on the resolvents

[TABLE]

which after a few manipulations rewrites

[TABLE]

Note the interesting structural replacement of $\overline{W}_{2}(x,x)$ appearing in (23) by $\overline{W}_{3}(x,x,x)$ and the appearance of a derivative term. Then we know from (41), (42) that $\overline{W}_{3}(x,x,x)=N^{3}W_{0,1}(x)^{3}+O(N)$ and $W_{1}(x)=NW_{0,1}(x)+O(1/N)$. Therefore we obtain the equation on $W_{0,1}(x)$

[TABLE]

This last equation relates to the equation satisfied by the generating function $G(u)$ of particular Fuss-Catalan numbers [Fus91, Mlo10, Riv18], $uG(u)^{3}-G(u)+1=0$ through the change of variables $W_{0,1}(x)=\frac{1}{x}G(1/x)$. Consequently we have

[TABLE]

where $C_{p}[D]$ are the Fuss-Catalan numbers of order $D$, the usual Catalan numbers $C_{p}$ being the Fuss-Catalan numbers of order $2$, that is $C_{p}=C_{p}[2]$, and have the binomial coefficient form

[TABLE]

An explicit form of $W_{0,1}(x)$ can be written as follows. First define

[TABLE]

then $G(u)$ writes

[TABLE]

Finally one has

[TABLE]

We study the solutions and the structure of (53) from a geometric perspective in the next sections.

Remark 4**.**

Though in principle we need to first focus on the cut structure of $W_{0,1}$ to use the arguments that follow, we will in this remark content ourselves with a formal computation. Starting from equation (53) we can also obtain a polynomial equation satisfied by the corresponding density by using the $\delta,s$ operators along the cut. Indeed with a similar method to that in Remark 3 we have the equalities

[TABLE]

This leads, using the same previously used notations, to the system

[TABLE]

which can be solved and leads to

[TABLE]

which is supported on $(0,27/4]$, see the plot of the distribution on Fig. 3. Notice that this result can also be obtained by computing the free multiplicative product of two Marc̆enko-Pastur distribution of parameters $c_{1,2}=1$. A functional form equivalent to (63) is given in [PŻ11].

Equation (52) possesses a $\frac{1}{N}$ expansion. This expansion results in a set of relations between $W_{g,1}(x)$, $W_{g^{\prime},2}(x,x)$ and $W_{g^{\prime\prime},3}(x,x,x)$. Indeed we have

[TABLE]

By collecting the coefficient of $N^{3-2g}$, we obtain the following tower of equations

[TABLE]

where we have $P_{g,1}(x)=\delta_{g,0}$. In particular, the coefficient of $N^{3}$ of equation (65) produces equation (53). The coefficient of $N$ produces an equation on the next-to-leading order $W_{1,1}(x)$ also involving $W_{0,1}(x)$ and $W_{0,2}(x,x)$

[TABLE]

More generally, the coefficient of $N^{3-2g}$ for a fixed value of $g$ produces the equation for $W_{g,1}(x)$ in terms of the functions $W_{g^{\prime},n^{\prime}}$ such that $2-2g-1<2-2g^{\prime}-n^{\prime}$ and $n^{\prime}\leq 3$.

3.3 Equation for $W_{2}(x_{1},x_{2})$

In this section we use Schwinger-Dyson equation techniques to obtain a loop equation for $W_{2}(x_{1},x_{2})$. We start with slightly different identities that involve an additional trace insertion $\mathrm{Tr}(S_{2}^{q})$. This allows us to access relations between more general moments.

Schwinger-Dyson equations and loop equation for $W_{2}(x_{1},x_{2})$ and $W_{0,2}(x_{1},x_{2})$. Consider the vanishing integrals of total derivatives

[TABLE]

and the higher derivative one

[TABLE]

where repeated indices are summed. After evaluating explicitly the derivatives, the two first equations (67) and (68) lead to

[TABLE]

where the first term of both equations (67) and (68) is obtained from the action of the derivative operator on the non-traced monomial. The second term is obtained via the action of the derivative operator on the traced monomial term $\mathrm{Tr}(S_{2}^{q})$. The third term comes from the action of the derivative operator on the exponential factor. These two equations involve mixed terms and cannot be written solely in terms of the moments of $S_{2}$. Meanwhile, the higher derivative equation (69) leads to

[TABLE]

where the two first terms come from the action of both derivatives operators on the non-traced monomial. Each term of the second line comes from the action of one of the derivative on the exponential factor and of the other on the non-traced monomial. The first term of the third line of (72) comes from the action of both derivatives on the exponential factor. The second term of the third line is obtained as a sum of the action of the $X_{1}^{\dagger}$ (resp. $X_{2}^{\dagger}$) derivative on the non-traced monomial and the action of the $X_{2}^{\dagger}$ (resp. $X_{1}^{\dagger}$) derivative on the traced monomial $\mathrm{Tr}(S_{2}^{q})$. The last term of the third line is obtained from the action of both derivative operators on the traced monomial. Finally the two terms of the fourth line of (72) are obtained by the action of $\partial_{X_{1,ab}^{\dagger}}$ (resp. $\partial_{X_{2,bc}^{\dagger}}$) on the traced monomial and $\partial_{X_{2,bc}^{\dagger}}$ (resp. $\partial_{X_{1,ab}^{\dagger}}$) on the exponential factor. Combining equations (70), (71) and (72), rewriting some of the sums in a nicer way and using our moments notation we obtain

[TABLE]

After performing the shift $k\rightarrow k-1$ in (73), we multiply (73) by $\frac{1}{x_{1}^{k+1}x_{2}^{q+1}}$, and sum over $k,q\geq 0$. Doing so we obtain the equation

[TABLE]

with $A_{2}(x_{1},x_{2})=-\frac{N}{x_{1}^{2}}W_{1}(x_{2})$. We re-express this equation in terms of the connected resolvents to obtain

[TABLE]

where we used the fact that the terms factoring in front of $W_{1}(x_{2})$ form the first loop equation (52). From this equation we can get an equation on $W_{0,2}$ by inserting the $1/N$ expansion of the resolvents appearing in (76) and collecting the coefficients of $N^{2}$. This equation involves only already computed quantities and can be re-expressed as

[TABLE]

First few relations between $c^{[0]}_{k},c^{[0]}_{k_{1},k_{2}}$. One can extract relations between the $c^{[0]}_{k},c^{[0]}_{k_{1},k_{2}}$ from equation (3.3). These relations are obtained by expanding the equation at $x_{1},x_{2}=\infty$. The first few examples are

[TABLE]

These relations allow to obtain the $c^{[0]}_{k_{1},k_{2}}$ recursively knowing that $c^{[0]}_{0},\ c^{[0]}_{1}=1$. We can check these first few relations combinatorially. For illustrative purposes we display the combinatorial maps interpretation of $3c^{[0]}_{0}c^{[0]}_{1}-c^{[0]}_{1,1}=0$

[TABLE]

More generally, one has

[TABLE]

3.4 General loop equations

In this section we describe the general loop equations for $W_{n}(x_{1},\ldots,x_{n})$. Because of the use of higher derivatives for Schwinger-Dyson equations, the case of $W_{3}(x_{1},x_{2},x_{3})$ is still special compared to the cases $W_{n<3}$. We thus give the corresponding Schwinger-Dyson equations in details before stating the corresponding loop equations. For the $W_{n>3}$ cases, the situation is very similar to the $W_{3}$ case. Therefore we refrain from presenting the detailed derivation, and only state the corresponding loop equations.

Loop and Schwinger-Dyson equations for $W_{3}(x_{1},x_{2},x_{3})$. We have to consider the equalities,

[TABLE]

The inspection of these Schwinger-Dyson equations reveals that the only type of terms that we have not already faced are obtained when both derivatives $\partial_{X_{1,ab}^{\dagger}}$, $\partial_{X_{2,bc}^{\dagger}}$ distribute over the two traced monomials $\mathrm{Tr}(S_{2}^{q_{1}})$, $\mathrm{Tr}(S_{2}^{q_{2}})$. The distributed action of derivatives on the traced monomial leads to the term

[TABLE]

The generating function of this term appearing in the corresponding loop equation will be

[TABLE]

where $\Delta(\{x_{1},x_{2},x_{3}\})=(x_{3}-x_{2})(x_{3}-x_{1})(x_{2}-x_{1})$ is the Vandermonde determinant of the family of variables $\{x_{1},x_{2},x_{3}\}$. The remaining terms of the loop equations can be inferred by realizing that for all terms involved in either (83), (84), (85), one of the two traced monomials plays a spectator role for the action of the derivatives. Consequently, one obtains the loop equation,

[TABLE]

where we have set $A_{3}(x_{1},x_{2},x_{3})=-\frac{N}{x_{1}^{2}}\overline{W}_{2}(x_{2},x_{3})$.

We now introduce some notations in order to shorten expressions. We denote

[TABLE]

The notation $\mu\vdash[x_{1},x_{1},x_{1}]$ needs to be explained. The summation runs over the partitions $\mu$ of the list $[x_{1},x_{1},x_{1}]$ in the following sense. Firstly, in our notation the object $[x_{a},x_{b},x_{c},\ldots]$ is a list of elements, that is an ordered multi-set. More concretely the order of appearance of the elements in the list is important and so for example the instances $[x_{1},x_{2},x_{1},x_{1},x_{4}]$, $[x_{1},x_{1},x_{1},x_{2},x_{4}]$ of lists are different (though they are the same multi-sets). We now come to explain what we mean by partitions of lists. A (denumerable444we will of course consider only the denumerable case since our lists are finite.) list of elements can be represented as a set in the following way. We send a list to the set of pairs $\{(\textrm{element},\textrm{position in the list})\}$. For instance, the list $[x_{1},x_{2},x_{1},x_{1},x_{4}]\mapsto\{(x_{1},1),(x_{2},2),(x_{1},3),(x_{1},4),(x_{4},5)\}$ while the second list $[x_{1},x_{1},x_{1},x_{2},x_{4}]\mapsto\{(x_{1},1),(x_{1},2),(x_{1},3),(x_{2},4),(x_{4},5)\}$ which are indeed two different sets. The partitions of the list $\mu$ are the partitions of the corresponding set of pairs $(\textrm{element},\textrm{position in the list})$. However, note that the elements of the partitions forget about the position in the list and thanks to the symmetry of the functions $W_{n}$ functions should be seen as subsets of the corresponding multi-set. For instance, due to the fact that $\mu$ is really a partition of a list, the partition $\mu=\{\{x_{1},x_{1}\},\{x_{1}\}\}$ with $\mu_{1}=\{x_{1},x_{1}\},\ \mu_{2}=\{x_{1}\}$ of the list $[x_{1},x_{1},x_{1}]$ appears three times in the sum.

Some further notations are also required. The sum over $\bigsqcup_{i=1}^{|\mu|}J_{i}=\{x_{2},\ldots,x_{n}\}$ means that we sum over the decompositions into $|\mu|$ (possibly empty) subsets $J_{i}$ of the set $\{x_{2},\ldots,x_{n}\}$. For instance, in the case $n=3$, one can consider the term indexed by the partition $\mu=\{\{x_{1},x_{1}\},\{x_{1}\}\}$ and the decomposition $J_{1}=\emptyset,\ J_{2}=\{x_{2},x_{3}\}$, which correspond to a term of the form $W_{2}(x_{1},x_{1})W_{3}(x_{1},x_{2},x_{3})$ in the sum. Note that these definitions are very similar to the ones appearing in [BE13, Definition 4]. We also introduce the notation

[TABLE]

Using these notations the corresponding equation for connected resolvents writes

[TABLE]

We can now extract the corresponding equation of order $g$ (that is the coefficient of $N^{-1-2g}$ in the expansion of (3.4)). The corresponding family of equations on $W_{g,3}$ can then be solved recursively provided that we know the $W_{g^{\prime},n^{\prime}}$ of lower orders,

[TABLE]

We now state in full generality the loop equations.

General loop equations. We obtain the higher order loop equations in full generality by starting with Schwinger-Dyson equalities of the same type than (83), (84), (85), but we now insert more traces of monomials of the matrix $S_{2}$. Doing so we obtain more relations between moments, and those relations can be translated into relations involving $W_{n}$ with higher values of $n$. As before, this first set of relations cannot be used to compute the $W_{n}$ as it does not close. To solve this problem we perform the $1/N$ expansion which leads to a closed set of equations on $W_{g,n}$. We display both the equations on $W_{n}$ and the equations on $W_{g,n}$ for $(g,n)$ such that $2g-2+n>0$. With $I_{ij}=\{x_{1},\ldots,x_{n}\}\backslash\{x_{i},x_{j}\}$,

[TABLE]

For the equations on $W_{g,n}$, write

[TABLE]

Using the family of equations (96) one can recursively compute any $W_{g,n}$ knowing the initial conditions $W_{0,1}(x)$ and $W_{0,2}(x_{1},x_{2})$. Moreover, starting from these equations it should be possible to obtain a topological recursion like formula. Such a recursion formula certainly looks like the Bouchard-Eynard topological recursion formula introduced in [BHL*+*14, BE13]. Establishing such a formula strongly depends on the analytic properties of the $W_{g,n}$ as well as the geometric information contained in $W_{0,1}$ and $W_{0,2}$. Thus in the next section we try to make explicit some of these properties. We first focus on the geometry underlying the equation satisfied by $W_{0,1}$, and then describe the analytic properties of the higher order terms, by: 1. doing explicit computations and 2. studying the structure of the loop equations. A more detailed and systematic study of the analytical properties of the loop equations is postponed to further work on the product of $p$ rectangular Ginibre matrices.

4 Spectral curve geometry

Before computing the first few solutions of the loop equations, we focus on studying the equation (53) on $W_{0,1}$. Indeed, this equation defines an affine algebraic curve $\mathcal{C}$, called the spectral curve, where by affine algebraic curve we mean the locus of zero in $(x,y)\in\hat{\mathbb{C}}^{2}=\left(\mathbb{C}\cup\{\infty\}\right)^{2}$ of the polynomial

[TABLE]

This set of zeros of $P$ in $\mathbb{C}^{2}$ is generically a (complex) codimension $1$ subset of $\mathbb{C}^{2}$. In particular it can be given the structure of a Riemann surface. Computing the solutions $W_{0,1}(x)$ of (53) gives a parametrization of the curve away from the ramification points. One of the goals of this section is to introduce a global, nicer parametrization called rational parametrization of the curve. Using this parametrization allows us to simplify the resulting expressions of the solutions. Indeed in the original $x$ variables, the solutions of (53) are multi-valued. However one can fix that by promoting these solutions to meromorphic functions on the full affine curve defined by equation (53), the curve being the Riemann surface of $W_{0,1}(x)$.

4.1 Basic properties of the curve

There are two finite ramification points in the $x$-plane, one at $(x_{r_{1}},y_{r_{1}})=(27/4,2/9)$, which is a simple ramification point and one at $(x_{r_{2}},y_{r_{2}})=(0,\infty)$ which is a double ramification point. There is also one ramification point at infinity $x_{r_{\infty}}=\infty$ which is a simple ramification point. These ramifications are found from the condition that $P(x,y)=0$ and $\partial_{y}P(x,y)=0$. We display the ramification profile in Fig. 4.

The cut structure is readily described in [FLZJ15, Section 2.1 & 2.2]. It is pictured in Fig. 5, where the lowest sheet of the figure corresponds to the physical sheet that is corresponding to the solution analytic at infinity, whose coefficients of the Laurent expansion are the moments of $S_{2}$. The other two sheets correspond to the two other solutions of (53) that are not analytic at infinity. Indeed they have a simple ramification point at infinity. From the Fig. 5 we can infer that the monodromy group is generated by the transposition $\tau_{1}=(12)$ (obtained by going around $x_{r_{1}}$ in the physical sheet) and $\tau_{2}=(132)$ (going around $x_{r_{2}}$). These permutations are represented using colors on Fig. 4.

The genus of the curve $\mathcal{C}$ can be obtained by considering the Newton polygon of the curve. The number of interior lattice points of the polygon drawn on Fig. 6 corresponds to the generic genus of the curve, that is the genus of the curve for generic enough coefficients of the polynomial $P$. However by fine tuning the coefficients of the polynomial one could in principle obtain a curve with smaller genus. The generic genus is the maximal genus the curve can have. In our case, $P(x,y)=x^{2}y^{3}-xy+1$, the number of lattice points in the Newton polygon is zero, thus the genus of the curve is zero. Since the genus of the curve is zero, there exists a rational parametrization. That is there exists two rational functions

[TABLE]

such that

[TABLE]

These two functions can be found by solving the following system on the coefficients of $Q_{x}(z),Q_{y}(z)$ and $P_{x}(z),P_{y}(z)$,

[TABLE]

where $Q_{x}(z),Q_{y}(z)$ and $P_{x}(z),P_{y}(z)$ are set to be polynomials of degree high enough for a solution to exist. Then one obtains explicitly one possible parametrization

[TABLE]

Note that from this point of view, $y(z)$ is the analytic continuation of $W_{0,1}(x(z))$. The function $x$ can be seen as a cover $x:\mathcal{C}\rightarrow\hat{\mathbb{C}}$ of generic degree $3$ (that is there are generically three values of $z$ corresponding to the same value of $x$). As such, the zeroes of $\mathrm{d}x$ corrrespond to the ramifications point of the cover. One can then check that $\mathrm{d}x=0$ at $z_{r_{1}}=0$ and $z_{r_{2}}=-3/2$, corresponding to the values $x(0)=0$ and $x(-3/2)=27/4$. One also notices that the zero of $\mathrm{d}x$ at $z=0$ is a double zero, thus confirming the fact that $x_{r_{1}}$ is a double ramification point. Finally since $x=27/4$ is a simple ramification point, there is another pre-image of $27/4$ in $z$ variable, that is we have $x(3)=27/4$. This leads to the ramification profile shown on Fig.4.

4.2 Computation of $w_{0,1}$ and $w_{0,2}$

Using this parametrization we compute the functions

[TABLE]

We also denote $\tilde{w}_{0,2}(z_{1},z_{2})=W_{0,2}(x(z_{1}),x(z_{2}))x^{\prime}(z_{1})x^{\prime}(z_{2})$. $w_{g,n}$ functions are meromorphic functions on $\mathcal{C}$, as such they are rational functions of their variables $z_{i}$. Consequently, they are much easier to manipulate than $W_{g,n}$ and their analytic properties are more transparent. For $w_{0,1}(z)$ we already know that $y(z)=W_{0,1}(x(z))$, thus

[TABLE]

The original functions $W_{g,n}$ can be recovered using the inverse function

[TABLE]

Indeed one has,

[TABLE]

Note also that the corresponding coefficients of the expansion of $W_{g,n}$ at infinity, that is the $c^{[g]}_{k_{1},\ldots,k_{n}}$, can be obtained by computing residues

[TABLE]

It is also true that the residue in $z$ variables can equivalently be computed at infinity. The passage from the $W_{g,n}$ to the $w_{g,n}$ functions takes into account the Jacobian of the change of variables.

For future convenience, we define

[TABLE]

where $\sigma$ relates to $\partial_{y}P$ since $\sigma(z)=\frac{1}{x(z)^{2}}\partial_{y}P(x(z),y(z))$. So in particular $\sigma$ vanishes at the ramification point $(x_{r_{1}},y_{r_{1}})=(27/4,2/9)$ and $x(z)^{2}\sigma(z)$ has a zero of order $2$ at $(x_{r_{2}},y_{r_{2}})=(0,\infty)$.

Expression of $\tilde{w}_{0,2}$. We have after multiplying (3.3) by $x^{\prime}(z_{1})x^{\prime}(z_{2})$ and performing a few additional manipulations

[TABLE]

From this equation $\tilde{w}_{0,2}(z_{1},z_{2})$ can be computed in the variables $z_{1},z_{2}$, so that one obtains

[TABLE]

From this expression we can recover the limiting cumulants of the product of traces,

[TABLE]

We provide the reader with the first few orders on Table 1. These numbers can be obtained easily via symbolic computation softwares.

Remark 5**.**

Using a table of coefficients $c_{i,j}^{[0]}$ for $i,j$ running from $1$ to $20$ it is possible to make an experimental guess for the explicit form of these coefficients. This is

[TABLE]

In particular we have checked that these numbers satisfy the recurrence equation (82) for the first few orders. It would be interesting to prove or disprove this guess via, for instance, combinatorial means.

Universality for $w_{0,2}$. In this paragraph we explain in detail and *a posteriori555Since they can already easily be inferred from the explicit result of equation (113). * the analytic properties of $\tilde{w}_{0,2}$ and $w_{0,2}$. We first argue that $\tilde{w}_{0,2}$ does not have poles at the ramification points that is $z=-3/2,0$. We then consider the situation when $x(z_{1})\rightarrow x(z_{2})$. First starting from the above remark that $x(z)^{2}\sigma(z)=\partial_{y}P(x(z),y(z))$, we know that $x(z)^{2}\sigma(z)$ has a double zero at $z=0$ and a simple zero at $z=-3/2$, which makes it a source of poles as this factor appears in the denominator in front of the two terms of (116), see below

[TABLE]

We start by focusing on poles at the simple ramification point $z=-3/2$. We remind ourselves that $\mathrm{d}x$ vanishes at the ramification points, and so $x^{\prime}(z)$ has a simple zero at $z=-3/2$. Therefore $\frac{x^{\prime}(z_{1})}{x(z_{1})^{2}\sigma(z_{1})}$ is holomorphic at $z_{1}=-3/2$. Moreover, $x(z_{1})$ and $y(z_{1})$ are holomorphic at $z_{1}=-3/2$. As a consequence $\tilde{w}_{0,2}$ is holomorphic at $z=-3/2$ in both $z_{1}$ and $z_{2}$ (thanks to the symmetry $z_{1}\leftrightarrow z_{2}$).

We now come back to the ratio $\frac{x^{\prime}(z_{1})}{x(z_{1})^{2}\sigma(z_{1})}$ for $z_{1}=0$. A similar argument is valid at $z_{1}=0$. Indeed $x^{\prime}(z_{1})$ has a double zero at $z_{1}=0$ and this cancels the double zero of $x(z_{1})^{2}\sigma(z_{1})$ at $z_{1}=0$. In fact one can explicitly compute the ratio and find

[TABLE]

which confirms our argument. $x(z)$ is holomorphic at $z=0$, but $y(z)$ is not, indeed it has a double pole at $z=0$. So the terms $x(z_{1})x(z_{2})y(z_{1})^{2}$ could bring a simple pole at $z_{1}=0$. However, using the fact that $\tilde{w}_{0,2}(z_{1},z_{2})$ is symmetric in its arguments, if such a simple pole exists at $z_{1}=0$ then one should have a simple pole at $z_{2}=0$. Using the fact that $x(z_{2})$ has a third order zero at $z_{2}=0$, and $y(z_{2})$ has a double pole at $z_{2}=0$ one can show that $\tilde{w}_{0,2}(z_{1},z_{2})$ is holomorphic at $z_{2}=0$, therefore the apparent singularity at $z_{1}=0$ is a removable singularity. Consequently, we have just shown that $\tilde{w}_{0,2}(z_{1},z_{2})$ is holomorphic at the ramification points $z=-3/2,0$ in both its variables.

Other possible singularities may occur at the singularities of $x(z)$ which possesses a simple pole at $z=-1$ and when $x(z_{1})\rightarrow x(z_{2})$. First note that

[TABLE]

has a double zero when $z_{1}\rightarrow-1$, thus

[TABLE]

is holomorphic when $z_{1}\rightarrow-1$ since $\frac{x(z_{1})x^{\prime}(z_{1})}{x(z_{1})^{2}\sigma(z_{1})}$ has a double pole at $z_{1}=-1$. A similar argument applies to the term

[TABLE]

thus showing that $\tilde{w}_{0,2}(z_{1},z_{2})$ is holomorphic at $z_{1}=-1$, and by symmetry at $z_{2}=-1$.

We are now left with the situation $x(z_{1})\rightarrow x(z_{2})$. A first possibility is $z_{1}\rightarrow z_{2}$. In this case both ratios

[TABLE]

are holomorphic since the denominators and numerators have simultaneous simple zeroes. So $\tilde{w}_{0,2}(z_{1},z_{2})$ is holomorphic when $z_{1}\rightarrow z_{2}$. However, since $x(z)$ is a covering of degree three, there exists two (not globally defined) functions, $d_{1}(z),d_{2}(z)$ that leaves $x$ invariant, that is $x\circ d_{i}=x,\,i\in\{1,2\}$. These functions are the (non-trivial) solutions of the equation

[TABLE]

This leads to the expressions

[TABLE]

One can check that $x(d_{1}(z))=x(d_{2}(z))=x(z)$. In order to understand the pole structure of $\tilde{w}_{0,2}(z_{1},z_{2})$, one also needs to know how does $y(z)$ changes when composed with one of the $d_{i}$. One has the simple identities for $i\in\{1,2\}$

[TABLE]

Using these identities, one expects poles when $z_{1}\rightarrow d_{1,2}(z_{2})$. Indeed, in this limit the numerators of (121) does not have zeroes anymore, while the denominators have simple zeroes. Thus $\tilde{w}_{0,2}(z_{1},z_{2})$ should have double poles when $z_{1}\rightarrow d_{1,2}(z_{2})$. This is indeed what we find by requiring that the denominator of (113) vanishes.

Remark 6**.**

*The functions $d_{i}$ have interesting properties. Indeed they permute the sheets of the covering $x:\mathcal{C}\rightarrow\hat{\mathbb{C}}$. Their behavior in a small neighborhood around a ramification point relates to the local deck transformation group of the cover.

Let us first focus on the double ramification point $z=0$. It is a fixed point of both $d_{1}$ and $d_{2}$ and around $z=0$, we have $d_{1}(z)\sim_{0}e^{-\frac{2i\pi}{3}}z$ and $d_{2}(z)\sim_{0}e^{\frac{2i\pi}{3}}z$ thus they are inverse of each other locally, and generate the cyclic group $\mathbb{Z}_{3}$. This cyclic group is the group generated by the permutation of the sheets $\tau_{2}=(132)$. This group is the local deck transformation group around the ramification point at $z=0$.

We now consider the behavior of $d_{1},d_{2}$ at $z=-3/2$. In this case, only $d_{1}$ fixes $z=-3/2$, while $d_{2}(-3/2)=3,\,d_{2}(3)=-3/2$, that is $d_{2}$ exchanges the ramification point with the point above it (see Fig. 4). Note however that one has $d_{1}(3)=d_{2}(3)=-3/2$ as the two solutions $d_{1},d_{2}$ of equation (122) merge at $z=3$ (as they also do at $z=1$). This merging has the following interpretation. At $z=-3/2$ two of the three sheets of the covering coincide. Therefore, there remains effectively only two sheets to be permuted, that is why $d_{1}$ fixes $z=-3/2$ while $d_{2}$ permutes $z=-3/2$ with $z=3$. The action of the local deck transformation group at $z=-3/2$ relates to the action of $d_{1}$ in a small neighborhood of $z=-3/2$. Since $d_{1}(-3/2+\epsilon)-d_{1}(-3/2)\sim_{0}-\epsilon$, $d_{1}$ locally generates the cyclic group $\mathbb{Z}_{2}$ corresponding to the group generated by the permutation $\tau_{1}=(12)$. Similar arguments can be used to describe the local deck transformation group at the ramification point $z=\infty$.*

We now come to the universality statement. Indeed, we expect that a slightly different object than $\tilde{w}_{0,2}(z_{1},z_{2})$ takes a universal form. This is the reason for the shift introduced in (105). The statement is that $w_{0,2}(z_{1},z_{2})$ should have a universal form, that is it should be the unique meromorphic function on the sphere with a double pole of order $2$ on the diagonal with coefficient $1$ and otherwise regular. Indeed if we compute $w_{0,2}(z_{1},z_{2})$ we obtain

[TABLE]

We find exactly the expected universal form for a genus zero spectral curve.

Comment on probabilistic interpretation of $W_{0,1}(x)$, $W_{1,1}(x)$ and $W_{0,2}(x_{1},x_{2})$. As stated earlier, $W_{1}(x)$ is the Stieltjes transform of the eigenvalues density of the matrix $S_{2}$, that is

[TABLE]

In particular in the large $N$ limit we have that

[TABLE]

and the computation of $W_{0,1}(x)$ uniquely determines $\rho_{0,1}(x)$. The same property is also true for the exact density, i.e. $W_{1}(x)$ uniquely determines $\rho_{1}(x)$. This can be traced back to the Carlemann condition [Akh65]. Indeed the Stieltjes transform $W_{1}(x)$, (resp. $W_{0,1}(x)$) contains the information on the whole moment sequence of $\rho_{1}(x)$ (resp. $\rho_{0,1}(x)$). The sequence of moments of both distributions can be shown to satisfy the Carlemann condition, and thus one expects that the knowledge of the Stieltjes transform is sufficient to reconstruct the densities $\rho_{1}(x)$, $\rho_{0,1}(x)$. However it is known [FFG06] that in general the truncation of the $1/N$ expansion of the resolvent does not determine a unique truncated density. Indeed, there exists, a priori, multiple densities truncated at order $p$, $\rho^{(p)}_{1}(x)=\sum_{g\geq 0}^{p}N^{-2g}\rho_{g,1}(x)$ with the same truncated resolvent

[TABLE]

That is the computation of the corrections to $W_{0,1}(x)$ only determines Stieltjes class of densities666Though this is not a rigorous justification, one can look at the truncated Carlemann criterion, for instance in the GUE case, and see that the Carlemann criterion is indeed not satisfied order-by-order in $1/N$. Only the large $N$ and the exact criterion are satisfied., often referred to as a smoothed density. This is sufficient however to compute the corrections to the average $\mathbb{E}(\phi(x))$ where $\phi(x)$ is any function analytic on the support of $\rho_{0,1}(x)$. In particular, our later computation of the first few corrections to the large $N$ resolvent does not determine corrections $\rho_{1,1}(x),\rho_{2,1}(x),\ldots$

The probabilistic interpretation of $W_{0,2}$ goes as follows. $W_{2}$ is the Stieltjes transform of the connected part of the eigenvalue correlation function

[TABLE]

and

[TABLE]

where the $\lambda_{i}$ are the eigenvalues of the matrix $S_{2}$. In the large $N$ limit, the centered random vector whose components are the traces of successive powers of the matrix $S_{2}$, $\left(\mathrm{Tr}(S_{2}^{i})-\mathbb{E}(\mathrm{Tr}(S_{2}^{i}))\right)_{i=1}^{k}$ converges to a normal random vector of zero mean and variance $\textrm{Var}_{m,n}$

[TABLE]

where the normality of this centered random vector at large $N$ follows from the fact that $W_{n}(x_{1},\ldots,x_{n})=O(1/N^{n-2})$, that is the higher cumulants of the limiting distribution of the family $\{\mathrm{Tr}(S_{2}^{i})\}$ vanish at large $N$. This statement extends to the large $N$ limit of any linear statistics $A$ of the eigenvalues of the form

[TABLE]

where $a$ is a sufficiently smooth function (analytic for instance), as we have

[TABLE]

with $\Gamma$ a contour encircling the cut $(0,27/4]$ of $W_{0,1}(x)$.

4.3 Computation of $w_{1,1}$ and higher correlation functions.

From these data one can access the first correction to the resolvent which allows in turn to access a first correction to the large $N$ density. The equation for $w_{1,1}(z)$ can be easily obtained from the equation (66) on $W_{1,1}(x)$. It reads

[TABLE]

This leads to the result of the next paragraph.

Expression of $w_{1,1}(z)$ and analytic properties of (135). We obtain,

[TABLE]

We notice that the poles are located at $z=0$ and $z=-3/2$, which are the zeroes of $\mathrm{d}x$. However, starting from (135) one can only infer that the poles of $w_{1,1}(z)$ can be located at $z=0,-3/2,-1$. Indeed, one can easily obtain from the analytic properties of $x(z),y(z)$ and $\tilde{w}_{0,2}(z,z)$ that the first term of the right hand side of (135) can have poles only at $z=0,-3/2$, and rule out singularities at $z=-1,\infty$. However when considering the derivatives term, that is the second term of equation (135), one can not rule out poles at $z=-1$. The explicit computation shows that the coefficient of these poles is zero.

Remark 7**.**

Note that we can also produce a guess for the coefficients $c^{[1]}_{n}$. We need however to prove our first guess of Remark 5 for $c_{i,j}^{[0]}$ to be able to prove this guess using the Schwinger-Dyson equations. We provide our guess for purely informative purposes,

[TABLE]

Expression for higher correlations. Using the loop equations (96) we can compute any $n$-point resolvents recursively at any order. We illustrate this claim by providing the first few resolvents of higher order.

One point case.

[TABLE]

Two points case.

[TABLE]

with $pol(z_{1},z_{2})$ a symmetric polynomial of $z_{1},z_{2}$ of degree $12$,

[TABLE]

Three points case.

[TABLE]

For all these computed $w_{g,n}$, $(g,n)\neq(0,1),(0,2)$ the poles are located at $z=0$ and $z=-3/2$. Therefore we can expect that the poles of $w_{g,n}$, for $2g-2+n>0$, are always located at $z=0$ and $z=-3/2$, however this remains to be proven.

Remark 8**.**

The computed $w_{g,n}$ are rational functions of the $z_{i}$. We notice that the numerator of these rational functions seems to be a polynomial with positive integer coefficients. If this property is true for every $w_{g,n}$, it would be interesting to understand if these positive integers have an enumerative (combinatorics or geometry) meaning.

5 Conclusion

In this first paper on loop equations for matrix product ensembles, we have shown how to obtain loop equations for any resolvents for a random matrix defined as a product of two square complex Ginibre matrices without resorting to an eigenvalues or singular values reformulation of the problem. Indeed, the eigenvalues reformulation is yet to access these observable quantities. We used these loop equations to compute several terms of the expansion of the any resolvents $W_{n}$. In particular we accessed $W_{0,2}$, giving us information on the fluctuations of linear statistics, as well as the first correction $W_{1,1}$ to $W_{0,1}$. We expect a similar technique to apply to the more general case of the product of $p\geq 2$ rectangular Ginibre (complex or real) as well as to some other product ensembles, for instance the ensembles introduced in [FIL18] that are closely related to the Hermite Muttalib-Borodin ensemble.

Several questions are suggested by this work. The most straightforward one concerns the establishment of a topological recursion formula for the $w_{g,n}$. In the present case this topological recursion formula is certainly similar to the one devised in [BHL*+*14, BE13] by Bouchard and al. and Bouchard and Eynard. We postpone the construction of such formula to further works. Another interesting question oriented towards enumerative geometry concerns the application of the same technical means to the matrix model introduced by Ambjørn and Chekhov in [AC14, AC18] which generates hypergeometric Hurwitz numbers. In these works the spectral curve is obtained, however this is done via a matrix-chain approach that requires $p-1$ of the $p$ matrices to be invertible, thus ruling out the fully general case of rectangular matrices. We hope this fully general case can be tackled using our higher derivatives technique.

Yet another related question is the following. Free probability provides us with tools to determine the equation satisfied by the large $N$ limit of the resolvent of a product of matrices knowing the large $N$ limit of the resolvents of the members of the product. These tools have been generalized to some extent to the $2$-point resolvent in the works of Collins and al. [CMSS07] in order to more systematically access the fluctuations of linear statistics. One question is then the following. Can we devise similar tools that would allow to construct the full set of loop equations for a product matrix knowing the loop equations satisfied by the member of the product (or, more realistically, the large $N$ sector of the loop equations)?

Finally, the loop equations can be interpreted as Tutte equations [Eyn16, Tut62, Tut68]. The loop equations described in this paper can also be interpreted combinatorially, and it would be interesting to understand the more general case of maps with an arbitrary number of black vertices in such a combinatorial setting. Moreover, one would also like to understand if it is possible to merge two sets of Tutte equations together for two independent sets of maps with one type of edge in order to obtain Tutte equations for maps with two types of edges. The combinatorial interpretation of the free multiplicative convolution described in [DLN18, section 3.3] may be a useful starting point.