Skew-orthogonal polynomials in the complex plane and their Bergman-like kernels

Non-Hermitian random matrices with symplectic symmetry provide examples for Pfaffian point processes in the complex plane. These point processes are characterised by a matrix valued kernel of skew-orthogonal polynomials. We develop their theory in providing an explicit construction of skew-orthogonal polynomials in terms of orthogonal polynomials that satisfy a three-term recurrence relation, for general weight functions in the complex plane. New examples for symplectic ensembles are provided, based on recent developments in orthogonal polynomials on planar domains or curves in the complex plane. Furthermore, Bergman-like kernels of skew-orthogonal Hermite and Laguerre polynomials are derived, from which the conjectured universality of the elliptic symplectic Ginibre ensemble and its chiral partner follow in the limit of strong non-Hermiticity at the origin. A Christoffel perturbation of skew-orthogonal polynomials as it appears in applications to quantum field theory is provided.


Introduction
The study of orthogonal and skew-orthogonal polynomials in the complex plane is closely related to the question of integrability of determinantal and Pfaffian point processes in the plane. Many of the known examples for such point processes can be realised as complex eigenvalues of non-Hermitian random matrices, such as the three Ginibre ensembles of Gaussian random matrices with real, complex or quaternion matrix elements [1,2,3], or their three chiral counterparts [4,5,6]. Compactly supported examples include the truncation of random orthogonal [7], unitary [8] and symplectic matrices [9], and we refer to [10] for a review on non-Hermitian random matrices. At the same time these point processes are examples for two-dimensional Coulomb gases in a confining potential with specific background charge, and we refer to [11] for details. A further example in this class is the circular quaternion ensemble of random matrices belonging to the symplectic group Sp(2N ), distributed according to Haar measure, see [12] for details. Its joint density of eigenvalues distributed on the unit circle represents a Pfaffian point process that is also determinantal.
Apart from this direct statistical mechanics interpretation as a Coulomb gas, further applications of these point processes include dissipative quantum maps [13], dynamical aspects of neural networks [14], properties of the quantum Hall effect [15] and quantum field theories with chemical potential [4,5,6]. In quantum optics the Bergman kernel of planar Hermite polynomials plays an important role in the construction of coherent and squeezed states [16]. Notably, in the symplectic symmetry class, that will be our focus here, there is a map from the symplectic Ginibre ensemble to disordered non-Hermitian Hamiltonians in an imaginary magnetic field [17]. The circular quaternion ensemble finds applications in the computation of thermal conduction in superconducting quantum dots [18]. The predictions of the non-Hermitian symplectic chiral symmetry class [5] were successfully compared with data from lattice simulations in Quantum Chromodynamics with two colours at non-vanishing chemical potential [19]. Here, the insertion of quark flavours, given in terms of characteristic polynomials in the random matrix setting, play an important role and were tested in [19]. We will investigate the effect of these insertions on the underlying skew-orthogonal polynomials (SOP) also known as Christoffel perturbation.
The real and symplectic Ginibre ensemble differ from the complex ensemble in the following way. First, in the two former ensembles eigenvalues come in complex conjugated pairs. Second, the real ensemble shows an accumulation of eigenvalues on the real axis, while the symplectic ensemble shows a depletion, as the probability to have real eigenvalues is zero. Let us briefly recall what is known about the symplectic Ginibre ensemble in the limit of large matrix size (or number of particles). It is not surprising that the gap probability [20] and all correlation functions at the origin [20,21] differ from the complex Ginibre ensemble. The same statement holds for the density at weak non-Hermiticity [17]. Below, we will present a proof that in the limit of strong non-Hermiticity the correlation functions at the origin of the symplectic Ginibre ensemble are universal, in the sense that they hold for the elliptic ensemble beyond the rotationally invariant case. This was conjectured in [17], and the same conjecture [5] for the chiral ensemble will be shown as well.
In contrast, the distribution of the largest eigenvalue in radius [22] as well as the local radial density in the bulk away from the real line agree with those of the complex Ginibre ensemble [23]. Without integration over the angles this agreement has been shown more recently [24] for all correlation functions (marginals) in the bulk of the spectrum of the symplectic Ginibre ensemble, away from the real axis. For the same agreement between the real and complex Ginibre ensemble we refer to [25]. This is in strong contrast with the Hermitian ensembles of random matrices, corresponding to a socalled Dyson gas on the real line at inverse temperature β = 1, 2, 4. The local statistics for the latter three ensembles differs everywhere in the spectrum, see e.g. [20,11] for a summary of results. It is one of the goals of this article to develop the theory of SOP in the complex plane. This will allow us to construct further examples of Pfaffian point processes with symplectic symmetry that are integrable, where this extended universality in the complex plane can be studied.
What is known for the construction of SOP for general ensembles with symplectic symmetry? In [21] it was shown that both polynomials of odd and even degree enjoy a Heine-like representation. It is given as the expectation value of a characteristic polynomial (times a trace in the case of polynomials of odd degree), that is a multiple integral representation of the order of the degree of the SOP. Only in a limited number of cases have these been used for an explicit construction, using Schur polynomials [26] or Grassmann integrals [6]. A second construction sets up a Gram-Schmidt skew-orthogonalisation procedure [27] that we will recall below. However, also this is of limited use for an explicit construction. Current explicit examples for SOP include Hermite [21] and Laguerre polynomials [5] for the elliptic symplectic Ginibre ensemble and its chiral counterpart, respectively.
In this article we will exploit the particularly simple structure of the skew-product in symplectic ensembles. It allows to relate the skew-product to the multiplication acting on the standard Hermitian inner product. It is well known that on subsets of the real line orthogonal polynomials (OP) always satisfy a three-term recurrence relation with respect to multiplication. In the plane this is no longer true, and on bounded domains we may not expect any finite-term recurrence in general [28]. It was shown much later on bounded domains with flat measure that if such a recurrence exists, the domain is an elliptic disc and the depth of recurrence is three, see [29] and references therein. However, in the weighted case the ellipse is no longer special, as also here OP without recurrence exist [30]. On the other hand, if we can promote classical OP from the real line to the complex plane, then the existence of such a recurrence is always guaranteed.
The remainder of this article is organised as follows. In the next Section 2 we define the class of ensembles of complex eigenvalues with symplectic symmetry we consider here, including their skewproduct. We recall the Gram-Schmidt skew-orthogonalisation that leads to the reproducing polynomial kernel, in terms of which all k-point correlation functions or marginals are given. In Section 3 an explicit construction is provided for SOP for a general class of weight functions. They are characterised by OP that satisfy a three-term recurrence relation with real coefficients, orthogonal with respect to the same real valued weight function. Several new examples of resulting planar SOP are presented, including a weight of Mittag-Leffler type in the plane and weights on an elliptic disc that lead to planar Gegenbauer SOP, as well as a subfamily of non-symmetric Jacobi polynomials that include Chebyshev polynomials. A realisation of the Chebyshev polynomials as Szegő SOP on an ellipse is given, too. In Appendix A we recover some known planar SOP from our construction. Section 4 is devoted to the derivation of the Bergman-like kernel for Hermite and Laguerre type SOP. These formulas imply the large-N limit at the origin of the spectrum for all correlation functions at strong non-Hermiticity in the corresponding ensembles. This confirms their conjectured universality within these two elliptic classes. In Section 5 a Christoffel perturbation of SOP is provided for general weight functions. As a corollary we obtain that a Christoffel perturbation does not preserve the three-term recurrence relation for measures on C. The corresponding Fourier coefficients are provided in Appendix B. Appendix C contains a collection of integrals needed throughout the article.

Acknowledgments
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -SFB 1283/2 2021 -317210226 "Taming uncertainty and profiting from randomness and low regularity in analysis, stochastics and their applications" (G.A. and I.P.) and IRTG2235 "Searching for the regular in the irregular: Analysis of singular and random systems" (M.E.), and by the grants DAAD-CONICYT/Becas Chile, 2016/91609937 and FWO research grant G.0910.20 (I.P.). We thank Boris Khoruzhenko for useful discussions and Sung-Soo Byun for discussions and many useful comments on this manuscript.

Symplectic ensembles and skew-orthogonal polynomials
The point processes on a subset of the complex plane with symplectic symmetry considered in this paper are defined by the following joint probability density Here, µ is a positive Borel measure on the domain D in the complex plane. Further conditions on the measure will be specified at the beginning of Subsection 2.1. The normalization constant Z N (partition function) is given by 2) The joint density (2.1) may result for example from the distribution of the 2N eigenvalues (z 1 , z 1 , . . . , z N , z N ) of an N × N dimensional quaternionic non-Hermitian random matrix (or its 2N dimensional complex representation). Examples for such random matrix realisations include the elliptic quaternionic Ginibre ensemble [21] and its chiral counterpart [5]. When the measure µ is supported on the unit circle, further representatives include the circular quaternion ensemble with a flat measure, see [12,Thm. 3.1]. It is not difficult to see, that for z k = e iθ k the first product in (2.1) leads to the Vandermonde determinant squared in the variables x k = cos θ k , and thus to a determinantal point process. The second product in (2.1), N k=1 4(1 − x 2 k ), is then taken as part of the weight function on (−1, 1). In the same way also the circular real ensembles of the Haar distributed groups SO(N ) for even and odd N provide examples for this Pfaffian point process [12,Thm. 3.1], with varying weight functions though. All of these are also determinantal.
We will be more general here in taking the complex eigenvalue model (2.1) as a starting point. In general, the ensemble defined in (2.1) is a Pfaffian point process, cf. [20,21]. Defining the k-point correlation functions (or marginal measures) as they can be shown to take the form [21] dR N, Here, K N is the 2 × 2 matrix-valued kernel of the corresponding point process. Equivalently, when µ has the density w with respect to the volume element dv on D, we have where the 2 × 2 matrix-valued kernel K N is defined as In that case we have K N (z, u) = w(z)w(u) K N (z, u). The function σ N (z, u) is called the pre-kernel or polynomial kernel, to be defined in (2.17) in terms of SOP. Compared with [21] we have taken the pre-factors (z − z)(u − u) out of the Pfaffian, to avoid cuts for the square root of these factors. As an example we get for the one-point and two-point function The pre-kernel can be expressed in terms of SOP to be defined below, being a particularly simple choice in a more general construction, see [21]. There, it is given in terms of the inverse moment matrix, where the de Bruijn integral formula is applied together with the fact that the joint density (2.1) is proportional to the product N j=1 (z j − z j ) times the Vandermonde determinant of size 2N of all eigenvalues and their complex conjugates.

Skew-orthogonal polynomials
From now on let µ be a positive Borel measure on C, with an infinite number of points in its support D, and such that |z| m dµ(z) < ∞ for all non-negative integers m. Further we require that D is symmetric about the real axis, i.e. z ∈ D if and only if z ∈ D. For any f, g ∈ C[z], we define the following skew-symmetric form (2.7) Equivalently, ·, · s is an alternating form. In particular, when the polynomials f, g have real coefficients, ·, · s is also a skew-Hermitian form.
Definition 2.1. A family of polynomials (q n ) n∈N with deg q n = n is called skew-orthogonal corresponding to µ, if they satisfy for all non-negative integers k, l ∈ N: with r k > 0 being their skew-norms.
The SOP are called monic if their leading coefficient is unity, i.e. q n (z) = z n + O(z n−1 ). Notice that this choice of the leading coefficient does not make the SOP unique. The reason is that the odd polynomials q 2n+1 can be modified by adding an arbitrary multiple of the even polynomial q 2n , q 2n+1 (z) = q 2n+1 (z) + d n q 2n (z), without changing the skew-orthogonality relations (2.8) and (2.9). Is it not difficult to see, using the Heine-like representation of the SOPs given in terms of 2n−fold integrals in [21], that different transformations do not preserve the skew orthogonality conditions (2.8) and (2.9). Hence, by imposing the SOP to be monic and the coefficient of z 2n in q 2n+1 (z) to be zero, we fix this ambiguity and make them unique, see Lemma 2.2 below. However, sometimes it is convenient to choose the coefficient of z 2n to be non-vanishing, in order to obtain closed formulas for the odd SOP, see Section 3 for various examples. The existence of the SOP is guaranteed by Gram-Schmidt skew-orthogonalisation to be provided below in Theorem 2.4. As a consequence of Gram-Schmidt skew-orthogonalisation and the definition of the skew symmetric form ·, · s , the SOP belong to the polynomial ring R[z]. As usual, we denote by R n [z] the space of polynomials with real coefficients and degree at most n . Lemma 2.2. Let (q n ) n∈N be a sequence of monic skew orthogonal polynomials, such that the coefficient z 2n in q 2n+1 (z) is zero. Then, the sequence of SOPs in Definition 2.1 is unique.

Remark 2.3.
When the measure µ has density w on some domain D, then we talk about SOP with respect to the weight function w.
A closed expression can be obtained for these SOP in terms of the skew-complex moments of µ and the real skew-symmetric Gram matrix (2.12) De Bruijn's integration formula implies that the Pfaffian of the skew-symmetric Gram matrix satisfies with Z k defined in (2.2). In analogy with the OP, the skew-complex moments already determine the SOP. In terms of these one can give explicit Pfaffian formulae for the orthonormalised polynomials.

14)
and for the odd degreeq 2k+1 , by replacing the (2k + 2)nd row and column by powers of z except z 2k , as well as shifting the index of the Gram matrix in the (2k + 1)st row and column by one, In particular this implies that the SOP have real coefficients. The proof of this theorem involves properties of Pfaffians. It was presented first in [31] for the case of SOP over the real line 1 , and extended to planar SOP in [27]. Here, we emphasise that this construction holds without specifying the support of the measure µ. Remark 2.5. We stress that the SOP always exist, under the minimal assumptions on the measure that we stated at the beginning of this subsection, even without the support being symmetric about the real axis. Also note that, by using de Bruijn's formula, we can easily express the partition function (2.2) in terms of the skew-norms from (2.9): see also [20,Sec. 15.2] for the symplectic Ginibre case. Given the positivity of Z N for all N this automatically implies the positivity of the skew-norms r k . Conversely, all examples in Section 3 where we determine these skew-norms explicitly for a given measure provide instances of Selberg-like integrals.
The representation of the SOP in terms of Gram matrices may be important from a theoretical point of view, but is not very useful for the actual computation of the SOP since it involves the evaluation of Pfaffians. In Section 3 we will propose a more explicit construction, given an orthonormal system in L 2 (dµ) that satisfies certain properties.
Let P 2n be the space of analytic polynomials of degree at most 2n − 1 (i.e. p ∈ P 2n with ∂ ∂z p = 0), and we equip P 2n with the skew-product ·, · s . Lemma 2.6. The polynomial kernel σ n (u, v) defined as follows is the reproducing skew-kernel of (P 2n , ·, · s ).
Proof. Let f ∈ P 2n and σ v (u) := σ n (u, v). Without restriction we can assume that f = q 2m or f = q 2m+1 with m ≤ n − 1, thus Remark 2.7. Lemma 2.6 tells us that σ n reproduces itself, σ v , σ u s = σ v (u) = σ n (u, v), and that σ n is a skew-symmetric function, σ n (u, v) = −σ n (v, u). These two properties are in complete analogy to the Hermitian inner product space L 2 (dµ) with scalar product ·, · , where the polynomial kernel is expressed in terms of the orthonormal polynomials P n corresponding to µ. It satisfies the reproducing property K v , K u = K n (u, v) and is a Hermitian function, i.e. K n (v, u) = K n (u, v).
The polynomial kernel σ n (u, v) is not affected by the non-uniqueness of the odd SOP. It remains unchanged if we redefineq 2n+1 (z) = q 2n+1 (z) + d n q 2n (z) as the latter terms drop out in (2.17), due to anti-symmetry. Remark 2.8. As was shown in [21] the pre-kernel σ n (u, v) is the main input of the 2×2 matrix-valued kernel K n in (2.5).

Remark 2.9.
Let c(u, v) = g(u)g(v), with g a continuous unimodular function, such that g(u) = 1/g(u). Then, the Pfaffian of a 2k × 2k skew-symmetric matrix (a i,j ) 1≤i,j≤2k remains unchanged if we multiply each a i,j by c(u i , u j ), where the u 1 , . . . , u 2k come in complex conjugate pairs, u k+i = u i . In particular, σ(u, v) = g(u)g(v) σ(u, v) is another pre-kernel giving rise to the same correlation functions R N,k . Two such pre-kernels as well as their corresponding kernels are then called equivalent.

Construction of skew-orthogonal polynomials
In this section we reduce the construction of the SOP corresponding to µ to that of OP corresponding to the same measure µ, that satisfy a (suitable) three-term recurrence relation. Let us emphasise that, in general, the construction of an orthonormal system in L 2 (C, dµ) is straightforward using Gram-Schmidt orthogonalisation, see [32]. However, generic OP do not satisfy such a recurrence, see [28] for a recent discussion. On the other hand all (polynomial) orthogonal systems in L 2 (R, dν) do satisfy such a three-term recurrence with real coefficients, and thus provide many potential candidates. If we can thus find a measure µ with support on the complex plane, such that the same orthogonal system in L 2 (R, dν) gives rise to an orthogonal system in L 2 (C, dµ), our Theorem 3.1 below immediately leads to SOP corresponding to the same measure µ. Below we will give several examples of such a situation.
Let µ be a measure that satisfies the properties from the previous Subsection 2.1, then by the Gram-Schmidt process one can construct a unique sequence of polynomials that form an orthogonal system in L 2 (dµ), Due to the linearity of the Hermitian form ·, · over R we can assume that γ n ≡ 1. In this case we say that the sequence (p n ) n∈N of OP is chosen in monic normalisation (i.e. p n (z) = z n + . . . ). These OP (corresponding to µ) satisfy a three-term recurrence relation with real coefficients, if The condition for b k and c k to be real is obviously equivalent to the condition that the planar OP have real coefficients, p n (z) = p n (z) for all n ∈ N.
Conversely, if the sequence (p n ) n∈N of monic OP satisfies a three-term recurrence relation (3.3), then, the sequence of monic polynomials defined in (3.4) forms a SOP system and the µ n,j 's in (3.4) can be explicitly computed: Proof. Assuming that the sequence in (3.4) satisfies the skew-orthogonality conditions (2.8) and (2.9), we want to show the three-term recurrence relation for the OP. Therefore consider the following Fourier expansion, starting with the odd polynomials and expressing zp 2m+1 (z) in terms of the p l 's, where we have to show that is the null polynomial. From p 2n+1 = q 2n+1 we have for all n, m Choosing in particular n < m, the first term vanishes due to orthogonality, zp 2n+1 , p 2m+1 = 0. From Consequently the odd coefficients vanish, due to h k = 0 for all k. Next, let us show that all even Fourier coefficients vanish too, a 2m+1,2l = 0 for 0 ≤ l ≤ m − 1. For this, we use that 0 = q 2l , q 2m+1 s for 0 ≤ l ≤ m − 1, and the fact that q 2l is a linear combination of even polynomials p 2j , see (3.4). Then, by inspection starting with l = 0, this implies that 0 = p 2l , p 2m+1 s for 0 ≤ l ≤ m − 1, and therefore for any linear combination of the p 2l up to degree 2m − 2 as in f 2m−1 .
The first scalar product is trivially zero due to orthogonality, and the second term yields 0 = f 2m−1 2 , the desired result. Now we turn to the Fourier expansion of the even polynomials zp 2m (z) given by By assumption µ m,j = λ m−1 µ m−1,j and that λ m−1 is independent of j, this implies that q 2m (z) = p 2m (z)+λ m−1 q 2(m−1) (z). Using 0 = q 2n , q 2m s for all n, m, we obtain 0 = q 2n , p 2m + λ m−1 q 2(m−1) s = q 2n , p 2m s . Starting with n = 0, by inspection we have 0 = p 2n , p 2m s . As before (choosing in particular n < m) this implies that all even Fourier coefficients a 2m,2n vanish for n = 0, . . . , m − 1.
Next we apply the skew-orthogonality 0 = q 2n+1 , q 2m s for n < m. In particular for 0 ≤ l ≤ m − 2, we have 0 = q 2l+1 , q 2m s = q 2l+1 p 2m + λ m−1 q 2(m−1) s = p 2l+1 , p 2m s . Therefore, any linear combination of the p 2l+1 up to degree 2m−3 (as in f 2m−2 ) will give zero in the skew-product. Applying the same argument as for the case of odd polynomials, we conclude that 0 = f 2m−2 2 . This concludes the proof of the first part of the theorem.
Conversely, let us assume that the OP satisfy a three-term recurrence relation, and define the (q n ) n by (3.4). For any polynomials f, g with real coefficients, the evaluation of the skew-product (2.7) can be reduced to the evaluation of the scalar product (3.2), f, g s = 2 Re[ zf , g − f, zg ]. Thus, for the polynomials with odd indices we can write To calculate the two scalar products we use the three-term recurrence relation and the orthogonality of the polynomials p n . This leads to zp 2k+1 , p 2l+1 = p 2k+1 , zp 2l+1 = b 2k+1 h 2k+1 δ k,l . Hence, q 2k+1 , q 2l+1 s = 0 for all k and l. Similarly, p 2i , zp 2j = zp 2i , p 2j = b 2i h 2i δ i,j , and after expanding the q 2n 's in terms of the p 2n 's, we obtain q 2k , q 2l s = 0 ∀k, l.
For the combination of odd and even indices in (2.9) we need to evaluate Depending on the values of k and l we need to distinguish between 3 cases: 1. Case k < l: Here, we have δ j,l = 0 and δ j,l+1 = 0 for all j ≤ k, and therefore q 2k , q 2l+1 s = 0 as claimed.

Case k = l:
It holds δ j,l = 1 only for j = l = k, and δ j,l+1 = 0 for all j ≤ k, therefore

Case k > l:
We have δ j,l = 1 only for j = l < k, and δ j,l+1 = 1 only for j = l + 1 ≤ k. Then, we want that is zero. This equation is fulfilled if and only if the µ k,l satisfy Thanks to our monic choice of the leading coefficient µ k,k = 1, this recursion can be solved to obtain the explicit form as claimed In particular, note that the µ k,j satisfy the assumption made in the first part of the theorem.
Remark 3.2. The recurrence coefficients b k (although non-zero in general) are not needed to express µ k,j (3.6) and r k (3.5). Therefore, in the examples below, we will only give formulas for c k and h k . Also note that we can invert (3.4) to get the following representation of the OP in terms of the SOP for all k ∈ N: For radially symmetric weight functions w(z) = w(|z|), the OP are given by monomials. Due to c k = 0 (= b k ) for all k ∈ N in that case, the above formulas simplify considerably (this result is known, see [23,Eqs. (31), (32)]): If the weight function of the planar OP in Theorem 3.1 is radially symmetric, w(z) = w(|z|), we have for the OP and their squared norms: where χ Dr denotes the characteristic function on D in radial direction. For the planar SOP we obtain: Let us give some examples for planar SOP (and SOP over a weighted analytic Jordan curve) that are constructed using Theorem 3.1 or Corollary 3.3. These immediately lead to examples of Pfaffian point processes (2.1) that are integrable, in the sense that the SOP and thus the corresponding (pre-) kernel is known explicitly. We will only state the SOP in what follows, and not the pre-kernel.

Example 3.4 (Gegenbauer ensemble).
Consider the measure µ supported in the interior of the standard ellipse ∂E, with semi-axes a > b > 0, such that µ has density function w with respect to planar Lebesgue measure: As shown in [30,Thm. 3.1], the Gegenbauer polynomials (also called ultraspherical) form an orthonormal basis on the Bergman space with respect to this weight function. The monic OP, recurrence coefficients and norms read (3.20) From Theorem 3.1 we obtain

23)
with coefficients and skew-norms In the special case when α = 0 the weight function is constant and the Gegenbauer polynomials reduce to the Chebyshev polynomials of the second kind, C (1) n (x) = U n (x). Apart from the example with Gegenbauer polynomials, which are symmetric Jacobi polynomials, for non-symmetric Jacobi polynomials with parameters (α + 1/2, ±1/2) the SOP can also be constructed. In particular they include the Chebyshev polynomials of the first up to fourth kind, see [30] and references therein.
The SOP coefficients and skew-norms take the form: . (3.26) They also can be obtained directly from our Corollary 3.3, since in this limit the weight is rotationally invariant. This weight describes the truncated unitary ensemble [8]: for integer values of α = N − M − 1 it appears in the eigenvalue distribution of Haar distributed unitary matrices of size N truncated to the upper sub-block of size M < N . The SOP can be used for the corresponding symplectic point process, see [9] for further details.
Also, when the weight function includes an extra charge insertion at the origin, that is, w(z) = |z| 2c (1 − |z| 2 ) α , α > −1, c > −1, the µ j,k and r k can be readily obtained. Example 3.6 (Mittag-Leffler ensemble). The weight function, that leads to a Bergman kernel in terms of the two-parametric Mittag-Leffler function [34,Thm 1.6], is given by including a charge insertion at the origin. When setting λ = 1, it reduces to the induced Ginibre ensemble, see [35] for the OP, and the SOP were determined in [5]. Setting also c = 0, it reduces to the standard Ginibre ensemble, see [20,21]. The squared norms h n of the monic OP p n (z) = z n with respect to the weight (3.27) can be easily obtained: From Corollary 3.3 we can read off the monic planar SOP with respect to (3.27) and their skew-norms. For λ = 1 and c = 0 this leads to the same pre-kernel as in the symplectic Ginibre ensemble [20].
Example 3.7 (Chebyshev on an elliptic contour). Consider the measure µ supported on the standard ellipse ∂E, with semi-axes a > b > 0, such that µ has the density function: Chebyshev polynomials of the third kind V n satisfy an orthogonality relation on this contour [36] (here |dz| stands for the arc length measure): Similar relations hold for Chebyshev polynomials of first, second and fourth kind for different weight functions [36]. From Theorem 3.1 we get the following Szegő SOP

32)
and for the skew-norms we obtain Note that in the limit b → a = 1 this example simplifies to the uniform weight on the unit circle, for which the Pfaffian point process is already well known [12]: it describes the eigenvalues of the circular quaternion ensemble (matrices of Sp(2N ) distributed according to Haar-measure). As mentioned already, this point process is determinantal and can be analysed via OP.
Further examples -including rotationally invariant weights on C -will be provided in Section 4 and Appendix A. The case resulting from products of M rectangular Ginibre matrices is deferred to the Appendix A as the resulting planar SOP are known [23].

Bergman-like kernel of skew-orthogonal polynomials
In this section we will derive the Bergman-like kernel for planar skew-orthogonal Hermite and Laguerre polynomials in separate subsections. They are given by the limiting pre-kernel, the sum over orthonormal SOP, hence the terminology. In both cases the corresponding weights are given by a one parameter family of elliptic ensembles, see e.g. (4.1) below for Hermite and Appendix A. The proofs will be based on the rotationally invariant limit, that is when the parameter is chosen such that the underlying domain and weight function has rotational symmetry. In that case the corresponding limiting pre-kernels are known, see [21] and [5] respectively. As a consequence, at the end of each subsection we will present the universality of a one-parameter family of kernels in such symplectic elliptic Ginibre ensembles, hinting at a much larger universality.

Bergman-like kernel of skew-orthogonal Hermite polynomials
In random matrix theory the planar Hermite polynomials appear in the solution of the elliptic complex Ginibre ensemble [37], with the one-parameter complex potential The monic polynomials (here with their third recurrence coefficient c n ), Here, dA(z) denotes the area measure on the complex plane, divided by π √ 1 − τ 2 to provide a probability measure, i.e. with h 0 = 1. We refer to Example A.2 for details, including a two-parameter complex normal distribution and the matching planar Hermite polynomials. The corresponding Bergman kernel [38] is the standard one of Hermite polynomials on R [39, 18.18.28], after continuation in the arguments ζ, η ∈ C, 0 ≤ τ < 1: This identity is also called Mehler-Hermite formula or Poisson kernel. The Hermite SOP and skew-norms with respect to the weight (4.1) are known [21] and recollected in Example A.2, as following from Theorem 3.1. We only give the resulting pre-kernel: (4.5) Here, the area measure is normalised by 1−τ 2 , to achieve r 0 = 1 in analogy to above. Both expressions (4.4) and (4.5) are the sum over (skew-)orthonormal polynomials. We are thus led to consider the following limit which is the first main result of this section. From now on we use the following notation f N (z) ⇒ f (z) to express that the sequence of functions (f N ) N converges to f uniformly on any compact subset of C as N → ∞. Theorem 4.1 (Bergman-like Hermite kernel). Let 0 < τ < 1 and z, u ∈ C, then we have that σ τ,N (z, u) ⇒ S τ (z, u), given by (4.6) We use the standard notation for the error function, in the form erf(x) = 2x √ π 1 0 e −x 2 s 2 ds, that can be continued to x ∈ C. The result (4.6) was already given in [10,Eq. (18.6.53)], without providing any details. In analogy to (4.4) being the infinite sum over orthonormalised OP, we call the corresponding sum over orthonormal SOP (4.6) Bergman-like Hermite kernel.
The proof of Theorem 4.1 will be in two steps. First, we recall the rotationally invariant case τ = 0, see Remark 4.3. In the second step we derive (4.6) by using an integral representation for the Hermite polynomials together with Lemma 4.2.
Then, as N → ∞ we have g N (u, v) ⇒ g(u, v), where the limiting function is given by Proof. Let u, v ∈ B(0, r), where B(a, r) denotes the disk of center a and radius r. Then each summand in (4.7) is bounded by re r 2 r 2k /(k!). Hence, by the Weierstraß M-test, its sum is an analytic function of u and v, being absolutely and uniformly convergent in each compact subset of the plane. Thus, the limiting function, as N → ∞, is given by the power series Now, we take the derivative term-wise to obtain: A convenient initial value of this linear inhomogeneous ordinary differential equation is u 0 = 0, then g(u 0 , v) = 0 for every v ∈ C. The solution of the initial value problem is given by (see [40]) (4.10) In the last step we have used the definition of the error function provided below (4.6), as well as its anti-symmetry.
This is a well-known limiting pre-kernel, first found by Mehta in [20] and later calculated in [21]. As can be seen from (4.1), the parameter τ controls the degree of Hermiticity of the elliptic ensemble. The case τ = 0 here corresponds to maximal non-Hermiticity, when real and imaginary part share the same variance. In Theorem 4.1 the limit N → ∞ is taken at fixed τ , which we call strong non-Hermiticity. The case when τ → 1 at a rate depending on N called weak non-Hermiticity will not be discussed further, and we refer to [37,21] for details. The kernel in (4.11) is the symplectic analogon of its Hermitian partner exp [uv] in the holomorphic Fock-space [41] in L 2 (e −|z| 2 dA(z)), as in Remark 2.7.
which we use to derive Finally we can complete the proof of our main result.
Proof of Theorem 4.1. As for (4.16) we can obtain a simple estimate for the sum in (4.6) with ζ, η ∈ B(0, r). Thus, the sum converges absolutely and uniformly on B(0, r), and therefore on each compact subset of C. The absolute convergence allows us to rearrange the summands and, using Lemma 4.4 with ϕ = φ = τ , we obtain: Remark 4.5. We note that the Bergman kernel K τ in (4.4) of the analytic space in L 2 (e −Qτ dA), and the Bergman-like kernel S τ in (4.6), both restricted to R × R, are related by the following differential equation

Universality of the symplectic elliptic Ginibre kernel
In this subsection we will prove the universality of all k-point correlation functions (2.4) in the symplectic elliptic Ginibre ensemble, in the large-N limit at strong non-Hermiticity close to the origin. Throughout this article we have considered weight functions that are N -independent. Therefore, at large-N the eigenvalues condense on a droplet with N -dependent radius, the support of the equilibrium measure. We refer to [42] for a discussion of the equilibrium problem for a general potential. In the case of our elliptic weight (4.1), we have that ([43, Thm. 2.1] applied to the result of [44]) yields The fact that in Theorem 4.1 we take the limit of the kernel at arguments independent of N , implies that we consider the vicinity of the origin. In contrast, to investigate the neighbourhood of a bulk (or edge) point, we would have to centre the arguments around this point z = √ 2N z 0 , with |z 0 | of order unity, and rescale accordingly. Furthermore, because we keep τ fixed when N → ∞, this implies that we consider strong non-Hermiticity. Let us first quote the known result at maximal non-Hermiticity τ = 0 from [21] 2 , the symplectic Ginibre kernel at the origin. In that case, from Remark 4.3 we can read off the matrix elements of the limiting kernel of K N (z, u) in (2.5), times the normalisation from the area measure As a corollary from Theorem 4.1, we can now prove the following universality statement. Proof. In order to compare the limits of the pre-kernels (4.5) and (4.21) and their pre-factors, we have to map to the same equilibrium measure (also called macroscopic density) first, implying the same mean level spacing. This is called unfolding and we refer to [13] for a standard reference that includes complex spectra, see also [42] for details about recentering and rescaling. Because the symplectic Ginibre kernel in (4.21) is given with respect to the limiting density 1 2π , in this case to unfold we simply have to rescale all arguments z → √ 1 − τ 2 z, in view of (4.20). Therefore we consider the limit The additional factor (1 − τ 2 ) 3 2 in the first line originates from the rescaling of the measure dA(z) and the factors (z − z) in the Pfaffian representation of the correlation functions, eq. (2.4). In the first line we also multiply with the τ -dependent normalisation of the area measure. After inserting (4.6) in the second line, we arrive at (4.21), apart from two pre-factors. These can be disregarded as they satisfy the condition under complex conjugation explained in Remark 2.9 to establish an equivalent kernel. Thus the universality of the kernel (4.21) holds for arbitrary fixed τ , with 0 ≤ τ < 1.
After completing this work, the universality we have found at the origin at strong non-Hermiticity has been extended to the entire bulk (and edge) along the real line [45]. The universality of the elliptic Ginibre ensemble, a one-parameter family of Gaussian random matrices, strongly suggests a more general universality to hold, when allowing for a larger class of potentials Q in the weight function (4.1). Despite our progress in Section 3 in the construction of SOP for more general weight functions based on OP, this is beyond the scope of the current article.
The main result of this subsection is the following limit.
Theorem 4.7 (Bergman-like Laguerre kernel). Let 0 < τ < 1 and z, u ∈ C, then for N → ∞ we have (4.28) The proof of the above theorem will proceed in two steps. First, we treat the rotationally invariant case τ = 0, then we apply the integral representation of Laguerre polynomials to establish a lemma analogous to Lemma 4.4.
Then, as N → ∞ we have g N (u, v) ⇒ g (u, v), where the limiting function is given by Proof. Let u, v ∈ B(0, r), then the same upper bound as in Lemma 4.2 serves as a summable upper bound for (4.29), and it only depends r. Hence, by the Weierstraß M-test, its sum is an analytic function of u and v. The convergence is absolute and uniform in each compact subset of the plane. The limiting function -as N → ∞ -is given by the power series: Next we derive a differential equation for this limit. Following the same steps as in [5, App. B.1], where a differential equation was derived for g(u, v) − g(v, u), we can be brief. We obtain (4.32) In the last step we have taken out the l = k term of the inner sum, before shifting the index k → k + 1.
The first sum that we are left with is equal to ug (u, v). In the last sum we can simplify the denominator, by using the doubling formula for the gamma function and (2k − 1)!!(2k)!! = (2k)!, leading to .
In summary, g(u, v) satisfies the following second order linear inhomogeneous differential equation To solve this equation, we may again use u 0 = 0 as an initial value, because g(u 0 , v) = 0 for every v ∈ C (and we will see later that a second initial condition is not needed for the uniqueness of the solution). We will solve this initial value problem in three steps: First, we find two linearly independent solutions γ homA (u) and γ homB (u) for the homogeneous equation, then we construct a solution γ inhom (u) for the inhomogeneous equation, finally we set g(u, v) = aγ homA (u) + bγ homB (u) + γ inhom (u) where a, b ∈ C are determined by the initial condition.
(4.39) Using the series representations (4.34) of the Bessel functions J ν and I ν , one can see that γ homA (u) and γ inhom (u) are continuous functions in a neighbourhood of u = 0. Since g(u, v) is continuous around 0, with g(0, v) = 0, we can set b = 0. A simple calculation using (4.34) for γ homA (u) and (C.4) together with [39, 10.32.2] for γ inhom (u) gives and thus By replacing the Bessel-I functions with the integral representations [39, 10.22.52] and [39, 10.43.24], we can match this term to one of the double-integrals in γ inhom (u), then (4.42) In the last step we used Fubini's theorem, and then we switched the variable names q and p. Finally, the claimed formula for g(u, v) follows from (C.6).
This is the known limiting kernel found in [5,Appendix B] and is expected to be universal. As in Remark 4.3, τ = 0 considered here corresponds to maximal non-Hermiticity of the underlying ensemble.
Applying the substitution p = qr, r ∈ [0, 1], we can calculate the q-integral with (C.5). In the last step we bring the two resulting one-dimensional integrals into the form in (4.45). This is achieved by following changes of variables (for the first and second integral respectively) We can now complete the proof of our main result.
Proof of Theorem 4.7. Lemma 4.8 together with Lemma 4.10 shows that the sum in (4.44) converges absolutely and uniformly in each ball of center 0 and radius r, and therefore in each compact subset of C. Due to the absolute convergence of the series we can rearrange it and, in particular when ϑ = θ = τ , we obtain (4.51)

Universality of the symplectic chiral elliptic Ginibre kernel
In this subsection we will prove the universality of all k-point correlation functions (2.4) in the symplectic chiral elliptic Ginibre ensemble, in the large-N limit at strong non-Hermiticity close to the origin. The comments from Subsection 4.1.1 about scaling and N -dependence of the weight apply here as well. To derive the macroscopic density and the droplet for our weight (4.23), we exploit the Bessel asymptotic K ν (z) ∼ π/(2z)e −z to construct the limiting potential.
For the behaviour at the origin let us first quote the known result at maximal non-Hermiticity τ = 0 from [5]. In that case we can read off from Remark 4.9 the matrix elements of the limiting kernel of K N (z, u) from (2.5), times the normalisation from the area measure (4.53) We call this end result the limiting symplectic chiral Ginibre kernel at the origin, after removing the first factor due to Remark 2.9. We can now prove the following universality statement.  (4.53) in the sense of Remark 2.9 for general values 0 < τ < 1 and thus universal.
Proof. As explained already in the proof of Corollary 4.6 we have to unfold. In this case we have to rescale all arguments z → (1 − τ 2 )z. Therefore, we take the limit (4.54) The pre-factor (1 − τ 2 ) 3 originates from the rescaling of the arguments and the factors (z − z) in (2.4), times the normalisation from the area measure. Furthermore, we have multiplied with the τ -dependent factor from the area measure. After inserting (4.28) in the second line we arrive at (4.53), apart from the two pre-factors which lead to an equivalent kernel, cf. Remark 2.9. Thus the universality of the kernel (4.53) holds.
Once again we expect the universality found for the chiral elliptic Ginibre ensemble to hold for a more general class of weight functions, that share the same singularity of the weight (4.23) at the origin.

Christoffel perturbation for skew-orthogonal polynomials
In this section we will relate the SOP q n with respect to the weight function w(z) to those q (1) n skew-orthogonal with respect to the weight function w (1) (z) = |z − m| 2 w(z). For OP on subsets of the real line such a relation between OP with respect to weights w(x) and (x − m)w(x) (or in fact P (x)w(x) for a polynomial P (x)) is well known under the name of Christoffel perturbation, and determinantal formulas exist, compare [49]. Such a transformation, including multiplication of the measure by a rational function, is closely related to the Darboux transformation of integrable systems. In the complex plane we consider quadratic factors |z − m| 2 , in order to preserve the non-negativity of the resulting weight 4 . For planar OP (and also weighted Szegő polynomials) such a Christoffel perturbation has already been studied for w (M ) (z) = M l=1 |z − m l | 2 w(z) in [50], from which we borrow the notation. There, determinantal formulas similar to those in [49] have been derived for arbitrary M . For SOP, no such formulas were know. Only the polynomial kernel σ (M ) n of w (M ) was given in terms of the Pfaffian determinant of the polynomial kernel σ n and odd SOP q n of w, see [51]. We will use these expressions to provide an explicit representation of the perturbed SOPs q (1) n in the following theorem.
Theorem 5.1. Let (q n ) n∈N be the family of monic SOP with respect to the weight function w(z), with norms r n and pre-kernel σ n (z, u). Then, the following expressions hold for the monic SOP q (1) n (z), their norms r (1) n and kernel σ (1) n (z, u) with respect to the perturbed weight w (1) 2n (z), where d n ∈ R is an arbitrary constant. Furthermore, it holds Notice that for z = m, in eqs. (5.1) and (5.2) both numerator and denominator vanish, leading to a finite expression after applying the rule of l'Hôpital. For comparison we state here the corresponding result for OP from [50], where the same statement about z = m applies.
Remark 5.2 (Perturbed OP). Let us assume that µ has density function w on some domain ⊆ C, m ∈ C. Then it follows from [50], that one can express the sequence (p (1) n ) n of OP in L 2 (| · −m| 2 w(·)) in terms of the sequence (p n ) n of OP in L 2 (w(·)) as: where K n+1 (z, u) is the polynomial kernel constituted by the partial sum -up to n -of the orthonormal polynomials p k / √ h k . For the norms h (1) n and polynomial kernel K The proof of this remark can be found as a special case in [50,Section 3], where both p as a ratio of Pfaffians following [51], but for the proof of the above theorem we will only consider the simplest case M = 1.
We note that the simple relationship between the odd SOP and odd OP found in Theorem 3.1 breaks down for q (1) 2k+1 (z) and p (1) 2k+1 (z), even if the initial OP p n (z) were to satisfy a three-term recurrence relation. For a particular case of p (1) n (z) given in terms of Gegenbauer polynomials p n (z), this was indeed shown in [30,Section 5].
Let us present the proof now for Theorem 5.1.
Proof. We start with eq. (5.2). In [51,Eq. (2.14)] the density R N,1 (z), defined as in (2.3), was expressed in terms of a ratio of Pfaffian determinants for an arbitrary product of M characteristic polynomials, which for M = 1 reads: In the first line we started with eq. (2.6) that holds for arbitrary weight functions. For the second step we used the result in [51,Eq. (2.14)], to establish the claim in eq. (5.2). From this equation we will recover the skew-norms r (1) n and the even and odd SOP q (1) n , by taking appropriate limits.
First, we determine the skew-norms. From the definition (2.17), together with the fact that the SOP are monic q n (z) ∼ z n , with ∼ meaning that lim |z|→∞ q n (z)z −n = 1, we can read off the asymptotic for the pre-kernel with both arguments being large, |z|, |u| ≫ 1: Similarly, we obtain the asymptotic for a single argument being large, |z| ≫ 1: Next, we insert this into (5.2) to determine the leading order expansion of σ (1) n+1 (z, u) for both arguments being large: In the second step we have only kept the leading order, and in the last step we used that (5.6) also holds for σ (1) n+1 . Comparing the last two expressions, this leads to r (1) n as in the last equation of (5.1). Likewise, we can use (5.7) together with (5.2) to read off the even SOP: where we have inserted r (1) n as well, to express the right hand side in terms of unperturbed quantities only. This agrees with the first equation in (5.1), upon using the anti-symmetry of the pre-kernel.
For the odd polynomials q (1) 2n+1 (z) we have to go beyond the leading order and use that they are only determined up to an arbitrary constant times the even SOP of one degree less, q (1) 2n (z). For that purpose we label the next-to-leading order coefficient in the SOP as follows q n (u) = u n + k n u n−1 + O(u n−2 ). (5.10) Next, we expand the definition of the pre-kernel (2.17) to next-to-leading order for one argument being large, |u| ≫ 1: It follows that the odd polynomials can be obtained as and likewise for the perturbed pre-kernel and SOP. Inserting the known expressions for the skew-norm, pre-kernel and even SOP of the perturbed weight on the right hand side, we thus obtain Recognising that the last term in the last line is just the perturbed even SOP q (1) 2n (z), this yields the expression for the odd perturbed SOP in the second equation of (5.1), with the constant given by d n = k (1) 2n+1 − k 2n+2 − m here. Because this constant is arbitrary, we didn't specify it in (5.1).
In principle, both even and odd perturbed SOP q (1) n (z) could be expanded in the basis of the perturbed OP p (1) n (z). However, these are not as simple as in Theorem 5.1 and include the full sum of even and odd polynomials down to lowest order. We refer to Appendix B for details.

A Recollection of known planar OP and SOP
In this appendix we collect more planar OP and SOP. In part they are already known, but for completeness (and because we use some of them in the main text) we state them here in as much generality as possible. In particular, we can rederive the planar SOP from Theorem 3.1. In contrast to the main text, in this appendix we will consider the flat Lebesgue measure d 2 z = dx dy for z = x + iy.
Example A.1 (Product of Ginibre matrices). When taking the product of M complex Ginibre matrices we obtain for the weight function and norms [52] Here, G M,0 0,M is the Meijer G-function, see [39,Chapter 16.17] for the definition, and we have slightly extended [52] by the insertion of a point charge c > −1 at the origin. From Corollary 3.3 we have for the SOP and their skew-norms (also stated in [23]) c,c |z| 2 = 2|z| 2c K 0 (2|z|), given in terms of the modified Bessel-function of the second kind K ν . At c = 0 this corresponds to the weight of the chiral symplectic Ginibre ensemble at maximal non-Hermiticity, compare [5] where the corresponding SOP were constructed. The monic OP p n (z), the recurrence coefficients c n and the squared norms h n are given by

From Theorem 3.1 we get for the SOP
For the skew-norms we obtain (A.7) These polynomials reduce to the OP from [37] and to the SOP from [21] when choosing A = 1/(1 − τ 2 ) and B = τ /(1 − τ 2 ) for 0 ≤ τ < 1. The orthogonality relation (4.3) was proven first in [38] and independently in [15]. For self-consistency we present a proof for the well-known form of Hermite polynomials depending on two parameters (A.5) following [10].
Proof. Inserting the following integral representation [39, into the orthogonality relation, where the contour integral is around the origin in positive direction, we obtain Because the integrals are absolutely convergent, the order of integration can be interchanged, and using (C.1) we have performed the two real integrations. Cauchy's integral theorem for derivatives leads to a single integral that gives δ n,m . Making the Hermite polynomials H n (x) = 2 n x 2 + O(z n−1 ) monic, and using the known recurrence relation for Hermite we arrive at (A.5).

B Fourier coefficients of the perturbed SOP
In this appendix we compute the expansion of the perturbed SOP q (1) n from Theorem 5.1, in the basis of the perturbed OP p (1) n from Remark 5.2, which are skew-orthogonal respectively orthogonal with respect to the perturbed weight w (1) (z) = |z −m| 2 w(z). Furthermore, we assume that the unperturbed OP p n obey a three-term recurrence relation, and thus Theorem 3.1 applies to determine the q n . As a result we will see that both even and odd polynomials q (1) n have Fourier coefficients in even and odd degree of p (1) n , down to the lowest degree. We begin with the polynomials of odd degree, defining the coefficients β 2k+1,j as It follows from the fact that both perturbed SOP and OP are monic, that β 2k+1,2k+1 = 1. Following the definition we have for the remaining coefficients, with l < 2k + 1, that and zero otherwise. In the final result (B.2) the last term is non-vanishing only when l = 2L + 1 is odd, whereas the previous term is also present for l = 2L even (the sum runs to ⌊l/2⌋ = L in both cases). Consequently, all even and odd Fourier coefficients β 2k+1,l are non-vanishing in general, down to the lowest degree l = 0, in contrast to Theorem 3.1.
Let us move to the Fourier coefficients of the even polynomials, defined as As for the odd polynomials, we have from the monic property of the two sets of polynomials that α 2k,2k = 1. For the remaining coefficients with l < 2k we obtain where we have followed the same procedure as before, and swapped the summation in the double sum in the last step, to facilitate the integration. Let us distinguish even and odd indices l now. For even l = 2L, with L < k we obtain Once again all even and odd coefficients α 2k,l are non-vanishing, down to the lowest degree l = 0.

C Some useful integrals
For completeness we collect a few simple Gaussian integrals that will be useful in several places throughout the main part.
1. For all α > 0 and β ∈ C it holds: In Subsection 4.2 we encounter the following integrals involving Bessel-J and Bessel-I functions.