Scaling Limits of Planar Symplectic Ensembles

We consider various asymptotic scaling limits $N\to\infty$ for the $2N$ complex eigenvalues of non-Hermitian random matrices in the symmetry class of the symplectic Ginibre ensemble. These are known to be integrable, forming Pfaffian point processes, and we obtain limiting expressions for the corresponding kernel for different potentials. The first part is devoted to the symplectic Ginibre ensemble with the Gaussian potential. We obtain the asymptotic at the edge of the spectrum in the vicinity of the real line. The unifying form of the kernel allows us to make contact with the bulk scaling along the real line and with the edge scaling away from the real line, where we recover the known determinantal process of the complex Ginibre ensemble. Part two covers ensembles of Mittag-Leffler type with a singularity at the origin. For potentials $Q(\zeta)=|\zeta|^{2\lambda}-(2c/N)\log|\zeta|$, with $\lambda>0$ and $c>-1$, the limiting kernel obeys a linear differential equation of fractional order $1/\lambda$ at the origin. For integer $m=1/\lambda$ it can be solved in terms of Mittag-Leffler functions. In the last part, we derive Ward's equation for planar symplectic ensembles for a general class of potentials. It serves as a tool to investigate the Gaussian and singular Mittag-Leffler universality class. This allows us to determine the functional form of all possible limiting kernels (if they exist) that are translation invariant, up to their integration domain.


Introduction
In the pioneering work of Ginibre [31], it was first discovered that the complex eigenvalues of Gaussian random matrices with quaternion elements (also known as the symplectic Ginibre ensemble) behave like equally charged particles with complex conjugation symmetry. They interact via the two-dimensional Coulomb repulsion, subject to a confining Gaussian potential Q(ζ) = |ζ| 2 . Below we will consider more general potentials Q, where the joint probability distribution P N of the point process ζ = (ζ 1 , . . . , ζ N ) ∈ C N is given by dA(ζ j ). (1.1) Here the Hamiltonian H N is given by Compared to the eigenvalues of the complex Ginibre or more general random normal matrix ensembles, which correspond to a genuine two-dimensional Coulomb gas at specific inverse temperature β = 2 without further symmetries, one of the most distinguished features of the symplectic ensemble is the existence of a local repulsion from the real axis, which follows from the term − log ζ j −ζ j 2 in (1.2). To be more precise, this local repulsion originates from complex conjugate eigenvalue pairs which repel each other when approaching the real axis. For illustration, see Figure 1(a) for some random samplings of eigenvalues from the symplectic Ginibre ensemble.  In (a), the local repulsion from the real axis is clearly visible. In (b), the rescaled process at the right endpoint of the spectrum clearly displays the complex conjugation symmetry.
Why are symplectic ensembles interesting, apart from their statistical mechanics interpretation? Together with its real and complex counterparts, the symplectic Ginibre ensemble represents one of the few examples which is integrable, in the sense that it constitutes a Pfaffian point process where the matrix-valued kernel can be explicitly constructed. This fact will be recalled in more detail in the next section. It makes the asymptotic analysis of the kernel in various scaling limits possible, as will be the main topic of this work. The situation is much more difficult for Coulomb gases at general values of β = 2, see [9,24,47] and references therein for recent developments.
Moreover, it has been found rather recently that non-Hermitian random matrices enjoy a much wider class of universality than their Hermitian counterparts, in the sense that away from the real axis the limiting complex eigenvalue correlation functions of all three Ginibre ensembles agree [7,20,27,46]. The same phenomenon has been observed very recently when comparing the spectra of truncated unitary and symplectic random matrices [37,49], respectively.
While random matrices with complex eigenvalues have many applications in physics in general, we provide some examples where the symplectic ensemble yields unique predictions and differs from other non-Hermitian symmetry classes, notably at the origin. These include disordered non-Hermitian Hamiltonians with an imaginary magnetic field [40], thermal conduction in superconducting quantum dots [26] in the related circular quaternion ensemble, or the spectrum of the Dirac operator in quantum chromodynamics with two colours at non-vanishing chemical potential, for the symplectic ensemble with additional chiral symmetry [2]. In the latter the complex eigenvalues of the Dirac operator in the vicinity of the origin are of particular importance because of their role in chiral symmetry breaking. The limiting random matrix predictions have been confirmed from field theory in [4].
It is the goal of this article to derive further asymptotic results which are specific for the symplectic ensembles, notably at the edge of the real axis and in the presence of singularities at the origin. The representation for the limiting edge kernel that we will derive provides a unifying picture that allows us to relate results at different parts of the spectrum. Before we present our main results in the next section, let us summarise what was previously known about symplectic ensembles in various parts of the spectrum. Here we include the elliptic symplectic Ginibre ensemble with the potential Q(ζ) = 1 1−τ 2 |ζ| 2 − τ Re ζ 2 , where the parameter τ ∈ [0, 1) controls the degree of non-Hermiticity. Its joint probability distribution (1.1) and kernel at finite-N were derived by Kanzieper [35].
In the asymptotic analysis of the kernel at various points of the spectrum one has to distinguish local or microscopic from global or macroscopic scales. Since the above-mentioned repulsion from the real line affects only the microscopic scale when the complex conjugate eigenvalues become close, it could be expected that the leading form of the macroscopic eigenvalue density is the same as that of the two-dimensional Coulomb gas in the symmetry class of the complex Ginibre ensemble, with external potential Q/2, see, e.g., [29]. Indeed, for a general Q satisfying complex conjugation symmetry Q(ζ) = Q ζ and suitable potential theoretic assumptions, it was shown by Benaych-Georges and Chapon that as N → ∞ the empirical measure of ζ converges to Frostman's equilibrium measure associated with Q/2, see [17,Theorem 3.1]. In particular, the density of ζ tends to 1 2 ∆Q(ζ) · 1 S (ζ), where S is a certain compact set called the droplet. For the Ginibre ensemble this is the wellknown circular (or elliptic) law.
In the local scaling limit at the origin at maximal non-Hermiticity (τ = 0), the limiting kernel of the symplectic Ginibre was derived in [35,44]. At weak non-Hermiticity, when τ scales as 1 − τ = O N −1 , a different limiting kernel was found at the origin by Kanzieper [35]. It interpolates between the former at τ = 0 and the sine kernel of the Gaussian symplectic ensemble in the limit τ → 1. Both limiting kernels are invariant under translations along the real line, and we will come back to this feature below.
At the edge of the spectrum, it was shown by Rider [46] and spelled out by Dubach [27] that the maximal modulus fluctuations of complex and symplectic Ginibre ensemble agree. The agreement between the two ensembles was also studied in [7]. It was shown that in the bulk away from the real axis both ensembles yield the same determinantal point processes. In this work we will first focus on the edge on the real line. The local statistics along the real line for the elliptic symplectic Ginibre ensemble is found in [21].
Secondly, we investigate what happens when a non-Gaussian potential develops a singularity by inserting a point charge at the origin. In the generic case of a potential of Mittag-Leffler type we will be able to provide an explicit expression for the limiting origin kernel, that differs from the Ginibre universality class. Our findings can be thought of as the counterparts for previous results in random normal matrix ensembles [15,25]. (See also [16,18,42] for extensive studies on the orthogonal polynomials associated with Mittag-Leffler type potentials.) As the third issue, we study the universality for kernels that are translation-invariant along the real line. In setting up Ward's equation for the symplectic ensemble -an identity satisfied by limits of the rescaled one-point functions, under some natural assumptions, we can completely characterise the class of all such possible limiting kernels for general potentials by an integral representation. It is unique (if it exists) up to the integration domain, which is a connected interval symmetric around the origin, see [13,14] for analogous works in random normal matrix ensembles.

Main results
Let us now come to the main objects in this work. We denote by D(η, r) the disc with centre η ∈ C and radius r. For a given sequence of points p N , the positive number We drop the subscript and write p ≡ p N if the sequence does not depend on N . By (1.3), the micro-scale r N corresponds to the mean eigenvalue spacing in radial distance of the ensemble (1.1) at the point p ∈ S (cf. [21] for a situation where p is outside the droplet). We define the rescaled process z = {z j } N j=1 as follows: for all j, distinguishing the interior and the boundary of the droplet S. Here the angle θ ≡ θ N ∈ R is chosen so that e iθ is outer normal to ∂S at p, see Figure 1(b) inset for an illustration of the rescaled process. The k-point correlation function of the rescaled process z is defined by In (3.2) a more standard definition of the k-point correlation function before rescaling denoted by R N,k (ζ 1 , . . . , ζ k ) is given. Throughout the paper we distinguish these and all other objects (kernels, Hamiltonian, joint distribution) before rescaling by bold symbols. For the precise relation on the level of the kernels see (3.6). It results into the following relation between the k-point correlation functions It is well known that before and after rescaling the set z forms a Pfaffian point process (see [35,44]), i.e., R N,k is expressed in terms of a certain 2 × 2 matrix-valued kernel K N as where the kernel K N is of the form Here Pf denotes a Pfaffian of the 2k × 2k skew-symmetric matrix and where θ is given as in (2.2). The arguments (2.4) are given according to the rescaling (2.2). For the spectral density at k = 1, let us denote R N ≡ R N,1 . The primary goal of this work is to derive the large-N limit of the kernel K N for various potentials Q. For the Gaussian potential Q(ζ) = |ζ| 2 , where the associated ensemble (1.1) corresponds to the symplectic Ginibre ensemble, it follows from the circular law that S = D 0, √ 2 . In [35], Kanzieper studied the elliptic potential Q(ζ) = 1 1−τ 2 |ζ| 2 − Re ζ 2 and derived the scaling limit for the pre-kernel at the origin p = 0 in the almost-Hermitian regime when 1 − τ = O 1 N . Also at p = 0 and for maximally non-Hermiticity at τ = 0, he showed that the associated ∞-point process {z j } ∞ j=1 of the symplectic Ginibre ensemble has the correlation kernel where the pre-kernel κ R bulk is given by We also refer to [44] for an alternative derivation of (2.6). Furthermore, it was shown in [21,38] that the kernel (2.5) also appears when p ∈ int(S)∩R in the symplectic elliptic Ginibre ensemble and in that sense is universal (in [7, Appendix B] a different strategy at τ = 0 was mentioned). Our first main result Theorem 2.1 below provides the boundary scaling limit when p ∈ ∂S ∩ R = ± √ 2 , see Figure 2 for the graphs of R N , respectively.
and write W (f, g) for the Wronskian of two functions f , g to formulate our first main result.
Theorem 2.1. Let Q(ζ) = |ζ| 2 and p = ± √ 2. Then R N,k (z 1 , . . . , z k ) converges uniformly for z 1 , . . . , z k in compact subsets of C to We remark that the pre-kernel κ R edge has the following alternative representation As a consequence of Theorem 2.1, we obtain the following corollary.
. Moreover, as t → ∞, we have where the o(1)-term is uniform on z 1 , . . . , z k in compact subsets of C and K C edge (z, w) := e −|z| 2 −|w| 2 +2zw 1 2 erfc(z +w). (2.11) We emphasise that the kernel (2.11) corresponds to that of the complex Ginibre ensemble at the edge of the spectrum (up to a trivial rescaling), which is known to form a determinantal point process, see, e.g., [20,30,36]. Therefore the convergence (2.10) implies that away from the real axis, the local edge statistics of the symplectic and complex Ginibre ensemble are equivalent in the large-N limit. We also refer the reader to [27,46] and [7] for the equivalence of the symplectic and complex Ginibre ensemble in the context of the scaled maximal modulus and the local bulk statistics away from the real axis, respectively.
Before moving on to our next topic, let us give some further remarks on Theorem 2.1. These complete the asymptotic study of the symplectic Ginibre ensemble with the Gaussian potential. (i) Recently, the scaling limit of the symplectic Ginibre ensemble at the edge of the spectrum has been obtained independently by Khoruzhenko and Lysychkin [38], see also [43] for more details. Contrary to our approach using a differential equation satisfied by the prekernel, their methods are based on contour integral representations, which yields the same result.
(ii) For p ∈ − √ 2, √ 2 , the pre-kernel (2.6) can be obtained from the same method used in the proof of Theorem 2.1. Indeed, the pre-kernel (2.6) can be represented by (2.8) with E = (−∞, ∞) using the well-known integral The same integral identity holds for the complementary error function, replacing erf → erfc.
(iii) Theorem 2.1 can be generalised to a moving boundary point p N . More precisely, if we set and rescale the process as (2.2), then the limiting correlation kernel is of the form (2.8) with E = (−∞, a). When a 1 we recover the bulk limit due to (ii).
(iv) It is well known that the local bulk/edge correlation kernel K C of the random normal matrix ensemble is given by where E = (−∞, ∞) for the bulk case and E = (−∞, 0) for the edge case, see [11,33]. Thus one can interpret the integral representation (2.8) as a symplectic analogue of the above expression.
(v) In [8], the authors studied the ensemble (1.1) with the elliptic potential For this model, the local correlation kernel at the right/left endpoint of the spectrum was derived in the regime of weak non-Hermiticity when 1 − τ = O N −1/3 . The kernel (2.8) was derived there as well from the large argument limit of such an intermediate process, see [8, equation (4.28)].
(vi) It was proved in [5] that the limiting local kernel of the chiral complex Ginibre ensemble at multi-criticality is given by the edge kernel of the complex Ginibre ensembles in squared variables. A quaternion version of such a result will appear in future work.
Next, we investigate ensembles containing certain types of singularities at the origin. More precisely, we consider the Mittag-Leffler ensemble whose associated potential Q is of the form Here the condition c > −1 is required to guarantee Z N < ∞, where Z N is the partition function in (1.1). Note that when c = 0, a "conical singularity" at the origin arises from an insertion of a point charge c. On the other hand, since the limiting global density of the ensemble is ∆ 1 2 Q(ζ) = λ 2 2 |ζ| 2λ−2 , when λ > 1 (resp., λ < 1) it vanishes (resp., diverges) at the origin. We aim to discover the local statistics of the symplectic Mittag-Leffler ensemble, which provides "non-standard" or multi-critical universality classes (2.14) beyond (2.6). Away from the origin we expect to be back in the Gaussian universality class of the Ginibre ensemble. Note that the micro-scale at the origin is given here by For each λ > 0 and c > −1, the rescaled limiting local kernel K λ,c at p = 0 is of the form where κ λ,c is a holomorphic function in both variables (with branch cuts). Since the correlation functions do not depend on the specific choice of these branches, we will ignore this issue in the sequel, such a monodromy of the pre-kernel can also be interpreted as a cocycle. (Cf. Section 3.2 for the definition of a cocycle.) In general, we derive a certain (1/λ)-order fractional differential equation for κ λ,c (Proposition 4.1), which can be recognised as a version of the Christoffel-Darboux formula for the kernel (2.14). Moreover, for λ = 1/m (m ∈ N), we obtain an explicit formula for κ λ,c by solving the associated differential equation of order m.
To describe the local kernel of the Mittag-Leffer ensemble, let us first recall the definitions of the Mittag-Leffler functions. By definition, the two-parametric Mittag-Leffler function E a,b is given by , (2.15) and the three-parametric (Kilbas-Saigo) Mittag-Leffler function E α,m,l is given by and W j,m (z) := W (g 1,m , . . . , g j−1,m , g j+1,m , . . . , g m,m ), (2.19) where W is the Wronskian.
where for m = 1, and for m > 1, Here G(m + 1) := m−1 j=1 j! is the Barnes G-function. Remark 2.5. For any m ≥ 1, one can write κ λ,c in a unified way as Theorem 2.4 is our second main result, and we shall present some examples for m = 1, 2.
Example 2.6 (λ = 1). We first discuss the case λ = 1, where the Mittag-Leffler ensemble corresponds to the eigenvalue statistics of the so-called induced symplectic Ginibre ensemble. This name was proposed for the matrix representation of real and complex Ginibre ensembles in the presence of zero eigenvalues [28]. We refer to [2] for the complex eigenvalue correlation functions in the symplectic Ginibre ensemble in the presence of zero eigenvalues in terms of skew-orthogonal polynomials at finite-N . For m = 1, we have g 1,1 (z) = e z 2 /2 . This immediately gives Note that by (2.15), we have is the (regularised) incomplete Gamma function. Using this, we have an alternative representation for c ≥ 0, Furthermore, for a non-negative integer c one can express the pre-kernel (2.20) in terms of error functions. For this we denote by (a) n Pochhammer's symbol: (a) 0 = 1, (a) n = a(a + 1) · · · (a + n − 1).
Then for an even integer c = 2n (n ∈ N) we have and for an odd integer c = 2n − 1 (n ∈ N), we have Here we use the convention that the summation with an empty index equals zero. Both (2.21) and (2.22) can be obtained by straightforward computations using see, e.g., [45, equation (8.4.9)].
In [3], the k-point correlation function R N,k for an integer-valued point charge c is presented in a different way as the ratio of Pfaffians of the correlation kernels of c = 0.  2 ). For m = 2, we have where I ν is the modified Bessel function of the first kind [45,Chapter 10]: .
Then using (2.23), we have In particular, for c = 0, − 1 2 , the pre-kernel (2.24) can also be expressed as (2.25) The expression (2.25) follows from [2, Appendix B] due to the relation between the Mittag-Leffler potential at m = 2, Q(ζ) = |ζ| − 2c N log |ζ| for c = 0, − 1 2 , and the potential of the chiral symplectic Ginibre ensemble at maximal non-Hermiticity τ = 0, Q(ζ) = − 1 N log |ζ| 2ν K 2ν (N |ζ|) (corresponding to µ = 1 and the change of variables ζ = z 2 therein). Namely, at ν = ± 1 4 the modified Bessel-function of the second kind simplifies, K ± 1 2 (x) = π 2x e −x , matching the two cases for c = 0, − 1 2 up to an additive constant. In [2] the relation between the chiral symplectic Ginibre ensemble in the origin scaling limit and quantum chromodynamics with two colours and chemical potential was pointed out, cf. [4]. We also refer to the recent work [6] where the pre-kernel at the origin was determined in the more general case of the chiral elliptic potential, with τ ∈ [0, 1).
Notice, however, that the equivalence between (2.24) and (2.25) is far from being obvious. One can verify this by showing that both of these expressions satisfy the same differential equation of second order (Proposition 4.1) with the same initial conditions, which uniquely determine the solution. In particular, the initial conditions can be easily checked using the integral representation [45, equation (11.5.6)] of the modified Struve function L ν .
In the third part, let us turn to the symplectic ensemble with general external potential Q. We present two important functional equations satisfied by the correlation kernel, the mass-one and Ward's equation.
First, we define the Berezin kernel is the 1-point function of the N -point process, conditioned to contain the prescribed point z. See Figure 5 for graphs of the Berezin kernel.
By the definition of R N,k in (3.2), one can easily see that the mass-one equation holds for finite-N . For a Pfaffian ∞-point process, let us define the associated Berezin kernel as where R k := lim N →∞ R N,k and R = R 1 . Here and in the sequel, a point process is called Pfaffian ∞-point process if its correlation functions R k are expressed in terms of the Pfaffian of a certain satisfies the mass-one equation. Here a ∈ R and f z is given in (2.7).
For a correlation kernel K with (associated) pre-kernel κ of the form the mass-one equation is equivalent to  In the opposite direction, for a general potential Q, suppose that the associated limiting Pfaffian point process z satisfies the mass-one equation (2.31). From an analytic point of view, this equation can be regarded as an integral equation satisfied by a pre-kernel κ. If one can characterise the solution to (2.31) under some conditions derived from intrinsic properties of Q, and the nature of the rescaling point p (e.g., whether p ∈ int(S) ∩ R or p ∈ ∂S ∩ R), this provides an alternative way to obtain all possible candidates for the scaling limit of the symplectic ensemble associated with potential Q. This would show local universality as this approach does not depend on the specific choice of Q. (We refer to [13,14] for extensive studies on the universality for random normal matrix ensembles based on this approach.) However, the massone equation (5.6) may have further solutions beyond (2.6) and (2.8). Therefore, one needs additional information about the kernel to determine the solution uniquely. This calls for an investigation of Ward's equation.
Remark 2.9 (Ward's equation and universality for the random normal matrix ensemble with translation invariant kernel). We pause here to briefly introduce recent developments in Ward's equation and its significance in the context of universality for the random normal matrix ensemble. Hopefully this gives more intuition as to why we aim to examine Ward's equation for symplectic ensembles as well.
For random normal matrix ensembles, the rescaled version of Ward's equation was studied in [13]. An important feature of this equation is that in the large-N limit, it does not depend on the choice of potential Q if p is regular, in that sense that 0 < ∆Q(p) < ∞ is non-vanishing and bounded. Due to this property, Ward's equation has been utilised to show the local universality conjectures in various situations, see, e.g., [10,14,15] and references therein. To be more precise, the overall strategy for the universality proof using Ward's equation is as follows.
• Derivation of Ward's equation (cf. [ The first step is to show that for a C 2 -smooth potential Q and a regular point p, the following form of Ward's equation holds: where R ≡ R 1 and R 2 are the (limiting) rescaled 1-and 2-point functions of the random normal matrix ensemble with potential Q, i.e., .
Here R N,k denotes the unscaled k-point functions. We emphasise that due to the rescaling factor N ∆Q(p) chosen according to the mean eigenvalue spacing of the random normal matrix model, the equation (2.32) does not depend on the choice of Q.
• Structure of the correlation kernel (cf. [13, Theorem 1.1], [15,Theorem 1.3]). The next step is to show that for a general potential Q and a regular point p, the limiting correlation kernel K exists and is of the form where Ψ is Hermitian, i.e., Ψ(z, w) = Ψ(w, z) and Ψ is entire as a function of z andw.
• Characterisation of translation-invariant solution to Ward's equation (cf. [14,Theorem 4]). The third step is to show that the only non-trivial horizontally translation-invariant solutions R (i.e., R(z + t) = R(z) for t ∈ R along the real axis) to Ward's equation (2.32) are given by where I is a connected interval. Due to this step, the possible translation-invariant scaling limit is determined only by one interval I.  , a) for some a > 0 if one considers the almost-Hermitian limit [10] or p is close to a cusp type singularity [14]. Determining the interval in each situation requires a separate analysis.
• Translation invariance of the correlation functions. The translation invariance is easy to check if Q is radially symmetric. As a consequence, for a radially symmetric potential, edge universality was shown in [13, Theorem 1.8].
We remark that for a more general class of non-radially symmetric potentials, edge universality was recently shown by Hedenmalm and Wennman in [33]. The authors used a different approach based on the theory of quasi-orthogonal polynomials. However, this theory is not directly applicable to the symplectic ensemble since it is far from obvious how to construct skew-orthogonal polynomials using (quasi-)orthogonal polynomials in general (see however [6]).
In contrast, the overall strategy using Ward's equation described above can be applied in parallel to the symplectic ensemble. This is the primary purpose of our remaining discussion. In the rest of this section, as an analogue of [14], we aim to particularly address the characterisation of the translation-invariant solution to Ward's equation for the symplectic ensemble.
Our next result is the resulting Ward's equation for symplectic ensembles. For this, let Then, we obtain the following form of Ward's equation at finite-N for general Q.
Proposition 2.10. Suppose that Q is C 2 -smooth and p ∈ R. Then for each N, we havē (2.33) Let us discuss the large-N limit of Ward's equation. Since the micro-scale r N is given by (2.1), if p is regular, we have (1)). Therefore if we formally take the large-N limit of Ward's equation (2.33) and if R is non-trivial, for a general potential Q and p regular, we arrive at One may compare Ward's equation (2.35) for the symplectic ensemble with that for the random normal matrix ensemble (2.32).
Remark 2.11 (limiting Ward's equation at a singular point of Mittag-Leffler type). As an analogue of [15], we briefly discuss the limiting form of Ward's equation at a singular point of Mittag-Leffler type. For λ > 0, let Q be a potential satisfying Q(ζ) = |ζ| 2λ + o |ζ| 2λ as ζ → 0. Then we consider a potential Q of the form By (2.13), in the sense of distributions, we have where δ 0 is the Dirac delta at the origin. Taking the limit N → ∞ of the equation (2.33), at least formally we arrive at the distributional Ward's equation For the symplectic ensemble with a general potential Q and for a regular point p, suppose that the limiting correlation kernel exists and is of the form (2.30). To obtain the Gaussian factor e −|z| 2 −|w| 2 in (2.30), notice the Taylor series expansion: Then it follows from (2.34) and (3.4) that the second line on the right-hand side contributes to the Gaussian factor, whereas the first line contributes to the pre-kernel. We refer to [13, Section 3.5] for a similar computation.
In the spirit of [14], we aim to characterise (horizontally) translation-invariant kernels, that is those invariant under translations along the real axis. Moving away from the real axis will change the universality class.
It is not difficult to see, that a pre-kernel κ of the form .
Therefore it is clear that the horizontal translation invariance of the k-point function R k holds along the real axis, i.e., for t ∈ R, R k (z 1 + t, . . . , z k + t) = R k (z 1 , . . . , z k ).
We write Ψ as Fourier's inversion form for some odd function J. Here and in the sequel, let us assume that J ∈ L 1 ∩ L 2 . Ward's equation of the form (2.35) will be used to show our final result valid for general potentials. To be more precise, we have discussed that for a general C 2 -smooth potential Q and a regular point p in the real bulk, the limiting point process should (at least intuitively) satisfy the mass-one equation (2.27) and Ward's equation (2.35). The steps that lack in the proof are existence of the limiting correlation kernel, and the one of taking the limit N → ∞ of (2.26) and (2.33). These would require a separate analysis, see [13] for the random normal matrix ensemble. As these are beyond the scope of this paper, we so far are able to characterise the translation-invariant solution of the limiting mass-one equation ( where E = (−a, a) for some a > 0.
The overall proof of this theorem is parallel to that of [14,Theorem 4].

Remark 2.13.
(i) Note that if a = ∞, the kernel (2.39) corresponds to the symplectic Ginibre kernel in the bulk along the real line (2.6). On the other hand, if a < ∞ is fixed, then it corresponds to the kernel in the almost-Hermitian limit at the origin [35] after an appropriate rescaling of the eigenvalues. (This result has been extended to the entire bulk along the real line in [22].) In particular, note that for the density we have see Figure 6. For random normal matrix ensembles in the almost-Hermitian regime, a way to characterise the precise interval E was presented in a recent work [10].
(ii) Similar to random normal matrix ensembles discussed in [13,14], the translation-invariant scaling limit (5.4) enjoys a Gaussian convolution structure. To be precise, let us write γ(z) := 1 √ π e −z 2 for the Gaussian kernel and define where ϕ is a suitable tempered distribution on R. Then we have Ψ(z) = γ * ϕ(z), where ϕ satisfies

Scaling limits of the Ginibre ensemble
In the first part of this section, we review the canonical representation of the correlation kernel in terms of skew-orthogonal polynomials. In the second part, we prove Theorem 2.1 and Corollary 2.2.

Skew-orthogonal polynomials
For a given potential Q, the anti-symmetric scalar product ·|· is given by By definition, a family {q k } ∞ k=0 of polynomials is called skew-orthogonal polynomials associated with Q if it satisfies where δ kl is the Kronecker delta and s k is a positive number that depends only on k.
For a radially symmetric potential Q, let Then by [6, Theorem 3.1] (see also [34, p. 7]) form a family of skew-orthogonal polynomials associated with the potential Q. Here we have s k = 2h 2k+1 .
Recall that the particle system (1.1) forms a Pfaffian point process, namely, the k-point correlation function is expressed as where the 2 × 2 matrix-valued kernel K N is of the form In particular, we write R N ≡ R N,1 for the 1-point function. It is well known that the prekernel κ κ κ N takes the form see [35], where the quaternionic determinant Q det is used instead of the Pfaffian Pf. We also refer to [44] for the relation between Q det and Pf.
The relation between the pre-kernels κ κ κ N and κ N for the change of variables (2.4) at point p ∈ S is given as Thus the rescaled process z retains its Pfaffian structure in terms of the rescaled matrix-valued kernel cf. (2.4) for the unscaled variables ζ and η inside the arguments of the potential Q. We remark that in the factor r 3 N in (3.6), one factor r N comes from the term k j=1 whereas the other two factors originate from the change of variables for the pre-kernel. This leads precisely to (2.3) as follows

Boundary Ginibre point processes
For the Gaussian case Q(ζ) = |ζ| 2 , by (3.1), we have Therefore it follows from (3.5) that the associated kernel κ κ κ N has an expression This recovers the expression obtained in Mehta's book [44] in a different approach, compared to the skew-orthogonal polynomials we use here from [35]. By definition, a function c(z, w) is called a cocycle if there exists a (continuous) unimodular function h (i.e., h(z) = h(1/z)) such that c(z, w) = h(z)h(w). Since cocycles cancel out when multiplying the pre-kernel and taking the Pfaffian, if (c N ) ∞ N =1 is a sequence of cocycles, then c N κ N is an equivalent realisation of the pre-kernel κ N for the same point process z. Let us denote by c N · K N the correlation kernel associated with c N κ N , i.e.,

Now we prove Theorem 2.1.
Proof of Theorem 2.1. We prove the theorem for p = √ 2 only. The other case p = − √ 2 can be proved in the same way. Note that for Q = |ζ| 2 , we have r N = 2/N , see (2.1).
By (3.7), we have Here it follows from (3.6) and (3.8) that Let us write Then we have where the cocycle c N is given by c N (z, w) := e 2i √ N Im(z+w) . We claim that κ := lim N →∞ κ N satisfies the differential equation Here the uniform convergence of κ N on compact subsets of C follows from the Weierstrass M -test, see [6,Section 4] for similar computations. Differentiating (3.9), we have Rearranging the terms, we rewrite this equation as Similarly, Combining all of the above equations with (3.10), we obtain We now derive (3.11) from (3.12). By combining Stirling's formula with the elementary asymptotic e −2 where the o(1)-terms are uniform on compact subsets of C.
Let us define Then integration by parts gives that which leads to Therefore κ(z, w) := e (z−w) 2 f (z, w) is an anti-symmetric solution to (3.11). Since (3.11) is a first-order differential equation, the anti-symmetry plays the role of the initial value, which determines the solution uniquely. This completes the proof.
We now prove Corollary 2.2.
Proof of Corollary 2.2. By (2.2), for p N = e iθ N with θ N = t √ N , we have Therefore the first assertion (2.9) follows along the same lines of the proof of Theorem 2.1 by replacing (z, w) → (z + it, w + it). We now prove (2.10). By (2.9), we have (3.15) as t → ∞. Here f (x) ∼ g(x) means that f (x)/g(x) tends to 1 as x → ∞. Using the following asymptotic of the error function (see, e.g., [45, equation (7.12.1)]) we have that as t → ∞, Therefore we obtain that as t → ∞, Here c(z, w) = e −(z+z)it+(w+w)it is a cocycle for the determinantal point process. On the other hand, it follows from similar computations that as t → ∞, as t → ∞. The remaining argument is similar to that used in [7]. By a proper reordering, we have

Now it follows from the identity
and (3.17) that This completes the proof.

Christoffel-Darboux type formula
For the Mittag-Leffler potentials (2.12), the orthogonal norm h k is given by Then by (3.1), we have and To describe the Christoffel-Darboux type formula, let us recall that the Caputo fractional derivative D ν z is given by see [39,Section 2.4]. Here f (n) is the usual n:th derivative. Recall that p = 0 in the rescaling (2.4) here, and K λ,c is given by (2.14).
uniformly for z, w in compact subsets of C, where Here κ λ,c is given in terms of a function G as κ λ,c (z, w) = G(z, w) − G(w, z) and the following fractional differential equations hold: We call the equation (4.4) (generalised) Christoffel-Darboux formula. Such a differential equation satisfied by the correlation kernel, that makes the asymptotic analysis possible, is broadly called Christoffel-Darboux type formula, see, e.g., [19,23,41].
We also remark that the inhomogeneous term (zw) corresponds to the kernel of the complex Mittag-Leffler ensemble, see [15]. Such a relation has been observed for other models as well, which include the elliptic Ginibre ensemble [21] and its chiral counter part [2]. See also [1,48] for a similar relation for Hermitian ensembles, which gives an expression of the kernel of the symplectic ensemble in terms of a small number of orthogonal polynomials.
Proof of Proposition 4.1. Let us define Here and henceforth we use the convention that if l = 0, then l j=1 ϕ = 1. By (4.1), (4.2) and (2.13), the kernel K N is given by By definition of the Caputo derivative, we have that if k > ν compact subsets of C. The existence of the limit follows from the Weierstrass M -test. By direct computations using (4.5), we have Note in particular that as z → 0, Applying the operator D 1 λ z to the identity (4.6), we obtain

4.2
Bulk singularities of the order λ = 1 m with m ∈ N A point p in which the eigenvalue density 1 2 ∆Q(p) vanishes or diverges is called bulk singularity, see [15] and references therein. In this subsection, we shall consider the case λ = 1 m , where m is a positive integer and prove Theorem 2.4 for general c > −1.
Proof of Theorem 2.4. By (4.3) and (4.7), the function G is a unique solution to the m-th order linear (ordinary) differential equation which satisfies the asymptotic behaviour Here the function F m,c is given by (2.17). Note that the m initial conditions are provided by differentiating (4.11). We aim to solve the initial value problem (4.10) and (4.11). It is well known that the functions g j,m given by (2.18) provide the solutions to the homogeneous equation of (4.10), see, e.g., [32,Theorem 5.18]. By (2.16), we have that Using this, it is easy to observe that where G is the Barnes G-function [45,Section 5.17]. Due to Abel's identity, for m > 1, this leads to On the other hand for m = 1, we have obviously W (g 1 )(z) = g 1 (z). For each j ∈ {1, . . . , m}, let us write Note here that by (2.19), we have W j,m (z, w) = (−1) m−j F m,c (zw)W j,m (z). (4.14) From now on, we shall only consider the case m > 1. The other case m = 1 follows along the same lines. In this case, one can also easily solve the associated linear non-homogenous ODE of degree 1 employing an integrating factor.
By virtue of the method of variation of parameters and (4.13), one can write for some functions C j,m . We shall show that in the above expression G satisfies the asymptotic behaviour (4.11) if and only if C j,m ≡ 0 for every j.
By (4.12) and (4.13), direct computations of the determinant give that as z → 0, On the other hand, by (2.15), as z → 0, Combining these two equations, we have that as z → 0, This leads to Next, we show the identity By the binomial theorem, we have Integrating this identity, we have where B is the beta function. This gives (4.15). Using (4.15), we have Thus by (4.11), we obtain that C j,m ≡ 0 for all j ∈ {1, . . . , m}. Therefore by (4.14), we conclude that G(z, w) = 1 G(m + 1)

Mass-one equation
In this subsection, we consider the Pfaffian ∞-point process with the correlation kernel (2.28).
By (2.3), we have and Therefore the Berezin kernel is written as We now prove Proposition 2.8.
Proof of Proposition 2.8. By (5.1) and conjugation symmetry, the mass-one equation (2.27) is equivalent to Note that the pre-kernel κ ≡ κ R a given in (2.29) is written as where E ≡ E a = (−∞, a). Using this expression, we rewrite the right-hand side of the equation (5.2) as We now compute Due to the translation invariance in x, it suffices to consider the case v = 0. Then we have g(u, 0) = 2 π R 2 ye −4y 2 e −2(x+iy−u) 2 e −2(x−iy) 2 dy dx = 2 π R e −4x 2 +4xu−2u 2 R ye 4iuy dy dx.
Note that in the sense of distributions, we have

Thus it follows from
This leads to where Θ denotes the Heaviside theta function. Combining the above equations, interchanging the integrals, and using integration by parts, we obtain which completes the proof.

Ward's equation at finite-N
In this subsection, we consider the symplectic ensemble (1.1) with general external potential Q and prove Proposition 2.10. In the sequel, we assume that θ = 0 without loss of generality.
Proof of Proposition 2.10. We denote by E N the expectation with respect to the Gibbs measure (1.1). Let ψ N be a fixed test function. Integration by parts gives that, for each j, where H N is the Hamiltonian given by (1.2). Summing in j, we obtain Thus we have Due to the conjugation symmetry Q(ζ) = Q ζ , one can easily observe that for any j, Then it follows from the definition (3.2) that R N,2 (ζ, η) = R N,2 (ζ,η). Using the change of variable η →η, we have which leads to Also it follows from the definition and integration by parts that To describe the rescaled version of Ward's equation, let ψ(z) := ψ N (ζ), where ζ = p + r N z. By the relation R N,k (z 1 , . . . , z k ) = r 2k N R N,k (ζ 1 , . . . , ζ k ), ζ j = p + r N z j and the change of variables ζ → p + r N z, η → p + r N w, we have Similarly, Therefore by (5.3), we obtain Dividing by R N (z), this identity is rewritten in terms of Berezin kernel B N as By taking∂ derivative on both sides of this equation, Proposition 2.10 follows.
In the next two subsections, we characterise translation-invariant scaling limits of symplectic ensembles, first, by solving the mass-one equation, and second, by the solution of the limiting Ward's equation. This leads to the proof of Theorem 2.12. By (2.37) together with (2.38), it suffices to show that where E is of the form (−a, a) with a ∈ [0, ∞].

Translation-invariant solutions to the mass-one equation
In this subsection, we prove the following proposition. It follows from this proposition that the specified form of the pre-kernel (2.39) with any symmetric Borel set satisfies the mass-one equation (2.27). However, in the following subsection, it will be shown that the if E is not a single interval, then the pre-kernel (2.39) does not satisfy Ward's equation. (An analogous result for the random normal matrix ensemble was obtained in [14].) Proof of Proposition 5.1. In terms of Ψ, the mass-one equation (5.2) is written as By translation invariance, it is equivalent to the fact that for any t ∈ R, Ψ(2it) = 4 πi R 2 y e −4y 2 |Ψ(x + i(y − t))| 2 dx dy. (5.5) In terms of J, the equation (5.5) is written as the above equation is simplified to π R e −2ut u e u 2 /4 J(u) 2 du.
Therefore we can readily see that the kernel (2.39) satisfies the mass-one equation.
Since J ∈ L 1 ∩ L 2 , we can extend the domain of f , g to the complex plane. Then by (5.6) we have f (z) = g(z) for z ∈ C, which in particular leads to f (it) = g(it) for t ∈ R, i.e., By Fourier inversion, this is equivalent to It gives rise to a.e.
for some Borel set E, which is symmetric with respect to the origin since J is odd. Here the argument u/2 is merely used for convenience.

Translation-invariant solution of Ward's equation
In this subsection, we characterise the translation-invariant solutions to Ward's equation (2.35) and complete the proof of Theorem 2.12.
Recall that by (2.37), we have Let us first observe the following.
Proof . Note that Notice that since Therefore by (5.10), we obtain d dt Thus Ward's equation is equivalent to d dt We are now ready to prove Theorem 2.12. It results from this part of the proof that E = (−a, a) for some a > 0. On the other hand, Combining the above equations, we obtain that Also since for some constants c 1 , c 2 . This gives that E is connected (up to a null set). Now the proof is complete.