{'fraction_non_alphanumeric': 0.05293244480939898, 'fraction_numerical': 0.018033667056185452, 'mean_word_length': 4.103803997421019, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Cauchy two-matrix model

We introduce a new class of two(multi)-matrix models of positive Hermitean matrices coupled in a chain; the coupling is related to the Cauchy kernel and differs from the exponential coupling more commonly used in similar models. The correlation functions are expressed entirely in terms of certain biorthogonal polynomials and solutions of appropriate Riemann-Hilbert problems, thus paving the way to a steepest descent analysis and universality results. The interpretation of the formal expansion of the partition function in terms of multicolored ribbon-graphs is provided and a connection to the O(1) model. A steepest descent analysis of the partition function reveals that the model is related to a trigonal curve (three-sheeted covering of the plane) much in the same way as the Hermitean matrix model is related to a hyperelliptic curve.


Introduction
In the last two decades or so, the theory of matrix models have been an incredibly fertile ground for a fruitful interaction between theoretical physics, statistics, analysis, number theory and dynamical systems (see the classical [35] and references therein).
The interplay with analysis has been particularly beneficial for the Hermitean matrix model [14], largely due to the realization that the matrix model could be "solved" in terms of orthogonal polynomials. This allowed to translate questions about the spectrum of large random matrices into questions about the asymptotics of orthogonal polynomials, connecting the former with a very well developed area of analysis. The availability of Riemann-Hilbert methods [28,13] was the crucial ingredient in addressing questions of universality in the bulk and at the edge of the spectrum.
It should be mentioned that this fortunate symbiosis relies on the following features: • the possibility of rewriting the matrix integral in terms of eigenvalues and the correlation functions in terms of suitable orthogonal polynomials; • the Riemann-Hilbert characterization of orthogonal polynomials [28]; • the Christoffel-Darboux formula, allowing to rewrite the kernel of the correlation functions in terms of only two orthogonal polynomials of degree N ; • the nonlinear steepest descent method [13] applied to the RH problem for orthogonal polynomials.
There are by now several matrix models based on various ensembles of matrices. For some (Symplectic, Orthogonal), a connection with (skew-)orthogonal polynomials can be established [35,37] and some of the steps in the above list have been performed, but typically only for a certain restricted class of "potentials".
For others (like O(n) models, rectangular models) the methods used rely on symmetries of the integral and dynamical system approaches (the "loop equations") [18,16,17,20,1] There are also the so-called multi-matrix-models which involve ensembles of several matrices (typically Hermitean) and one of the most studied among them is the Itzykson-Zuber-Harish-Chandra (IZHC) chain of matrices. In its simplest form, the two-matrix model, it amounts to the study of the spectral properties of a pair of Hermitean matrices of size N × N with a measure dµ(M 1 , M 2 ) = dM 1 dM 2 e −N Tr(V1(M1)+V2(M2)−M1M2)  where the interaction term is e N TrM1M2 . The first three items of the bulleted list above can be implemented at least for potentials whose derivative is a rational function [19,6,4,5] while the last item is still not fully under control. One of the main reasons for this difficulty is that the size of the RHP depends on the potentials: for example if V j are polynomials of degree d j then there are two relevant RHPs of the same size. Clearly the question arises as to whether a Riemann-Hilbert method can be utilized for a more general class of potentials, for example real-analytic, like in the case of the single-matrix model.
The present paper is a part of a larger project initiated in [8], which puts forward a new multimatrix model (at the moment primarily two matrices) that can be completely solved along the lines described above.
While the new model is very close in spirit to the IZHC model, the different interaction links the model to a new class of biorthogonal polynomials which were studied in [8] and termed "Cauchy biorthogonal polynomials". These polynomials have several of the desirable features of classical orthogonal polynomials.
We will fill in more details in Sect. 3 but here we just indicate that the model we want to study is (1-4) whose dependence on the measures α, β carries all relevant information about the model.
We will show that this model is related to Cauchy biorthogonal polynomials R+ R+ p n (x)q m (y) α(x)dx β(y)dy x + y = δ mn (1)(2)(3)(4)(5) defined and studied in [8] in relation with the spectral theory of the cubic string and the Degasperis-Procesi wave equation (see also [2,3,34,33]). Contrary to the biorthogonal polynomials of [19,6,7] the properties of these polynomials do not depend on the measures α, β and the main highlights are 1. they solve a four-term recurrence relation; 2. their zeroes are positive and simple; 3. their zeroes have the interlacing property; 4. they possess Christoffel-Darboux identities relevant to matrix models; 5. they can be characterized by a (pair of) Riemann-Hilbert problem(s) of size 3 × 3; 6. the steepest descent method is fully applicable.
Our recent paper [8] contains a detailed discussion of all these points except the last one which will be addressed in a forthcoming publication. the solution of the Riemann-Hilbert problems; this will be used in later section to describe the spectral statistics of the two matrices. In Sect. 3 we introduce in detail the two-matrix model that we have outlined above and show how to reduce the study of its spectral statistics to Cauchy biorthogonal polynomials using a formula appeared in [25]. We also indicate (without details) how to deal with a similar model where R matrices are linked in a chain (Sect. 3.2).
In Sect. 4 we show how a formal treatment of the partition function of the model for large sizes of the matrices can be used to extract combinatorial information for certain bi-colored ribbon graphs (as it has been done for the Hermitean one-matrix model in [11,32] and for the multi-matrix models with exponential coupling in [16]).
In Sect. 5 we show that a saddle-point treatment of the partition function (see for example [15] for the one-matrix case) leads to a three-sheeted covering of the spectral plane (a trigonal curve); this (pseudo) algebraic curve plays the same rôle as the hyperelliptic curve in the Hermitean matrix model.

Cauchy biorthogonal polynomials
Let K(x, y) = 1 x+y be the Cauchy kernel on R + × R + . It is known that it is totally positive in the sense of the classical definition [22,31]: is a function such that for all n ∈ N and ordered n-tuples x 1 < x 2 < · · · < x n , y 1 < y 2 < · · · < y n we have the strict inequality It was shown in [8] that for any totally positive kernel K(x, y) (hence also for the Cauchy kernel) and for any pair of measures α(x)dx, β(y)dy supported in I, J ⊂ R + respectively, the matrix of bimoments is a totally positive matrix, that is, every square submatrix has a positive determinant. This guarantees the existence of biorthogonal polynomials {p j (x), q j (y)} of exact degree j such that The polynomials p j (x), q j (y) are defined uniquely up to a C × action p j → λ j p j (x), q j → 1 λj q j and the ambiguity can be disposed of by requiring that their leading coefficients are the same and positive. With this understanding we have p n (x) = 1 cn x n + . . . and q n = 1 cn y n + . . . , with the positive constant c n given by [8] c n = D n+1 D n , D n := det I jk 0≤j,k≤n−1 , I jk := The determinantal expressions for the BOPs in terms of bimoments can be obtained using Cramer's rule [8].
Proposition 2.1 (Thm. 4.5 in [8]) For any totally positive kernel K(x, y) the zeroes of the biorthogonal polynomials p j 's (q j 's) are simple, real and contained in the convex hull of the support of the measure α (β respectively).
In the case K(x, y) = 1 x+y we named the corresponding polynomials {p j (x), q j (y)} Cauchy BOPs and proved, in addition, Proposition 2.2 (Thm. 5.2 in [8]) The roots of adjacent polynomials in the sequences {p j (x)}, {q j (y)} are interlaced.
Cauchy BOPs enjoy more structure than the BOPs associated to a generic kernel: they solve a four term recurrence relation of the form As proved in [8] π n , η n are strictly positive. The other coefficients appearing in the recurrence relations are described in loc. cit.
After lengthy manipulations one obtains several Christoffel-Darboux-like identities (CDIs) for Cauchy BOPs which play a crucial rôle in what follows and hence will be carefully described.

Riemann-Hilbert characterization of the integrable kernels
The sums appearing on the left hand side in Prop. 2.3 are all examples a general framework of "integrable kernels" that were studied in great generality in [30,27,29,26,24].

Remark 2.2
The quantities referring to the letter q on the right hand side in Prop. 2.3 can be extracted from the RHP involving Γ. The next proposition will achieve the same goal for the remaining quantities.

Γ(z) satisfies the jump conditions
3. its asymptotic behavior at z = ∞ (z) = 0 is Then such a Γ(z) exists and is unique. Moreover Γ(z) can equivalently be written as: or also where the normalization constants η n , c n have been introduced in (2-4, 2-6).
n (z) were introduced in [8] independently from a RHP formulation, but for all practical purposes the formulation above is sufficiently explicit.
Finally the RHPs above allow us to reconstruct the ratio of two consecutive principal minors of the bimoment matrix (this will become relevant when discussing the partition function of the matrix model) Proposition 2.6 (Corollary 8.1 in [8]) If D n is a leading principal n × n minor of the bimoment The CDIs in Prop. 2.3 can be written in a very simple form in terms of the solutions of the Riemann-Hilbert problems. To show this let us momentarily denote by Y (w) and Y (z) the matrices with columns given by the q Now the perfect duality of Cor. 2.1 implies that Note that the solutions of the RHPs Γ(w), Γ(z) differ from Y (w), Y (z) only by some constant invertible left multipliers, and hence We collect this into the following proposition for later reference is given in terms of the solution of the Riemann-Hilbert problem (2-24) -  as The relevance of the matrix kernel H n (z, w) will become clear when we will discuss the spectral statistics of the matrix model.

Matrix Models
Consider the vector space H N of Hermitean matrices of size N endowed with the U (N )-invariant Lebesgue Inside H N we consider the convex cone of positive definite matrices with the induced measure. Let α, β be two positive densities of finite mass supported on the positive real axis: α(M ) will simply mean the product measure on the eigenvalues of M .
Define now the measure on The resulting random matrix model falls into the general class of two-matrix-models although the "coupling" term is not the most common one (which is e T r(M1M2) ). However these and much more general couplings have been considered relatively recently in [25].
As customary, we re-express the Lebesgue measures dM i in terms of the normalized Haar measure dU i of the unitary group U (N ) and the Lebesgue measure on the cone R N At this point we need to compute the integral over U (N ); we can use the result on pages 23-24 of [25] for the special example (A-28) is the Pochammer symbol.
Multiplying both sides by det(X) −r we then obtain The case of main relevance to us is r = N , which yields (C N, (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13) This shows that the measure dµ(M 1 , M 2 ) can be reduced to a measure on the spectra of the two matrices (we use the same symbol for the measure on the eigenvalues) In fact we could have used the general formula (3-8) for any r := N + h; note, however, that for h integer and less than −1 Harnad-Orlov's formula in the form presented above cannot be used, since the determinant in the numerator vanishes (for N ≥ −h) and so do some denominators in the definition of the constants C N,r . In other words, one should take an appropriate limit and use de l'Hôpital's rule. For any h ≥ 0 or any h / ∈ −N, tracing the steps above one could obtain general models where the reduced measures on the spectra have a form (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) corresponding to the unreduced measure (3-16) It will be important for us in what follows that for any value of h ∈ −N the kernel K h is totally positive or at least sign regular on R + × R + . (Sign-regularity means that determinants in Def. 2.1 are nonzero and their sign depends only on their size.) Proof. To prove this assertion define for 0 < y 0 < · · · < y n functions k i (x) = 1 (x+yi) s , i = 0, . . . , n, s = 1 + h . Denote for j = 1, 2, ..., c j (s) = s(s + 1) . . . (s + j − 1) and c 0 (s) = 1. Then the Wronskian Thus, W (k 0 , . . . , k n )(x) is a nonzero multiple of (−1) n(n−1) 2 ∆(z 0 , . . . , z n ), where z i = 1 (x+yi) and ∆(z 0 , . . . , z n ) is the Vandermonde determinant constructed out of the z i 's. Since z 0 > . . . > z n , we see that is positive for any n if s > 0 since all the c j (s) are positive numbers; if s < 0 then -denoting by [s] the greatest integer less than s (hence negative) -we see that the sign of W is for any n ≥ 0 and x ≥ 0. Together these observations imply (Thm 2.3, ch. 2 of [31]), that K h is totally positive for s > 0 on R + × R + or at least sign regular for s ∈ R − \ −N. Q.E.D.
In App. A we will show that for h integer or half-integer, the model (3-16) is the reduction of a 3-matrix model, where a "ghost" gaussian matrix A has been integrated out. Depending on the range h < N or h > N the matrix A consists of ordinary variables (bosons) or Grassmann variables.

Remark 3.1
The fact that the measure dµ h depends on N (in the determinant in the denominator) is a feature of the model rather than a problem; indeed, in studying the large N limit it is natural to make the strength of the interaction increase at the same rate as the size of the matrices. For example, in the "standard" two-matrix model one considers the interaction e cN T r(M1M2)
The computations of these correlation functions in term of biorthogonal functions reported below follows the general approach in [19] and [23].
Using the well-known formula for the Cauchy determinant we obtain the measure We will be using the correlation functions only to compute expectations of spectral functions, namely functions of X, Y which are separately symmetric in the permutations of the x j 's or y j 's.
Proof. Under the assumption F (X) = F (X σ ) for any σ ∈ S N we have Q.E.D.
Corollary 3.1 Let F (X, Y ) be symmetric with respect to either set of variables X or Y . Then Since we are interested in the unordered spectrum of the matrices M 1 , M 2 , in view of the above Cor. 3.1 we will focus henceforth on the following unnormalized measure From the properties of the Vandermonde determinant we can write where p i , q i are any monic polynomials of exact degree i in the respective variables (i = 0, . . . , N − 1).
It is therefore natural [19] to choose the sets of monic polynomials { p j (x), q j (y)} j∈N to be biorthogonal with respect to the bi-measure α(x)β(y) x+y dxdy. This means that The constant c k was defined in (2-4, with K(x, y) = 1/(x + y)). The normalization constant for d ν is If we introduce the orthonormal biorthogonal polynomials p n := 1 c n p n , q n := 1 c n q n (3-28) then we can write the normalized measure as Notice now that the product of determinants in (3-29) is a determinant of the product of the two The kernel K N (x, y) is a "reproducing" kernel which follows immediately from the biorthogonality. In addition we have Summarizing, in the new notation we have The probability measure on X × Y induced by (3-3) is given by: while the correlation functions are: Example 3.1 Consider r = 1, s = 0. In this case, with the help of Lemma 3.2, we obtain:

Correlation functions in terms of biorthogonal polynomials
In the paper by Eynard and Mehta [19] (generalized in [23]) they never used any specific information about the model they were considering (with the Itzykson-Zuber interaction) but only the fact that the matrices were coupled in a chain. We recall the relevant result here 4 . Define Since H 00 is a reproducing kernel and similarly Integrating these two equations against (x + y ) −1 α(x)dx (the first) or (x + y ) −1 β(y )dy (the second) we find also H 11 (y, x)H 10 (y, y )β(y)dy = H 11 (x, y)H 01 (y, y )α(y)dy ≡ 0 (3-41) The correlation functions for r eigenvalues x 1 , . . . , x r of M 1 and s eigenvalues y 1 , . . . , y s of M 2 were computed in [19] and are given by where 0 ≤ r ≤ N , 0 ≤ s ≤ N , 1 ≤ r + s, with the understanding that if either r or s is 0 then the corresponding blocks labeled by r, s respectively, are absent. Thus, for example R (1,0) This is a nontrivial result and perhaps one can get a bit of an insight by considering the special case r = s = N . We know that in this case so we have the following not at all obvious identity: Proof It suffices to observe that the 2N × 2N matrix on the right hand side has the following block structure: where P 0 (x) = [p i−1 (x j )] 1≤i,j≤N , P 1 (y) = [p (1) i−1 (y j )] 1≤i,j≤N and, after exchanging p i s with q i s, we define the remaining symbols (see equations (2)(3)(4)(5)(6)(7)(8)(9)(10)(11) and (2-10) for definitions) accordingly. Moreover, K(x, y) = [K(x i , y j ] 1≤i,j≤N . With this notation in place it is clear that K admits the following (Bruhat) decomposition: . (3-47) A comparison of the definitions (3-38) with the entries of the Christoffel-Darboux identities of Prop.

shows that they are intimately related, in fact
and we recognize that (up to some signs) these are precisely the entries of the H n kernel of Prop. 2.7: More explicitly we have The potential relevance of these formulae is that it allows one to compute the large N asymptotic behaviour of the correlation functions in terms of the asymptotic behaviour of only 3 consecutive biorthogonal polynomials and auxiliary functions associated with them that enter in the Riemann-Hilbert formulation of Props. 2.4, 2.5.
The steepest-descent analysis of these problems which will appear in a forthcoming paper.
We expect that, after a complete description of the asymptotics of the BOPs is obtained, these formulae can be used to address the issue of universality for this matrix model via the Riemann-Hilbert approach.
Remark 3.2 Note also that the formula above is not symmetric inasmuch as the BOPs p, q play different roles in the Christoffell-Darboux theorem; we could however rewrite the theorem as where B(z) was defined in (2-9). The auxiliary vectorq 0,N enters in a similar Riemann-Hilbert formulation with the rôles of the densities α and β interchanged.

A multi-matrix model
It is possible to extend the model (3-3) to a chain of matrices with the nearest neighbor interaction The "strength" of the interaction (i.e. h j 's) may depend on a site along the chain. Specifically, consider the space H + Using the Harnad-Orlov formula one we obtain, up to normalization constant, the following measure which we denote by the same symbol Following the same steps that led to the expression (3-24) as positive density for (S N ) 2 -invariant observables, we obtain the following reduced measure on the spectra (up to a normalization depending only on h 's but not on the measures) In the case where all interactions are the same h = 0 we have which can be seen as a generalization of . This model too can be treated with the aid of biorthogonal polynomials that will now satisfy a R + 3 recurrence relation; note that the length of the recurrence relations does not depend on the "potentials", or densities, α j . This is in sharp contrast with the usual multi-matrix model with interaction e cTr(Mj Mj+1) [19,6]. To characterize these biorthogonal polynomials, an (R + 2) × (R + 2) Riemann-Hilbert problem can be set up and the strong asymptotics can be dealt with, but the complexity of this problem definitely warrants a separate paper.

Diagrammatic expansion
In parallel with the 1 N 2 -expansion for the Hermitean matrix model and the IZHC two-matrix model, we would like to sketch the similar formal expansion of the model in terms of colored ribbon graphs. The The procedure amounts to a perturbative Taylor-expansion around a Gaussian integral. In other words we will be considering a partition function in the form where U p , V p are the perturbations of the Gaussian (including a quadratic term as well) which, for convenience, we parametrize as the following formal series We thus have where the average is taken w.r.t. the underlying (uncoupled!) Gaussian measure Note that, due to the shifts in u j , v j , the = 1 term in the third sum above cancels exactly against the shifts in the first two sums; similarly, only the term ζη 2 Tr(M 1 M 2 ) remains in the quadratic part: Using Wick's theorem for the evaluation of Gaussian integrals and the frequently used combinatorial interpretation (see [15] for an excellent introduction) one sees that the partition function is the sum over all possible bi-colored ribbon graphs respecting rules which we now specify • There are two colors for vertices/edges (say, red/blue); vertices are distinguished in monochromatic and bichromatic.
• Each monochromatic vertex of valency j enters with a weight N u j /j (blue) or N v j /j (red).
• Each edge (red or blue) enters with a weight 1 N .
• Each bichromatic vertex of valency enters with a weight ζ k η −k  Summing up over all possible labeling of the bicolored fat-graphs Γ leaves a factor |Aut(Γ)| in the partition function (see [15]) and the result is hence where n j = n j (Γ) is the number of blue vertices with valency j, m j = m j (Γ) the number of j-valent red vertices and r k = r k (Γ) is the number of -valent bichromatic vertices with k red (and hence − k blue) legs.

Large N behaviour
Consider two densities α(x), β(y) of the form Here we have introduced a dependence on the small parameter on the measures (but we will not emphasize this dependence in the notation). From the experience amassed in the literature on the ordinary orthogonal polynomials we start by considering the heuristic "saddle-point" for the partition function Z N , namely the total mass of the reduced measure . The heuristics calls for a scaling approach where we send the size of the matrices N to infinity and the scaling parameter to zero as O(1/N ), namely that we vary the densities α, β by raising them to the power 1 . We thus have Using 3-22 we have Taking 1 N 2 times the logarithm of the integrand we have the expression It is convenient -in order to deal with a more standard potential-theoretic problem-to map Y → −Y and define V (y) = V (−y) so that we can rewrite the action as which describes the energy of a gas of 2N particles of charge +1 (for the x j 's) and −1 (for the y j 's) separated by an impenetrable, electrically neutral, partition confining the positively charged particles to the positive real axis under the potential U and the negatively charged particles to the negative axis under the potential V .
The usual argument is that the configuration of the minimum contribute the most to the integral and -under suitable assumptions of regularity for the potential-one wants to show that the sequence of these minima configurations tends in-measure to some probability distributions.
We define and rewrite the action as

Continuum version: cubic spectral curve and solution of the potential problem
We immediately rephrase the above minimization problem in a continuum version. In order to emphasize the symmetry of the problem it is convenient to denote U (x) by V 1 (x) and V (y) = V (−y) by V 2 (y) and denote the corresponding equilibrium densities by ρ 1 (x), ρ 2 (y).
The main point of this section is Thm. 5.1 which states that the resolvents (Markov functions) of the equilibrium distributions are solutions of a cubic equation that defines a trigonal curve 5 ; this result is the analogue of the better-known result for the equilibrium measure appearing in the one-matrix model ( [12] and references therein). The derivation that we present here is of a formal heuristic nature inasmuch as we discount several important issues about the regularity of the equilibrium measures. However this sort of manipulations is quite common and they can be obtained also from the loop equations as done in [18] for the O(n) model.
Rewriting the functional in the new notation for the potentials we find We will be studying the measures ρ 1 , ρ 2 that minimize the above functional; here we will assume their regularity (which shall be proved in a separate paper). Also, the potentials V 1 , V 2 will be assumed real analytic.
More precisely, we make the following assumption on the nature of the potentials • the derivatives V j (x) have at most finitely many poles in a strip of finite width around R; • finally, V 1 (x) > C 1 | ln(x)| for some constant C 1 > 0 and x > 0 and V 2 (y) > C 2 | ln |y| | for y < 0 and some constant C 2 > 0.
It will be shown in a separate paper that this assumption is sufficient to guarantee that the equilibrium densities exist, have compact support, are regular and their supports do not include the origin. In order to enforce the normalization of ρ 1 , ρ 2 -as customary-we will introduce in the action two suitable Lagrange multipliers S := S[γ 1 , γ 2 , ρ 1 , ρ 2 ] := S[ρ 1 , The vanishing of the first variation of the functional S yields the equations Differentiating w.r.t. z yields where P.V. indicates the Cauchy principal value. These equations are best written in terms of the resolvents (or Weyl functions) Indeed, with the help of the Sokhotskyi-Plemelj formula, the equations above take a simpler form Note that in the RHS of the system of equations, the resolvent of the other measure is evaluated at a regular point due to the fact that the two supports are disjoint.
We define, partly motivated by [21], the shifted resolvents . (5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) Observe that, in view of Assumption (5.1) , {Y 1 , Y 2 } have the same analytic structure as {W 1 , W 2 }, while the equations describing the jumps simplify to: Using these equations one obtains by direct inspection that has no jumps on either supports. Multiplying eq. (5-18) on both sides by Y 1 − Y 2 one obtains which can be rewritten as (5-20) The first expression may a priori have at most jumps across the support of ρ 1 , while the second may have jumps only across the support of ρ 2 : since the two supports are disjoint we conclude that D(x) is a regular function on the real axis.
Proof. The statement about the support of the measures follows from the argument used in [12] that that we sketch here. By results of Chapter 13 of [36]), Assumption 5.1 implies that the supports of the equilibrium measures are compact. Moreover, it follows that the branchpoints (and branchcuts) of the functions Y 1,2,0 coincide with the supports of ρ 1 , ρ 2 . Since they solve a cubic algebraic equation, the branchpoints are determined as zeroes of the discriminant ∆ := 4R 3 − 27D 2 of (5-22).
As suggested by the previous discussion, even without an explicit expression for R, D, it follows from Morera's theorem that R, D are also real-analytic and so is ∆. Thus ∆ cannot have infinitely many zeroes in a compact domain and so the number of endpoints of the supports is a-priori finite.
The only part that is left to be proven are the formulae for R, D.
The second formula for D(x) can be obtained by noticing that D(x) = Y 0 Y 1 Y 2 . Explicitly, this reads (after rearranging the terms and using the definitions for the shifted resolvents (5-15) and the U j 's (5-26) Substituting in the last two terms the expressions (5-31, 5-32) we obtain Using W 1 W 2 = W 12 + W 21 (5-31, 5-32) we obtain Now we can rewrite so that we have obtained We want to show that R ≡ 0: indeed it is clear that the only part of D(x) which may have jumpdiscontinuities is R, but we know that D(x) has no such discontinuities.
Hence R has no discontinuities; on the other hand it is clear from its definition that R cannot have any other singularities and hence it is an entire function. Inspection shows that R(z) → 0 as z → ∞ and hence R must be identically vanishing by Liouville's theorem. This concludes the proof of the theorem. Q.E.D.
The presence of a "spectral curve" will be one of the crucial ingredients for the large-degree asymptotic analysis of the biorthogonal polynomials using the Riemann-Hilbert formulation given in [8] and the Deift-Zhou nonlinear steepest descent method. Indeed we will show in a separate publication that the OPs are modeled by (spinorial) Baker-Akhiezer vectors (similarly to [9]) that naturally live on the threesheeted covering specified by (5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22).  h) × N, C) we will use a shorthand notation for the volume element. Let, as above, M 1 , M 2 ∈ H + N and consider the following normalizable measure on the space Since the measure is Gaussian in A one immediately sees that (up to inessential proportionality constants) where h ∈ 1 2 Z and h < N obtaining then Of course, if h is actually an integer then this case reduces to the previous one.

A.3 Integer and half integer h > N
For the sake of completeness we note that when h > N the determinant in the reduced measure (3-16) is a positive power in the numerator; this can be obtained from a Gaussian integral over Grassmann variables which can be obtained from "complex" anticommuting variables (for h integer) or "real" (for h ∈ 1 2 Z) much along the line of the previous (commuting variable) case.
Clearly, the plethora of models is endless and each value of h can be studied with the aid of a different set of biorthogonal polynomial.
We singled out the case h = 0, since it corresponds to the Cauchy BOPs that have many features in common with the much better-known orthogonal polynomials and yet seem to have a very rich but still tractable asymptotic theory.

B Relation to the O(1)-model
The O(n) matrix model [18]  Since we are dealing with the case N = 2k even, the adjoint of a skew-symmetric matrix is skewsymmetric as well and hence the second term has to vanish. Thus we have det(I) N = det M N = (P f (M N )) 2 (B-7) We see that up to a multiplicative factor 1 We denote by Z the vector of size N with components Z = (z 1 , z 2 , . . . , z N ) = (x 1 , y 1 , x 2 , y 2 , . . . , x k , y k ) . The rational function Γ(Z) is invariant under permutations of the variables, whereas R(Z) is not; however R(Z) has the same total degree as ∆(Z), in particular its degree in y 1 is N − 1.
Symmetrizing under the integral sign the variables gives: (with the coefficient 1), and has the same total degree in all variables; since the antisymmetrization must be divisible by ∆(Z) the assertion follows. We thus have (B-12) We see then tha up to a proportionality constant Pf M N is the total integral of the measure (B-2) and thus the proposition is proved. Q.E.D.
The relationship between the two model does not seem to go much further in the sense that there is no direct and simple relationship between the correlation functions of the two models. It seems, however, that some connection should be present and is worth exploring. We leave it as an open problem to establish a connection between these two models on the level of the correlation functions.