# Winding Number Statistics for Chiral Random Matrices: Averaging Ratios of
Parametric Determinants in the Orthogonal Case

Nico Hahn1∗, Mario Kieburg2, Omri Gat3 and Thomas Guhr1 1 Fakultät für Physik,
Universität Duisburg–Essen, Duisburg, Germany
2 School of Mathematics and Statistics, The University of Melbourne,
Melbourne, Australia
3 Racah Institute of Physics, The Hebrew University, Jerusalem, Israel
<EMAIL_ADDRESS>
###### Abstract

We extend our recent study of winding number density statistics in Gaussian
random matrix ensembles of the chiral unitary (AIII) and chiral symplectic
(CII) classes. Here, we consider the chiral orthogonal (BDI) case which is the
mathematically most demanding one. The key observation is that we can map the
topological problem on a spectral one, rendering the toolbox of random matrix
theory applicable. In particular, we employ a technique that exploits
supersymmetry structures without reformulating the problem in superspace.

††preprint: AIP/123-QED

##  I Introduction

The topological classification of spectrally gapped Hamiltonians plays a major
role in contemporary condensed matter physics. This is due to the possible
emergence of localized states at the boundary of the system, which arrive with
a certain robustness against perturbations [1, 2, 3, 4, 5]. The number of
these edge modes is determined by a topological invariant [6, 7, 8], that is a
property of the eigenstates of the Bloch Hamiltonian, which parametrically
depends on the (quasi-)momenta in the Brillouin zone. In one-dimensional
chiral symmetric systems this topological invariant is the winding number.

In disordered systems the winding number can become random, which motivates a
statistical analysis. We do this by setting up a parametric random matrix
model for the Hamiltonian for which we then calculate the winding number
statistics. Random matrix theory is known to successfully model universal
properties in systems with a sufficient degree of complexity in the limit of
large matrix dimensions [9, 10]. Unlike earlier results on parametric
Hamiltonians using similar random matrix models [11, 12] we are not
investigating spectral properties, but topological ones.

When additionally considering time reversal invariance next to chiral
symmetry, one can distinguish between three different symmetry classes: the
chiral unitary, chiral orthogonal and chiral symplectic class (AIII, BDI and
CII respectively in the tenfold way classification [13, 14, 15, 16]). In
previous works, we considered our random matrix model in the classes AIII and
CII [17, 18]. Here, we want to catch up on the chiral orthogonal class BDI.

Our goal is to calculate averages over ratios of determinants with parametric
dependence. These will serve as generators for the winding number statistics.
The problem is related to averaging ratios of characteristic polynomials,
which is a common task in random matrix theory [19, 20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30], because they entail information about the spectral
statistics. Correspondingly, powerful methods, like the supersymmetry method
[31], have been developed to treat such ensemble averages. In this way we are
mapping a topological problem to a spectral one.

In this study, just as we did in Ref. 18, we employ a technique introduced by
two of the present authors [21, 22]. It establishes exact expressions for the
ensemble averages while also revealing supersymmetric structures without
actually mapping to superspace. For this reason it has been coined
supersymmetry without supersymmetry. The resulting expressions are
determinants in the unitary case and Pfaffians in the orthogonal and
symplectic case containing simplified averages of only two determinants. The
remaining task is to calculate these simplified averages over the
corresponding random matrix ensemble, which is the step we show in this paper.

The outline of the paper is as follows: in Sec. II we set up the random matrix
model and define the problem we want to solve. In Sec. III we present our
results, while Sec. IV contains the derivations. We relegated some details to
the appendix. In Sec. V we conclude our study.

##  II Posing the Problem

Formally, chiral symmetry of a Hamiltonian $H$ can be expressed as

$\\{\mathcal{C},H\\}=0\quad{\rm with}\quad\mathcal{C}^{2}=\mathds{1},$ (1)

where $\\{,\\}$ denotes the anticommutator and $\mathcal{C}$ is the chiral
operator. In its eigenbasis it reads as

$\mathcal{C}=\begin{bmatrix}\mathds{1}&0\\\ 0&-\mathds{1}\end{bmatrix}$ (2)

and the Hamiltonian takes a block off-diagonal form

$H=\begin{bmatrix}0&K\\\ K^{\dagger}&0\end{bmatrix},$ (3)

where the matrix $K$ and its Hermitian adjoint $K^{\dagger}$ have dimension
$N\times N$. We want to consider the chiral orthogonal symmetry class, labeled
BDI in the tenfold way. In this class the time reversal operator $\mathcal{T}$
is given by a complex conjugation, such that $\mathcal{T}^{2}=+1$. This
renders the Hamiltonian, fulfilling the equation

$\mathcal{T}H\mathcal{T}^{-1}=H^{*}=H,$ (4)

as well as the matrix $K$ as real. In random matrix theory this case is
usually labeled as $\beta=1$, where the Dyson index $\beta$ is the dimension
of the underlying number field. Correspondingly, class AIII, featuring complex
entries and class CII, featuring quaternion entries are labeled by $\beta=2$
and $\beta=4$, respectively.

In the following we choose the dimension $N$ of the off-diagonal blocks as
even. This is due to technical reasons, that exist exclusively for the real
case. However, this does not interfere with our goal to uncover universal
behaviour in the large $N$ limit as is outlined in more detail further below.

Drawing all independent elements of $H$ from a Gaussian probability
distribution leads to the chiral Gaussian orthogonal ensemble (chGOE), while
at the same time $K$ can be viewed as a member of the real Ginibre ensemble
[32]. The matrix probability densities of these ensembles are invariant under
orthogonal transformations.

To examine topological properties we endow the Hamiltonian with a parametric
dependence, $H=H(p)$ and thus $K=K(p)$, that we assume to be periodic
$H(p+2\pi)=H(p)$. Physically, we interpret $H(p)$ as the Bloch Hamiltonian and
$p$ as a (quasi-)momentum in the Brillouin zone. Therefore upon time reversal
$K(p)$ behaves as

$\mathcal{T}K(p)\mathcal{T}^{-1}=K^{*}(p)=K(-p).$ (5)

Consequently, $K(p)$ will be complex in general and real only at the time
reversal invariant momenta $p=0$ and $p=\pi$. Then also $H(p)$ will be non-
real, despite being in the class BDI. However, our results remain valid also
in the case that condition (5) is relaxed and the parameter $p$ is interpreted
otherwise.

Figure 1: A realization of a BDI Hamiltonian with
$K(p)=\cos(p)K_{1}+i\sin(p)K_{2}$ and some fixed $4\times 4$ real matrices
$K_{1}$ and $K_{2}$. The top left plot shows the real eigenvalues of $H(p)$,
the top right one shows the generically complex eigenvalues of $K(p)$, and the
bottom left plot depicts the determinant $\det K(p)$. The bottom right plot
shows the determinant of a different random matrix field
$K(p)=(a_{1}+a_{2}e^{ip}+a_{3}e^{2ip}+a_{4}e^{3ip})K_{1}+(b_{1}+b_{2}e^{ip}+b_{3}e^{2ip}+b_{4}e^{3ip})K_{2}$
with fixed $4\times 4$ real matrices $K_{1}$ and $K_{2}$ and real coefficients
$a_{j}$ and $b_{j}$. All plots show the parametric dependence in
$p\in[0,2\pi]$ where we have employed the step size $2\pi/100$ and a B-Spline
to obtain the curves. In the parametric plots the starting points $p=0$ are
marked by black points and the directions are marked by a color gradient resp.
arrows.

The pertinent topological invariant in the chiral classes AIII, BDI and CII is
the winding number of $\det\,K(p)$ around the origin. It can be expressed as
the integral[33, 34]

$W=\frac{1}{2\pi i}\int\limits_{0}^{2\pi}dp\ w(p),$ (6)

where the logarithmic derivative

$w(p)=\frac{d}{dp}\ln\det K(p)=\frac{1}{\det K(p)}\frac{d}{dp}\det K(p)$ (7)

is the winding number density. This definition corresponds directly to the
usual definition in complex analysis by Cauchy’s argument principle when
writing the integral as a contour integral in $e^{ip}$. The real- and
imaginary parts of $w(p)$ have different interpretations. While the real part
describes the radial part of $\det K(p)$ and always integrates to zero over a
closed contour, the imaginary part describes its angle. Therefore the latter
would be sufficient to define the winding number.

In Fig. 1 we illustrate the spectral flow of $K(p)=\cos(p)K_{1}+i\sin(p)K_{2}$
with generic real $K_{1},K_{2}\in\mathbb{R}^{4\times 4}$. We plotted the eight
real eigenvalues of the Hamiltonian $H(p)$, showing a reflection symmetry at
$E=0$ due to chiral symmetry, and the four complex eigenvalues of $K(p)$ as
well as the determinant $\det\,K(p)$. Due to level repulsion no degeneracies
occur in the spectrum of $H(p)$, which implies that each eigenvalue is a
$2\pi$-periodic function. This is not the case for the complex eigenvalues of
$K(p)$, which may experience a permutation when running from $p=0$ to
$p=2\pi$, as can be observed in Fig. 1. However, the determinant of $K(p)$ is
always a $2\pi$-periodic function and therefore suitable to define a winding
number. For our specific choice of the parametric dependence we have
$K(p+\pi)=-K(p)$, which restricts the amount of times $\det\,K(p)$ winds
around the origin to be an even resp. odd number for even resp. odd matrix
dimensions. This symmetry does not exist in general as we also illustrated in
Fig. 1 by means of another random matrix field
$K(p)=(a_{1}+a_{2}e^{ip}+a_{3}e^{2ip}+a_{4}e^{3ip})K_{1}+(b_{1}+b_{2}e^{ip}+b_{3}e^{2ip}+b_{4}e^{3ip})K_{2}$
with generic real $K_{1},K_{2}\in\mathbb{R}^{4\times 4}$ and real coefficients
$a_{j}$ and $b_{j}$. In this example we find an odd winding number.

The constraint (5) is also reflected in the spectral flow of the Hamiltonian
when considering the symmetry class CII, like we did in Ref. 18. The Bloch
Hamiltonian becomes self-dual at the time reversal invariant momenta and
therefore its eigenvalues show Kramers’ degeneracy at these points on the
parameter manifold.

The simplest choice for a non-trivial parametric dependence on a one-
dimensional manifold is the random matrix field

$K(p)=a(p)K_{1}+b(p)K_{2},$ (8)

where $a(p)$ and $b(p)$ are smooth $2\pi$-periodic and not simultaneously
vanishing functions and $K_{1}$ and $K_{2}$ are real Ginibre matrices. For
notational reasons we gather the coefficient functions in a complex two-
dimensional vector

$v(p)=(a(p),b(p))\in\mathbb{C}^{2}\setminus\\{0\\}.$ (9)

In our case, where time reversal invariance is present, condition (5)
translates to

$v^{*}(p)=v(-p).$ (10)

We denote an average over the Ginibre matrices $K_{1}$ and $K_{2}$ as

$\langle F\rangle=\int d[K_{1},K_{2}]P(K_{1})P(K_{2})F(K_{1},K_{2}).$ (11)

Our goal is to compute the ensemble average over a ratio of determinants

$Z^{(1,N)}_{k|l}(q,p)=\left\langle\frac{\prod_{j=1}^{l}\det
K(p_{j})}{\prod_{j=1}^{k}\det K(q_{j})}\right\rangle.$ (12)

The superscript refers to the Dyson index $\beta=1$, i.e. it defines the
symmetry class, and the dimension $N$ of the block matrix $K(p)$. Although we
are only interested in the case $k=l$ we employ the more general definition
with $k\neq l$ as this function can be traced back to averages of only two
determinants for which $k+l=2$. This is shown in Sec. IV, see also Refs. 21,
22 for more details and the derivation of this result. This quantity is the
key to the winding number statistics in form of the $k$-point correlation of
winding number densities

$C^{(1,N)}_{k}(p_{1},\ldots,p_{k})=\langle w(p_{1})\cdots w(p_{k})\rangle,$
(13)

generated by taking the derivative of (12)

$C^{(1,N)}_{k}(p_{1},\ldots,p_{k})=\frac{\partial^{k}}{\prod_{j=1}^{k}\partial
p_{j}}Z^{(1,N)}_{k|k}(q,p)\Bigg{|}_{q=p}.$ (14)

We relegate the evaluation of these derivatives and taking the large $N$ limit
in order to reveal universal properties of the statistics to a subsequent
paper.

##  III Result

It is well known that the orthogonal case (as well as the symplectic one) is
described by a Pfaffian point process [9]. This means that the eigenvalue
correlations can be written as a Pfaffian of an ensemble dependent kernel
function. This will also be reflected in our results. Via the methods
developed in Refs. 21, 22 the ensemble average (12) for the two matrix model
(8) evaluates to a Pfaffian

$Z^{(1,N)}_{k|k}(q,p)=\frac{1}{\det\left[\displaystyle\frac{1}{iv^{T}(q_{m})\tau_{2}v(p_{n})}\right]_{1\leq
m,n\leq
k}}{\text{Pf}}\left[\begin{array}[]{c|c}\displaystyle\underset{}{\widehat{\rm
K}_{1}(p_{m},p_{n})}&\displaystyle\underset{}{\widehat{\rm
K}_{2}(p_{m},q_{n})}\\\ \hline\cr\displaystyle\overset{}{-\widehat{\rm
K}_{2}(p_{n},q_{m})}&\displaystyle\overset{}{\widehat{\rm
K}_{3}(q_{m},q_{n})}\end{array}\right]_{1\leq m,n\leq k}$ (15)

of a $2\times 2$-matrix with $k\times k$ blocks. We chose the sign of the
Pfaffian as

${\text{Pf}}[i\tau_{2},i\tau_{2},\ldots,i\tau_{2}]=1,$ (16)

where $\tau_{2}$ is the second Pauli matrix. The blocks are the three kernel
functions

$\begin{split}\widehat{\rm
K}_{1}(p_{m},p_{n})=&\frac{N(N-1)}{4\pi}\,iv^{T}(p_{n})\tau_{2}v(p_{m})\left[v^{T}(p_{m})v(p_{n})\right]^{N-2},\\\
\widehat{\rm
K}_{2}(p_{n},q_{m})=&\frac{1}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\frac{N(N-1)}{4\pi}\left(\frac{v^{T}(p_{n})v(p_{n})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{N}\\\
&\times\int\limits_{\mathbb{C}^{2}}d[z]\,z_{1}^{N-2}g^{(1,N)}(z_{1},z_{2})\left(z_{2}+\frac{v^{T}(p_{n})v(q_{m})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{-1},\\\
\widehat{\rm
K}_{3}(q_{m},q_{n})=&\frac{1}{[b(q_{m})b(q_{n})]^{N}}\Bigg{[}\int\limits_{\mathbb{C}^{2}}d[z]\frac{g^{(1,N)}(z_{1},z_{2})}{(a(q_{m})+b(q_{m})\,z_{1})(a(q_{n})+b(q_{n})\,z_{2})}\\\
&-\frac{N(N-1)}{4\pi}\int\limits_{\mathbb{C}^{4}}d[z]\frac{(z_{3}-z_{1})(z_{1}z_{3}+1)^{N-2}}{(a(q_{m})+b(q_{m})\,z_{2})(a(q_{n})+b(q_{n})\,z_{4})}g^{(1,N)}(z_{1},z_{2})g^{(1,N)}(z_{3},z_{4})\Bigg{]}\end{split}$
(17)

which are calculated in Sec. IV. We made use of the notation (9), which,
except for the third kernel function, nicely reflects a symmetry of our random
matrix field, that we will illuminate more closely in a moment. The
antisymmetric function

$g^{(1,N)}(z_{1},z_{2})=\frac{\left|z_{2}-z_{1}\right|}{z_{2}-z_{1}}\frac{B(1/2,(N+1)/2)\delta(y_{1})\delta(y_{2})+2\delta(x_{1}-x_{2})\delta(y_{1}+y_{2})Q(z_{1},z_{1}^{*})}{\left[(1+z_{1}^{2})(1+z_{2}^{2})\right]^{(N+1)/2}}$
(18)

directly emerges from our random matrix problem as part of a joint eigenvalue
probability density. For the complex integration variables we use the notation
$z_{j}=x_{j}+iy_{j}$ with $x_{j},y_{j}\in\mathbb{R}$. Furthermore, a lower
incomplete Beta function appears in our results

$\begin{split}Q(z,z^{*})&=B\left(\frac{4y^{2}}{\left|1+z^{2}\right|^{2}+4y^{2}};1/2,(N+1)/2\right),\\\
B(x;a,b)&=\int\limits_{x}^{1}dt\,t^{a-1}(1-t)^{b-1}=2\int\limits_{\sqrt{\frac{x}{1-x}}}^{\infty}dt\frac{t^{2a-1}}{(1+t^{2})^{a+b}},\text{
with }x\in[0,1],\\\
B(0;a,b)&=B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}.\end{split}$ (19)

Inserting (18) into the second and third kernels yields the full expressions

$\begin{split}\widehat{\rm
K}_{2}(p_{n},q_{m})=&\frac{1}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\frac{N(N-1)}{2\pi}\left(\frac{v^{T}(p_{n})v(p_{n})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{N}\\\
&\times\Bigg{[}\frac{(-1)^{N/2}\pi\,B\left(1/2,\frac{N+1}{2}\right)}{N-1}{\frac{N-1}{2}\choose
N}{{}_{2}F_{1}}\left(1,\frac{N+1}{2};N+1;1+\left(\frac{v^{T}(p_{n})v(q_{m})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{2}\right)\\\
&+i\int\limits_{\mathbb{C}}d[z]\frac{z^{N-2}\,\text{sgn}(\text{Im}\,z)Q(z,z^{*})}{\left|1+z^{2}\right|^{N+1}}\left(z^{*}+\frac{v^{T}(p_{n})v(q_{m})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{-1}\Bigg{]}\end{split}$
(20)

with the hypergeometric function ${{}_{2}F_{1}}$ and

$\begin{split}\widehat{\rm
K}_{3}(q_{m},q_{n})=&\frac{1}{[b(q_{m})b(q_{n})]^{N}}\left[\int\limits_{\mathbb{R}}dx\frac{r(x,v(q_{n}))}{(a(q_{m})+b(q_{m})x)}+\int\limits_{\mathbb{C}}d[z]\frac{s(z,z^{*},v(q_{n}))}{(a(q_{m})+b(q_{m})z)}\right]\\\
&-\frac{N(N-1)}{2\pi[b(q_{m})b(q_{n})]^{N-1}}\Bigg{[}\int\limits_{\mathbb{R}^{2}}d[x]r(x_{1},v(q_{m}))r(x_{2},v(q_{n}))(x_{2}-x_{1})(x_{1}x_{2}+1)^{N-2}\\\
&+\int\limits_{\mathbb{R}}dx\int\limits_{\mathbb{C}}d[z]\det\begin{bmatrix}r(x,v(q_{m}))&s(z,z^{*},v(q_{m}))\\\
r(x,v(q_{n}))&s(z,z^{*},v(q_{n}))\end{bmatrix}(z-x)(zx+1)^{N-2}\\\
&+\int\limits_{\mathbb{C}^{2}}d[z]s(z_{1},z_{1}^{*},v(q_{m}))s(z_{2},z_{2}^{*},v(q_{n}))(z_{2}-z_{1})(z_{1}z_{2}+1)^{N-2}\Bigg{]}\end{split}$
(21)

with the functions

$\begin{split}r(x,v(q))&=B(1/2,(N+1)/2)\int\limits_{\mathbb{R}}dx^{\prime}\frac{\text{sgn}(x^{\prime}-x)}{(a(q)+b(q)x^{\prime})[(1+x^{2})(1+x^{\prime
2})]^{(N+1)/2}},\\\
s(z,z^{*},v(q))&=\frac{2i\,\text{sgn}(\text{Im}\,z)Q(z,z^{*})}{(a(q)+b(q)z^{*})\left|1+z^{2}\right|^{N+1}}.\end{split}$
(22)

These expressions still contain integrals over up to four real dimensions. A
further analytic computation of them is difficult due to the singularities
appearing at various points in the complex plane. However, an evaluation in
the large $N$ limit, which is our ultimate goal for the correlation function
(13), should be feasible.

##  IV Derivation

We start by highlighting a symmetry of the parametric matrix model (8), that
we will use to simplify the necessary integrals. We continue by tracing the
averages over the Ginibre ensembles back to an average over another random
matrix ensemble. Then we give the expressions for the three kernel functions,
as obtained from the supersymmetry without supersymmetry technique and
calculate those one by one.

The matrices $K_{1}$ and $K_{2}$ are drawn from a real Ginibre distribution,
i.e. their joint matrix distribution is

$P(K_{1})P(K_{2})=(2\pi)^{-N^{2}}\exp\left[-\frac{1}{2}\text{tr}\,\left(K_{1}K_{1}^{T}+K_{2}K_{2}^{T}\right)\right].$
(23)

This function is invariant under O(2) transformations

$\widehat{K}=\begin{bmatrix}K_{1}\\\
K_{2}\end{bmatrix}\to\left[U\otimes\mathds{1}_{N}\right]\begin{bmatrix}K_{1}\\\
K_{2}\end{bmatrix},$ (24)

where $U\in\text{O}(2)$ acts on the two component matrix valued vector
$\widehat{K}$. This invariance is inherited by the vectors $v(p)$

$K(p)=a(p)K_{1}+b(p)K_{2}=v^{T}(p)\widehat{K}.$ (25)

Assuming for the moment $b(p)\neq 0$, we may rephrase the random matrix field

$K(p)=a(p)K_{1}+b(p)K_{2}=b(p)K_{1}\left(\kappa(p)\mathds{1}_{N}+K_{1}^{-1}K_{2}\right),\
{\rm with}\ \kappa(p)=\frac{a(p)}{b(p)},$ (26)

such that the generating function becomes

$Z^{(1,N)}_{k|k}(q,p)=\left(\prod_{j=1}^{k}\frac{b(p_{j})}{b(q_{j})}\right)^{N}\left\langle\prod_{j=1}^{k}\frac{\det(\kappa(p_{j})\mathds{1}_{N}+K_{1}^{-1}K_{2})}{\det(\kappa(q_{j})\mathds{1}_{N}+K_{1}^{-1}K_{2})}\right\rangle.$
(27)

The random matrix $Y=K_{1}^{-1}K_{2}$ defines the real spherical ensemble,
which was studied in Ref. 35. The spherical ensemble is also well studied for
the complex and the quaternion case, see Refs. 36, 37, which contain the
results we used for our studies in the cases AIII and CII[17, 18]. Its matrix
probability density function in the real case is

$\widetilde{G}^{(1,N)}(Y)=\pi^{-N^{2}/2}\prod_{j=1}^{N}\frac{\Gamma((N+j)/2)}{\Gamma(j/2)}\frac{1}{\det^{N}\left(\mathds{1}_{N}+YY^{T}\right)}.$
(28)

Unfortunately, in the literature the joint eigenvalue probability density
exists only in a form, that is impractical for our purpose. Therefore we
derive in Appendix A the following expression for $N$ even

$\begin{split}G^{(1,N)}(z)&=C^{(N)}\Delta_{N}(z)\prod_{j=1}^{N/2}g^{(1,N)}(z_{2j-1},z_{2j}),\\\
C^{(N)}&=\frac{2^{-N/2}\pi^{-N/4}}{(N/2)!}\prod_{j=1}^{N}\frac{\Gamma((N+j)/2)}{\Gamma(j/2)}\prod_{j=1}^{N/2}\frac{\Gamma\left((N+1)/2\right)\Gamma\left(N/2+1\right)}{\Gamma\left(N+1/2-j\right)\Gamma\left(N+1-j\right)}\end{split}$
(29)

with the antisymmetric function (18). The eigenvalues of a real matrix are
either real or come in complex conjugated pairs, which is what (18) describes.
Although the part of this function concerning the complex eigenvalues is
complex valued, the full joint eigenvalue probability density will be real as
the corresponding phases drop out when considering the product with the
Vandermonde determinant.

Only in the odd dimensional case we have to supplement this with an additional
function for the unpaired eigenvalue, that must be real. This would complicate
our already cumbersome expressions even more, which is why we decided to limit
ourselves to the even case. This restriction is not an obstacle to our plan to
deduce universal statistics in the large $N$ limit as the number of expected
real eigenvalues is asymptotically given by $\sqrt{\pi N/2}$, regardless of
the parity of $N$ [35].

In the supersymmetry without supersymmetry technique [21, 22] one starts from
integrals in a joint eigenvalue density. These can be reformulated by
introducing superspace Jacobians, which in the present case are mixtures of
Vandermonde and Cauchy determinants. This allows the application of a
generalized version of de Bruijn’s integration theorem [38, 22] resulting in a
high dimensional block structured Pfaffian. After reducing the dimensionality
by working out the Schur complement, the blocks of the Pfaffian can be
identified as simplified averages over only two determinants.

Applying this method to (27) yields

$Z^{(1,N)}_{k|k}(q,p)=\frac{1}{\sqrt{\text{Ber}_{k|k}^{(2)}(\kappa)}}{\text{Pf}}\left[\begin{array}[]{c|c}\displaystyle\underset{}{{\rm
K}_{1}^{(1)}(p_{m},p_{n})}&\displaystyle\underset{}{{\rm
K}_{2}^{(1)}(p_{m},q_{n})}\\\ \hline\cr\displaystyle\overset{}{-{\rm
K}_{2}^{(1)}(p_{n},q_{m})}&\displaystyle\overset{}{{\rm
K}_{3}^{(1)}(q_{m},q_{n})}\end{array}\right]_{1\leq m,n\leq k},$ (30)

where the kernel functions are the simplified averages

$\begin{split}{\rm
K}_{1}^{(1)}(p_{m},p_{n})&=[\kappa(p_{n})-\kappa(p_{m})][b(p_{m})b(p_{n})]^{N}\widetilde{Z}_{0|2}^{(1,N-2)}(p_{m},p_{n}),\\\
{\rm
K}_{2}^{(1)}(p_{n},q_{m})&=\left(\frac{b(p_{n})}{b(q_{m})}\right)^{N}\frac{\widetilde{Z}_{1|1}^{(1,N)}(q_{m},p_{n})}{\kappa(q_{m})-\kappa(p_{n})}=\frac{Z_{1|1}^{(1,N)}(q_{m},p_{n})}{\kappa(q_{m})-\kappa(p_{n})},\\\
{\rm
K}_{3}^{(1)}(q_{m},q_{n})&=\frac{\kappa(q_{n})-\kappa(q_{m})}{[b(q_{m})b(q_{n})]^{N}}\widetilde{Z}_{2|0}^{(1,N+2)}(q_{m},q_{n}).\end{split}$
(31)

We define for $2M+l-k\leq N$

$\begin{split}\widetilde{Z}^{(1,2M)}_{k|l}(q,p)=&\frac{1}{M!\,{\text{Pf}}\,D^{(2M+l-k)}}\\\
&\times\int\limits_{\mathbb{C}^{2M}}d[z]\Delta_{2M}(z)\prod_{j=1}^{M}g^{(1,N)}(z_{2j-1},z_{2j})\prod_{j=1}^{2M}\frac{\prod_{n=1}^{l}(\kappa(p_{n})+z_{j})}{\prod_{m=1}^{k}(\kappa(q_{m})+z_{j})}\end{split}$
(32)

with the moment matrix

$\begin{split}D^{(d)}&=\left[D_{ab}\right]_{1\leq a,b\leq d},\\\ D_{ab}&=\int
d[z]g^{(1,N)}(z_{1},z_{2})\,(z_{1}^{a-1}z_{2}^{b-1}-z_{1}^{b-1}z_{2}^{a-1})=2\int
d[z]g^{(1,N)}(z_{1},z_{2})\,z_{1}^{a-1}z_{2}^{b-1}.\end{split}$ (33)

The Berezinian is a superspace analogue to the Jacobian. In this case it is
given by a Cauchy determinant[39, 21]

$\sqrt{\text{Ber}_{k|k}^{(2)}(\kappa)}=\det\left[\displaystyle\frac{1}{\kappa(q_{m})-\kappa(p_{n})}\right]_{1\leq
a,b\leq k}.$ (34)

In the following subsections we calculate more explicit expressions of the
kernels (31). Furthermore, to arrive at the results (15) and (17) one has to
rearrange some factors in order to highlight the symmetry (24) of our problem.

###  IV.1 The Kernel ${\rm K}_{1}^{(1)}$

The first kernel is determined by
$\widetilde{Z}_{0|2}^{(1,N-2)}(p_{m},p_{n})$, which is related to
$Z_{0|2}^{(1,N-2)}(p_{m},p_{n})$. For the latter function we can use the
O(2)-symmetry to reduce the amount of determinants in the average by one. We
obtain the proper normalization via the limits

$\lim_{\kappa(p)\to\infty}\frac{\widetilde{Z}_{0|2}^{(1,N-2)}(p_{m},p_{n})}{[\kappa(p_{m})\kappa(p_{n})]^{N-2}}=\frac{{\text{Pf}}\,D^{(N-2)}}{{\text{Pf}}\,D^{(N)}},\qquad\qquad\lim_{a(p)\to\infty}\frac{Z_{0|2}^{(1,N-2)}(p_{m},p_{n})}{[a(p_{m})a(p_{n})]^{N-2}}=\left\langle\det\nolimits^{2}K_{1}\right\rangle.$
(35)

These limits are equivalent because $\kappa(p)$ is directly proportional to
$a(p)$, see (26). Therefore

$\widetilde{Z}_{0|2}^{(1,N-2)}(p_{m},p_{n})=\frac{{\text{Pf}}\,D^{(N-2)}}{{\text{Pf}}\,D^{(N)}}\frac{1}{\left\langle\det^{2}K_{1}\right\rangle}\frac{Z_{0|2}^{(1,N-2)}(p_{m},p_{n})}{\left[b(p_{m})b(p_{n})\right]^{N-2}}.$
(36)

Now let $a_{1},b_{1},a_{2},b_{2}\in\mathbb{R}$ and

$\Xi_{1}=\frac{\left\langle\det\left(a_{1}K_{1}+b_{1}K_{2}\right)\det\left(a_{2}K_{1}+b_{2}K_{2}\right)\right\rangle}{\left\langle\det^{2}K_{1}\right\rangle}.$
(37)

This function is a polynomial in $a_{1},b_{1},a_{2},b_{2}$ and can thus be
analytically continued to the partition function
$Z_{0|2}^{(1,N-2)}(p_{m},p_{n})$, that has, in general, complex coefficients.
We are now finally in the position to exploit the O(2)-symmetry as we are
allowed to rotate only real vectors. We use

$U=\frac{1}{\sqrt{a_{1}^{2}+b_{1}^{2}}}\begin{bmatrix}a_{1}&-b_{1}\\\
b_{1}&a_{1}\end{bmatrix}\in\text{SO}(2)$ (38)

yielding

$\Xi_{1}=\frac{(a_{1}b_{2}-b_{1}a_{2})^{N-2}}{\left\langle\det^{2}K_{1}\right\rangle}\int\limits
d[K_{1}]d[K_{2}]P(K_{1})P(K_{2})\det\nolimits^{2}K_{1}\det\left(\frac{a_{1}a_{2}+b_{1}b_{2}}{a_{1}b_{2}-b_{1}a_{2}}+K_{1}^{-1}K_{2}\right).$
(39)

Up to the factor $\det^{2}K_{1}$ appearing in the integrand this is almost an
integral over the spherical ensemble. In fact, this type of measure leads to a
generalized form of the spherical ensemble of matrices $K_{1}^{-1}K_{2}$,
where one or both of $K_{1}$ and $K_{2}$ are drawn from a deformed Ginibre
ensemble

$P_{\mu}(K)=2^{-N\mu}(2\pi)^{-N^{2}/2}\prod_{j=1}^{N}\frac{\Gamma(j/2)}{\Gamma(\mu+j/2)}\exp\left(-\frac{1}{2}\text{tr}\,\,KK^{T}\right)\det\nolimits^{\mu}KK^{T}.$
(40)

Matrix ensembles of such kind are called induced spherical ensembles. They are
well studied for $\beta=1,2,4$ [40, 41, 42]. In appendix B we show that the
characteristic polynomial of the real spherical ensemble, induced or not, is
given by a monomial. Thus, we find

$\Xi_{1}=(a_{1}a_{2}+b_{1}b_{2})^{N-2}$ (41)

and the partition function $Z_{0|2}^{(1,N-2)}(p_{m},p_{n})$ has the form

$\frac{Z_{0|2}^{(1,N-2)}(p_{m},p_{n})}{\left\langle\det^{2}K_{1}\right\rangle}=[a(p_{m})a(p_{n})+b(p_{m})b(p_{n})]^{N-2}=v^{T}(p_{m})v(p_{n}),$
(42)

which is indeed O(2)-invariant. The Pfaffians of the moment matrix can be
related to the normalization of an induced spherical ensemble. This is covered
in appendix A. We find

$\frac{{\text{Pf}}\,D^{(N-2)}}{{\text{Pf}}\,D^{(N)}}=\frac{N(N-1)}{8\pi}$ (43)

and altogether for the first kernel

$\begin{split}{\rm
K}_{1}^{(1)}(p_{m},p_{n})=&\frac{N(N-1)}{8\pi}b(p_{m})b(p_{n})[a(p_{n})b(p_{m})-a(p_{m})b(p_{n})]\\\
&\times[a(p_{m})a(p_{n})+b(p_{m})b(p_{n})]^{N-2}\\\
=&\frac{N(N-1)}{8\pi}b(p_{m})b(p_{n})\,iv^{T}(p_{n})\tau_{2}v(p_{m})[v^{T}(p_{m})v(p_{n})]^{N-2},\end{split}$
(44)

which by itself is not O(2)-invariant, just as the other kernels (31) will not
be. However, the final expression for $Z_{k|k}^{(1,N)}(q,p)$ will be, because
the non-invariant prefactor will be compensated.

###  IV.2 The Kernel ${\rm K}_{2}^{(1)}$

The second kernel is essentially the generator $Z_{1|1}^{(1,N)}(q_{m},p_{n})$
of the one-point function. Once again we want to apply the O(2)-symmetry to
eliminate one of the determinants, for this purpose we consider

$\Xi_{2}=\left\langle\frac{\det\left(a_{1}K_{1}+b_{1}K_{2}\right)}{\det\left(a_{2}K_{1}+b_{2}K_{2}\right)}\right\rangle.$
(45)

This function is a polynomial in $a_{1},b_{1}$ and non-holomorphic in
$a_{2},b_{2}$. Therefore we may apply our arguments only to the numerator and
have to set $a_{1},b_{1}\in\mathbb{R}$ and $a_{2},b_{2}\in\mathbb{C}$. We
rotate with the same matrix (38) as for the first kernel and obtain an
integral over the spherical ensemble

$\Xi_{2}=\left(\frac{a_{1}^{2}+b_{1}^{2}}{a_{1}b_{2}-b_{1}a_{2}}\right)^{N}\int\limits
d[Y]\widetilde{G}^{(1,N)}(Y)\frac{1}{\det\left(\widehat{\kappa}+Y\right)}$
(46)

with

$\widehat{\kappa}=\frac{a_{1}a_{2}+b_{1}b_{2}}{a_{1}b_{2}-b_{1}a_{2}}.$ (47)

We want to perform the integral in the joint probability density of the
eigenvalues

$\Xi_{2}=C^{(N)}\left(\frac{a_{1}^{2}+b_{1}^{2}}{a_{1}b_{2}-b_{1}a_{2}}\right)^{N}\int\limits_{\mathbb{C}^{N}}d[z]\Delta_{N}(z)\prod_{j=1}^{N/2}g^{(1,N)}(z_{2j-1},z_{2j})\prod_{j=1}^{N}\frac{1}{\widehat{\kappa}+z_{j}}.$
(48)

First, we identify part of the integrand as a Berezinian

$\Delta_{N}(z)\prod_{j=1}^{N}\frac{1}{\widehat{\kappa}+z_{j}}=\sqrt{\text{Ber}^{(2)}_{N|1}\left(z;-\widehat{\kappa}\right)}=(-1)^{N+1}\det\left[\begin{array}[]{c|c}\displaystyle
z_{a}^{b-1}&\displaystyle\frac{1}{z_{a}+\widehat{\kappa}}\end{array}\right]_{\begin{subarray}{c}1\leq
a\leq N\\\ 1\leq b\leq N-1\end{subarray}},$ (49)

which in turn can be written as a determinant [21]. We expand this determinant
in the last column. Using its skew-symmetry under row permutation and the
antisymmetry of $g^{(1,N)}(z_{2j-1},z_{2j})$ we find that each terms yields
the same contribution, which is

$\begin{split}\Xi_{2}&=(-1)^{N+1}N\,C^{(N)}\left(\frac{a_{1}^{2}+b_{1}^{2}}{a_{1}b_{2}-b_{1}a_{2}}\right)^{N}\int\limits_{\mathbb{C}^{N}}d[z]\frac{\Delta_{N-1}(z_{1},\ldots,z_{N-1})}{z_{N}+\widehat{\kappa}}\prod_{j=1}^{N/2}g^{(1,N)}(z_{2j-1},z_{2j})\\\
&=-\frac{N(N-1)}{4\pi}\left(\frac{a_{1}^{2}+b_{1}^{2}}{a_{1}b_{2}-b_{1}a_{2}}\right)^{N}\int\limits_{\mathbb{C}^{2}}d[z]\frac{z_{1}^{N-2}}{z_{2}+\widehat{\kappa}}g^{(1,N)}(z_{1},z_{2}).\end{split}$
(50)

The integral over $z_{1},\ldots,z_{N-2}$ can be conceived as a characteristic
polynomial of an induced spherical ensemble in the variable $z_{N-1}$, which
gives up to a prefactor a monomial. The details of this calculation are
covered in appendix B. The integral over the remaining two eigenvalues yields
two contributions, one for the case in which they are real and one for the
case in which they are a complex conjugate pair

$\Xi_{2}=-\frac{N(N-1)}{4\pi}\left(\frac{a_{1}^{2}+b_{1}^{2}}{a_{1}b_{2}-b_{1}a_{2}}\right)^{N}[I_{R}(\widehat{\kappa})+I_{C}(\widehat{\kappa})].$
(51)

The contribution of the real eigenvalues is

$\begin{split}I_{R}(\widehat{\kappa})&=B(1/2,(N+1)/2)\int\limits_{\mathbb{R}^{2}}d[x]\frac{\text{sgn}(x_{2}-x_{1})}{\left[(1+x_{1}^{2})(1+x_{2}^{2})\right]^{(N+1)/2}}\frac{x_{1}^{N-2}}{x_{2}+\widehat{\kappa}}\\\
&=B(1/2,(N+1)/2)\int\limits_{x_{2}>x_{1}}d[x]\frac{1}{\left[(1+x_{1}^{2})(1+x_{2}^{2})\right]^{(N+1)/2}}\left(\frac{x_{1}^{N-2}}{x_{2}+\widehat{\kappa}}-\frac{x_{2}^{N-2}}{x_{1}+\widehat{\kappa}}\right),\end{split}$
(52)

where we treat the sign function by splitting the integral into two terms,
that we calculate separately. In the first term we integrate $x_{1}$ over
$(-\infty,x_{2}]$ yielding

$\begin{split}&\int\limits_{-\infty}^{\infty}dx_{2}\int\limits_{-\infty}^{x_{2}}dx_{1}\frac{1}{\left[(1+x_{1}^{2})(1+x_{2}^{2})\right]^{(N+1)/2}}\frac{x_{1}^{N-2}}{x_{2}+\widehat{\kappa}}\\\
=&\int\limits_{-\infty}^{\infty}dx_{2}\frac{1}{(1+x_{2}^{2})^{(N+1)/2}}\frac{1}{x_{2}+\widehat{\kappa}}\frac{(-1)^{N}+x_{2}^{N-1}(1+x_{2}^{2})^{(1-N)/2}}{N-1}\end{split}$
(53)

and in the second term we integrate $x_{2}$ over $[x_{1},\infty)$

$\begin{split}&\int\limits_{-\infty}^{\infty}dx_{1}\int\limits_{x_{1}}^{\infty}dx_{2}\frac{1}{\left[(1+x_{1}^{2})(1+x_{2}^{2})\right]^{(N+1)/2}}\frac{x_{2}^{N-2}}{x_{1}+\widehat{\kappa}}\\\
=&\int\limits_{-\infty}^{\infty}dx_{1}\frac{1}{(1+x_{1}^{2})^{(N+1)/2}}\frac{1}{x_{1}+\widehat{\kappa}}\frac{1-x_{1}^{N-1}(1+x_{1}^{2})^{(1-N)/2}}{N-1}.\end{split}$
(54)

For both integrals the antiderivative is

$\int
dx\frac{x^{N-2}}{(1+x^{2})^{(N+1)/2}}=\frac{1}{N-1}\frac{x^{N-1}}{(1+x^{2})^{(N-1)/2}}.$
(55)

Considering $N$ even one of the terms drops out and (52) results in an
integral over the real axis

$I_{R}(\widehat{\kappa})=\frac{2\,B(1/2,(N+1)/2)}{N-1}\int\limits_{\mathbb{R}}dx\frac{x^{N-1}}{(1+x^{2})^{N}}\frac{1}{x+\widehat{\kappa}}.$
(56)

In the case that $\widehat{\kappa}$ is real it has to be understood as a
principal value integral. We continue by symmetrizing the integrand and
substituting $t=x^{2}$

$\begin{split}I_{R}(\widehat{\kappa})&=\frac{2\,B(1/2,(N+1)/2)}{N-1}\int\limits_{0}^{\infty}dx\frac{x^{N-1}}{(1+x^{2})^{N}}\left(\frac{1}{x+\widehat{\kappa}}+\frac{1}{x-\widehat{\kappa}}\right)\\\
&=\frac{2\,B(1/2,(N+1)/2)}{N-1}\int\limits_{0}^{\infty}dt\frac{t^{(N-1)/2}}{(1+t)^{N}}\frac{1}{t-\widehat{\kappa}^{2}}.\end{split}$
(57)

The last integrand has a branch cut along the negative real axis. It is
equivalent to an integral of the function

$f(z)=\frac{z^{(N-1)/2}}{(1-z)^{N}}\frac{-1}{z+\widehat{\kappa}^{2}}$ (58)

over a keyhole contour around the negative real axis. Applying the residue
theorem yields a truncated binomial series

$\begin{split}I_{R}(\widehat{\kappa})&=\frac{(-1)^{N/2}2\pi\,B(1/2,(N+1)/2)}{N-1}\left[\sum_{l=0}^{N-1}{(N-1)/2\choose
l}\frac{(-1)^{N-1-l}}{(1+\widehat{\kappa}^{2})^{N-l}}+\frac{(-\widehat{\kappa}^{2})^{(N-1)/2}}{(1+\widehat{\kappa}^{2})^{N}}\right]\\\
&=\frac{(-1)^{N/2}2\pi\,B(1/2,(N+1)/2)}{N-1}{(N-1)/2\choose
N}{{}_{2}F_{1}}\left(1,(N+1)/2;N+1;1+\widehat{\kappa}^{2}\right).\end{split}$
(59)

The second representation involves the hypergeometric function [43]. The
contribution of the complex conjugated eigenvalues is

$I_{C}(\widehat{\kappa})=2i\int\limits_{\mathbb{C}}d[z]\text{sgn}(\text{Im}\,z)\frac{z^{N-2}}{z^{*}+\widehat{\kappa}}\frac{Q(z,z^{*})}{\left|1+z^{2}\right|^{N+1}}.$
(60)

Therefore, we find for the second kernel

$\begin{split}{\rm
K}_{2}^{(1)}(p_{n},q_{m})=&\frac{N(N-1)}{2\pi}\frac{b(p_{n})b(q_{m})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\left(\frac{v^{T}(p_{n})v(p_{n})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{N}\\\
&\times\Bigg{[}\frac{(-1)^{N/2}\pi
B\left(1/2,\frac{N+1}{2}\right)}{N-1}{\frac{N-1}{2}\choose
N}{{}_{2}F_{1}}\left(1,\frac{N+1}{2};N+1;1+\left(\frac{v^{T}(p_{n})v(q_{m})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{2}\right)\\\
&+i\int\limits_{\mathbb{C}}d[z]\text{sgn}(\text{Im}\,z)\frac{z^{N-2}Q(z,z^{*})}{\left|1+z^{2}\right|^{N+1}}\left(z^{*}+\frac{v^{T}(p_{n})v(q_{m})}{iv^{T}(p_{n})\tau_{2}v(q_{m})}\right)^{-1}\Bigg{]}\end{split}$
(61)

Once again we note the O(2)-invariance of this function in the vectors
$v(p)=(a(p),b(p))$ up to a prefactor.

###  IV.3 The Kernel ${\rm K}_{3}^{(1)}$

The third kernel is determined by
$\widetilde{Z}^{(1,N+2)}_{2|0}(q_{m},q_{n})$. Therefore we need to evaluate
the following integral

$\Xi_{3}=\frac{1}{(N/2+1)!\,{\text{Pf}}\,D^{(N)}}\int\limits_{\mathbb{C}^{N+2}}d[z]\Delta_{N+2}(z)\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})\prod_{j=1}^{N+2}\frac{1}{(\kappa_{1}+z_{j})(\kappa_{2}+z_{j})}.$
(62)

This time we cannot use the O(2)-symmetry to simplify the integral, because it
cannot be traced back to the ensemble average (12). Instead, we proceed by
applying the identity [21]

$\begin{split}\Delta_{N+2}(z)\prod_{j=1}^{N+2}\frac{1}{(\kappa_{1}+z_{j})(\kappa_{2}+z_{j})}&=\frac{1}{\kappa_{1}-\kappa_{2}}\sqrt{\text{Ber}^{(2)}_{N+2|2}\left(z;-\kappa\right)}\\\
&=\frac{1}{\kappa_{2}-\kappa_{1}}\det\left[\begin{array}[]{c|c|c}\displaystyle
z_{a}^{b-1}&\displaystyle\frac{1}{z_{a}+\kappa_{1}}&\displaystyle\frac{1}{z_{a}+\kappa_{2}}\end{array}\right]_{\begin{subarray}{c}1\leq
a\leq N+2\\\ 1\leq b\leq N\end{subarray}}\end{split}$ (63)

and expand the determinant in the last two columns. Similar to (50) we use the
antisymmetry of $g^{(1,N)}(z_{2j-1},z_{2j})$ to reduce the amount of terms
appearing in the integral

$\begin{split}\Xi_{3}=&\frac{2\,C^{(N)}}{(N+2)(\kappa_{2}-\kappa_{1})}\int\limits_{\mathbb{C}^{N+2}}d[z]\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})\det\left[\begin{array}[]{c|c|c}\displaystyle
z_{a}^{b-1}&\displaystyle\frac{1}{z_{a}+\kappa_{1}}&\displaystyle\frac{1}{z_{a}+\kappa_{2}}\end{array}\right]_{\begin{subarray}{c}1\leq
a\leq N+2\\\ 1\leq b\leq N\end{subarray}}\\\
=&\frac{2\,C^{(N)}}{\kappa_{2}-\kappa_{1}}\Bigg{[}\int\limits_{\mathbb{C}^{N+2}}d[z]\frac{\Delta_{N}(z)}{(z_{N+1}+\kappa_{1})(z_{N+2}+\kappa_{2})}\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})\\\
&-N\int\limits_{\mathbb{C}^{N+2}}d[z]\frac{\Delta_{N}(z_{1},\ldots,z_{N-1},z_{N+1})}{(z_{N}+\kappa_{1})(z_{N+2}+\kappa_{2})}\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})\Bigg{]}.\end{split}$
(64)

In the first term the pair of eigenvalues $z_{N+1},z_{N+2}$ is decoupled from
$z_{1},\ldots,z_{N}$. Integration over the latter yields only a constant

$C^{(N)}\int\limits_{\mathbb{C}^{N+2}}d[z]\frac{\Delta_{N}(z)}{(z_{N+1}+\kappa_{1})(z_{N+2}+\kappa_{2})}\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})=\int\limits_{\mathbb{C}^{2}}d[z]\frac{g^{(1,N)}(z_{1},z_{2})}{(z_{1}+\kappa_{1})(z_{2}+\kappa_{2})}.$
(65)

In the second term we expand the Vandermonde determinant in the last two
variables

$\begin{split}&\phantom{=}C^{(N)}\int\limits_{\mathbb{C}^{N+2}}d[z]\frac{\Delta_{N}(z_{1},\ldots,z_{N-1},z_{N+1})}{(z_{N}+\kappa_{1})(z_{N+2}+\kappa_{2})}\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})\\\
&=C^{(N)}\int\limits_{\mathbb{C}^{N+2}}d[z]\Delta_{N-2}(z)\frac{z_{N+1}-z_{N-1}}{(z_{N}+\kappa_{1})(z_{N+2}+\kappa_{2})}\prod_{j=1}^{N-2}(z_{N+1}-z_{j})(z_{N-1}-z_{j})\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j}).\end{split}$
(66)

This allows us to identify the integral over $z_{1},\ldots z_{N-2}$ as the
function $\Xi_{1}$, see Eq. (37), that we calculated for the first kernel,
resulting in

$\begin{split}&\phantom{=}\frac{N-1}{4\pi}\int\limits_{\mathbb{C}^{4}}d[z]\frac{z_{3}-z_{1}}{(z_{2}+\kappa_{1})(z_{4}+\kappa_{2})}g^{(1,N)}(z_{1},z_{2})g^{(1,N)}(z_{3},z_{4})\frac{\left\langle\det\left(z_{1}K_{1}-K_{2}\right)\det\left(z_{3}K_{1}-K_{2}\right)\right\rangle}{\left\langle\det^{2}K_{1}\right\rangle}\\\
&=\frac{N-1}{4\pi}\int\limits_{\mathbb{C}^{4}}d[z]\frac{z_{3}-z_{1}}{(z_{2}+\kappa_{1})(z_{4}+\kappa_{2})}g^{(1,N)}(z_{1},z_{2})g^{(1,N)}(z_{3},z_{4})\,(z_{1}z_{3}+1)^{N-2}.\end{split}$
(67)

Altogether the third kernel function is given by

$\begin{split}{\rm
K}_{3}^{(1)}(q_{m},q_{n})=&\frac{2}{[b(q_{m})b(q_{n})]^{N-1}}\Bigg{[}\int\limits_{\mathbb{C}^{2}}d[z]\frac{g^{(1,N)}(z_{1},z_{2})}{(a(q_{m})+b(q_{m})z_{1})(a(q_{n})+b(q_{n})z_{2})}\\\
&-\frac{N(N-1)}{4\pi}\int\limits_{\mathbb{C}^{4}}d[z]\frac{(z_{3}-z_{1})(z_{1}z_{3}+1)^{N-2}}{(a(q_{m})+b(q_{m})z_{2})(a(q_{n})+b(q_{n})z_{4})}g^{(1,N)}(z_{1},z_{2})g^{(1,N)}(z_{3},z_{4})\Bigg{]}.\end{split}$
(68)

Unlike for the first and the second kernel in this result the O(2)-symmetry is
not visible. However, our end result has to entail this symmetry. We define
the functions

$\begin{split}r(x,v(q))&=B(1/2,(N+1)/2)\int\limits_{\mathbb{R}}dx^{\prime}\frac{\text{sgn}(x^{\prime}-x)}{(a(q)+b(q)x^{\prime})\,[(1+x^{2})(1+x^{\prime
2})]^{(N+1)/2}},\\\
s(z,z^{*},v(q))&=\frac{2i\,\text{sgn}(\text{Im}\,z)Q(z,z^{*})}{(a(q)+b(q)z^{*})\left|1+z^{2}\right|^{N+1}}\end{split}$
(69)

Inserting (18) yields the cumbersome expression

$\begin{split}{\rm
K}_{3}^{(1)}(q_{m},q_{n})=&\frac{2}{[b(q_{m})b(q_{n})]^{N-1}}\left[\int\limits_{\mathbb{R}}dx\frac{r(x,v(q_{n}))}{(a(q_{m})+b(q_{m})x)}+\int\limits_{\mathbb{C}}d[z]\frac{s(z,z^{*},v(q_{n}))}{(a(q_{m})+b(q_{m})z)}\right]\\\
&-\frac{N(N-1)}{2\pi[b(q_{m})b(q_{n})]^{N-1}}\Bigg{[}\int\limits_{\mathbb{R}^{2}}d[x]r(x_{1},v(q_{m}))r(x_{2},v(q_{n}))(x_{2}-x_{1})(x_{1}x_{2}+1)^{N-2}\\\
&+\int\limits_{\mathbb{R}}dx\int\limits_{\mathbb{C}}d[z]\det\begin{bmatrix}r(x,v(q_{m}))&s(z,z^{*},v(q_{m}))\\\
r(x,v(q_{n}))&s(z,z^{*},v(q_{n}))\end{bmatrix}(z-x)(zx+1)^{N-2}\\\
&+\int\limits_{\mathbb{C}^{2}}d[z]s(z_{1},z_{1}^{*},v(q_{m}))s(z_{2},z_{2}^{*},v(q_{n}))(z_{2}-z_{1})(z_{1}z_{2}+1)^{N-2}\Bigg{]},\end{split}$
(70)

which contains six integrals over at most four real variables.

##  V Conclusions

The winding number is the pertinent topological invariant for chiral symmetric
Bloch Hamiltonians over a one-dimensional parameter manifold. We studied its
statistics in the chiral orthogonal class BDI by calculating the ensemble
averages for ratios of parametric determinants. This is the key quantity for
the winding number statistics as it generates the correlation functions of the
winding number densities. We did this by first mapping our problem to a
spectral one, namely finding the ratios of characteristic polynomials in the
spherical ensemble. Then we applied the supersymmetry without supersymmetry
technique to write the ensemble average as a Pfaffian of three kernel
functions, that are simplified averages over only two determinants. We could
simplify the kernel function even further by employing a symmetry of our
random matrix model. This complements our recent studies on the chiral unitary
and chiral symplectic class [17, 18].

For technical reasons we considered only the case of even matrix dimensions.
This does not interfere with our goal to study the universal behaviour in the
limit of large matrix dimensions as the even and odd case behave equally in
this limit. Universality is reached in the limit of large matrix dimension
when also considering the parametric dependence on proper scales. In Ref. 17
we additionally addressed this question for the correlation functions of
winding number densities in the chiral unitary class AIII. Obtaining the
universal limits for the correlation functions involves taking the derivatives
of the here computed ensemble averages. Although the resulting expressions for
the ensemble averages are quite cumbersome, we are positive that in the
universality limit the correlation functions are feasible. We want to address
this issue further in a future work.

###### Acknowledgements.

We thank Boris Gutkin for fruitful discussions. This work was funded by the
German–Israeli Foundation within the project Statistical Topology of Complex
Quantum Systems, grant number GIF I-1499-303.7/2019 (N.H.,O.G. and T.G.).
Furthermore, M.K. acknowledges support by the Australian Research Council via
Discovery Project grant DP210102887.

##  Appendix A The joint eigenvalue probability density of the real induced
spherical ensemble

We consider the real induced spherical ensemble with the matrix probability
density

$\widetilde{G}^{(1,N)}_{\mu,\nu}(Y)=\pi^{-N^{2}/2}\prod_{j=1}^{N}\frac{\Gamma(j/2)\Gamma(\mu+\nu+(N+j)/2)}{\Gamma(\nu+j/2)\Gamma(\mu+j/2)}\frac{\det^{2\nu}Y}{\det^{N+\mu+\nu}\left(\mathds{1}_{N}+YY^{T}\right)}.$
(71)

It is the ensemble of matrices $Y=K_{1}^{-1}K_{2}$, where $K_{1}$ and $K_{2}$
are distributed according to deformed Gaussians

$P_{\mu}(K)=2^{-N\mu}(2\pi)^{-N^{2}/2}\prod_{j=1}^{N}\frac{\Gamma(j/2)}{\Gamma(\mu+j/2)}\exp\left(-\frac{1}{2}\text{tr}\,\,KK^{T}\right)\det\nolimits^{\mu}KK^{T}$
(72)

for the parameters $\mu$ resp. $\nu$. In the body of the text we omitted these
indices, implying that they are zero, e.g.
$\widetilde{G}^{(1,N)}_{0,0}(Y)=\widetilde{G}^{(1,N)}(Y)$, which leads to the
ordinary spherical ensemble. The induced spherical ensemble is well-studied
for $\beta=1,2,4$ [40, 41, 42]. Unfortunately, in the real case the known
joint eigenvalue probability density is available only in a form, that is
impractical for our purpose. For this reason we will reproduce this result
once again. In doing so, we adhere to methods employed in the aforementioned
works.

Generally for real matrices, we have to distinguish between even
($\widetilde{N}=N/2$) and odd ($\widetilde{N}=(N-1)/2$) matrix dimensions. We
apply a real Schur decomposition

$Y=U(D+T)U^{T},$ (73)

where $D$ is a diagonal matrix of $\widetilde{N}$ real $2\times 2$-blocks
$D_{j}$. In the odd case these are complemented by a real $1\times 1$-block
$D_{\widetilde{N}+1}$. The matrix $T$ is a real strict upper triangular matrix
and $U$ is orthogonal, $U\in\text{O}(N)/\text{O}(2)^{\widetilde{N}}$ resp.
$U\in\text{O}(N)/(\text{O}(2)^{\widetilde{N}}\times\text{O}(1))$ for even
resp. odd $N$. Upon this transformation the measure changes according to [44,
45]

$d[Y]=\Delta_{N}(z)\prod_{j=1}^{\widetilde{N}}\frac{1}{z_{-j}-z_{+j}}d[D]d[T]d\mu(U),$
(74)

where $z_{+j}$ and $z_{-j}$ are the eigenvalues of $D_{j}$ and thus also of
$Y$. First, one integrates over the upper triangular matrix $T$. Only the
denominator of the matrix density (71) depends on $T$. The determinant is
expanded and the free parameters are integrated out column wise. In the even
case this yields [35, 42]

$\begin{split}\int
d[T]\frac{\det^{2\nu}Y}{\det^{N+\mu+\nu}\left(\mathds{1}_{N}+YY^{T}\right)}=&\prod_{j=1}^{\widetilde{N}}\pi^{N-2j}\frac{\Gamma\left(N/2+\mu+\nu+1/2\right)\Gamma\left(N/2+\mu+\nu+1\right)}{\Gamma\left(N+\mu+\nu+1/2-j\right)\Gamma\left(N+\mu+\nu+1-j\right)}\\\
&\times\prod_{j=1}^{\widetilde{N}}\frac{\det^{2\nu}D_{j}}{\det^{N/2+\mu+\nu+1}(\mathds{1}_{N}+D_{j}D_{j}^{T})}\end{split}$
(75)

and in the odd case

$\begin{split}\int
d[T]\frac{\det^{2\nu}Y}{\det^{N+\mu+\nu}\left(\mathds{1}_{N}+YY^{T}\right)}=&\prod_{j=1}^{\widetilde{N}}\pi^{N-2j}\frac{\Gamma\left(N/2+\mu+\nu+1/2\right)\Gamma\left(N/2+\mu+\nu+1\right)}{\Gamma\left(N+\mu+\nu+1/2-j\right)\Gamma\left(N+\mu+\nu+1-j\right)}\\\
&\times
h_{\mu,\nu}^{(N)}(z_{2\widetilde{N}+1})\prod_{j=1}^{\widetilde{N}}\frac{\det^{2\nu}D_{j}}{\det^{N/2+\mu+\nu+1}(\mathds{1}_{N}+D_{j}D_{j}^{T})}.\end{split}$
(76)

Here, we obtain an additional factor

$h_{\mu,\nu}^{(N)}(z)=\frac{z^{2\nu}}{\left(1+z^{2}\right)^{N/2+\mu+\nu+1/2}}\delta\left(y\right)$
(77)

which effectively renders the variable $z_{2\widetilde{N}+1}$ to be real. We
use the notation $z=x+iy$ with $x,y\in\mathbb{R}$ for the generally complex
eigenvalues. Next, we work on the $2\times 2$-blocks. Gathering all factors
depending on the eigenvalues of $D_{j}$, except those in the Vandermonde
determinant, we obtain

$A_{j}=\frac{1}{z_{-j}-z_{+j}}\frac{\det^{2\nu}D_{j}}{\det^{N/2+\mu+\nu+1}(\mathds{1}_{N}+D_{j}D_{j}^{T})}.$
(78)

Following Ref. 46, we start by diagonalizing the symmetric part

$O_{j}^{T}D_{j}O_{J}=\begin{bmatrix}\lambda_{1j}&\rho_{j}\\\
-\rho_{j}&\lambda_{2j}\end{bmatrix}\qquad
O_{j}=\begin{bmatrix}\cos\varphi_{j}&\sin\varphi_{j}\\\
-\sin\varphi_{j}&\cos\varphi_{j}\end{bmatrix}\in\text{SO}(2).$ (79)

The range of parameters is $\rho_{j}\in\mathbb{R}$, $\varphi_{j}\in[0,\pi)$
and $\lambda_{1j},\lambda_{2j}\in\mathbb{R}$ with
$\lambda_{1j}\geq\lambda_{2j}$. The flat measure of the four real independent
variables transforms as

$d[D_{j}]=2(\lambda_{1j}-\lambda_{2j})d\varphi_{j}d\rho_{j}d\lambda_{1j}d\lambda_{2j}.$
(80)

Next, we want to change coordinates from the eigenvalues of the symmetric part
$\lambda_{1j},\lambda_{2j}$ to the eigenvalues $z_{\pm j}$ of the full matrix.
They are related via

$z_{\pm
j}=\frac{\lambda_{1j}+\lambda_{2j}}{2}\pm\sqrt{\left(\frac{\lambda_{1j}-\lambda_{2j}}{2}\right)^{2}-\rho_{j}^{2}}\qquad\lambda_{1,2j}=\frac{z_{+j}+z_{-j}}{2}\pm\sqrt{\left(\frac{z_{+j}-z_{-j}}{2}\right)^{2}+\rho_{j}^{2}}.$
(81)

According to the ordering of $\lambda_{1}$ and $\lambda_{2}$, the eigenvalues
$z_{\pm j}$ also have to be ordered

$\begin{split}z_{\pm j}\in\mathbb{R}&:x_{+j}\geq x_{-j}\\\ z_{\pm
j}=z^{*}_{\mp j}&:y_{+j}=-y_{-j}>0,\end{split}$ (82)

where we distinguished between the cases of real and complex conjugated
eigenvalues. However, in the following we will disregard this ordering, which
can be compensated by an overall factor $1/2$. We obtain the Jacobian

$\left|\det\frac{\partial(z_{+j},z_{-j})}{\partial(\lambda_{1j},\lambda_{2j})}\right|=\frac{\lambda_{1j}-\lambda_{2j}}{\left|z_{-j}-z_{+j}\right|}.$
(83)

The determinants in the new coordinates are

$\det\left(\mathds{1}_{2}+D_{j}D_{j}^{T}\right)=4\rho_{j}^{2}+(1+z_{+j}^{2})(1+z_{-j}^{2})\qquad\det\,D_{j}=z_{+j}z_{-j}.$
(84)

We also need to consider that in the case, where $z_{\pm j}$ are complex
conjugates, the integration regime of $\rho_{j}$ is confined to

$\rho_{j}^{2}\geq-\left(\frac{z_{+j}-z_{-j}}{2}\right)^{2}=y^{2}_{\pm
j}\Rightarrow\left|\rho_{j}\right|\geq\left|y_{\pm j}\right|.$ (85)

This is due to the condition $\lambda_{1,2j}\in\mathbb{R}$. Gathering
everything and integrating over $\varphi_{j}$ and $\rho_{j}$ we obtain for
(78)

$\begin{split}\phantom{=}&2(\lambda_{1j}-\lambda_{2j})d\lambda_{1j}d\lambda_{2j}\int\limits_{\mathbb{R}}d\rho_{j}\int\limits_{0}^{\pi}d\varphi_{j}\,A_{j}\\\
&=\frac{\pi}{2}\frac{\left(z_{+j}z_{-j}\right)^{2\nu}}{\left[(1+z_{+}^{2})(1+z_{-}^{2})\right]^{N/2+\nu+\mu+1/2}}\frac{\left|z_{-j}-z_{+j}\right|}{z_{-j}-z_{+j}}\\\
&\times\bigg{[}B(1/2,N/2+\mu+\nu+1/2)\Theta(x_{+j}-x_{-j})\delta(y_{-j})\delta(y_{+j})\\\
&+2\,\Theta(y_{+j})\delta(x_{+j}-x_{-j})\delta(y_{+j}+y_{-j})Q_{\mu,\nu}^{(N)}(z_{+j},z_{-j})\bigg{]}dx_{+j}dx_{-j}dy_{+j}dy_{-j},\end{split}$
(86)

where the function

$\begin{split}Q_{\mu,\nu}^{(N)}(z_{+j},z_{-j})&=2\int\limits_{\frac{2\left|y_{\pm
j}\right|}{\left|1+z_{\pm
j}^{2}\right|}}^{\infty}\frac{d\rho_{j}}{\left(1+\rho_{j}^{2}\right)^{N/2+\mu+\nu+1}}\\\
&=B\left(\frac{4y^{2}}{\left|1+z^{2}\right|^{2}+4y^{2}};1/2,N/2+\mu+\nu+1/2\right)\end{split}$
(87)

emerges. Now only the integral over the Haar measure of the respective cosets
remains to be computed. It yields the constants

$\begin{split}\text{Vol}\left(\text{O}(N)/\text{O}(2)^{\widetilde{N}}\right)&=\frac{\text{Vol}(\text{O}(N))}{(4\pi)^{N/2}}=\frac{1}{(4\pi)^{N/2}}\prod_{j=1}^{N}\frac{2\pi^{j/2}}{\Gamma(j/2)},\\\
\text{Vol}\left(\text{O}(N)/(\text{O}(2)^{\widetilde{N}}\times\text{O}(1))\right)&=\frac{\text{Vol}(\text{O}(N))}{2(4\pi)^{(N-1)/2}}=\frac{1}{2(4\pi)^{(N-1)/2}}\prod_{j=1}^{N}\frac{2\pi^{j/2}}{\Gamma(j/2)}\end{split}$
(88)

for $N$ even resp. odd. We summarize this under

$\int
d\mu(U)=\frac{1}{(1+N-2\widetilde{N})(4\pi)^{\widetilde{N}}}\prod_{j=1}^{N}\frac{2\pi^{j/2}}{\Gamma(j/2)}.$
(89)

As we not only disregard the ordering of the eigenvalues $z_{\pm j}$ in each
$2\times 2$-block, but also the ordering of the $2\times 2$-blocks themselves,
all of this has to be supplemented with a factor $1/\widetilde{N}!$.

Finally bringing everything together we obtain the correctly normalized joint
eigenvalue density

$\begin{split}G^{(1,N)}_{\mu,\nu}(z)&=C_{\mu,\nu}^{(N)}\Delta_{N}(z)\prod_{j=1}^{\widetilde{N}}g_{\mu,\nu}^{(1,N)}(z_{2j-1},z_{2j}),\\\
G^{(1,N)}_{\mu,\nu}(z)&=C_{\mu,\nu}^{(N)}\Delta_{N}(z)h_{\mu,\nu}^{(N)}(z_{N})\prod_{j=1}^{\widetilde{N}}g_{\mu,\nu}^{(1,N)}(z_{2j-1},z_{2j})\end{split}$
(90)

for $N$ even resp. odd with the antisymmetric function

$\begin{split}g_{\mu,\nu}^{(1,N)}(z_{1},z_{2})=&\left(z_{1}z_{2}\right)^{2\nu}\frac{\left|z_{2}-z_{1}\right|}{z_{2}-z_{1}}\\\
&\times\frac{B(1/2,N/2+\mu+\nu+1/2)\delta(y_{1})\delta(y_{2})+2\,\delta(x_{1}-x_{2})\delta(y_{1}+y_{2})Q_{\mu,\nu}^{(N)}(z_{1},z_{1}^{*})}{\left[(1+z_{1}^{2})(1+z_{2}^{2})\right]^{N/2+\mu+\nu+1/2}}.\end{split}$
(91)

The antisymmetry of this function is due to us disregarding an ordering (82).
The normalization is

$\begin{split}C_{\mu,\nu}^{(N)}=&\frac{2^{N-3\widetilde{N}}\pi^{-\widetilde{N}^{2}+\widetilde{N}(N-1)-N(N-1)/4}}{(1+N-2\widetilde{N})\widetilde{N}!}\prod_{j=1}^{N}\frac{\Gamma(\mu+\nu+(N+j)/2)}{\Gamma(\nu+j/2)\Gamma(\mu+j/2)}\\\
&\times\prod_{j=1}^{\widetilde{N}}\frac{\Gamma\left(N/2+\mu+\nu+1/2\right)\Gamma\left(N/2+\mu+\nu+1\right)}{\Gamma\left(N+\mu+\nu+1/2-j\right)\Gamma\left(N+\mu+\nu+1-j\right)}.\end{split}$
(92)

Setting $\mu=\nu=0$ one immediately finds the joint probability density of
eigenvalues (29) of the ordinary real spherical ensemble.

The Pfaffian of the moment matrix (33) can be related to this constant. One
finds for the even case [46]

$\int\limits_{\mathbb{C}^{2M}}d[z]\Delta_{2M}(z)\prod_{j=1}^{M}g^{(1,N)}(z_{2j-1},z_{2j})=M!\,{\text{Pf}}\,D^{(2M)},$
(93)

where $M\leq N/2$. On the other hand by equations (90) and (91) this integral
is

$\int\limits_{\mathbb{C}^{2M}}d[z]\Delta_{2M}(z)\prod_{j=1}^{M}g_{(N-2M)/2,0}^{(1,2M)}(z_{2j-1},z_{2j})=\frac{1}{C^{(2M)}_{(N-2M)/2,0}}$
(94)

and therefore

${\text{Pf}}\,D^{(2M)}=\frac{1}{M!}\frac{1}{C^{(2M)}_{(N-2M)/2,0}}.$ (95)

We use this result at multiple points in the body of the text to evaluate the
prefactors of our integrals.

##  Appendix B Characteristic polynomials of the real induced spherical
ensemble

We are interested in integrals of the type

$I_{M}=\int\limits_{\mathbb{C}^{2M}}d[z]\Delta_{2M}(z)\prod_{j=1}^{M}g^{(1,N)}(z_{2j-1},z_{2j})\prod_{j=1}^{2M}(x-z_{j}),$
(96)

where $M\leq N/2$. Applying equations (90) we may map them to integrals over a
matrix density

$\begin{split}I_{M}=&\int\limits_{\mathbb{C}^{2M}}d[z]\Delta_{2M}(z)\prod_{j=1}^{M}g_{(N-2M)/2,0}^{(1,2M)}(z_{2j-1},z_{2j})\prod_{j=1}^{2M}(x-z_{j})\\\
=&\frac{1}{C^{(2M)}_{(N-2M)/2,0}}\int
d[Y]\widetilde{G}_{(N-2M)/2,0}^{(1,2M)}(Y)\det(x-Y).\end{split}$ (97)

The matrix density $\widetilde{G}^{(1,N)}_{\mu,\nu}(Y)$ is invariant under
left and right actions of the orthogonal group, see (71),

$Y\to O_{1}YO_{2}\qquad\text{with }O_{1},O_{2}\in\text{O}(N).$ (98)

Let us choose $O_{1}=\mathds{1}_{N}$ and $O_{2}$ and as diagonal with entries
$\pm 1$. This effectively changes the sign of $Y$ column-wise, $Y_{jl}\to\pm
Y_{jl}$ for all $l$. Expanding the determinant we see that only one term
survives the average

$\begin{split}\left\langle\det\left(x-Y\right)\right\rangle=\sum_{\sigma\in\mathbb{S}_{N}}\text{sgn}\,\sigma\left\langle\prod_{j=1}^{N}\left(x\delta_{j\sigma(j)}-Y_{j\sigma(j)}\right)\right\rangle=\sum_{\sigma\in\mathbb{S}_{N}}\text{sgn}\,\sigma\,x^{N}\prod_{j=1}^{N}\delta_{j\sigma(j)}=x^{N},\end{split}$
(99)

which is a monomial of power $N$. Therefore we obtain

$I_{M}=\frac{x^{2M}}{C^{(2M)}_{(N-2M)/2,0}}.$ (100)

This result corresponds to the skew-orthogonal polynomial of even degree [18,
47].

##  Appendix C Alternative expression of ${\rm K}_{3}^{(1)}$

Unlike in (64) we expand the determinant only in the last column

$\Xi_{3}=\frac{2\,C^{(N)}}{\kappa_{2}-\kappa_{1}}\int\limits_{\mathbb{C}^{N+2}}d[z]\frac{1}{z_{N+2}+\kappa_{2}}\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})\det\left[\begin{array}[]{c|c}\displaystyle
z_{a}^{b-1}&\displaystyle\frac{1}{z_{a}+\kappa_{1}}\end{array}\right]_{\begin{subarray}{c}1\leq
a\leq N+1\\\ 1\leq b\leq N\end{subarray}}$ (101)

and apply the identity (see (49))

$\begin{split}\det\left[\begin{array}[]{c|c}\displaystyle
z_{a}^{b-1}&\displaystyle\frac{1}{z_{a}+\kappa_{1}}\end{array}\right]_{\begin{subarray}{c}1\leq
a\leq N+1\\\ 1\leq b\leq
N\end{subarray}}&=(-1)^{N+2}\Delta_{N+1}(z)\prod_{j=1}^{N+1}\frac{1}{z_{j}+\kappa_{1}}\\\
&=(-1)^{N+2}\frac{\Delta_{N}(z)}{z_{N+1}+\kappa_{1}}\prod_{j=1}^{N}\frac{z_{N+1}-z_{j}}{z_{j}+\kappa_{1}}.\end{split}$
(102)

This allows us to identify the integral over $z_{1},\ldots,z_{N}$ as the
function $\Xi_{2}$, see Eq. (45), that we calculated for the second kernel,
resulting in

$\begin{split}\Xi_{3}&=\frac{2\,C^{(N)}}{\kappa_{2}-\kappa_{1}}\int\limits_{\mathbb{C}^{N+2}}d[z]\frac{\Delta_{N}(z)}{(z_{N+1}+\kappa_{1})(z_{N+2}+\kappa_{2})}\prod_{j=1}^{N}\frac{z_{N+1}-z_{j}}{z_{j}+\kappa_{1}}\prod_{j=1}^{N/2+1}g^{(1,N)}(z_{2j-1},z_{2j})\\\
&=\frac{2}{\kappa_{2}-\kappa_{1}}\int\limits_{\mathbb{C}^{2}}d[z]\frac{g^{(1,N)}(z_{1},z_{2})}{(z_{1}+\kappa_{1})(z_{2}+\kappa_{2})}\left\langle\frac{\det\left(-z_{1}K_{1}+K_{2}\right)}{\det\left(\kappa_{1}K_{1}+K_{2}\right)}\right\rangle.\end{split}$
(103)

Inserting our result (51) for the ensemble average we find the following
expression for the third kernel

$\begin{split}{\rm
K}_{3}^{(1)}(q_{m},q_{n})=&\frac{2}{[b(q_{m})b(q_{n})]^{N-1}}\int\limits_{\mathbb{C}^{2}}d[z]\frac{g^{(1,N)}(z_{1},z_{2})}{(a(q_{m})+b(q_{m})z_{1})(a(q_{n})+b(q_{n})z_{2})}\left\langle\frac{\det(-z_{1}K_{1}+K_{2})}{\det(\kappa(q_{m})K_{1}+K_{2})}\right\rangle\\\
=&\frac{-N(N-1)}{2\pi[b(q_{m})b(q_{n})]^{N-1}}\int\limits_{\mathbb{C}^{2}}d[z]\frac{g^{(1,N)}(z_{1},z_{2})}{(a(q_{m})+b(q_{m})z_{1})(a(q_{n})+b(q_{n})z_{2})}\\\
&\times\Bigg{[}\frac{(-1)^{N/2}2\pi B(1/2,(N+1)/2)}{N-1}{(N-1)/2\choose N}\\\
&\times{{}_{2}F_{1}}\left(1,(N+1)/2;N+1;1+\left(\frac{\kappa(q_{m})-z_{1}}{\kappa(q_{m})+z_{1}}\right)^{2}\right)\\\
&+2i\int\limits_{\mathbb{C}}d[z]\frac{z^{N-2}\text{sgn}(\text{Im}\,z)Q(z,z^{*})}{\left|1+z^{2}\right|^{N+1}}\left(z^{*}+\frac{\kappa(q_{m})-z_{1}}{\kappa(q_{m})+z_{1}}\right)^{-1}\Bigg{]}.\end{split}$
(104)

Furthermore inserting (18) yields

$\begin{split}{\rm K}_{3}^{(1)}(q_{m},q_{n})=&\frac{-N(N-1)b(q_{m})}{2\pi
b^{N-1}(q_{n})}\Bigg{(}\int\limits_{\mathbb{R}}dx\frac{r(x,v(q_{n}))}{a(q_{m})+b(q_{m})x}\left(\frac{x^{2}+1}{a(q_{m})+b(q_{m})x}\right)^{N}\\\
&\times\Bigg{[}\frac{(-1)^{N/2}2\pi B(1/2,(N+1)/2)}{N-1}{(N-1)/2\choose N}\\\
&\times{{}_{2}F_{1}}\left(1,(N+1)/2;N+1;1+\left(\frac{a(q_{m})-b(q_{m})x}{a(q_{m})+b(q_{m})x}\right)^{2}\right)\\\
&+\int\limits_{\mathbb{C}}d[z]s\left(z,z^{*},\left(\frac{a(q_{m})-b(q_{m})x}{a(q_{m})+b(q_{m})x}\right)\right)z^{N-2}\Bigg{]}\\\
&+\int\limits_{\mathbb{C}}d[z_{1}]\frac{s(z_{1},z_{1}^{*},v(q_{n}))}{a(q_{m})+b(q_{m})z_{1}}\left(\frac{z_{1}^{2}+1}{a(q_{m})+b(q_{m})z_{1}}\right)^{N}\\\
&\times\Bigg{[}\frac{(-1)^{N/2}2\pi B(1/2,(N+1)/2)}{N-1}{(N-1)/2\choose N}\\\
&\times{{}_{2}F_{1}}\left(1,(N+1)/2;N+1;1+\left(\frac{a(q_{m})-b(q_{m})z_{1}}{a(q_{m})+b(q_{m})z_{1}}\right)^{2}\right)\\\
&+\int\limits_{\mathbb{C}}d[z_{2}]s\left(z_{2},z_{2}^{*},\left(\frac{a(q_{m})-b(q_{m})z_{1}}{a(q_{m})+b(q_{m})z_{1}}\right)\right)z_{2}^{N-2}\Bigg{]}\Bigg{)}\end{split}$
(105)

with the functions (22). In this representation also a dependence solely in
O(2)-invariants cannot be identified.

## References

  * Bravyi and Hastings [2011] S. Bravyi and M. B. Hastings, “A short proof of stability of topological order under local perturbations,” Commun. Math. Phys. 307, 609 (2011), arXiv:1001.4363 .
  * Oliviero  _et al._ [2022] S. F. E. Oliviero, L. Leone, Y. Zhou, and A. Hamma, “Stability of topological purity under random local unitaries,” SciPost Phys. 12, 096 (2022), arXiv:2106.04600 .
  * Clark  _et al._ [2022] O. J. Clark, F. Freyse, L. V. Yashina, O. Rader, and J. Sánchez-Barriga, “Robust behavior and spin-texture stability of the topological surface state in Bi2Se3 upon deposition of gold,” npj Quantum Materials 7, 36 (2022).
  * Chaunsali  _et al._ [2021] R. Chaunsali, H. Xu, J. Yang, P. G. Kevrekidis, and G. Theocharis, “Stability of topological edge states under strong nonlinear effects,” Phys. Rev. B 103, 024106 (2021), arXiv:2010.13142 .
  * Bisianov  _et al._ [2019] A. Bisianov, M. Wimmer, U. Peschel, and O. A. Egorov, “Stability of topologically protected edge states in nonlinear fiber loops,” Phys. Rev. A 100, 063830 (2019).
  * Prodan and Schulz-Baldes [2016] E. Prodan and H. Schulz-Baldes,  _Bulk and Boundary Invariants for Complex Topological Insulators: From K-Theory to Physics_, Mathematical Physics Studies (Springer International Publishing, Cham, 2016).
  * Chen and Chiou [2020] B.-H. Chen and D.-W. Chiou, “An elementary rigorous proof of bulk-boundary correspondence in the generalized Su-Schrieffer-Heeger model,” Phys. Lett. A 384, 126168 (2020).
  * Shapiro [2020] J. Shapiro, “The bulk-edge correspondence in three simple cases,” Rev. Math. Phys. 32, 2030003 (2020), https://doi.org/10.1142/S0129055X20300034 .
  * Mehta [2004] M. L. Mehta,  _Random Matrices_ (Academic Press, 2004).
  * Guhr, Müller-Groeling, and Weidenmüller [1998] T. Guhr, A. Müller-Groeling, and H. A. Weidenmüller, “Random-matrix theories in quantum physics: common concepts,” Phys. Rep. 299, 189–425 (1998).
  * Simons and Altshuler [1993a] B. D. Simons and B. L. Altshuler, “Universal velocity correlations in disordered and chaotic systems,” Phys. Rev. Lett. 70, 4063–4066 (1993a).
  * Simons and Altshuler [1993b] B. D. Simons and B. L. Altshuler, “Universalities in the spectra of disordered and chaotic systems,” Phys. Rev. B 48, 5422–5438 (1993b).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, “Classification of topological quantum matter with symmetries,” Rev. Mod. Phys. 88, 035005 (2016).
  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, “Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures,” Phys. Rev. B 55, 1142–1161 (1997).
  * Heinzner, Huckleberry, and Zirnbauer [2005] P. Heinzner, A. Huckleberry, and M. R. Zirnbauer, “Symmetry classes of disordered fermions,” Commun. Math. Phys. 257, 725–771 (2005).
  * Kitaev [2009] A. Kitaev, “Periodic table for topological insulators and superconductors,” AIP Conf. Proc. 1134, 22–30 (2009).
  * Braun  _et al._ [2022] P. Braun, N. Hahn, D. Waltner, O. Gat, and T. Guhr, “Winding number statistics of a parametric chiral unitary random matrix ensemble,” J. Phys. A: Math. Theor. 55, 224011 (2022).
  * Hahn  _et al._ [2023] N. Hahn, M. Kieburg, O. Gat, and T. Guhr, “Winding number statistics for chiral random matrices: Averaging ratios of determinants with parametric dependence,” Journal of Mathematical Physics 64 (2023), 10.1063/5.0112423, 021901, https://pubs.aip.org/aip/jmp/article-pdf/doi/10.1063/5.0112423/16724006/021901_1_online.pdf .
  * Fyodorov and Strahov [2003] Y. V. Fyodorov and E. Strahov, “An exact formula for general spectral correlation function of random hermitian matrices,” J. Phys. A: Math. Gen. 36, 3203–3213 (2003).
  * Borodin and Strahov [2006] A. Borodin and E. Strahov, “Averages of characteristic polynomials in random matrix theory,” Commun. Pure Appl. Math. 59, 161 (2006).
  * Kieburg and Guhr [2010a] M. Kieburg and T. Guhr, “Derivation of determinantal structures for random matrix ensembles in a new way,” J. Phys. A: Math. Theor. 43, 075201 (2010a).
  * Kieburg and Guhr [2010b] M. Kieburg and T. Guhr, “A new approach to derive pfaffian structures for random matrix ensembles,” J. Phys. A: Math. Theor. 43, 135204 (2010b).
  * Akemann, Strahov, and Würfel [2020] G. Akemann, E. Strahov, and T. R. Würfel, “Averages of products and ratios of characteristic polynomials in polynomial ensembles,” Ann. Henri Poincare 21, 3973–4002 (2020), arXiv:2003.08128 .
  * Ipsen and Forrester [2018] J. R. Ipsen and P. J. Forrester, “Kac-rice fixed point analysis for single- and multi-layered complex systems,” J. Phys. A: Math. Theor. 51, 474003 (2018), arXiv:1807.05790 .
  * Webb [2015] C. Webb, “The characteristic polynomial of a random unitary matrix and gaussian multiplicative chaos - the $l^{2}$-phase,” Electron. J. Probab. 20, 1–21 (2015), arXiv:1410.0939 .
  * Mehta and Normand [2001] M. L. Mehta and J.-M. Normand, “Complexity of random energy landscapes, glass transition and absolute value of spectral determinant of random matrices,” J. Phys. A: Math. Gen. 34, 4627 (2001), arXiv:cond-mat/0101469 .
  * Conrey  _et al._ [2003] J. B. Conrey, D. W. Farmer, J. P. Keating, M. O. Rubinstein, and N. C. Snaith, “Autocorrelation of random matrix polynomials,” Commun. Math. Phys. 237, 365–395 (2003), arXiv:math-ph/0208007 .
  * Bourgade  _et al._ [2008] P. Bourgade, C. P. Hughes, A. Nikeghbali, and M. Yor, “The characteristic polynomial of a random unitary matrix: A probabilistic approach,” Duke Math. J. 145, 45 – 69 (2008).
  * Eberhard [2022] S. Eberhard, “The characteristic polynomial of a random matrix,” Combinatorica 42, 491 (2022), arXiv:2008.01223 .
  * Fyodorov [2004] Y. V. Fyodorov, “Complexity of random energy landscapes, glass transition and absolute value of spectral determinant of random matrices,” Phys. Rev. Lett. 92, 240601 (2004), arXiv:cond-mat/0401287 .
  * Efetov [1983] K. B. Efetov, “Supersymmetry and theory of disordered metals,” Adv. Phys. 32, 53–127 (1983).
  * Ginibre [1965] J. Ginibre, “Statistical ensembles of complex, quaternion, and real matrices,” J. Math. Phys. 6, 440–449 (1965).
  * Maffei  _et al._ [2018] M. Maffei, A. Dauphin, F. Cardano, M. Lewenstein, and P. Massignan, “Topological characterization of chiral models through their long time dynamics,” New J. Phys. 20, 013023 (2018).
  * Asbóth, Oroszlány, and Pályi [2016] J. K. Asbóth, L. Oroszlány, and A. Pályi,  _A Short Course on Topological Insulators_ (Springer International Publishing, 2016).
  * Forrester and Mays [2012] P. J. Forrester and A. Mays, “Pfaffian point process for the Gaussian real generalised eigenvalue problem,” Probab. Theory Relat. Fields 154, 1–47 (2012).
  * Forrester and Krishnapur [2009] P. J. Forrester and M. Krishnapur, “Derivation of an eigenvalue probability density function relating to the poincaré disk,” J. Phys. A: Math. Theor. 42, 385204 (2009).
  * Mays [2013] A. Mays, “A real quaternion spherical ensemble of random matrices,” J. Stat. Phys. 153, 48–69 (2013).
  * Bruijn, de [1955] N. Bruijn, de, “On some multiple integrals involving determinants,” J. Indian Math. Soc. 19, 133–151 (1955).
  * Guhr [1991] T. Guhr, “Dyson’s correlation functions and graded symmetry,” J. Math. Phys. 32, 336–347 (1991).
  * Fischmann and Forrester [2011] J. Fischmann and P. J. Forrester, “One-component plasma on a spherical annulus and a random matrix ensemble,” Journal of Statistical Mechanics: Theory and Experiment 2011, P10003 (2011).
  * Mays and Ponsaing [2017] A. Mays and A. Ponsaing, “An induced real quaternion spherical ensemble of random matrices,” Random Matrices: Theory and Applications 06, 1750001 (2017), https://doi.org/10.1142/S2010326317500010 .
  * Fischmann [2015] J. A. Fischmann,  _Eigenvalue distributions on a single ring_ , Ph.D. thesis, Queen Mary University of London (2015).
  * [43] DLMF, “NIST Digital Library of Mathematical Functions,” http://dlmf.nist.gov/, Release 1.1.3 of 2021-09-15, edited by F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain.
  * Ipsen and Kieburg [2014] J. R. Ipsen and M. Kieburg, “Weak commutation relations and eigenvalue statistics for products of rectangular random matrices,” Phys. Rev. E 89, 032106 (2014).
  * Edelman [1997] A. Edelman, “The probability that a random real gaussian matrix haskreal eigenvalues, related distributions, and the circular law,” Journal of Multivariate Analysis 60, 203–232 (1997).
  * Sommers and Wieczorek [2008] H.-J. Sommers and W. Wieczorek, “General eigenvalue correlations for the real ginibre ensemble,” Journal of Physics A: Mathematical and Theoretical 41, 405003 (2008).
  * Akemann, Kieburg, and Phillips [2010] G. Akemann, M. Kieburg, and M. J. Phillips, “Skew-orthogonal laguerre polynomials for chiral real asymmetric random matrices,” J. Phys. A: Math. Gen. 43, 375207 (2010), arXiv:1005.2983 .