# A Riemann-Hilbert Approach to Skew-Orthogonal Polynomials of Symplectic Type

Alex Little

► To cite this version:

Alex Little. A Riemann-Hilbert Approach to Skew-Orthogonal Polynomials of Symplectic Type. Symmetry, Integrability and Geometry: Methods and Applications, 2024, 20 (076), 32 p. 10.3842/SIGMA.2024.076. [hal-04871503]

HAL Id: hal-04871503

https://hal.science/hal-04871503v1

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A Riemann–Hilbert Approach to Skew-Orthogonal Polynomials of Symplectic Type

Alex LITTLE

Unité de Mathématiques Pures et Appliquées, ENS de Lyon, France

E-mail: alexander.little@ens-lyon.fr

Received December 27, 2023, in final form August 06, 2024; Published online August 16, 2024

https://doi.org/10.3842/SIGMA.2024.076

**Abstract.** We present a representation of skew-orthogonal polynomials of symplectic type ($\beta = 4$) in terms of a matrix Riemann–Hilbert problem, for weights of the form $e^{-V(z)}$ where $V$ is a polynomial of even degree and positive leading coefficient. This is done by representing skew-orthogonality as a kind of multiple orthogonality. From this, we derive a $\beta = 4$ analogue of the Christoffel–Darboux formula. Finally, our Riemann–Hilbert representation allows us to derive a Lax pair whose compatibility condition may be viewed as a $\beta = 4$ analogue of the Toda lattice.

**Key words:** Riemann–Hilbert problem; skew-orthogonal polynomials; random matrices

2020 Mathematics Subject Classification: 60B20; 33C45; 30E15; 30E25

## 1 Introduction

In this paper we develop a Riemann–Hilbert approach to the theory of skew-orthogonal polynomials of symplectic type. In order to motivate our work, let us recall the basic theory of orthogonal polynomials. One begins with a measurable function $w: \mathbb{R} \to [0, +\infty]$ with finite moments to all orders, that is $\int_{\mathbb{R}} |x|^k w(x)dx < +\infty$ for all $k \in \mathbb{N}$. Then one constructs a Hermitian form

$$ \langle P, Q \rangle_2 = \int_{\mathbb{R}} P(x)Q(x)w(x)dx $$

for (real) polynomials $P$ and $Q$. We would like that $\langle \cdot, \cdot \rangle_2$ be nondegenerate on the linear space of real polynomials so that we can make it into an inner product, which would be true if $w$ is continuous and has positive mass $\int_{\mathbb{R}} wdx > 0$. Given such an inner product there is a unique sequence of polynomials $\{P_n\}_{n \in \mathbb{N}}$ such that $P_n(x) = x^n + O(x^{n-1})$ and $\langle P_n, P_m \rangle_2 = 0$ for $n \neq m$, constructable by a Gram–Schmidt procedure. $w$ is referred to as the weight function and $P_n$ is called the $n$-th monic orthogonal polynomial with respect to $w$.

The applications of orthogonal polynomials are extremely diverse, however in this work we shall be interested solely in their application to random matrix theory (see, e.g., [14]). Let $M$ be an $n \times n$ Hermitian random matrix distributed according to a probability measure

$$ \frac{1}{Z_n} e^{-\text{tr} V(M)} dM, \qquad (1.1) $$

where $dM = \prod_{k=1}^n dM_{kk} \prod_{1 \le i < j \le n} d\text{Re}(M_{ij})d\text{Im}(M_{ij})$, $V$ is a real-valued analytic function on $\mathbb{R}$ and $Z_n > 0$ is a constant of normalisation. For this measure to be normalisable, $V(x) \to +\infty$

<footnote>This paper is a contribution to the Special Issue on Evolution Equations, Exactly Solvable Models and Random Matrices in honor of Alexander Its’ 70th birthday. The full collection is available at https://www.emis.de/journals/SIGMA/Its.html</footnote>
as $x \to \pm\infty$ sufficiently rapidly, more precisely, $\liminf_{x \to \pm\infty} \frac{V(x)}{\log|x|} = +\infty$. This condition would be satisfied if $V$ was a polynomial of even degree and positive leading coefficient, as will be the case in this paper.

A remarkable calculation (see [35, Chapter 5] or [28, Chapter 5]) then shows that the correlation functions of the eigenvalues of such an ensemble are expressible entirely in terms of the orthogonal polynomials $P_n$ and $P_{n-1}$ and their respective norms, $\|P_n\|$ and $\|P_{n-1}\|$, with respect to the weight function $w(x) = e^{-V(x)}$. Thus to make an asymptotic analysis of the correlation functions as $n \to \infty$ it is necessary to make an asymptotic analysis of orthogonal polynomials with respect to more or less arbitrary weight functions. For certain “classical” weight functions, this can be done by classical steepest descent techniques, however the general case required the discovery of a Riemann-Hilbert representation of orthogonal polynomials [27] and the development of the nonlinear steepest descent technique [23]. A Riemann-Hilbert representation involves representing an object of interest in terms of the solution to a boundary value problem in the complex plane with prescribed behaviour at infinity, and can be regarded as a generalisation of the notion of a contour integral [30]. This led to a proof of universality of local correlations in the bulk and at the edge of the spectrum for ensembles of type (1.1) by means of the so-called Deift-McLaughlin-Kriecherbauer-Venakides-Zhou (DMKVZ) scheme [20, 21, 22], building off earlier work of Pastur-Shcherbina [37] and Bleher-Its [6].

Motivated by, amongst others, applications to quantum mechanical systems with time-reversal invariance (see [35, Chapter 2]), random matrix theorists have also been interested for a long time in studying not only Hermitian but also real symmetric and quaternion self-dual random matrices with probability measures of the form of (1.1). These three cases, real symmetric, complex Hermitian and quaternion self-dual, correspond to the three finite-dimensional, associative division algebras over the real numbers, $\mathbb{R}$, $\mathbb{C}$ and $\mathbb{H}$, respectively, [24]. These are often labelled by the “Dyson index” $\beta = 1, 2, 4$, respectively.

In this work, we shall be interested in the case of $\beta = 4$, the so-called symplectic case, since the probability measure is invariant under conjugation by elements of the compact symplectic group, and since we prefer to think of an $n \times n$ quaternion matrix as a $2n \times 2n$ matrix where the quaternions are represented by $2 \times 2$ blocks. For such a $2n \times 2n$ matrix, there are, with probability 1, only $n$ distinct eigenvalues since every eigenvalue occurs with multiplicity 2 (known in the physics literature as Kramers’ degeneracy). These $n$ eigenvalues have the joint PDF

$$ \rho_n(x_1, \dots, x_n) = \frac{1}{\tilde{Z}_n} \prod_{1 \le i < j \le n} |x_i - x_j|^4 \prod_{k=1}^n e^{-2V(x_k)}, \quad (1.2) $$

where $\tilde{Z}_n$ is a normalisation constant. We shall furthermore restrict to the case of $V$ being a real polynomial of even degree $D$ and leading coefficient $\gamma > 0$,

$$ V(x) = \gamma x^D + \mathcal{O}(x^{D-1}), \quad x \to \infty. $$

Here the classical approach (see [35, Chapter 5], [28, Chapter 6]) involves finding an analogue of orthogonal polynomials and this is supplied by skew-orthogonal polynomials. Here the inner product $\langle \cdot, \cdot \rangle_2$ is replaced by a skew-symmetric nondegenerate bilinear form $\langle \cdot, \cdot \rangle_4$, where

$$ \langle P, Q \rangle_4 \stackrel{\text{def}}{=} \frac{1}{2} \int_{\mathbb{R}} (P(x)Q'(x) - P'(x)Q(x))e^{-2V(x)}dx \quad (1.3) $$

for $P$ and $Q$ real polynomials. By “non-degenerate”, we mean the following. Let $\mathcal{P}_k^{\mathbb{R}}$ be the space of real polynomials of degree $\le k$. Then we say $\langle \cdot, \cdot \rangle_4$ is non-degenerate if for all $n \in \mathbb{N}$, we have the implication

$$ \forall P \in \mathcal{P}_{2n+1}^{\mathbb{R}} (\langle P, Q \rangle_4 = 0, \forall Q \in \mathcal{P}_{2n+1}^{\mathbb{R}} \Rightarrow P = 0). $$
Then a sequence of monic polynomials {$P_n$}$_{n \in \mathbb{N}}$, $P_n(x) = x^n + \mathcal{O}(x^{n-1})$, is said to be skew-orthogonal with respect to $\langle \cdot, \cdot \rangle_4$ if the following holds:

$$
(1) \quad \langle P_{2n}, P_{2m} \rangle_4 = \langle P_{2n+1}, P_{2m+1} \rangle_4 = 0 \text{ for all } n, m \in \mathbb{N}.
$$

$$
(2) \quad \langle P_{2n}, P_{2m+1} \rangle_4 = 0 \text{ for all } n, m \in \mathbb{N}, n \neq m.
$$

The existence of such a sequence is guaranteed by non-degeneracy (it can be constructed by
a procedure analogous to the Gram–Schmidt method), however these properties do not uniquely
define the sequence because one can replace $P_{2n+1} \mapsto P_{2n+1} + \alpha P_{2n}$ for any constant $\alpha \in \mathbb{R}$.
This degree of freedom may be fixed by demanding that the next-to-leading coefficient of $P_{2n+1}$
is equal to a specified value (it is common to choose this value to be zero however this will not
be the most convenient choice for us). It can then be verified that with this extra constraint the
sequence is unique.

The importance of skew-orthogonal polynomials for random matrix theory can be seen by the
following. Let $h_k = \langle P_{2k}, P_{2k+1} \rangle_4$ be the “skew norm”, where by the non-degeneracy of $\langle \cdot, \cdot \rangle_4$ we
necessarily have $h_k \neq 0$. Then (see [28, Propositions 6.1.6 and 6.1.7]) all eigenvalue correlation
functions associated to the point process (1.2) may be expressed in terms of the so-called pre-
kernel,

$$
S_n(x, y) \stackrel{\text{def}}{=} \sum_{k=0}^{n-1} \frac{P_{2k}(x) e^{-V(x)} \frac{d}{dy} (P_{2k+1}(y) e^{-V(y)}) - P_{2k+1}(x) e^{-V(x)} \frac{d}{dy} (P_{2k}(y) e^{-V(y)})}{2h_k}. \quad (1.4)
$$

Precisely, from $S_n$ one constructs a 2 × 2 matrix “kernel”, and then the $k$-point correla-
tion function of the ensemble is written as a $2k \times 2k$ Pfaffian, or equivalently, quaternion
determinant (see [28, Chapter 6]). In particular, the 1-point correlation function is given
by $\rho_n^{(1)}(x) = S_n(x, x)$. A similar treatment in terms of skew-orthogonal polynomials, though
with a different skew-inner product than (1.3), can be done in the case of real symmetric ma-
trices $(\beta = 1)$.

The importance of understanding skew-orthogonal polynomials of real and symplectic type ($\beta = 1, 4$) was stressed by Dyson [25] when he wrote:

We need to develop the theory of skew-orthogonal polynomials until it becomes
a working tool as handy as the existing theory of orthogonal polynomials.

Unfortunately, this theory remains comparatively underdeveloped relative to that of orthogonal
polynomials (β = 2). Nevertheless let us state some known results. Adler, Horozov and van Mo-
erbeke [2] give Pfaffian formulas for the even and odd skew-orthogonal polynomials for β = 1, 4.
Eynard [26] derives Heine-type multiple integral formulas for even and odd skew-orthogonal
polynomials for β = 1, 4 and he uses a non-rigorous but intriguing resolvent method to ob-
tain n → ∞ asymptotics for the case of a one-cut polynomial potential. These multiple-integral
formulas are morally equivalent to the Pfaffian formulas since the two may be related by de
Bruijn identities [13]. In the case of V' being rational, Adler and van Moerbeke [3] find a nat-
ural basis of orthogonal polynomials in which to express the corresponding skew-orthogonal
polynomials and Adler, Forrester, Nagao and van Moerbeke [1] use this to write simple ex-
pressions for the skew-orthogonal polynomials in the classical (Hermite, Laguerre and Jacobi)
cases.

A key inspiration for the present paper is the paper of Pierce [38] in which a Riemann–Hilbert representation for skew-orthogonal polynomials of real type (β = 1) with a polynomial V is derived; unfortunately this work has not received nearly the same attention as the corresponding problem for orthogonal polynomials because there is no known Christoffel–Darboux-type formula to express the pre-kernel in terms of the Riemann–Hilbert problem, and even if there was, the complexity of the jump matrix makes it unclear how to perform a nonlinear steepest descent
analysis. Pierce's paper is similar to the present one in that it considers polynomial potentials and this allows one to represent ($\beta = 1$) skew-orthogonality as a kind of multiple-orthogonality, and this naturally leads to a Riemann-Hilbert problem.

Although skew-orthogonal polynomials on the real line have received relatively little recent attention, there has been much recent work on planar symplectic ensembles and their associated skew-orthogonal polynomials. Akemann, Ebke and Parra [4] give a construction of skew-orthogonal polynomials in the plane and relate these to a corresponding set of orthogonal polynomials; while several recent works [8, 9, 10] derive a “generalised Christoffel-Darboux formula” for several different ensembles, which is a differential equation relating the kernel for the symplectic ensemble to that of the corresponding complex ensemble. These works are in a similar spirit to the method of Widom [44] in that they establish relations between the symplectic ensemble and the corresponding complex ensemble, the latter of which is often more tractable. This should be contrasted with the present work where our Riemann-Hilbert formulation does not make any direct connection with the corresponding random matrix ensemble for $\beta = 2$.

The difficulty of working with skew-orthogonal polynomials has led some to believe that they are the wrong way to approach $\beta = 1, 4$ ensembles and to develop a different formalism that avoids skew-orthogonal polynomials entirely. This is the Widom formalism [44], presented and developed in the excellent book of Deift and Gioev [18], and it is this formalism that has yielded proofs of universality of local correlations for $\beta = 1, 4$ [16, 17, 19, 31, 40, 39, 41, 42]. One of the aims of this paper is to show that an approach involving skew-orthogonal polynomials can nevertheless be highly fruitful and that it is possible to treat $\beta = 4$ ensembles in a self-contained way, without making a reduction to $\beta = 2$.

Let us now state our results. Given the sequence of skew-orthogonal polynomials {$P_n$}$_{n \in \mathbb{N}}$, we define the “dual” polynomials as $\Psi_n(x) \stackrel{\text{def}}{=} -\frac{1}{\gamma D} e^{V(x)} \frac{d}{dx} (P_n(x)e^{-V(x)})$.

(1) In Section 2, we show how $\Psi_{2n}$ and $\Psi_{2n+1}$ may be thought of as “multiple-orthogonal polynomials of Type II” (see Corollary 2.23). This yields two families of Riemann-Hilbert Problems 2.7 and 2.13 ($A_n$ and $B_n$, respectively), where $\Psi_{2n}$ and $\Psi_{2n+1}$ appear as the (1, 1) matrix element of $A_n$ and $B_n$, respectively. $A_n$ has dimension $D+1$ and $B_n$ has dimension $D+2$. Furthermore, we find that $P_{2n}$ is the (2, 2) matrix element of $A_n^{-1}$, however, somewhat contrary to what one might expect, neither $P_{2n+1}$ nor any multiple of it appears as a matrix element of $B_n^{-1}$. This seems to be related to the fact that $P_{2n}$ can be thought of as a multiple-orthogonal polynomial of Type I (see Corollary 2.25), but we are not aware of a way to represent $P_{2n+1}$ in such a way.

(2) In Section 3, we prove an analogue of the Christoffel-Darboux formula, namely

$$S_n(x, y) = -\frac{1}{4\pi i} e^{-V(x)-V(y)} \frac{(A_n^{-1}(x)A_n(y))_{21}}{x-y}, \quad (1.5)$$

where $S_n$ is the pre-kernel (1.4) and $A_n$ is the solution of Riemann-Hilbert Problem 2.7. If a nonlinear steepest descent analysis could be carried out for $A_n$, we would thus obtain a new proof of universality for $\beta = 4$, perhaps easier than existing proofs. It is also interesting to remark that the pre-kernel is expressible entirely in terms of the “even” Riemann-Hilbert problem $A_n$, i.e., from the point of view of random matrix theory the odd problem $B_n$ can be ignored. This is not the first instance that an analogue of the Christoffel-Darboux formula for skew-orthogonal polynomials has appeared in the literature; in [29], Ghosh derives a formula for $(x-y)S_n(x, y)$ in terms of something which, for polynomial potentials, has finite rank. However, no connection to a Riemann-Hilbert problem was made so an asymptotic analysis of this formula was limited to special cases.

(3) In Section 4, we show how the $A_n$ problem may be transformed to have constant jump matrices. Then allowing the potential $V$ to depend on a parameter $t \in \mathbb{R}$ by $V(z, t) =$
$V_0(z)+zt$, we derive an $(n, t)$ Lax pair. The linear system in $n$ could be regarded as a $\beta = 4$ analogue of the “three term recurrence” which appears in the theory of orthogonal polynomials (see Remark 4.7). The compatibility condition of this Lax pair yields a dynamical system which one could regard as a $\beta = 4$ analogue of the Toda lattice (see Theorem 4.1). In the $D = 2$ case an exact solution can be given (see Example 4.4). We should remark that an analogue of the Toda lattice for $\beta = 1, 4$, with connections to skew-orthogonal polynomials, already exists in the literature under the name of the “Pfaff lattice” [2]. The relationship between the Pfaff lattice and our dynamical system is left as a topic for future investigation.

As already remarked, if a nonlinear steepest descent analysis of $A_n$ could be carried out then (1.5) would give access to the asymptotic behaviour of many observables that were either out of reach of previous techniques or could only be approached by cumbersome methods. For example, it could be used to probe non-standard universality classes, analogous to those that have been studied for $\beta = 2$ when two cuts in the support of the equilibrium measure merge [7, 11, 12] and also in planar ensembles when one has a merging singularity [32]. Such non-standard universality classes are much less well-understood in the $\beta = 1, 4$ cases as compared to $\beta = 2$.

## 1.1 Notation

(1) We shall use the notation $X^k$ to stand for the rational function $X^k(x) = x^k$ for $k \in \mathbb{Z}$.

(2) For $a, b \in \mathbb{Z}$, we let $\llbracket a, b \rrbracket$ be the integer interval from $a$ to $b$, $\llbracket a, b \rrbracket = [a, b] \cap \mathbb{Z}$.

(3) For $k \in \mathbb{N}$, define the space of real polynomials of degree at most $k$ as $\mathcal{P}_k^{\mathbb{R}}$, and let $\mathcal{P}^{\mathbb{R}} = \bigcup_{k \ge 0} \mathcal{P}_k^{\mathbb{R}}$ be the space of all real polynomials. Similarly, let $\mathcal{P}_k^{\mathbb{C}}$ be the space of polynomials with complex coefficients of degree $\le k$ and $\mathcal{P}^{\mathbb{C}} = \bigcup_{k \ge 0} \mathcal{P}_k^{\mathbb{C}}$.

(4) For $m$ a positive integer and $\mathbb{F}$ a field (in our case always $\mathbb{R}$ or $\mathbb{C}$), we let $M_m(\mathbb{F})$ be the space of $m \times m$ matrices with entries in $\mathbb{F}$.

(5) We use the convention that $0 \in \mathbb{N}$.

(6) We denote the 1-torus as $\mathbb{T} = \{z \in \mathbb{C} : |z| = 1\}$.

(7) We let $\chi_A$ be the indicator function of the set $A$, i.e., $\chi_A(x) = 1$ if $x \in A$ and $\chi_A(x) = 0$ otherwise.

(8) Let $\Gamma \subset \mathbb{C}$ be a piecewise smooth contour in $\mathbb{C}$. Then the *Cauchy transform* of the measurable function $f: \Gamma \to \mathbb{C}$ is defined as

$$C_{\Gamma}(f)(z) = \frac{1}{2\pi i} \int_{\Gamma} \frac{f(x)}{x-z}dx, \quad z \in \mathbb{C} \setminus \Gamma,$$

whenever this integral exists. In our case, we will always consider sufficiently “nice” functions $f$ and contours $\Gamma$ such that the Cauchy transform is well-defined and analytic for all $z \in \mathbb{C} \setminus \Gamma$. Furthermore, we shall always consider functions with finite moments $\int_{\Gamma} |x|^k |f(x)||dx| < +\infty$ for all $k \ge 0$ and for which $C_{\Gamma}(f)(z)$ admits an asymptotic expansion in $\frac{1}{z}$ as $z \to \infty$ in the sense of Poincaré. More precisely, for every $K \in \mathbb{N}$, we have

$$C_{\Gamma}(f)(z) = -\frac{1}{2\pi i} \sum_{k=0}^{K} z^{-k-1} \int_{\Gamma} x^k f(x)dx + \frac{z^{-K-1}}{2\pi i} \int_{\Gamma} \frac{x^{K+1}f(x)}{x-z}dx.$$

If $\text{dist}(z, \Gamma) \ge \epsilon > 0$, we find that $\int_{\Gamma} \frac{x^{K+1}f(x)}{x-z}dx$ is bounded. However, if $f$ is analytic we may deform the contour to achieve boundedness in a full neighbourhood of (complex)
infinity. This will be the case in all the Cauchy integrals that appear in this paper. The general theory of such Cauchy integrals for functions of varying degrees of decay and regularity is discussed in notes of Deift [15] and the book of Muskelishvilli [36].

## 2 Riemann–Hilbert representations

Before proceeding, let us demonstrate that the problem we are considering is well-posed. The following proposition is of course known to experts but we include a proof to keep the presentation complete.

**Proposition 2.1.** The skew-symmetric bilinear form defined by (1.3) is non-degenerate. Furthermore, the skew-norms are all positive, i.e., $h_n = \langle P_{2n}, P_{2n+1} \rangle_4 > 0$, for all $n \ge 0$.

**Proof.** By one of de Bruijn's identities [13], we have

$$ \mathrm{pf}(\langle X^{i-1}, X^{j-1} \rangle_4)^{2n}_{i,j=1} = \frac{1}{n!} \int_{\mathbb{R}^n} \prod_{1 \le i < j \le n} |x_i - x_j|^4 \prod_{k=1}^n e^{-2V(x_k)} dx_1 \cdots dx_n > 0, $$

where pf is the Pfaffian. Hence the matrix $(\langle X^{i-1}, X^{j-1} \rangle_4)^{2n}_{i,j=1}$ is invertible for all $n \ge 1$ and so defines a non-degenerate bilinear form. Furthermore, by performing row and column operations,

$$ 0 < \mathrm{pf}(\langle X^{i-1}, X^{j-1} \rangle_4)^{2n}_{i,j=1} = h_0 \cdots h_{n-1}, \quad \forall n \ge 1, $$

which shows that $h_k > 0$ for all $k \ge 0$.

<figure><img src="image_2.png" /></figure>

We begin by observing that the symplectic inner product (1.3) can be re-written as

$$ \langle P, Q \rangle_4 = \int_{\mathbb{R}} P(x) e^{-V(x)} \frac{d}{dx} (Q(x) e^{-V(x)}) dx, \quad P, Q \in \mathcal{P}^{\mathbb{R}}. $$

**Definition 2.2 (dual skew-orthogonal polynomials).** Given a sequence of skew-orthogonal polynomials $\{P_n\}_{n \in \mathbb{N}}$, let us introduce the sequence of “dual” polynomials $\{\Psi_n\}_{n \in \mathbb{N}}$, where

$$ \Psi_n(x) \stackrel{\mathrm{def}}{=} -\frac{1}{\gamma D} e^{V(x)} \frac{d}{dx} (e^{-V(x)} P_n(x)) = x^{n+D-1} + \mathcal{O}(x^{n+D-2}), \quad x \to \infty. $$

The monic skew-orthogonal polynomials can be reconstructed from the dual polynomials by integration,

$$ P_n(x) = \gamma D e^{V(x)} \int_x^{+\infty} e^{-V(y)} \Psi_n(y) dy. $$

In order to characterise $\Psi_n$, let us ask the following question. Given a polynomial $Q \in \mathcal{P}^{\mathbb{C}}$, under what conditions will

$$ e^{V(x)} \int_x^{+\infty} e^{-V(y)} Q(y) dy $$

also be a polynomial? Note that this is equivalent to the question of whether $Q$ can be written $Q(x) = e^{V(x)} \frac{d}{dx}(P(x)e^{-V(x)})$ for some polynomial $P$. Let us now proceed to answer this question.

**Remark 2.3.** Note that the map $P \mapsto e^{V(x)} \frac{d}{dx}(Pe^{-V})$ which maps $\mathcal{P}^{\mathbb{C}} \to \mathcal{P}^{\mathbb{C}}$ is linear and injective. Thus if $Q = e^{V(x)} \frac{d}{dx}(Pe^{-V})$ then this $P$ is unique.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** This figure depicts the contours $C_0, \dots, C_{D-1}$ in the case $D = 8$. The arrows denote the orientation of the contours. Note the outgoing orientation of $C_0$.</figcaption></figure>

Consider the contours

$$C_k \stackrel{\text{def}}{=} e^{2\pi i \frac{k}{D}} [0, +\infty), \quad k \in [0, D-1]. \qquad (2.1)$$

Let us orient these contours so $C_0 = [0, +\infty)$ has left to right orientation, whereas the remaining contours $C_1, \dots, C_{D-1}$ have orientation towards the origin. Then let

$$\Gamma_k \stackrel{\text{def}}{=} C_0 \cup C_k, \quad k \in [1, D-1]. \qquad (2.2)$$

For example, when $k = \frac{D}{2}$ (recall $D$ is even and positive), $\Gamma_k = \mathbb{R}$ with standard left to right orientation. These contours are depicted in Figure 1 for the case $D = 8$.

**Proposition 2.4.** Let $Q \in \mathcal{P}^C$ be a polynomial, then $Q$ can written as

$$Q(x) = e^{V(x)} \frac{d}{dx} (e^{-V(x)} P(x)) = P'(x) - V'(x) P(x)$$

for some polynomial $P \in \mathcal{P}^C$ if and only if the following integrals vanish:

$$\int_{\Gamma_k} Q(x) e^{-V(x)} dx = 0, \quad \forall k \in [1, D-1] \qquad (2.3)$$

**Proof.** “Only if” is trivial since for $z \in C_k \setminus \{0\}$, $z^D > 0$, hence $P(z)e^{-V(z)}$ is rapidly decreasing, hence

$$\int_{\Gamma_k} Q(x) e^{-V(x)} dz = \int_{\Gamma_k} \frac{d}{dx} (P(x) e^{-V(x)}) dx = 0, \quad \forall k \in [1, D-1].$$

For the converse claim, suppose the integrals (2.3) all vanish. Let us then divide $D$ into $\text{deg} Q + 1$ with remainder, so that $\text{deg} Q+1 = D\ell+r$, where $\ell, r \in \mathbb{N}$ and $r \le D-1$. We will construct a certain polynomial $R$ as follows. If $\text{deg} Q \le D-2$, then $R := Q$, so we may assume $\text{deg} Q \ge D-1$. We can then factorise

$$Q(x) = F_1(x)V'(x) + G_1(x),$$
where $F_1 \in \mathcal{P}^\mathbb{C}$ has degree $\deg Q - D + 1$ and $G_1 \in \mathcal{P}^\mathbb{C}$ has degree at most $D - 2$. We can then
repeat the procedure,

$$
\begin{align*}
F'_1(x) &= F_2(x)V'(x) + G_2(x), & F'_2(x) &= F_3(x)V'(x) + G_3(x), & \dots, \\
F'_{\ell-1}(x) &= F_\ell(x)V'(x) + G_\ell(x).
\end{align*}
$$

$F_j \in \mathcal{P}^\mathbb{C}$ has degree $\deg Q - Dj + 1$ and $G_j \in \mathcal{P}^\mathbb{C}$ has degree at most $D - 2$ for $j = 1, \dots, l$. In
particular $F_l$ has degree $r$ and so the process terminates at this stage. Let

$$
R(x) := G_1(x) + \dots + G_l(x) + F'_l(x).
$$

$R$ has degree at most $D - 2$.

By repeated integration by parts, we see that

$$
\int_{\Gamma_k} R(x) e^{-V(x)} dx = 0, \quad \forall k \in [1, D-1].
$$

We now claim that $R \equiv 0$. Consider the function

$$
C(z) \stackrel{\text{def}}{=} e^{V(z)} \int_z^{+\infty} R(x) e^{-V(x)} dx,
$$

where the integration path is from $z \in \mathbb{C}$ to 0 then from 0 to $+\infty$. $C$ is clearly an entire function.
We claim that as $|z| \to +\infty$, $C(z) \to 0$, uniformly in arg $z$. Hence by Liouville's theorem $C \equiv 0$,
and hence $R \equiv 0$.

From our assumptions,

$$
\int_z^{+\infty} R(x) e^{-V(x)} dx = \int_z^{(+\infty) \cdot e^{2\pi i \frac{k}{D}}} R(x) e^{-V(x)} dx, \quad \forall k \in [1, D-1].
$$

Now introduce the function $\varphi(z) = \tilde{V}(z)^{\frac{1}{D}}$ (where $\tilde{V}(z) = \frac{1}{\gamma}V(z)$) which is holomorphic and injective so long as $|z|$ is sufficiently large, i.e., on a neighbourhood of (complex) infinity. Let us choose $k \in \{0, \dots, D-1\}$, so that

$$
2\pi \frac{k - \frac{1}{2}}{D} \leq \arg \varphi(z) \leq 2\pi \frac{k + \frac{1}{2}}{D}
$$

(this $k$ may not be unique but this does not matter). Let us also define $\arg$ so that there is no
branch cut through this sector. Finally, let us choose the contour such that $\text{Re}(V(z))$ is always
non-decreasing along the contour. The existence of such a contour is proven in Lemma 2.5.
Integrating by parts, we find

$$
\int_z^{(+\infty) \cdot e^{2\pi i \frac{k}{D}}} R(x) e^{-V(x)} dx = \frac{R(z)}{V'(z)} e^{-V(z)} + \int_z^{(+\infty) \cdot e^{2\pi i \frac{k}{D}}} \frac{d}{dx} \left( \frac{R(x)}{V'(x)} \right) e^{-V(x)} dx.
$$

$\frac{d}{dx} \left( \frac{R(x)}{V'(x)} \right)$ is a rational function which scales like $x^{-2}$ and so for $|x|$ sufficiently large there is a constant $\kappa > 0$ such that $|\frac{d}{dx} \left( \frac{R(x)}{V'(x)} \right)| \le \kappa |x|^{-2}$. Thus

$$
\left| e^{V(z)} \int_z^{(+\infty) \cdot e^{2\pi i \frac{k}{D}}} R(x) e^{-V(x)} dx \right| \le \frac{|R(z)|}{|V'(z)|} + \kappa \int_z^{(+\infty) \cdot e^{2\pi i \frac{k}{D}}} |x|^{-2} dx \to 0
$$

as $z \to \infty$. This limit holds so long as $z$ is confined to the sector $2\pi \frac{k-1}{D} \le \arg \varphi(z) \le 2\pi \frac{k+1}{D}$,
however since there are only finitely many such sectors this limit thus holds uniformly in $\arg \varphi(z)$
(and so uniformly in $\arg z$). Thus $R \equiv 0$.
Finally, this shows that $e^{V(z)} \int_{z}^{+\infty} Q(x)e^{-V(x)}dx$ is a polynomial, since if we factorise $Q$ and integrate by parts, we find

$$e^{V(z)} \int_{z}^{+\infty} Q(x)e^{-V(x)}dx = F_{1}(z) + \cdots + F_{l}(z),$$

which is a polynomial of degree $\deg Q - D + 1$.

**Lemma 2.5.** Let $\tilde{V}(w) = \frac{1}{\gamma}V(w) = w^D + \mathcal{O}(w^{D-1})$ and let $\varphi(w) := \tilde{V}(w)^{\frac{1}{D}}$. There exists a $\varrho > 0$ sufficiently large such that $\varphi$ is analytic and injective on

$$U_{\varrho} := \{z \in \mathbb{C} : |z| > \varrho\}$$

and such that the following holds. For each $z \in U_{\varrho}$ such that $\frac{2\pi(k-\frac{1}{2})}{D} \le \arg \varphi(z) \le \frac{2\pi(k+\frac{1}{2})}{D}$, there is a of contour $\tau_z$ starting at $z$ and tending to $(+\infty) \cdot e^{2\pi i \frac{k}{D}}$ such that $\text{Re}(V(w))$ is non-decreasing along the contour and $\int_{\tau_z} |x|^{-2}|dx| \to 0$ as $|z| \to \infty$ uniformly in $\arg z$.

**Proof.** Let

$$\tilde{\tau}(s) = \begin{cases} |\varphi(z)| \exp \left( i(1-s) \arg \varphi(z) + i \frac{2\pi k}{D} \right), & s \in [0, 1], \\ \left( |\varphi(z)| + s - 1 \right) \exp \left( i \frac{2\pi k}{D} \right), & s \in [1, +\infty). \end{cases}$$

Then we choose our contour to be the preimage of $\tilde{\tau}$, that is $\tau = \varphi^{-1} \circ \tilde{\tau}$. From this, we see that $V(\tau(s)) = \gamma \tilde{\tau}(s)^D$. A short calculation then shows that for $s \in (0, 1)$

$$\frac{d}{ds} \text{Re}(V(\tau(s))) = \gamma |\varphi(z)|^D D \left( \arg \varphi(z) - \frac{2\pi k}{D} \right) \sin \left( D(1-s) \left( \arg \varphi(z) - \frac{2\pi k}{D} \right) \right) \ge 0.$$

It is also easy to show for $s \in (1, +\infty)$, $\frac{d}{ds}(\text{Re} \tilde{\tau}(s)^D) \ge 0$. Furthermore, $\varphi(w) = w(1 + o(1))$ and $\varphi'(w) = 1 + o(1)$ uniformly on $U_{\varrho}$ as $\varrho \to +\infty$, and so $\int_{\tau} |x|^{-2}|dx| \to 0$ as $|z| \to \infty$ uniformly in $\arg z$.

A central result in the theory of orthogonal polynomials is that their zeros are simple and lie in the support of the weight function. The following proposition gives a partial extension of this to dual skew-orthogonal polynomials.

**Proposition 2.6.** $\Psi_{2n}$ has at least $2n + 1$ real zeros (not counting multiplicity), of which at least $2n + 2 - \frac{D}{2}$ must be simple. $\Psi_{2n+1}$ has at least $2n$ real zeros (not counting multiplicity), of which at least $2n - \frac{D}{2}$ must be simple.

**Proof.** Let $x_1, \dots, x_\ell$ be the collection of real zeros of $\Psi_{2n}$ with odd multiplicity, i.e., exactly those points on $\mathbb{R}$ at which $\Psi_{2n}$ changes sign, where we do not repeat according to multiplicity. Then $\Psi_{2n}(x)(x - x_1) \cdots (x - x_\ell)$ has the same sign everywhere and vanishes only on a finite set, hence

$$\int_{\mathbb{R}} \Psi_{2n}(x)(x - x_1) \cdots (x - x_\ell) e^{-2V(x)} dx \ne 0.$$

However, if $\ell \le 2n$ then the above integral vanishes, thus $\ell \ge 2n + 1$. Furthermore, it is easily seen that at least $2n + 2 - \frac{D}{2}$ of these must be simple. A similar argument works for $\Psi_{2n+1}$. ■

We now arrive at a pair of Riemann-Hilbert problems, for the even and odd dual skew-orthogonal polynomials, respectively. In what follows, it is more convenient to ensure the uniqueness of the sequence by fixing the next-to-leading coefficient of $\Psi_{2n+1}$ to be zero rather than $P_{2n+1}$. We begin with the “even problem”, so called because the $(1, 1)$ matrix element of the solution will turn out to be $\Psi_{2n}$.
**Riemann–Hilbert Problem 2.7** (even problem). Let $V$ be a real polynomial of even degree $D \ge 2$ and positive leading coefficient, $\Gamma_j$ be defined by (2.1) and (2.2), and

$$ \Gamma = \bigcup_{j=1}^{D-1} \Gamma_j \qquad (2.4) $$

with the orientations of $\Gamma_j$ inherited by $\Gamma$ (see, e.g., Figure 1). We look for a $(D+1) \times (D+1)$ matrix valued function

$$ A_n: \mathbb{C} \setminus \Gamma \to \mathcal{M}_{D+1}(\mathbb{C}) $$

such that

(1) $A_n$ is analytic on $\mathbb{C} \setminus \Gamma$ (that is, its matrix elements are analytic on $\mathbb{C} \setminus \Gamma$).

(2) $A_n$ has continuous boundary values up to $\Gamma \setminus \{0\}$. That is, given $x \in \Gamma \setminus \{0\}$, as $\mathbb{C} \setminus \Gamma \ni z \to x$ from either the + (left) side or – (right) side non-tangentially, the limit of $A_n(z)$ exists. These limits are denoted

$$ A_n^+ : \Gamma \setminus \{0\} \to \mathcal{M}_{D+1}(\mathbb{C}), \quad A_n^- : \Gamma \setminus \{0\} \to \mathcal{M}_{D+1}(\mathbb{C}) $$

and are continuous functions. “Left” and “right” mean relative to the orientation of the contour. The boundary values are related by the “jump matrix”

$$ A_n^+(x) = A_n^-(x) \begin{pmatrix} 1 & e^{-2V} \chi_{\mathbb{R}} & e^{-V} \chi_{\Gamma_1} & \dots & e^{-V} \chi_{\Gamma_{D-1}} \\ & 1 & & & \\ & & 1 & & \\ & & & \ddots & \\ & & & & 1 \end{pmatrix}, \quad x \in \Gamma \setminus \{0\}. $$

(3) $A_n$ is bounded in a neighbourhood of 0.

(4) We have the asymptotic normalisation

$$ A_n(z) = (\mathbb{I} + \mathcal{O}(z^{-1})) \begin{pmatrix} z^{2n+D-1} & & & \\ & z^{-2n} & & \\ & & z^{-1} & \\ & & & \ddots \\ & & & & z^{-1} \end{pmatrix}, \quad z \to \infty. $$

**Lemma 2.8.** If a solution to Riemann–Hilbert Problem 2.7 exists, then it is unique.

**Proof.** This method is standard but we repeat it for completeness. We note that the jump matrix has determinant 1 so $\det A_n(z)$ has no jump across $\Gamma \setminus \{0\}$ and takes continuous boundary values, so by Morera's theorem is analytic everywhere except possibly at 0. However, the requirement of boundedness in a neighbourhood of 0 implies that the singularity at 0 is removable. Hence $z \mapsto \det A_n(z)$ is entire. Furthermore, we have the normalisation $\det A_n(z) = 1 + \mathcal{O}(z^{-1})$, as $z \to \infty$ and so $\det A_n(z) \equiv 1$. Then let $\tilde{A}_n$ be another solution. Then the function $z \mapsto \tilde{A}_n(z)A_n^{-1}(z)$ has no jump across $\Gamma \setminus \{0\}$ and obtains continuous boundary values, so by Morera's theorem is analytic everywhere except possibly at 0. Again by the requirement of boundedness of $\tilde{A}_n$ and $A_n$ in a neighbourhood of 0 and that $\det A_n(z) \equiv 1$, we see that the singularity at 0 is removable. Hence $z \mapsto \tilde{A}_n(z)A_n^{-1}(z)$ is actually entire. Finally, since

$$ \tilde{A}_n(z) A_n^{-1}(z) = \mathbb{I} + \mathcal{O}(z^{-1}), $$

we find that $\tilde{A}_n(z)A_n^{-1}(z) \equiv \mathbb{I}$.

<figure><img src="image_4.png" /></figure>
The jump condition implies that $(A_n)_{11}$ is an entire function which scales like $z^{2n+D-1}$ as $z \to \infty$. It is therefore a monic polynomial of degree $2n + D - 1$ by Liouville's theorem,

$$(A_n)_{11} := Q \in \mathcal{P}_{2n+D-1}^{\mathbb{C}}.$$

The Plemelj formula then implies that the first row of $A_n$ is

$$(Q, C_{\mathbb{R}}(e^{-2V}Q), C_{\Gamma_1}(e^{-V}Q), \dots, C_{\Gamma_{D-1}}(e^{-V}Q)).$$

From the asymptotic condition

$$C_{\Gamma_j}(e^{-V}Q)(z) = \mathcal{O}(z^{-2}), \quad \forall j \in [1, D-1],$$

we see, by Proposition 2.4, that there exists a unique $P \in \mathcal{P}^{\mathbb{C}}$ of degree exactly $2n$ such that $Q = e^V \frac{d}{dx}(Pe^{-V})$. Finally, the normalisation

$$C_{\mathbb{R}}(e^{-2V}Q)(z) = \mathcal{O}(z^{-2n-1}), \quad z \to \infty$$

implies that $\langle X^k, P \rangle_4 = 0$ for $k \in [0, 2n-1]$, and so $P$ is skew-orthogonal. Thus $Q = \Psi_{2n}$.

For the second row, we see that $(A_n)_{21} := \tilde{Q}$ is a polynomial of degree at most $2n + D - 2$. Following the same procedure as above, we see that the second row is

$$(\tilde{Q}, C_{\mathbb{R}}(e^{-2V}\tilde{Q}), C_{\Gamma_1}(e^{-V}\tilde{Q}), \dots, C_{\Gamma_{D-1}}(e^{-V}\tilde{Q})),$$

where now $\tilde{Q} = -\frac{2\pi i \gamma D}{h_{n-2}} \Psi_{2n-2}$. For the remaining rows, we see that $(A_n)_{j+2,1} := R_n^{(j)} \in \mathcal{P}_{2n+D-2}^{\mathbb{C}}$ must be a polynomial of degree at most $2n + D - 2$. To show there exists a unique solution to Riemann-Hilbert Problem 2.7, let us introduce the following linear functionals:

$$\alpha_j: \mathcal{P}_{2n+D-2}^{\mathbb{C}} \to \mathbb{C}, \quad \beta_j: \mathcal{P}_{2n+D-2}^{\mathbb{C}} \to \mathbb{C},$$

where

$$\alpha_j(P) \stackrel{\text{def}}{=} \int_{\mathbb{R}} e^{-2V(x)} x^j P(x) dx, \quad j \in [0, 2n-1], \qquad (2.5)$$

$$\beta_j(P) \stackrel{\text{def}}{=} -\frac{1}{2\pi i} \int_{\Gamma_j} e^{-V(x)} P(x) dx, \quad j \in [1, D-1]. \qquad (2.6)$$

**Lemma 2.9.** The linear functionals $\alpha_0, \dots, \alpha_{2n-1}, \beta_1, \dots, \beta_{D-1}$ form a linearly independent set (and thus span $(\mathcal{P}_{2n+D-2}^{\mathbb{C}})^*$).

**Proof.** To show this, we need only show that

$$\bigcap_{j=0}^{2n-1} \ker \alpha_j \cap \bigcap_{k=1}^{D-1} \ker \beta_k = 0.$$

Thus let $Q \in \mathcal{P}_{2n+D-2}^{\mathbb{C}}$ be such that $\alpha_0(Q) = \dots = \alpha_{2n-1}(Q) = \beta_1(Q) = \dots = \beta_{D-1}(Q) = 0$. Thus by Proposition 2.4, $Q = e^V \frac{d}{dx}(Pe^{-V})$ for some $P \in \mathcal{P}_{2n-1}^{\mathbb{C}}$. Furthermore, $\langle X^0, P \rangle_4 = \dots = \langle X^{2n-1}, P \rangle_4 = 0$. Hence $P$ is skew-orthogonal to all of $\mathcal{P}_{2n-1}^{\mathbb{C}}$. Thus by the non-degeneracy of the inner product $\langle \cdot, \cdot \rangle_4$, $P \equiv 0$, and hence $Q \equiv 0$. ■

Let $R_n^{(j)} \in \mathcal{P}_{2n+D-2}^{\mathbb{C}}$ be the unique polynomial such that $\alpha_0(R_n^{(j)}) = \dots = \alpha_{2n-1}(R_n^{(j)}) = 0$ and

$$\beta_k(R_n^{(j)}) = \delta_{kj}, \quad \forall k \in [1, D-1].$$

That such a polynomial exists and is unique follows from the linear independence of these functionals.
**Theorem 2.10.** The unique solution of the even Riemann-Hilbert Problem 2.7 is

$$A_n = \begin{pmatrix} \Psi_{2n} & C_{\mathbb{R}}(e^{-2V}\Psi_{2n}) & C_{\Gamma_1}(e^{-V}\Psi_{2n}) & \cdots & C_{\Gamma_{D-1}}(e^{-V}\Psi_{2n}) \\ -\frac{2\pi i \gamma D}{h_{n-1}}\Psi_{2n-2} & -\frac{2\pi i \gamma D}{h_{n-1}}C_{\mathbb{R}}(e^{-2V}\Psi_{2n-2}) & -\frac{2\pi i \gamma D}{h_{n-1}}C_{\Gamma_1}(e^{-V}\Psi_{2n-2}) & \cdots & -\frac{2\pi i \gamma D}{h_{n-1}}C_{\Gamma_{D-1}}(e^{-V}\Psi_{2n-2}) \\ R_n^{(1)} & C_{\mathbb{R}}(e^{-2V}R_n^{(1)}) & C_{\Gamma_1}(e^{-V}R_n^{(1)}) & \cdots & C_{\Gamma_{D-1}}(e^{-V}R_n^{(1)}) \\ \vdots & & & & \vdots \\ R_n^{(D-1)} & C_{\mathbb{R}}(e^{-2V}R_n^{(D-1)}) & C_{\Gamma_1}(e^{-V}R_n^{(D-1)}) & \cdots & C_{\Gamma_{D-1}}(e^{-V}R_n^{(D-1)}) \end{pmatrix}$$

**Lemma 2.11** (symmetry for even case). Let $P \in \mathcal{M}_{D-1}(\mathbb{R})$ be the permutation matrix on $\mathbb{R}^{D-1}$ for which $Pe_k = e_{D-k}$ for $k = 1, \dots, D-1$. Here $e_k = (0, \dots, 0, 1, 0, \dots, 0)^T$, where the 1 lies in the k-th entry. Let

$$Q := \operatorname{diag}(-1, 1) \oplus P = \begin{pmatrix} -1 & 0 \\ 0 & 1 \end{pmatrix} \in \mathcal{M}_{D+1}(\mathbb{R}).$$

Note that $Q = Q^T = Q^{-1}$ is an orthogonal matrix. Then we have the relation $A_n(z) = \overline{QA_n(z)Q}$.

**Proof.** Let $\tilde{A}_n(z) = \overline{QA_n(z)Q}$. We show that $\tilde{A}_n$ solves the even Riemann-Hilbert problem and thus by uniqueness (see Lemma 2.8) $\tilde{A}_n = A_n$. Clearly, $\tilde{A}_n$ is analytic on $\mathbb{C} \setminus \Gamma$ and it is easy to see that

$$\tilde{A}_n(z) = (\mathbb{I} + \mathcal{O}(z^{-1})) \begin{pmatrix} z^{2n+D-1} & & & \\ & z^{-2n} & & \\ & & z^{-1} & \\ & & & z^{-1} \end{pmatrix}, \quad z \to \infty.$$

Furthermore, $\tilde{A}_n$ is clearly bounded at 0 since $A_n$ is.

Because $V$ is a polynomial with real coefficients, $\overline{V(\bar{x})} = V(x)$. Furthermore, $\chi_{\mathbb{R}}(\bar{x}) = \chi_{\mathbb{R}}(x)$ and $\chi_{\Gamma_k}(\bar{x}) = \chi_{\Gamma_{D-k}}(x)$. Now for $x \in \Gamma \setminus \{0\}$, we have

$$\begin{align*} (\tilde{A}_n)_+(x) &= \overline{Q(A_n)_-(\bar{x})Q} \\ &= \overline{Q(A_n)_+(\bar{x})} \begin{pmatrix} 1 & -e^{-2V}\chi_{\mathbb{R}} & -e^{-V}\chi_{\Gamma_{D-1}} & \cdots & -e^{-V}\chi_{\Gamma_1} \\ & 1 & & & \\ & & 1 & & \\ & & & \ddots & \\ & & & & 1 \end{pmatrix} Q \\ &= (\tilde{A}_n)_-(x) \begin{pmatrix} 1 & e^{-2V}\chi_{\mathbb{R}} & e^{-V}\chi_{\Gamma_1} & \cdots & e^{-V}\chi_{\Gamma_{D-1}} \\ & 1 & & & \\ & & 1 & & \\ & & & \ddots & \\ & & & & 1 \end{pmatrix}. \end{align*}$$

Thus $\tilde{A}_n = A_n$.

<figure><img src="image_5.png" /></figure>

**Remark 2.12.** Our method assumes $V$ is a polynomial, however there are interesting cases, e.g., Jacobi weights, that do not fit into this class. How might our analysis be extended? The natural way would be to redefine $\Psi_n$ and $R_n^{(j)}$ so that these are no longer polynomials. For example, if $V'$ was rational with simple poles at $\lambda_1, \dots, \lambda_k$, then $\Psi_n$ would also be rational with simple poles, with the property that $\mathrm{Res}_{z=\lambda_i} \Psi_n(z)e^{-V(z)} = 0$ for all $i = 1, \dots, k$ and similarly for the auxiliary functions $R_n^{(j)}$. Investigation of this is left as a topic for future research.
We now move onto the “odd problem”, so called because it will later turn out that the (1, 1) matrix element of the solution is $\Psi_{2n+1}$.

**Riemann–Hilbert Problem 2.13** (odd problem). Let $V$ be a real polynomial of even degree $D \ge 2$ and positive leading coefficient, and let

$$ \Sigma = \mathbb{T} \cup \bigcup_{j=1}^{D-1} \Gamma_j, $$

where $\Gamma_j$ have the specified orientations and $\mathbb{T}$ is positively oriented (see, e.g., Figure 2). Let $\omega = e^{\frac{2\pi i}{D}}$ be a primitive root of unity and let

$$ S = \{0, 1, \omega, \omega^2, \dots, \omega^{D-1}\} $$

be the points of self-intersection of $\Sigma$. We look for a $(D+2) \times (D+2)$ matrix valued function

$$ B_n: \mathbb{C} \setminus \Sigma \to \mathcal{M}_{D+2}(\mathbb{C}) $$

such that

(1) $B_n$ is analytic on $\mathbb{C} \setminus \Sigma$ (that is, its matrix elements are analytic on $\mathbb{C} \setminus \Sigma$).

(2) $B_n$ has continuous boundary values up to $\Sigma \setminus S$. That is, given $x \in \Sigma \setminus S$, as $\mathbb{C} \setminus \Sigma \ni z \to x$ from either the + (left) side or – (right) side non-tangentially, the limit of $B_n(z)$ exists. These limits are denoted

$$ B_n^+ : \Sigma \setminus S \to \mathcal{M}_{D+2}(\mathbb{C}), \quad B_n^- : \Sigma \setminus S \to \mathcal{M}_{D+2}(\mathbb{C}) $$

and are continuous functions. “Left” and “right” mean relative to the orientation of the contour. The boundary values are related by the “jump matrix”

$$ B_n^+(x) = B_n^-(x) \begin{pmatrix} 1 & e^{-2V} \chi_{\mathbb{R}} & e^{-V} \chi_{\Gamma_1} & \dots & e^{-V} \chi_{\Gamma_{D-1}} & X^{-2n-D} \chi_{\mathbb{T}} \\ 1 & & 1 & & & \\ & & & \ddots & & \\ & & & & 1 & \\ & & & & & 1 \\ & & & & & \end{pmatrix}, $$ $$ x \in \Sigma \setminus S. $$

(3) $B_n$ is bounded in a neighbourhood of $S$.

(4) We have the asymptotic normalisation

$$ B_n(z) = (\mathbb{I} + \mathcal{O}(z^{-1})) \begin{pmatrix} z^{2n+D} & & & \\ & z^{-2n} & & \\ & & z^{-1} & \\ & & & \ddots \\ & & & & z^{-1} \end{pmatrix}, \quad z \to \infty. $$

The method of solution of this Riemann–Hilbert problem is very similar to the even case. The asymptotic condition implies that $(B_n)_{11} = Q$ is a monic polynomial of degree $2n + D$. From the jump condition and the Plemelj formula we find that the first row is

$$ (Q, C_{\mathbb{R}}(e^{-2V}Q), C_{\Gamma_1}(e^{-V}Q), \dots, C_{\Gamma_{D-1}}(e^{-V}Q), C_{\mathbb{T}}(X^{-2n-D}Q)). $$

From the normalisation $C_{\Gamma_j}(e^{-V}Q)(z) = \mathcal{O}(z^{-2})$ for all $j \in [1, D-1]$, we find that there exists a polynomial $P$ of degree $2n + 1$ such that $Q = e^V \frac{d}{dx}(P e^{-V})$.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** The contours of the odd Riemann–Hilbert problem, here shown for the case $D = 8$. Orientations of the contours are indicated by arrows.</figcaption></figure>

**Remark 2.14.** The scaling $C_T(X^{-2n-D}Q)(z) = O(z^{-2})$ implies that

$$Q(x) = x^{2n+D} + O(x^{2n+D-2}), \quad z \to \infty.$$

This fixes the next-to-leading coefficient of $P$ (where $Q = e^V \frac{d}{dx}(Pe^{-V})$) though not necessarily to 0. This explains the remark made in the introduction that fixing the next-to-leading coefficient of $P$ to 0 is not always the most convenient choice. The addition of the contour $\mathcal{T}$ could also be done for the Riemann–Hilbert problem of Pierce [38] in order to guarantee a unique solution.

Finally, the scaling $C_R(e^{-2V}Q)(z) = O(z^{-2n-1})$ implies that $\langle X^0, P \rangle_4 = \dots = \langle X^{2n-1}, P \rangle_4 = 0$, so $P$ is indeed skew-orthogonal. Thus $Q = \Psi_{2n+1}$.

In the same manner as for the even case, we find that the second row is

$$(\tilde{Q}, C_R(e^{-2V}\tilde{Q}), C_{\Gamma_1}(e^{-V}\tilde{Q}), \dots, C_{\Gamma_{D-1}}(e^{-V}\tilde{Q}), C_T(X^{-2n-D}\tilde{Q})),$$

where $\tilde{Q} = -\frac{2\pi i \gamma_D}{h_{n-1}} \Psi_{2n-2}$. Indeed, the remaining rows except the final one are the same as in the even case. To show a unique solution in all remaining rows, let us introduce the linear functional $\beta_D: \mathcal{P}_{2n+D-2}^{\mathcal{C}} \to \mathbb{C}$ defined by

$$\beta_D(P) \stackrel{\text{def}}{=} -\frac{1}{2\pi i} \oint_{\mathcal{T}} P(z) z^{-2n-D} dz.$$

Now let $\alpha_0, \dots, \alpha_{2n-1}$ and $\beta_1, \dots, \beta_{D-1}$ act on $\mathcal{P}_{2n+D-1}^{\mathcal{C}}$ according to the formulas given in (2.5) and (2.6). That the 2nd to $(D+2)$-th rows have a unique solution follows from the lemma below.

**Lemma 2.15.** Consider the linear functionals $\alpha_0, \dots, \alpha_{2n-1}$ and $\beta_1, \dots, \beta_D$ viewed as acting on $\mathcal{P}_{2n+D-1}^{\mathcal{C}}$ according to (2.5) and (2.6). Then

$$\bigcap_{j=0}^{2n-1} \ker \alpha_j \cap \bigcap_{k=1}^{D} \ker \beta_k = 0.$$

**Proof.** If $\beta_D(Q) = 0$, then $Q \in \mathcal{P}_{2n+D-2}^{\mathcal{C}}$. The proof then reduces to that of Lemma 2.9. ■
**Theorem 2.16.** The unique solution of the odd Riemann-Hilbert Problem 2.13 is then

$$B_n = \begin{pmatrix} \Psi_{2n+1} & C_{\mathbb{R}}(e^{-2V}\Psi_{2n+1}) & C_{\Gamma_1}(e^{-V}\Psi_{2n+1}) & \cdots & C_{\mathbb{T}}(X^{-2n-D}\Psi_{2n+1}) \\ -\frac{2\pi i \gamma_D}{h_{n-1}}\Psi_{2n-2} & -\frac{2\pi i \gamma_D}{h_{n-1}}C_{\mathbb{R}}(e^{-2V}\Psi_{2n-2}) & -\frac{2\pi i \gamma_D}{h_{n-1}}C_{\Gamma_1}(e^{-V}\Psi_{2n-2}) & \cdots & -\frac{2\pi i \gamma_D}{h_{n-1}}C_{\mathbb{T}}(X^{-2n-D}\Psi_{2n-2}) \\ R_n^{(1)} & C_{\mathbb{R}}(e^{-2V}R_n^{(1)}) & C_{\Gamma_1}(e^{-V}R_n^{(1)}) & \cdots & C_{\mathbb{T}}(X^{-2n-D}R_n^{(1)}) \\ \vdots & & & & \vdots \\ R_n^{(D-1)} & C_{\mathbb{R}}(e^{-2V}R_n^{(D-1)}) & C_{\Gamma_1}(e^{-V}R_n^{(D-1)}) & \cdots & C_{\mathbb{T}}(X^{-2n-D}R_n^{(D-1)}) \\ R_n^{(D)} & C_{\mathbb{R}}(e^{-2V}R_n^{(D)}) & C_{\Gamma_1}(e^{-V}R_n^{(D)}) & \cdots & C_{\mathbb{T}}(X^{-2n-D}R_n^{(D)}) \end{pmatrix}$$

Note that these polynomials $R_n^{(1)}, \dots, R_n^{(D-1)}$ are the same polynomials that appear in the even problem, so we have not made a confusion of notation. Note also that the integrals in the final column can be computed explicitly.

**Proposition 2.17.** $R_n^{(D)} = -\Psi_{2n}$.

**Proof.** From the normalisation, we see that $R_n^{(D)}$ is a polynomial of degree at most $2n + D - 1$, and from $C_{\mathbb{T}}(X^{-2n-D}R_n^{(D)})(z) = z^{-1} + \mathcal{O}(z^{-2})$, we find that $R_n^{(D)}(x) = -Q(x)$, where $Q$ is a monic polynomial of degree $2n + D - 1$. From the fact that

$$\beta_1(R_n^{(D)}) = \cdots = \beta_{D-1}(R_n^{(D)}) = 0,$$

we find that $Q = e^V \frac{d}{dx}(Pe^{-V})$ for some polynomial $P$ of degree exactly $2n$. Finally, the requirement that

$$\alpha_0(R_n^{(D)}) = \cdots = \alpha_{2n-1}(R_n^{(D)}) = 0$$

implies that $P$ is skew-orthogonal to all of $\mathcal{P}_{2n-1}^C$.

<figure><img src="image_7.png" /></figure>

From this, one might think the $B_n$ problem is strictly superior to the $A_n$ problem, since the former encodes both $\Psi_{2n+1}$ and $\Psi_{2n}$, however our Christoffel-Darboux type formula (3.1) implies that from the point of view of random matrix theory it is completely sufficient to study $A_n$.

**Lemma 2.18.** $2n \le \deg R_n^{(j)} \le 2n + D - 2$ for $j \in [1, D-1]$.

**Proof.** Suppose by way of contradiction that $\deg R_n^{(j)} \le 2n - 1$. Then by the conditions $\alpha_0(R_n^{(j)}) = \cdots = \alpha_{2n-1}(R_n^{(j)}) = 0$, we find that $R_n^{(j)} \equiv 0$. But this contradicts $\beta_j(R_n^{(j)}) = 1$. ■

**Remark 2.19.** For the case $n = 0$, the polynomials $R_0^{(1)}, \dots, R_0^{(D-1)}$ can be constructed as follows. Let $M$ be the matrix

$$M = (M_{ij})_{i,j=1}^{D-1} = \left( \int_{\Gamma_i} x^{j-1} e^{-V(x)} dx \right)_{i,j=1}^{D-1},$$

then we have

$$R_0^{(j)}(x) = -2\pi i \sum_{k=1}^{D-1} (M^{-1})_{kj} x^{k-1}, \quad j \in [1, D-1].$$

**Lemma 2.20 (symmetry for odd problem).** Let $P \in \mathcal{M}_{D-1}(\mathbb{R})$ be defined as in Lemma 2.11 and let

$$Q := \operatorname{diag}(-1, 1) \oplus P \oplus (-1) = \left( \begin{array}{c|c} -1 & 0 \\ 0 & 1 \\ \hline & P \\ \hline & -1 \end{array} \right) \in \mathcal{M}_{D+2}(\mathbb{R}).$$

Note that $Q = Q^T = Q^{-1}$ is an orthogonal matrix. Then we have the relation $B_n(z) = Q \overline{B_n(z)} Q$.
**Proof.** Morally, the same as for the even case.

■

Finally, our results allow us to characterise the “dual” skew-orthogonal polynomials $\Psi_{2n}$ and auxiliary polynomials $R_n^{(j)}$ in terms of the corresponding sequence of orthogonal polynomials. This is morally the same construction as in Remark 2.19.

**Corollary 2.21** (skew-orthogonal polynomials in terms of orthogonal polynomials). Let $\{H_{n,j} \}_{j \in \mathbb{N}}$ be the sequence of monic orthogonal polynomials corresponding to the weight $e^{-2V}$, where $V(x) = \gamma x^D + O(x^{D-1})$ is a real polynomial of even degree $D \ge 2$ and $\gamma > 0$. That is,

$$H_n(x) = x^n + O(x^{n-1})$$

and

$$\int_{\mathbb{R}} H_n(x) H_m(x) e^{-2V(x)} dx = 0, \quad \forall n, m \in \mathbb{N}, \quad n \neq m.$$

Let $M$ be the matrix

$$M = (M_{jk})_{j,k=1}^{D-1} = \left( \int_{\Gamma_j} H_{2n+k-1}(z) e^{-V(z)} dz \right)_{j,k=1}^{D-1},$$

and let $\lambda_n$ be the next-to-leading coefficient of $H_{n+D}$, i.e.,

$$H_{n+D}(x) = x^{n+D} + \lambda_n x^{n+D-1} + O(x^{n+D-2}).$$

Then the corresponding “dual” skew-orthogonal polynomials and auxiliary polynomials may be expressed

$$\Psi_{2n} = H_{2n+D-1} - \sum_{k,j=1}^{D-1} \left[ (M^{-1})_{kj} \int_{\Gamma_j} H_{2n+D-1}(z) e^{-V(z)} dz \right] H_{2n+k-1},$$

$$\Psi_{2n+1} = H_{2n+D} - \lambda_n H_{2n+D-1}$$

$$- \sum_{k,j=1}^{D-1} \left[ (M^{-1})_{kj} \int_{\Gamma_j} (H_{2n+D}(z) - \lambda_n H_{2n+D-1}(z)) e^{-V(z)} dz \right] H_{2n+k-1},$$

$$R_n^{(j)} = -2\pi i \sum_{k=1}^{D-1} (M^{-1})_{kj} H_{2n+k-1}.$$

Indeed, these $\{H_k\}_{k \in \mathbb{N}}$ are exactly the “natural basis” identified by Adler and van Moerbeke [3]. A similar construction to Corollary 2.21 in which orthogonal and skew-orthogonal polynomials were related was considered recently in the planar case by [4], however we note that in Corollary 2.21 the relation is really between orthogonal polynomials and dual skew-orthogonal polynomials.

**Example 2.22.** In the case $D = 2$, for the potential $V(z) = \frac{1}{2}z^2$, we can write the corresponding orthogonal polynomials explicitly as

$$H_n(z) = \frac{(-1)^n}{2^n} e^{z^2} \frac{d^n}{dz^n} e^{-z^2}.$$

Then using $\int_{\mathbb{R}} H_{2n}(x) e^{-\frac{1}{2}x^2} dx = \sqrt{2} (n - \frac{1}{2})!$, one finds that

$$\Psi_{2n} = H_{2n+1}, \quad \Psi_{2n+1} = H_{2n+2} - \left(n + \frac{1}{2}\right) H_{2n}, \quad P_{2n} = n! \sum_{k=0}^{n} \frac{1}{k!} H_{2k},$$

$$R_n^{(1)} = - \frac{2\pi i}{\sqrt{2} (n - \frac{1}{2})!} H_{2n}, \quad P_{2n+1} = H_{2n+1}.$$
These formulas for $P_{2n}$ and $P_{2n+1}$ already appeared in the paper of Adler, Forrester, Nagao and van Moerbeke [1].

What we have done in the preceding discussion amounts to representing $\Psi_{2n}$ and $\Psi_{2n+1}$ as “multiple-orthogonal polynomials of Type II” (see, e.g., [34] for a definition of this term). Let us summarise this with the following corollary.

**Corollary 2.23** (dual skew-orthogonal polynomials as Type II multiple-orthogonal polynomials). *Consider the following problem. Find a monic polynomial $Q$ of degree exactly $2n + D - 1$ such that*

$$ \int_{\mathbb{R}} x^k Q(x) e^{-2V(x)} dx = 0, \quad \forall k \in [0, 2n-1], $$

$$ \int_{\Gamma_j} Q(x) e^{-V(x)} dx = 0, \quad \forall j \in [1, D-1]. $$

Then this problem is uniquely solved by $Q = \Psi_{2n}$. Similarly, if we are asked to find a monic polynomial $Q$ of degree exactly $2n + D$ such that

$$ \int_{\mathbb{R}} x^k Q(x) e^{-2V(x)} dx = 0, \quad \forall k \in [0, 2n-1], $$

$$ \int_{\Gamma_j} Q(x) e^{-V(x)} dx = 0, \quad \forall j \in [1, D-1], $$

$$ \int_{\mathbb{T}} Q(x) x^{-2n-D} dx = 0, $$

then this problem is uniquely solved by $Q = \Psi_{2n+1}$. We should note that our case differs from usual multiple-orthogonality because our “measures” are not all supported on $\mathbb{R}$.

Following the standard theory of multiple-orthogonal polynomials [43], it is natural to consider the Riemann–Hilbert problems satisfied by $\tilde{A}_n := A_n^{-\mathrm{T}}$ and $\tilde{B}_n := B_n^{-\mathrm{T}}$ (inverse transpose).

**Riemann–Hilbert Problem 2.24** (dual even problem). *By virtue of the fact that $\det A_n \equiv 1$, $\tilde{A}_n \stackrel{\mathrm{def}}{=} A_n^{-\mathrm{T}}$ satisfies the following Riemann–Hilbert problem:*

(1) $\tilde{A}_n$ is analytic on $\mathbb{C} \setminus \Gamma$ (that is, its matrix elements are analytic on $\mathbb{C} \setminus \Gamma$).

(2) $\tilde{A}_n$ has continuous boundary values up to $\Gamma \setminus \{0\}$. That is, given $x \in \Gamma \setminus \{0\}$, as $\mathbb{C} \setminus \Gamma \ni z \to x$ from either the + (left) side or – (right) side non-tangentially, the limit of $\tilde{A}_n(z)$ exists. These limits are denoted

$$ \tilde{A}_n^+ : \Gamma \setminus \{0\} \to \mathcal{M}_{D+1}(\mathbb{C}), \quad \tilde{A}_n^- : \Gamma \setminus \{0\} \to \mathcal{M}_{D+1}(\mathbb{C}) $$

and are continuous functions. “Left” and “right” mean relative to the orientation of the contour. The boundary values are related by the “jump matrix”,

$$ \tilde{A}_n^+(x) = \tilde{A}_n^-(x) \begin{pmatrix} 1 & & & \\ -e^{-2V} \chi_{\mathbb{R}} & 1 & & \\ -e^{-V} \chi_{\Gamma_1} & & 1 & \\ \vdots & & \ddots & \\ -e^{-V} \chi_{\Gamma_{D-1}} & & & 1 \end{pmatrix}, \quad x \in \Gamma \setminus \{0\}. $$

(3) $\tilde{A}_n$ is bounded in a neighbourhood of 0.
(4) We have the asymptotic normalisation

$$ \widehat{A}_n(z) = (\mathbb{I} + \mathcal{O}(z^{-1})) \begin{pmatrix} z^{-2n-D+1} & & & \\ & z^{2n} & & \\ & & z & \\ & & & z \\ & & & z \end{pmatrix}, \quad z \to \infty. $$

In what follows we shall only be interested in computing the second column of $\widehat{A}_n$ (the motivation for this is that this column will appear in the upcoming Christoffel-Darboux-type formula), though it is an easy exercise for the reader to compute the remaining columns. By virtue of the jump condition, the 2nd to $(D+1)$-th columns all have no jump across $\Gamma \setminus \{0\}$ and by boundedness around 0 we see that the singularity at 0 is removable. Hence the 2nd to $(D+1)$-th columns are entire functions, indeed polynomials by the asymptotic condition,

$$ \begin{aligned} (\widehat{A}_n)_{i2}(z) &= z^{2n}\delta_{i,2} + \mathcal{O}(z^{2n-1}), && i \in [1, D+1], \\ (\widehat{A}_n)_{ij}(z) &= z\delta_{ij} + c_{ij}, && i \in [1, D+1] \quad \text{and} \quad j \in [3, D+1], \end{aligned} $$

where the $c_{ij}$ are constants. The jump condition on the first column may be solved by Plemelj's formula

$$ (\widehat{A}_n)_{i1} = -C_{\mathbb{R}}((\widehat{A}_n)_{i2}e^{-2V}) - \sum_{j=3}^{D+1} c_{ij}C_{\Gamma_{j-2}}(e^{-V}) - \chi_{i \ge 3}C_{\Gamma_{i-2}}(Xe^{-V}), \quad i \in [1, D+1]. $$

The asymptotic condition $(\widehat{A}_n)_{11}(z) = z^{-2n-D+1}(1 + \mathcal{O}(z^{-1}))$ yields the conditions

$$ \int_{\mathbb{R}} (\widehat{A}_n)_{12}(x) x^k e^{-2V(x)} dx + \sum_{j=3}^{D+1} c_{1j} \int_{\Gamma_{j-2}} x^k e^{-V(x)} dx = 0, \quad \forall k \in [0, 2n + D - 3]. $$

By taking linear combinations of these equations, we can replace $x^k$ with any polynomial of degree $\le 2n+D-3$; in particular, we can replace it with $\Psi_0, \dots, \Psi_{2n-2}$. This implies that $(\widehat{A}_n)_{12}$ is a scalar multiple of $P_{2n-2}$. This multiple is fixed by

$$ \int_{\mathbb{R}} (\widehat{A}_n)_{12}(x) \Psi_{2n-1}(x) e^{-2V(x)} dx = 2\pi i, $$

and so

$$ (\widehat{A}_n)_{12} = -\frac{2\pi i \gamma D}{h_{n-1}} P_{2n-2}. \qquad (2.7) $$

By a similar argument, we see that

$$ (\widehat{A}_n)_{22} = P_{2n}. \qquad (2.8) $$

The remaining entries in the second column are polynomials of degree at most $2n-1$. These satisfy the conditions

$$ \begin{aligned} & \int_{\mathbb{R}} (\widehat{A}_n)_{i2}(x) \Psi_k(x) e^{-2V(x)} dx + \int_{\Gamma_{i-2}} x \Psi_k(x) e^{-V(x)} dx = 0, \\ \forall k \in [0, 2n-1] \quad \text{and} \quad i \ge 3. \end{aligned} $$
Expanding this in the basis $P_0, \dots, P_{2n-1}$, we find

$$ (\widehat{A}_n)_{i2} = \gamma D \sum_{k=0}^{n-1} h_k^{-1} \left( P_{2k} \int_{\Gamma_{i-2}} x \Psi_{2k+1}(x) e^{-V(x)} dx - P_{2k+1} \int_{\Gamma_{i-2}} x \Psi_{2k}(x) e^{-V(x)} dx \right) \quad (2.9) \\ i \in [3, D+1]. $$

What the above calculation demonstrates is that $P_{2n}$ may be thought of as a Type I multiple-orthogonal polynomial (again see [34] for a definition of this term). Let us formalise this with the following corollary.

**Corollary 2.25** (even skew-orthogonal polynomials as Type I multiple-orthogonal polynomials). *Consider the following problem. Find a monic polynomial $P$ of degree exactly $2n$ and a collection of constants $c_1, \dots, c_{D-1}$ such that*

$$ \int_{\mathbb{R}} x^k P(x) e^{-2V(x)} dx + \sum_{j=1}^{D-1} c_j \int_{\Gamma_j} x^k e^{-V(x)} dx = 0, \quad \forall k \in [0, 2n + D - 2]. $$

Then this problem has a unique solution $(P, c_1, \dots, c_{D-1})$, where $P = P_{2n}$.

It is curious that Types I and II multiple orthogonality both make their appearance here (a similar situation, in which both Types I and II multiple orthogonality appear in the study of the same ensemble, was observed recently in certain planar ensembles [5, 33]). On the other hand, it is not clear if it is possible to think of $P_{2n+1}$ as a Type I multiple-orthogonal polynomial and, as we shall soon see, $P_{2n+1}$ (or any multiple of it) does not appear as a matrix element of $B_n^{-1}$. It is possible to characterise $P_{2n+1}$ by a combination of both Types I and II orthogonality conditions, as follows. We look for a monic polynomial $P$ of degree $2n+1$ and a collection of constants $c_1, \dots, c_{D-1}$ such that

$$ \int_{\mathbb{R}} P(x) x^k e^{-2V(x)} dx + \sum_{j=1}^{D-1} c_j \int_{\Gamma_j} x^k e^{-V(x)} dx = 0, \quad \forall k \in [0, 2n + D - 2] $$

and

$$ \int_{\mathbb{T}} P(x) e^{-V(x)} \frac{d}{dx} (x^{-2n-D} e^{V(x)}) dx = 0. $$

$(P, c_1, \dots, c_{D-1}) = (P_{2n+1}, \dots, \dots)$ is the unique solution of this problem. However, this construction seems somewhat ugly and we shall not investigate the implications of this here.

**Riemann–Hilbert Problem 2.26** (dual odd problem). By virtue of the fact that $\det B_n \equiv 1$, $\widehat{B}_n \stackrel{\text{def}}{=} B_n^{-\top}$ satisfies the following Riemann–Hilbert problem:

(1) $\widehat{B}_n$ is analytic on $\mathbb{C} \setminus \Sigma$ (that is, its matrix elements are analytic on $\mathbb{C} \setminus \Sigma$).

(2) $\widehat{B}_n$ has continuous boundary values up to $\Sigma \setminus S$. That is, given $x \in \Sigma \setminus S$, as $\mathbb{C} \setminus \Sigma \ni z \to x$ from either the + (left) side or – (right) side non-tangentially, the limit of $\widehat{B}_n(z)$ exists. These limits are denoted

$$ \widehat{B}_n^+ : \Sigma \setminus S \to \mathcal{M}_{D+2}(\mathbb{C}), \quad \widehat{B}_n^- : \Sigma \setminus S \to \mathcal{M}_{D+2}(\mathbb{C}) $$

and are continuous functions. “Left” and “right” mean relative to the orientation of the contour. The boundary values are related by the “jump matrix”

$$ \widehat{B}_n^+(x) = \widehat{B}_n^-(x) \begin{pmatrix} 1 & & & \\ -e^{-2V} \chi_{\mathbb{R}} & 1 & & \\ -e^{-V} \chi_{\Gamma_1} & & 1 & \\ \vdots & & \ddots & \\ -e^{-V} \chi_{\Gamma_{D-1}} & & & 1 \\ -X^{-2n-D} \chi_{\mathbb{T}} & & & 1 \end{pmatrix}, \quad x \in \Sigma \setminus S. $$
(3) $\widehat{B}_n$ is bounded in a neighbourhood of S.

(4) We have the asymptotic normalisation

$$ \widehat{B}_n(z) = (\mathbb{I} + \mathcal{O}(z^{-1})) \begin{pmatrix} z^{-2n-D} & & & \\ & z^{2n} & & \\ & & \ddots & \\ & & & z \end{pmatrix}, \quad z \to \infty. $$

Here again we restrict ourselves to computing the second column. We shall simply state the result since the calculation is practically the same as in the case of $\widehat{A}_n$,

$$ \begin{align*} (\widehat{B}_n)_{12} &= 0, & (\widehat{B}_n)_{22} &= P_{2n}, & (\widehat{B}_n)_{i2} &= (\widehat{A}_n)_{i2}, & i \in [3, D+1], \\ (\widehat{B}_n)_{D+2,2} &= \frac{2\pi i \gamma D}{h_{n-1}} P_{2n-2}. \end{align*} $$

### 3 Christoffel–Darboux type formula

In this section, we prove the following theorem.

**Theorem 3.1 (Christoffel–Darboux type formula).** Let the pre-kernel $S_n$ be defined by (1.4). Then

$$ S_n(x, y) = -\frac{1}{4\pi i} e^{-V(x)-V(y)} \frac{(A_n^{-1}(x)A_n(y))_{21}}{x-y}, \quad (3.1) $$

where $A_n$ is the unique solution of Riemann–Hilbert Problem 2.7. Note that we do not specify whether we take the + or – boundary value of $(A_n^{-1}(x)A_n(y))_{21}$ since this will turn out to be an entire function.

**Remark 3.2.** We call this a “Christoffel–Darboux type” formula because an analogous formula appears in the theory of orthogonal polynomials ($\beta = 2$), i.e., if one takes $A_n$ to be the $2 \times 2$ Fokas–Its–Kitaev Riemann–Hilbert problem [27], then (3.1) yields the classical Christoffel–Darboux formula, up to a factor of $\frac{1}{2}$.

**Remark 3.3.** We could, if we like, have expressed the Christoffel–Darboux type formula (3.1) in terms of the solution of the odd Riemann–Hilbert Problem 2.13, since we see from our solution to $B_n$ and $\widehat{B}_n$ that $(A_n^{-1}(x)A_n(y))_{21} = (B_n^{-1}(x)B_n(y))_{21}$.

To prove this result, we need to derive a variety of relations. First note that the polynomial $R_k^{(j)} - R_n^{(j)}$ for $k = 0, \dots, n-1$ is in the kernel of $\beta_j$ for all $j$ and has degree at most $2n + D - 2$. It is therefore in the span of $\Psi_0, \dots, \Psi_{2n-1}$. Thus write

$$ R_k^{(j)} - R_n^{(j)} = \sum_{i=0}^{n-1} \eta_{ik}^{(j)} \Psi_{2i} + \sum_{i=0}^{n-1} \xi_{ik}^{(j)} \Psi_{2i+1}. $$

We then have

$$ \eta_{ik}^{(j)} = \frac{\gamma D}{h_i} \int_{\mathbb{R}} R_k^{(j)} P_{2i+1} e^{-2V} dx, \quad \xi_{ik}^{(j)} = -\frac{\gamma D}{h_i} \int_{\mathbb{R}} R_k^{(j)} P_{2i} e^{-2V} dx. $$

Note that for $i \le k-1$, we have $\eta_{ik}^{(j)} = \xi_{ik}^{(j)} = 0$. We may thus write

$$ R_k^{(j)} = R_n^{(j)} + \sum_{i=k}^{n-1} \eta_{ik}^{(j)} \Psi_{2i} + \sum_{i=k}^{n-1} \xi_{ik}^{(j)} \Psi_{2i+1}. $$
Next consider $X\Psi_{2k+1} - \Psi_{2k+2}$. This has degree at most $2n + D$ and is orthogonal to $\mathcal{P}_{2k-2}^C$. Then $X\Psi_{2k+1} - \Psi_{2k+2} - \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k+1})R_k^{(j)}$ has degree at most $2n + D$ and is in the kernel of $\beta_j$ for all $j = 1, \dots, D-1$. It is thus a linear combination of $\Psi_{2k-2}$, $\Psi_{2k}$ and $\Psi_{2k+1}$. We thus find

$$ 
\begin{aligned}
X\Psi_{2k+1} &= \Psi_{2k+2} + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k+1})R_k^{(j)} + a_k\Psi_{2k+1} + b_k\Psi_{2k} + c_k\Psi_{2k-2} \\
&= \Psi_{2k+2} + a_k\Psi_{2k+1} + b_k\Psi_{2k} + c_k\Psi_{2k-2} + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k+1})R_n^{(j)} \\
&\quad + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k+1})\sum_{i=k}^{n-1} \eta_{ik}^{(j)}\Psi_{2i} + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k+1})\sum_{i=k}^{n-1} \xi_{ik}^{(j)}\Psi_{2i+1}.
\end{aligned}
$$

These coefficients $a_k$, $b_k$ and $c_k$ are unique as can be seen from the formulas

$$ a_k = -\frac{\gamma D}{h_k} \int_{\mathbb{R}} x P_{2k} \Psi_{2k+1} e^{-2V} dx - \sum_{j=1}^{D-1} \beta_j(X \Psi_{2k+1}) \xi_{kk}^{(j)}, $$

$$ b_k = \frac{\gamma D}{h_k} \int_{\mathbb{R}} x \Psi_{2k+1} P_{2k+1} e^{-2V} dx - \sum_{j=1}^{D-1} \beta_j(X \Psi_{2k+1}) \eta_{kk}^{(j)}, \quad c_k = -\frac{h_k}{h_{k-1}}. $$

Similarly, $X\Psi_{2k} - \Psi_{2k+1}$ has degree at most $2n + D - 1$ and is orthogonal to $\mathcal{P}_{2k-1}^C$. Thus $X\Psi_{2k} - \Psi_{2k+1} - \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k})R_k^{(j)}$ is in the kernel of all $\beta_j$ and is orthogonal to $\mathcal{P}_{2k-1}^C$. Thus it must be proportional to $\Psi_{2k}$. This gives

$$ 
\begin{aligned}
X\Psi_{2k} &= \Psi_{2k+1} + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k})R_k^{(j)} + \tilde{a}_k\Psi_{2k} \\
&= \Psi_{2k+1} + \tilde{a}_k\Psi_{2k} + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k})R_n^{(j)} + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k})\sum_{i=k}^{n-1} \eta_{ik}^{(j)}\Psi_{2i} \\
&\quad + \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k})\sum_{i=k}^{n-1} \xi_{ik}^{(j)}\Psi_{2i+1}.
\end{aligned}
\quad (3.2) $$

As before, $\tilde{a}_k$ is unique and may be written

$$ \tilde{a}_k = \frac{\gamma D}{h_k} \int_{\mathbb{R}} x \Psi_{2k} P_{2k+1} e^{-2V} dx - \sum_{j=1}^{D-1} \beta_j(X \Psi_{2k}) \eta_{kk}^{(j)}. $$

Equation (3.2) also implies the intriguing identity

$$ 2 = - \sum_{j=1}^{D-1} \beta_j(X\Psi_{2k})\xi_{kk}^{(j)}. \quad (3.3) $$

These formulas imply “dual” formulas for $P_{2k}$ and $P_{2k+1}$. Let us begin with the even case. $XP_{2k}$ is monic of degree $2k+1$. Thus consider $XP_{2k} - P_{2k+1}$ which has degree at most $2k$. Thus we may expand

$$ XP_{2k} - P_{2k+1} = \sum_{m=0}^{k} \mu_{mk} P_{2m} + \sum_{m=0}^{k-1} \lambda_{mk} P_{2m+1}. $$
We then find

$$ \mu_{mk} = -\frac{\gamma D}{h_m} \int_{\mathbb{R}} x P_{2k} \Psi_{2m+1} e^{-2V} dx = \frac{h_k}{h_m} \sum_{j=1}^{D-1} \beta_j(X \Psi_{2m+1}) \xi_{km}^{(j)} + a_k \delta_{km}, $$

$$ \lambda_{mk} = \frac{\gamma D}{h_m} \int_{\mathbb{R}} x P_{2k} \Psi_{2m} e^{-2V} dx = -\frac{h_k}{h_m} \sum_{j=1}^{D-1} \beta_j(X \Psi_{2m}) \xi_{km}^{(j)}. $$

Finally, consider $XP_{2k+1} - P_{2k+2}$. This has degree at most $2k+1$ and so may be expanded as follows:

$$ XP_{2k+1} - P_{2k+2} = \sum_{m=0}^{k} \tilde{\mu}_{mk} P_{2m} + \sum_{m=0}^{k} \tilde{\lambda}_{mk} P_{2m+1}. $$

Then after a calculation, we find

$$ \tilde{\mu}_{mk} = \delta_{k,m+1} c_k - \frac{h_k}{h_m} \sum_{j=1}^{D-1} \beta_j(X \Psi_{2m+1}) \eta_{km}^{(j)} - b_k \delta_{km}, $$

$$ \tilde{\lambda}_{mk} = \frac{h_k}{h_m} \sum_{j=1}^{D-1} \beta_j(X \Psi_{2m}) \eta_{km}^{(j)} + \tilde{a}_k \delta_{km}. $$

Let us now derive a Christoffel-Darboux type formula for the pre-kernel (1.4). Note that the pre-kernel may be written

$$ S_n(x, y) = -\frac{1}{2} \gamma D e^{-V(x)-V(y)} \sum_{k=0}^{n-1} \frac{P_{2k}(x) \Psi_{2k+1}(y) - P_{2k+1}(x) \Psi_{2k}(y)}{h_k}. $$

Then applying our recursion formulas, we find

$$ (x-y)S_n(x,y) = -\frac{1}{2}\gamma D e^{-V(x)-V(y)} \sum_{k=0}^{n-1} (\textcircled{1}_k - \textcircled{2}_k - \textcircled{3}_k + \textcircled{4}_k), $$

where

$$ \textcircled{1}_k = h_k^{-1} \left( P_{2k+1}(x) + \sum_{m=0}^{k} \mu_{mk} P_{2m}(x) + \sum_{m=0}^{k-1} \lambda_{mk} P_{2m+1}(x) \right) \Psi_{2k+1}(y), $$

$$ \textcircled{2}_k = h_k^{-1} P_{2k}(x) \left( \sum_{i=k}^{n-1} \sum_{j=1}^{D-1} \beta_j(X \Psi_{2k+1}) [\eta_{ik}^{(j)} \Psi_{2i}(y) + \xi_{ik}^{(j)} \Psi_{2i+1}(y)] \right. \\ \left. + \Psi_{2k+2}(y) + a_k \Psi_{2k+1}(y) + b_k \Psi_{2k}(y) + c_k \Psi_{2k-2}(y) + \sum_{j=1}^{D-1} \beta_j(X \Psi_{2k+1}) R_n^{(j)}(y) \right), $$

$$ \textcircled{3}_k = h_k^{-1} \left( P_{2k+2}(x) + \sum_{m=0}^{k} \tilde{\mu}_{mk} P_{2m}(x) + \sum_{m=0}^{k} \tilde{\lambda}_{mk} P_{2m+1}(x) \right) \Psi_{2k}(y), $$

$$ \textcircled{4}_k = h_k^{-1} P_{2k+1}(x) \left( \Psi_{2k+1}(y) + \tilde{a}_k \Psi_{2k}(y) + \sum_{j=1}^{D-1} \beta_j(X \Psi_{2k}) R_n^{(j)}(y) \right. \\ \left. + \sum_{i=k}^{n-1} \sum_{j=1}^{D-1} \beta_j(X \Psi_{2k}) [\eta_{ik}^{(j)} \Psi_{2i}(y) + \xi_{ik}^{(j)} \Psi_{2i+1}(y)] \right). $$
Let us then write

$$ \sum_{k=0}^{n-1} \binom{2}{k} = h_{n-1}^{-1} P_{2n-2}(x) \Psi_{2n}(y) + h_{n-1}^{-1} P_{2n}(x) \Psi_{2n-2}(y) \\ + \sum_{j=1}^{D-1} \sum_{k=0}^{n-1} \beta_j(X \Psi_{2k+1}) h_k^{-1} P_{2k}(x) R_n^{(j)}(y) + \sum_{k=0}^{n-1} \binom{5}{k}, $$

where

$$ \binom{5}{k} = h_{k-1}^{-1} P_{2k-2}(x) \Psi_{2k}(y) + h_k^{-1} P_{2k}(x) (a_k \Psi_{2k+1}(y) + b_k \Psi_{2k}(y)) - h_k^{-1} P_{2k+2}(x) \Psi_{2k}(y) \\ + \sum_{m=0}^{k} h_m^{-1} P_{2m}(x) \sum_{j=1}^{D-1} \beta_j(X \Psi_{2m+1}) [\eta_{km}^{(j)} \Psi_{2k}(y) + \xi_{km}^{(j)} \Psi_{2k+1}(y)], $$

and similarly we may write

$$ \sum_{k=0}^{n-1} \binom{4}{k} = \sum_{j=1}^{D-1} \sum_{k=0}^{n-1} h_k^{-1} P_{2k+1}(x) \beta_j(X \Psi_{2k}) R_n^{(j)}(y) + \sum_{k=0}^{n-1} \binom{6}{k}, $$

where

$$ \binom{6}{k} = h_k^{-1} P_{2k+1}(x) \Psi_{2k+1}(y) + \tilde{a}_k h_k^{-1} P_{2k+1}(x) \Psi_{2k}(y) \\ + \sum_{m=0}^{k} h_m^{-1} P_{2m+1}(x) \sum_{j=1}^{D-1} \beta_j(X \Psi_{2m}) [\eta_{km}^{(j)} \Psi_{2k}(y) + \xi_{km}^{(j)} \Psi_{2k+1}(y)]. $$

Thus

$$ (x - y)S_n(x, y) = \frac{\gamma D}{2} e^{-V(x)-V(y)} (h_{n-1}^{-1} P_{2n-2}(x) \Psi_{2n}(y) + h_{n-1}^{-1} P_{2n}(x) \Psi_{2n-2}(y) \\ + \sum_{j=1}^{D-1} \sum_{k=0}^{n-1} h_k^{-1} [\beta_j(X \Psi_{2k+1}) P_{2k}(x) - \beta_j(X \Psi_{2k}) P_{2k+1}(x)] R_n^{(j)}(y) \\ + \frac{\gamma D}{2} e^{-V(x)-V(y)} \sum_{k=0}^{n-1} (-\binom{1}{k} + \binom{5}{k} + \binom{3}{k} - \binom{6}{k}). $$

The reader may verify that a miraculous cancellation occurs so that

$$ -\binom{1}{k} + \binom{5}{k} + \binom{3}{k} - \binom{6}{k} = 0. $$

Finally, using our solutions of $A_n$ and $\tilde{A}_n$,

$$ (x - y)S_n(x, y) = \frac{\gamma D}{2} e^{-V(x)-V(y)} (h_{n-1}^{-1} P_{2n-2}(x) \Psi_{2n}(y) + h_{n-1}^{-1} P_{2n}(x) \Psi_{2n-2}(y) \\ + \sum_{j=1}^{D-1} \sum_{k=0}^{n-1} h_k^{-1} [\beta_j(X \Psi_{2k+1}) P_{2k}(x) - \beta_j(X \Psi_{2k}) P_{2k+1}(x)] R_n^{(j)}(y) \\ = - \frac{e^{-V(x)-V(y)}}{4\pi i} (A_n^{-1}(x) A_n(y))_{21}. $$
**Remark 3.4.** We refer to this cancellation as “miraculous”, however one could have predicted that some sort of cancellation must occur simply from the form of the pre-kernel (1.4). Introduce the integral operator

$$S_n: L^2(\mathbb{R}) \to L^2(\mathbb{R}), \quad (S_n\psi)(x) = \int_{\mathbb{R}} S_n(x, y)\psi(y)dy, \quad \psi \in L^2(\mathbb{R}).$$

It is then easily seen that if $\{H_k\}_{k \in \mathbb{N}}$ is the sequence of orthogonal polynomials with respect to the weight $w = e^{-2V}$, then $S_n(H_k e^{-V}) = \frac{1}{2} H_k e^{-V}$ for $k \le 2n - 1$ and $S_n(H_k e^{-V}) = 0$ for $k \ge 2n + D - 1$. It thus follows that

$$S_n = \frac{1}{2}K_n + \sum_{j=2n}^{2n+D-2} \frac{1}{\nu_j} S_n(H_j e^{-V}) \otimes H_j e^{-V},$$

where $K_n = \sum_{k=0}^{n-1} \frac{1}{\nu_k} H_k \otimes H_k$ for $\nu_k = \int_{\mathbb{R}} H_k^2 e^{-2V} dx$. By the classical Christoffel–Darboux formula for orthogonal polynomials and the three term recurrence relation, it follows that $(x - y)S_n(x, y)$ has rank at most $D + 1$. Indeed, this is one of the ideas behind Widom’s method.

# 4 A dynamical system

This paper is part of a broader project to develop structures for $\beta = 1, 4$ random matrix ensembles parallel to those of $\beta = 2$ ensembles. A central part of the $\beta = 2$ theory is the relation between orthogonal polynomials and the celebrated Toda lattice (see [14, Chapter 2]). One way to derive this relation is, starting with the Fokas–Its–Kitaev Riemann–Hilbert problem for orthogonal polynomials [27], to pass to a problem with constant jump matrices and then by standard techniques of varying parameters one derives a Lax pair whose compatibility condition gives the Toda lattice in Flaschka variables. In this section we repeat this method to see what analogous dynamical system arises for our Riemann–Hilbert Problem 2.7. The calculation is complicated and so we confine ourselves simply to deriving the compatibility condition, and leave the question of implications for the theory of integrable systems to future research.

We begin, analogously to the Toda lattice, by introducing a $t \in \mathbb{R}$ dependence into the potential as

$$V(z, t) := V_0(z) + tz. \quad (4.1)$$

Let $\hat{A}_n(\cdot, t)$ be the solution of Riemann–Hilbert Problem 2.24 with potential $V(\cdot, t)$ and let $\hat{A}_n(\cdot, t)$ have the following expansion at $z = \infty$:

$$\hat{A}_n(z,t) = (\mathbb{I} + z^{-1}\hat{A}_n^{(1)}(t) + z^{-2}\hat{A}_n^{(2)}(t) + \mathcal{O}(z^{-3})) \begin{pmatrix} z^{-2n-D+1} & & & \\ & z^{2n} & & \\ & & z & \\ & & \ddots & \\ & & & z \end{pmatrix}, \\ z \to \infty.$$

Introduce the variables $a_{ij}^{(1)}(n, t) \stackrel{\text{def}}{=} \hat{A}_n^{(1)}(t)_{ij}$ and $a_{ij}^{(2)}(n, t) \stackrel{\text{def}}{=} \hat{A}_n^{(2)}(t)_{ij}$ for $i, j \in [1, D+1]$. In what follows, we shall often not write the dependence on $t$ since this is always understood. Furthermore, introduce the variable

$$b_j(n, t) \stackrel{\text{def}}{=} \sum_{k=3}^{D+1} a_{2k}^{(1)}(n, t) a_{kj}^{(1)}(n, t), \quad j \in [3, D+1]. \quad (4.2)$$
**Theorem 4.1.** The matrix elements $a_{22}^{(1)}$, $a_{2i}^{(1)}$, $a_{i2}^{(1)}$, $a_{22}^{(2)}$, $a_{12}^{(2)}$, $a_{i2}^{(2)}$ and $b_i$, for $i = 3, \dots, D+1$, satisfy the following first-order dynamical system ($n \ge 0$, $t \in \mathbb{R}$):

$$ \frac{da_{22}^{(1)}(n)}{dt} = \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k2}^{(1)}(n), \quad \frac{da_{2i}^{(1)}(n)}{dt} = b_i(n), \quad i \in [3, D+1], $$

$$ \frac{da_{i2}^{(1)}(n)}{dt} = a_{i2}^{(2)}(n) - a_{i2}^{(1)}(n)a_{22}^{(1)}(n), \quad i \in [3, D+1], $$

$$ \frac{da_{22}^{(2)}(n)}{dt} = \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k2}^{(2)}(n), \quad \frac{d}{dt} \ln a_{12}^{(2)}(n+1) = 2a_{22}^{(1)}(n) - \sum_{j=3}^{D+1} a_{2j}^{(1)}(n) a_{j2}^{(2)}(n+1), $$

$$ \begin{aligned} \frac{d}{dt} a_{i2}^{(2)}(n+1) = & a_{i2}^{(2)}(n+1)a_{22}^{(1)}(n) + a_{i2}^{(1)}(n+1)a_{22}^{(2)}(n) + a_{i2}^{(1)}(n) - a_{i2}^{(1)}(n+1)a_{22}^{(2)}(n+1) \\ & + a_{i2}^{(1)}(n+2) \frac{a_{12}^{(2)}(n+1)}{a_{12}^{(2)}(n+2)} - a_{i2}^{(1)}(n+1)a_{22}^{(1)}(n)^2 \\ & - a_{i2}^{(1)}(n+1) \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k2}^{(1)}(n), \quad i \in [3, D+1], \end{aligned} $$

and for $n \ge 1$ and $t \in \mathbb{R}$,

$$ \begin{aligned} \frac{db_i(n)}{dt} = & \frac{a_{12}^{(2)}(n)}{a_{12}^{(2)}(n+1)} a_{2i}^{(1)}(n-1) + a_{22}^{(1)}(n+1) a_{22}^{(1)}(n) a_{2i}^{(1)}(n) + a_{2i}^{(1)}(n+1) \\ & - a_{22}^{(2)}(n+1) a_{2i}^{(1)}(n) + a_{22}^{(2)}(n) a_{2i}^{(1)}(n) - a_{2i}^{(1)}(n) (a_{22}^{(1)}(n))^2 - a_{22}^{(1)}(n) b_i(n) \\ & - 2a_{2i}^{(1)}(n) \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k2}^{(1)}(n) + a_{22}^{(1)}(n+1) b_i(n), \quad i \in [3, D+1]. \end{aligned} $$

Note that $a_{12}^{(2)}(n+1) \ne 0$ so the above is well defined. The above dynamics are furthermore constrained by the relation

$$ \sum_{j=3}^{D+1} a_{2j}^{(1)}(n) a_{j2}^{(1)}(n+1) = 2. \quad (4.3) $$

Our “initial data” with respect to $n$ must satisfy $a_{22}^{(1)}(0, t) = a_{i2}^{(1)}(0, t) = a_{22}^{(2)}(0, t) = a_{12}^{(2)}(0, t) = a_{i2}^{(2)}(0, t) = b_i(0, t) = 0$ for $i \in [3, D+1]$.

The identity (4.3) may be seen as analogous to (3.3). What follows in the rest of this section is a proof of Theorem 4.1.

**Remark 4.2.** Note that repeated differentiation of (4.3) yields an infinite number of identities, so that the $k$-th derivative of (4.3) yields an identity involving matrix elements at $k$ different values of $n$.

**Remark 4.3.** Some of the matrix elements not included in Theorem 4.1 may be related to those included by the following identities:

$$ a_{i1}^{(1)}(n) = a_{21}^{(2)}(n) a_{i2}^{(1)}(n+1), \quad i \in [3, D+1], \quad (4.4) $$

$$ a_{1i}^{(1)}(n+1) = a_{12}^{(2)}(n+1) a_{2i}^{(1)}(n), \quad i \in [3, D+1], \quad (4.5) $$

$$ \begin{aligned} 1 = a_{21}^{(2)}(n) a_{12}^{(2)}(n+1), \quad & 0 = a_{12}^{(1)}(n) = a_{21}^{(1)}(n) = a_{ij}^{(2)}(n), \\ i \in [1, D+1], \quad & j \in [3, D+1]. \end{aligned} $$
These identities can be seen directly from the solution, though (4.4) and (4.5) will be deduced in the course of proving Theorem 4.1. The matrix elements $a_{i1}^{(2)}$ for $i \neq 2$ are left undetermined by these relations. Furthermore, $a_{ij}^{(1)}(n, t)$, for $i, j \ge 3$ do not enter into the dynamical system on their own but only in the combination of the variable $b_i(n, t)$ (4.2).

**Example 4.4.** In the case $D = 2$, it is possible to give an exact solution. Without loss of generality (by shifting and rescaling), one can take $V_0(z) = \frac{1}{2}z^2$. With this initial condition one finds the following:

$$ 
\begin{aligned}
\hat{A}_n^{(1)}(t) &= \begin{pmatrix}
-t(2n+1) & 0 & \frac{2\pi i e^{-\frac{1}{2}t^2}}{\sqrt{2}(n-\frac{1}{2})!} \\
0 & 2nt & -\frac{e^{\frac{1}{2}t^2} n!}{\sqrt{2}} \\
\frac{e^{\frac{1}{2}t^2} 2\sqrt{2}(n+\frac{1}{2})!}{2\pi i} & -\frac{2\sqrt{2}e^{-\frac{1}{2}t^2}}{(n-1)!} & t
\end{pmatrix}, \\
\hat{A}_n^{(2)}(t) &= \begin{pmatrix}
* & -i\sqrt{\pi}\frac{2^{2n}e^{-t^2}}{(2n-1)!} & 0 \\
-\frac{(2n+1)!\sqrt{\pi}e^{t^2}}{2^{2n+1}2\pi i} & n(2n-1)(t^2-\frac{1}{2})+n & 0 \\
* & -\frac{2\sqrt{2}(2n-1)te^{-\frac{1}{2}t^2}}{(n-1)!} & 0
\end{pmatrix}.
\end{aligned}
$$

We find also $b_3(n, t) = -\frac{te^{\frac{1}{2}t^2}n!}{\sqrt{2}}$. The reader may verify that this is an exact solution of the dynamical system in Theorem 4.1 for the case $D = 2$.

Let us introduce the $(D+1) \times (D+1)$ matrices

$$ E_1 = \operatorname{diag}(1, \underbrace{0, \dots, 0}_{D \text{ times}}), \quad E_2 = \operatorname{diag}(0, 1, \underbrace{0, \dots, 0}_{D-1 \text{ times}}), \quad \sigma = E_1 - E_2. $$

Let us then define

$$ \tilde{A}_n(z, t) \stackrel{\text{def}}{=} \hat{A}_n(z, t) e^{V(z, t)\sigma}. $$

It is then easy to see that $\tilde{A}_n$ satisfies the following Riemann–Hilbert problem.

**Riemann–Hilbert Problem 4.5.** Let $\Gamma$ be as in (2.4). Then $\tilde{A}_n$ has the following properties:

(1) $z \mapsto \tilde{A}_n(z, t)$ is an analytic function on $\mathbb{C} \setminus \Gamma$.

(2) $z \mapsto \tilde{A}_n(z, t)$ has continuous left (+) and right (-) non-tangential boundary valuesx $\mapsto \tilde{A}_n^{\pm}(x, t)$ for $x \in \Gamma \setminus \{0\}$ which are related by the jump condition

$$ \tilde{A}_n^+(x, t) = \tilde{A}_n^-(x, t) \begin{pmatrix} 1 & & & \\ & -\chi_{\mathbb{R}} & 1 & \\ & -\chi_{\Gamma_1} & 1 & \\ & \vdots & \ddots & \\ & -\chi_{\Gamma_{D-1}} & & 1 \end{pmatrix}, \quad x \in \Gamma \setminus \{0\}. $$

(3) For each $t \in \mathbb{R}$, the function $z \mapsto \tilde{A}_n(z, t)$ is bounded in a neighbourhood of $0 \in \mathbb{C}$.

(4) For each $t \in \mathbb{R}$, $z \mapsto \tilde{A}_n(z, t)$ satisfies the asymptotic normalisation

$$ \tilde{A}_n(z, t) = (\mathbb{I} + \mathcal{O}(z^{-1})) \begin{pmatrix} z^{-2n-D+1} & & & \\ & z^{2n} & & \\ & & z & \\ & & \ddots & \\ & & & z \end{pmatrix} e^{V(z,t)\sigma}, \quad z \to \infty. $$
Because the jump matrix is independent of both $n$ and $t$, standard Riemann-Hilbert theory now enables us to derive a dynamical system. In particular $\tilde{A}_{n+1}(z,t)$ and $\frac{\partial \tilde{A}_n(z,t)}{\partial t}$ both satisfy the first to third properties, differing only on the fourth property. Because of the continuity of the boundary values and boundedness at 0, we find that

$$z \mapsto \tilde{A}_{n+1}(z,t) \tilde{A}_n(z,t)^{-1}, \quad z \mapsto \frac{\partial \tilde{A}_n(z,t)}{\partial t} \tilde{A}_n(z,t)^{-1}$$

extend to entire functions on $\mathbb{C}$. Let us write the asymptotic expansion $(z \to \infty)$

$$\tilde{A}_n(z,t) = (\mathbb{I} + z^{-1}\hat{A}_n^{(1)}(t) + z^{-2}\hat{A}_n^{(2)}(t) + \mathcal{O}(z^{-3})) \begin{pmatrix} z^{-2n-D+1} & & & \\ & z^{2n} & & \\ & & z & \\ & & & \vdots \\ & & & z \end{pmatrix} e^{V(z,t)\sigma}.$$

A short calculation yields the *Lax pair*,

$$\tilde{A}_{n+1}(z,t) = (z^2 E_2 + z \rho_n^{(1)}(t) + \rho_n^{(2)}(t)) \tilde{A}_n(z,t), \quad (4.6)$$

$$\frac{\partial \tilde{A}_n(z,t)}{\partial t} = (z\sigma + \kappa_n(t)) \tilde{A}_n(z,t), \quad (4.7)$$

where

$$\rho_n^{(1)}(t) = \hat{A}_{n+1}^{(1)}(t)E_2 - E_2\hat{A}_n^{(1)}(t),$$

$$\rho_n^{(2)}(t) = \mathbb{I} - E_1 - E_2 + \hat{A}_{n+1}^{(2)}(t)E_2 + E_2(\hat{A}_n^{(1)}(t)^2 - \hat{A}_n^{(2)}(t)) - \hat{A}_{n+1}^{(1)}(t)E_2\hat{A}_n^{(1)}(t),$$

$$\kappa_n(t) = [\hat{A}_n^{(1)}(t), \sigma].$$

**Remark 4.6.** The reader will observe that the jump matrices of $\tilde{A}_n$ are not only independent of $n$ and $t$ but also independent of $z$. Thus one could also derive $(n, z)$ and $(t, z)$ Lax pairs. Indeed, one could consider more general deformations than (4.1).

**Remark 4.7** (analogue of three term recurrence). We note that (4.6) implies a similar equation for $\hat{A}_n$, i.e.,

$$\hat{A}_{n+1}(z) = (z^2 E_2 + z \rho_n^{(1)} + \rho_n^{(2)}) \hat{A}_n(z).$$

Considering the second column of both sides yields a recursion relation, quadratic in $z$, which may be regarded as a $\beta = 4$ analogue of the three term recurrence for orthogonal polynomials (recall equations (2.7), (2.8), (2.9)).

Before proceeding, we must indicate a point of rigour. To compute this Lax pair we have exchanged a $z \to \infty$ limit with a derivative with respect to $t$. How can this be justified? The proof is elementary but tedious to write out and so we only sketch how this is done. The solution of our even problem $\hat{A}_n$ is expressible as rational combinations of integrals of the form

$$\int_{\Gamma_j} x^k e^{-V(x,t)} dx, \quad \int_{\mathbb{R}} x^k e^{-2V(x,t)} dx, \quad \int_{\Gamma_j} \frac{x^k}{x-z} e^{-V(x,t)} dx, \quad \int_{\mathbb{R}} \frac{x^k}{x-z} e^{-2V(x,t)} dx.$$

The proof reduces to showing that such integrals are $\mathscr{C}^\infty$ functions of $t$ and that we can differentiate under the integral sign. Both can be justified by a dominated convergence argument. This implies that $\hat{A}_n^{(1)}$ and $\hat{A}_n^{(2)}$ are $\mathscr{C}^\infty$ functions of $t$.
Moving forward, let us write

$$ \rho(z, t, n) = z^2 E_2 + z \rho_n^{(1)}(t) + \rho_n^{(2)}(t), \quad \lambda(z, t, n) = z\sigma + \kappa_n(t). $$

Then the compatibility condition (or zero curvature condition) of the above Lax pair (4.6) and (4.7) is

$$ \frac{\partial \rho(z, t, n)}{\partial t} + \rho(z, t, n)\lambda(z, t, n) - \lambda(z, t, n+1)\rho(z, t, n) = 0, \quad \forall z \in \mathbb{C}. \quad (4.8) $$

The left-hand side is at most a cubic polynomial in $z$, all of whose coefficients must vanish. In fact the coefficient of $z^3$ is $[E_2, \sigma] = 0$ which yields no information. The coefficient of $z^2$ yields the equation

$$ E_2\kappa_n(t) - \kappa_{n+1}(t)E_2 + [\rho_n^{(1)}(t), \sigma] = 0. $$

However, when we substitute in our formulas for $\rho_n^{(1)}$ and $\kappa_n$ we find that the left-hand side is identically zero. So the quadratic coefficient also yields no information. Thus the left-hand side of (4.8) is in fact a linear polynomial in $z$.

The coefficient of $z$ in (4.8) yields

$$ \frac{d\rho_n^{(1)}(t)}{dt} + \rho_n^{(1)}(t)\kappa_n(t) - \kappa_{n+1}(t)\rho_n^{(1)}(t) + [\rho_n^{(2)}(t), \sigma] = 0, \quad (4.9) $$

while the constant term gives

$$ \frac{d\rho_n^{(2)}(t)}{dt} + \rho_n^{(2)}(t)\kappa_n(t) - \kappa_{n+1}(t)\rho_n^{(2)}(t) = 0. \quad (4.10) $$

Equation (4.9), after substituting our expressions for $\rho_n^{(1)}(t)$ and $\rho_n^{(2)}(t)$, yields

$$ \begin{aligned} & -\hat{A}_{n+1}^{(2)} E_2 + \hat{A}_{n+1}^{(1)} E_2 \hat{A}_{n+1}^{(1)} E_2 - E_1 \hat{A}_{n+1}^{(2)} E_2 - E_2 \hat{A}_n^{(2)} + E_2 \hat{A}_n^{(2)} E_2 - E_2 \hat{A}_n^{(2)} E_1 + E_2 \hat{A}_{n+1}^{(2)} E_2 \\ & - E_2 \frac{d\hat{A}_n^{(1)}}{dt} - E_2 (\hat{A}_{n+1}^{(1)})^2 E_2 - E_2 \hat{A}_n^{(1)} E_2 \hat{A}_n^{(1)} + E_2 (\hat{A}_n^{(1)})^2 + \frac{d\hat{A}_{n+1}^{(1)}}{dt} E_2 + E_1 (\hat{A}_{n+1}^{(1)})^2 E_2 \\ & = 0. \end{aligned} \quad (4.11) $$

The (1,1) matrix element of (4.11) is 0. The (1,2) matrix element yields the identity

$$ 2a_{12}^{(2)}(n) = \sum_{k=3}^{D+1} a_{1k}^{(1)}(n) a_{k2}^{(1)}(n). \quad (4.12) $$

The (2,1) matrix element yields

$$ 2a_{21}^{(2)}(n) = \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k1}^{(1)}(n). \quad (4.13) $$

The (2,2) element gives that the quantity

$$ (n, t) \mapsto \frac{d}{dt} a_{22}^{(1)}(n, t) - \sum_{k=3}^{D+1} a_{2k}^{(1)}(n, t) a_{k2}^{(1)}(n, t) $$

is independent of $n$.
**Lemma 4.8.**

$$ \frac{d}{dt} a_{22}^{(1)}(n) = \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k2}^{(1)}(n). $$

**Proof.** To show this, we need only show that

$$ \frac{d}{dt} a_{22}^{(1)}(0, t) - \sum_{k=3}^{D+1} a_{2k}^{(1)}(0, t) a_{k2}^{(1)}(0, t) = 0. $$

However, if $n = 0$, then $(\hat{A}_n)_{22} = P_{2n} = P_0 = 1$, and hence $a_{22}^{(1)}(0, t) = 0$. Similarly, $(\hat{A}_n)_{k2} = 0$ for $k \ge 3$, and so $a_{k2}^{(1)}(0, t) = 0$. ■

The $(1, i)$ and $(i, 1)$ matrix elements for $i \ge 3$ give 0. The $(2, i)$ matrix element for $i \ge 3$ gives

$$ \frac{da_{2i}^{(1)}(n)}{dt} = \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{ki}^{(1)}(n), \quad i \in [3, D+1]. $$

Finally, the $(i, 2)$ matrix element for $i \ge 3$ gives the equation

$$ \frac{da_{i2}^{(1)}(n, t)}{dt} = a_{i2}^{(2)}(n, t) - a_{i2}^{(1)}(n, t) a_{22}^{(1)}(n, t), \quad i \in [3, D+1]. $$

Next let us move onto the constant term (4.10), which gives us the complicated expression

$$
\begin{align}
& E_2 \hat{A}_n^{(2)} E_1 \hat{A}_n^{(1)} - E_2 \hat{A}_n^{(2)} \hat{A}_n^{(1)} E_1 + E_2 \hat{A}_n^{(2)} \hat{A}_n^{(1)} E_2 - E_2 \hat{A}_{n+1}^{(1)} \hat{A}_{n+1}^{(2)} E_2 + E_2 \hat{A}_{n+1}^{(1)} E_2 \nonumber \\
& - E_2 (\hat{A}_n^{(1)})^3 E_2 - \hat{A}_{n+1}^{(1)} E_2 \hat{A}_n^{(1)} E_2 \hat{A}_n^{(1)} - \hat{A}_{n+1}^{(2)} E_2 \hat{A}_n^{(1)} E_2 + E_2 \hat{A}_{n+1}^{(1)} E_2 \hat{A}_n^{(2)} \nonumber \\
& + E_1 \hat{A}_{n+1}^{(1)} \hat{A}_{n+1}^{(2)} E_2 - \hat{A}_{n+1}^{(1)} E_2 (\hat{A}_n^{(1)})^2 E_1 + E_2 \hat{A}_n^{(1)} \frac{d\hat{A}_n^{(1)}}{dt} - E_1 (\hat{A}_{n+1}^{(1)})^2 E_2 \hat{A}_n^{(1)} \nonumber \\
& - E_2 \frac{d\hat{A}_n^{(2)}}{dt} - E_2 \hat{A}_{n+1}^{(1)} E_2 (\hat{A}_n^{(1)})^2 - E_1 \hat{A}_n^{(1)} E_1 + E_2 \frac{d\hat{A}_n^{(1)}}{dt} \hat{A}_n^{(1)} - \hat{A}_{n+1}^{(1)} E_2 \frac{d\hat{A}_n^{(1)}}{dt} \nonumber \\
& + \frac{d\hat{A}_{n+1}^{(2)}}{dt} E_2 - E_2 \hat{A}_{n+1}^{(1)} - E_2 (\hat{A}_n^{(1)})^2 E_1 \hat{A}_n^{(1)} + E_2 (\hat{A}_{n+1}^{(1)})^2 E_2 \hat{A}_n^{(1)} + \hat{A}_{n+1}^{(2)} E_2 \hat{A}_n^{(1)} \nonumber \\
& + \hat{A}_n^{(1)} E_1 - \hat{A}_{n+1}^{(1)} E_2 \hat{A}_n^{(2)} - \hat{A}_n^{(1)} E_2 + E_2 \hat{A}_n^{(1)} E_2 + E_2 (\hat{A}_n^{(1)})^2 E_2 \hat{A}_n^{(1)} \nonumber \\
& - \hat{A}_{n+1}^{(1)} E_2 \hat{A}_{n+1}^{(1)} E_2 \hat{A}_n^{(1)} - \frac{d\hat{A}_{n+1}^{(1)}}{dt} E_2 \hat{A}_n^{(1)} + \hat{A}_{n+1}^{(1)} E_2 \hat{A}_{n+1}^{(2)} E_2 + E_2 (\hat{A}_n^{(1)})^3 E_1 + E_1 \hat{A}_{n+1}^{(1)} \nonumber \\
& + \hat{A}_{n+1}^{(1)} E_2 (\hat{A}_n^{(1)})^2 E_2 + \hat{A}_{n+1}^{(1)} E_2 (\hat{A}_n^{(1)})^2 - E_2 \hat{A}_n^{(2)} E_2 \hat{A}_n^{(1)} - E_1 \hat{A}_{n+1}^{(1)} E_1 \nonumber \\
& - \hat{A}_{n+1}^{(1)} E_1 \hat{A}_{n+1}^{(2)} E_2 = 0. \tag{4.14}
\end{align}
$$

The $(1, 1)$ matrix element of (4.14) is 0. The $(2, 2)$ matrix element yields, after some simplification, the relation

$$ \frac{da_{22}^{(2)}(n+1)}{dt} - \sum_{k=3}^{D+1} a_{2k}^{(1)}(n+1) a_{k2}^{(2)}(n+1) = \frac{da_{22}^{(2)}(n)}{dt} - \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k2}^{(2)}(n). $$

**Lemma 4.9.**

$$ \frac{da_{22}^{(2)}(n)}{dt} = \sum_{k=3}^{D+1} a_{2k}^{(1)}(n) a_{k2}^{(2)}(n). $$
**Proof.** Follows by the same observations as in Lemma 4.8.

The (2,1) matrix element of (4.14) gives the following equation:

$$2(a_{22}^{(1)}(n) - a_{22}^{(1)}(n+1))a_{21}^{(2)}(n) + \frac{d}{dt}a_{21}^{(2)}(n) + \sum_{j=3}^{D+1} b_j(n)a_{j1}^{(1)}(n) = 0. \quad (4.15)$$

The (1,2) matrix element gives

$$\frac{d}{dt}a_{12}^{(2)}(n+1) + \sum_{k=3}^{D+1} a_{1k}^{(1)}(n+1)a_{k2}^{(2)}(n+1) - 2a_{22}^{(1)}(n)a_{12}^{(2)}(n+1) = 0. \quad (4.16)$$

The (i, 1) matrix element for $i \ge 3$ gives

$$a_{i1}^{(1)}(n) = a_{21}^{(2)}(n)a_{i2}^{(1)}(n+1), \quad i \in [3, D+1]. \quad (4.17)$$

The (1, i) matrix element gives

$$a_{1i}^{(1)}(n+1) = a_{12}^{(2)}(n+1)a_{2i}^{(1)}(n), \quad i \in [3, D+1]. \quad (4.18)$$

The (i, 2) matrix element yields after some simplification

$$\begin{aligned} \frac{d}{dt}a_{i2}^{(2)}(n+1) = & a_{i2}^{(2)}(n+1)a_{22}^{(1)}(n) + a_{i2}^{(1)}(n+1)a_{22}^{(2)}(n) + a_{i2}^{(1)}(n) - a_{i2}^{(1)}(n+1)a_{22}^{(2)}(n+1) \\ & - a_{i2}^{(1)}(n+1) \sum_{k=3}^{D+1} a_{2k}^{(1)}(n)a_{k2}^{(1)}(n) + a_{21}^{(2)}(n+1)a_{i2}^{(1)}(n+2)a_{12}^{(2)}(n+1) \\ & - a_{i2}^{(1)}(n+1)a_{22}^{(1)}(n)^2, \quad i \in [3, D+1]. \end{aligned}$$

The (2, i) matrix element gives

$$\begin{aligned} & -a_{21}^{(2)}(n)a_{12}^{(2)}(n)a_{2i}^{(1)}(n-1) - a_{22}^{(1)}(n+1)a_{22}^{(1)}(n)a_{2i}^{(1)}(n) - a_{2i}^{(1)}(n+1) + a_{22}^{(2)}(n+1)a_{2i}^{(1)}(n) \\ & + a_{2i}^{(1)}(n)(a_{22}^{(1)}(n))^2 + 2a_{2i}^{(1)}(n) \sum_{k=3}^{D+1} a_{2k}^{(1)}(n)a_{k2}^{(1)}(n) + \frac{db_i(n)}{dt} \\ & + (a_{22}^{(1)}(n) - a_{22}^{(1)}(n+1))b_i(n) - a_{22}^{(2)}(n)a_{2i}^{(1)}(n) = 0, \quad i \in [3, D+1]. \end{aligned}$$

Let us now combine our formulas. Substituting (4.17) into (4.13), we find

$$\sum_{k=3}^{D+1} a_{2k}^{(1)}(n)a_{k2}^{(1)}(n+1) = 2.$$

We would have found the same if we had substituted (4.18) into (4.12). Note that if we combine (4.15) and (4.16), we find

$$\begin{aligned} & \frac{d}{dt} \ln a_{21}^{(2)}(n) + \frac{d}{dt} \ln a_{12}^{(2)}(n+1) \\ & = 2a_{22}^{(1)}(n+1) - \sum_{k=3}^{D+1} b_k(n)a_{k2}^{(1)}(n+1) - \sum_{k=3}^{D+1} a_{2k}^{(1)}(n)a_{k2}^{(2)}(n+1) \\ & = -\frac{d}{dt} \underbrace{\sum_{k=3}^{D+1} a_{2k}^{(1)}(n)a_{k2}^{(1)}(n+1)}_{=2} = 0, \end{aligned}$$

which we already knew since $a_{21}^{(2)}(n)a_{12}^{(2)}(n+1) = 1$. Finally, the (i, j) matrix element of (4.14) for $i, j \ge 3$, after substituting the identities we have derived so far, is identically zero so yields no further information. This completes the proof of Theorem 4.1.
Acknowledgements

The author would like to express his gratitude to Thomas Bothner and Mattia Cafasso, whose comments and suggestions have been invaluable for this paper. The author would also like to express his gratitude to the anonymous referees whose comments have improved the presentation of the paper and have drawn attention to connections to recent work. This work was supported by the UK Engineering and Physical Sciences Research Council through grant EP/T013893/2 and by the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No. 884584).

References

[1] Adler M., Forrester P.J., Nagao T., van Moerbeke P., Classical skew orthogonal polynomials and random matrices, *J. Stat. Phys.* **99** (2000), 141–170, arXiv:solv-int/9907001.

[2] Adler M., Horozov E., van Moerbeke P., The Pfaff lattice and skew-orthogonal polynomials, *Int. Math. Res. Not.* **1999** (1999), 569–588, arXiv:solv-int/9903005.

[3] Adler M., van Moerbeke P., The Pfaff lattice, matrix integrals and a map from Toda to Pfaff, arXiv:solv-int/9912008.

[4] Akemann G., Ebke M., Parra I., Skew-orthogonal polynomials in the complex plane and their Bergman-like kernels, *Comm. Math. Phys.* **389** (2022), 621–659, arXiv:2103.12114.

[5] Berezin S., Kuijlaars A.B.J., Parra I., Planar orthogonal polynomials as type I multiple orthogonal polynomials, *SIGMA* **19** (2023), 020, 18 pages, arXiv:2212.06526.

[6] Bleher P., Its A., Semiclassical asymptotics of orthogonal polynomials, Riemann–Hilbert problem, and universality in the matrix model, *Ann. of Math.* **150** (1999), 185–266, arXiv:math-ph/9907025.

[7] Bleher P., Its A., Double scaling limit in the random matrix model: the Riemann–Hilbert approach, *Comm. Pure Appl. Math.* **56** (2003), 433–516, arXiv:math-ph/0201003.

[8] Byun S.-S., Ebke M., Universal scaling limits of the symplectic elliptic Ginibre ensemble, *Random Matrices Theory Appl.* **12** (2023), 2250047, 33 pages, arXiv:2108.05541.

[9] Byun S.-S., Ebke M., Seo S.-M., Wronskian structures of planar symplectic ensembles, *Nonlinearity* **36** (2023), 809–844, arXiv:2110.12196.

[10] Byun S.-S., Noda K., Scaling limits of complex and symplectic non-Hermitian Wishart ensembles, arXiv:2402.18257.

[11] Claeys T., Kuijlaars A.B.J., Universality of the double scaling limit in random matrix models, *Comm. Pure Appl. Math.* **59** (2006), 1573–1603, arXiv:math-ph/0501074.

[12] Claeys T., Kuijlaars A.B.J., Vanlessen M., Multi-critical unitary random matrix ensembles and the general Painlevé II equation, *Ann. of Math.* **168** (2008), 601–641, arXiv:math-ph/0508062.

[13] de Bruijn N.G., On some multiple integrals involving determinants, *J. Indian Math. Soc.* **19** (1955), 133–151.

[14] Deift P., Orthogonal polynomials and random matrices: a Riemann–Hilbert approach, *Courant Lect. Notes Math.*, Vol. 3, American Mathematical Society, Providence, RI, 1999.

[15] Deift P., Riemann–Hilbert problems, in *Random Matrices*, IAS/Park City Math. Ser., Vol. 26, American Mathematical Society, Providence, RI, 2019, 1–40, arXiv:1903.08304.

[16] Deift P., Gioev D., Universality at the edge of the spectrum for unitary, orthogonal, and symplectic ensembles of random matrices, *Comm. Pure Appl. Math.* **60** (2007), 867–910, arXiv:math-ph/0507023.

[17] Deift P., Gioev D., Universality in random matrix theory for orthogonal and symplectic ensembles, *Int. Math. Res. Pap. IMRP* **2007** (2007), 004, 116 pages, arXiv:math-ph/0411075.

[18] Deift P., Gioev D., Random matrix theory: invariant ensembles and universality, *Courant Lect. Notes Math.*, Vol. 18, American Mathematical Society, Providence, RI, 2009.

[19] Deift P., Gioev D., Kriecherbauer T., Vanlessen M., Universality for orthogonal and symplectic Laguerre-type ensembles, *J. Stat. Phys.* **129** (2007), 949–1053, arXiv:math-ph/0612007.

[20] Deift P., Kriecherbauer T., McLaughlin K.T.-R., Venakides S., Zhou X., Asymptotics for polynomials orthogonal with respect to varying exponential weights, *Int. Math. Res. Not.* **1997** (1997), 759–782.
[21] Deift P., Kriecherbauer T., McLaughlin K.T.-R., Venakides S., Zhou X., Strong asymptotics of orthogonal polynomials with respect to exponential weights, *Comm. Pure Appl. Math.* **52** (1999), 1491–1552.

[22] Deift P., Kriecherbauer T., McLaughlin K.T.-R., Venakides S., Zhou X., Uniform asymptotics for polynomials orthogonal with respect to varying exponential weights and applications to universality questions in random matrix theory, *Comm. Pure Appl. Math.* **52** (1999), 1335–1425.

[23] Deift P., Zhou X., A steepest descent method for oscillatory Riemann–Hilbert problems. Asymptotics for the MKdV equation, *Ann. of Math.* **137** (1993), 295–368.

[24] Dyson F.J., The threefold way. Algebraic structure of symmetry groups and ensembles in quantum mechanics, *J. Math. Phys.* **3** (1962), 1199–1215.

[25] Dyson F.J., A class of matrix ensembles, *J. Math. Phys.* **13** (1972), 90–97.

[26] Eynard B., Asymptotics of skew orthogonal polynomials, *J. Phys. A* **34** (2001), 7591–7605, arXiv:cond-mat/0012046.

[27] Fokas A.S., Its A.R., Kitaev A.V., Discrete Painlevé equations and their appearance in quantum gravity, *Comm. Math. Phys.* **142** (1991), 313–344.

[28] Forrester P.J., Log-gases and random matrices, *London Math. Soc. Monogr. Ser.*, Vol. 34, Princeton University Press, Princeton, NJ, 2010.

[29] Ghosh S., Generalized Christoffel–Darboux formula for skew-orthogonal polynomials and random matrix theory, *J. Phys. A* **39** (2006), 8775–8782, arXiv:math-ph/0603013.

[30] Its A.R., The Riemann–Hilbert problem and integrable systems, *Notices Amer. Math. Soc.* **50** (2003), 1389–1400.

[31] Kriecherbauer T., Shcherbina M., Fluctuations of eigenvalues of matrix models and their applications, arXiv:1003.6121.

[32] Krüger T., Lee S.-Y., Yang M., Local statistics in normal matrix models with merging singularity, arXiv:2306.12263.

[33] Lee S.-Y., Yang M., Planar orthogonal polynomials as type II multiple orthogonal polynomials, *J. Phys. A* **52** (2019), 275202, 14 pages, arXiv:1801.01084.

[34] Martínez-Finkelshtein A., Van Assche W., What is ... a multiple orthogonal polynomial?, *Notices Amer. Math. Soc.* **63** (2016), 1029–1031, arXiv:1707.09511v1.

[35] Mehta M.L., Random matrices, 3rd ed., *Pure Appl. Math. (Amst.)*, Vol. 142, Elsevier, Amsterdam, 2004.

[36] Muskelishvili N.I., Singular integral equations. Boundary problems of functions theory and their applications to mathematical physics, Springer, Dordrecht, 1958.

[37] Pastur L., Shcherbina M., Universality of the local eigenvalue statistics for a class of unitary invariant random matrix ensembles, *J. Stat. Phys.* **86** (1997), 109–147.

[38] Pierce V.U., A Riemann–Hilbert problem for skew-orthogonal polynomials, *J. Comput. Appl. Math.* **215** (2008), 230–241, arXiv:math/0610592.

[39] Shcherbina M., Edge universality for orthogonal ensembles of random matrices, *J. Stat. Phys.* **136** (2009), 35–50, arXiv:0812.3228.

[40] Shcherbina M., On universality for orthogonal ensembles of random matrices, *Comm. Math. Phys.* **285** (2009), 957–974, arXiv:math-ph/0701046.

[41] Shcherbina M., Orthogonal and symplectic matrix models: universality and other properties, *Comm. Math. Phys.* **307** (2011), 761–790, arXiv:1004.2765.

[42] Stojanovic A., Universality in orthogonal and symplectic invariant matrix models with quartic potential, *Math. Phys. Anal. Geom.* **3** (2000), 339–373.

[43] Van Assche W., Geronimo J.S., Kuijlaars A.B.J., Riemann–Hilbert problems for multiple orthogonal polynomials, in *Special Functions 2000: Current Perspective and Future Directions*, NATO Sci. Ser. II Math. Phys. Chem., Vol. 30, Kluwer, Dordrecht, 2001, 23–59.

[44] Widom H., On the relation between orthogonal, symplectic and unitary matrix ensembles, *J. Stat. Phys.* **94** (1999), 347–363, arXiv:solv-int/9804005.