# Moments of quantum purity and biorthogonal polynomial recurrence

Shi-Hao Li and Lu Wei  Shi-Hao Li, Department of Mathematics, Sichuan
University, Chengdu, 610064, China<EMAIL_ADDRESS><EMAIL_ADDRESS>Lu Wei, Department of Electrical and Computer
Engineering, University of Michigan, Dearborn, MI 48128, USA<EMAIL_ADDRESS>
###### Abstract.

The Bures-Hall ensemble is a unique measure of density matrices that satisfies
various distinguished properties in quantum information processing. In this
work, we study the statistical behavior of entanglement over the Bures-Hall
ensemble as measured by the simplest form of an entanglement entropy - the
quantum purity. The main results of this work are the exact second and third
moment expressions of quantum purity valid for any subsystem dimensions, where
the corresponding results in the literature are limited to the scenario of
equal subsystem dimensions. In obtaining the results, we have derived
recurrence relations of the underlying integrals over the Cauchy-Laguerre
biorthogonal polynomials that may be of independent interest.

###### Key words and phrases:

entanglement entropy, quantum purity, orthogonal polynomials, recurrence
relation

##  1\. Introduction and main results

Quantum information theory aims at studying the theoretical foundations of
quantum technologies including quantum computing and quantum communications.
Crucial to successful exploitation of the quantum revolutionary advances is
the understanding of the phenomenon of quantum entanglement. Entanglement is
the most fundamental characteristic trait of quantum mechanics, which is also
the resource that enables quantum technologies.

In this work, we study the statistical behavior of entanglement of quantum
bipartite systems over the Bures-Hall ensemble [1, 2, 3]. In particular, we
investigate the degree of entanglement as measured by quantum purity over such
an ensemble. Quantum purity measures how far a state is from a pure state and
is the simplest form of an entanglement entropy. We focus on finding the exact
moments of purity, which in practice can be utilized to construct finite-size
approximations to the distribution [4, 5]. The considered bipartite model is
useful in describing the entanglement between the two subsystems of various
real-world quantum systems, in which one subsystem represents a physical
object (such as a set of spins) and the other subsystem is the environment
(such as a heat bath). Existing results on the moments of purity are limited
to the special case of equal subsystem dimensions, where the exact first three
moments are known [6, 7]. In the general case of subsystems of possibly
unequal dimensions, the exact first moment formula has been recently obtained
in [8, 9]. The contribution of this work lies in the corresponding exact
second and third moment expressions of quantum purity.

The density matrix formalism [10] introduced by von Neumann provides a natural
framework to describe density matrices of quantum states. Under this
framework, the formulation that has led to the Bures-Hall ensemble [1, 2, 3]
is outlined as follows. Consider a bipartite system that consists of two
subsystems $A$ and $B$ of Hilbert space (complex vector space) dimensions $m$
and $n$, respectively. The Hilbert space $\mathcal{H}_{A+B}$ of the composite
system is given by the tensor product of the subsystems,
$\mathcal{H}_{A+B}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. A random pure state
of the composite system $\mathcal{H}_{A+B}$ is defined as a linear combination
of the coefficients $z_{i,j}$ and the complete bases
$\left\\{\Ket{i^{A}}\right\\}$ and $\left\\{\Ket{j^{B}}\right\\}$ of
$\mathcal{H}_{A}$ and $\mathcal{H}_{B}$,

$\Ket{\psi}=\sum_{i=1}^{m}\sum_{j=1}^{n}z_{i,j}\Ket{i^{A}}\otimes\Ket{j^{B}},$
(1.1)

where each $z_{i,j}$ follows the standard complex Gaussian distribution of
zero mean and unit variance with the probability constraint
$\sum_{i,j}|z_{i,j}|^{2}=1$. We now consider a superposition of the state
(1.1) of the form

$\Ket{\varphi}=\Ket{\psi}+\left(\mathbf{U}\otimes\mathbf{I}_{m}\right)\Ket{\psi},$
(1.2)

where $\mathbf{U}$ is an $m\times m$ random unitary matrix with the measure
proportional to $\det\left(\mathbf{I}_{m}+\mathbf{U}\right)^{2\alpha+1}$ [8].
The corresponding density matrix of the pure state (1.2) is

$\rho=\Ket{\varphi}\Bra{\varphi},$ (1.3)

which satisfies the natural probability constraint

${\mathrm{tr}}(\rho)=1.$ (1.4)

We assume without loss of generality that the dimension of subsystem $A$ is no
greater than that of the subsystem $B$, i.e., $m\leq n$. The reduced density
matrix $\rho_{A}$ of the smaller subsystem $A$ is computed by partial tracing
(a.k.a. purification) of the full density matrix (1.3) over the other
subsystem $B$, interpreted as the environment, as

$\rho_{A}={\mathrm{tr}}_{B}\rho.$ (1.5)

The resulting density of eigenvalues of $\rho_{A}$ is the
generalized111Hereafter, we refer to this generalized ensemble as the Bures-
Hall ensemble despite the term is often referred to the special case
$\alpha=-1/2$ in the literature. Bures-Hall ensemble [1, 2, 3, 8]

$f\left(\bm{\lambda}\right)=\frac{1}{c}~{}\delta\left(1-\sum_{i=1}^{m}\lambda_{i}\right)\prod_{1\leq
i<j\leq
m}\frac{\left(\lambda_{i}-\lambda_{j}\right)^{2}}{\lambda_{i}+\lambda_{j}}\prod_{i=1}^{m}\lambda_{i}^{\alpha}$
(1.6)

supported in the probability simplex

$\mathcal{D}=\Bigg{\\{}0\leq\lambda_{m}<\ldots<\lambda_{1}\leq
1,~{}~{}\sum_{i=1}^{m}\lambda_{i}=1\Bigg{\\}},$ (1.7)

where the parameter $\alpha$ takes half-integer values

$\alpha=n-m-\frac{1}{2}$ (1.8)

and the constant $c$ is

$c=\frac{2^{-m(m+2\alpha)}\pi^{m/2}}{\Gamma\left(m(m+2\alpha+1)/2\right)}\prod_{i=1}^{m}\frac{\Gamma(i+1)\Gamma(i+2\alpha+1)}{\Gamma(i+\alpha+1/2)}.$
(1.9)

A relatively detailed derivation of the Bures-Hall ensemble (1.6) can be
found, for example, in [7, Sec. 3]. The presence of the Dirac delta function
$\delta(\cdot)$ in (1.6) reflects the constraint (1.4). Note that another
approach to define the Bures-Hall ensemble (1.6) is by introducing a distance
metric, known as the Bures-Hall metric, over reduced density matrices [11].
The Bures-Hall ensemble satisfies several distinguished properties among the
measures of random density matrices. It is the only monotone metric that is
simultaneously Fisher adjusted and Fubini-Study adjusted. The Bures-Hall
metric, related to quantum distinguishability, is known to be the minimal
monotone metric [11]. It is also a function of fidelity [12], a key
performance indicator in quantum information processing. In addition, the
Bures-Hall measure enjoys the property that, without any prior knowledge on a
density matrix, the optimal way to estimate the density matrix is to generate
a state at random with respect to this measure. It is therefore often used as
a prior distribution, referred to as the Bures prior, in reconstructing
quantum states from measurements.

Entanglement serves as a measure of the non-classical correlation between the
subsystems $A$ and $B$. The degree of entanglement can be assessed by
entanglement entropies, which are functions of the eigenvalues (entanglement
spectrum) of a density matrix. Any function that satisfies a list of axioms
qualifies as an entanglement entropy [11]. In particular, an entropy should
monotonically change from the separable state

$\lambda_{1}=1,~{}~{}\lambda_{2}=\dots=\lambda_{m}=0$ (1.10)

to the maximally-entangled state

$\lambda_{1}=\lambda_{2}=\dots\lambda_{m}=\frac{1}{m}$ (1.11)

corresponding to boundaries of the support (2.82). In this work, we focus on
the study of quantum purity

$\text{P}=\sum_{i=1}^{m}\lambda_{i}^{2}$ (1.12)

supported in $\text{P}\in[1/m,1]$, which attains maximally-entangled state and
separable state when $\text{P}=1/m$ and $\text{P}=1$, respectively. Quantum
purity estimates how far a state is from a pure state that corresponds to
$\text{P}=1$, i.e., when the density matrix becomes a one dimensional
projector. Being a simple polynomial function, quantum purity is one of the
few entanglement entropies that can be measured experimentally [13].

Statistical information of an entanglement entropy such as the considered
quantum purity is encoded through its moments: the first moment (average
value) implies the typical behavior of entanglement, the second moment
(variance) specifies the fluctuation around the typical value, and the higher
order moments (such as skewness and kurtosis) describe the tails of the
distribution. Moreover, the moments can be also utilized to construct non-
asymptotic approximations to the distribution of purity. Due to the compact
support of purity, the resulting approximations become exact as the number of
moments involved increases as promised by the Weierstrass approximation
theorem. In the literature, the first moment of purity has been recently
obtained in [8, 4] as

$\mathbb{E}_{f}\\!\left[\text{P}\right]=\frac{m^{2}-2mn-4n^{2}-1}{2n\left(m^{2}-2mn-2\right)}.$
(1.13)

Closed-form expressions of the higher order moments remain open except for the
special case of equal subsystem dimensions $m=n$ (i.e., $\alpha=-1/2$), where
the exact second and third moments have been derived in [7] via an integrable
system method. In this work, we propose a framework that allows one to obtain
any higher order moment of purity of arbitrary subsystem dimensions in a
systematic manner. A key ingredient is the recurrence relations of the
underlying integrals of biorthogonal polynomials that we study in the first
place. In particular, using the proposed framework we derive the second and
third moments valid for any subsystem dimensions as summarized in the
following proposition, which contains the main results of this work.

Figure 1. The first three moments ($k=1,2,3$) of quantum purity
$\mathbb{E}_{f}\\!\left[\text{P}^{k}\right]$: analytical results versus
simulations. Solid lines represent the obtained analytical results (1.13),
(1.14), (1.15), and the diamond-shape scatters represent numerical
simulations. The left subplot and the right subplot refer to the cases of
equal subsystem dimensions $n=m$ and unequal subsystem dimensions $n=2m$,
respectively.

######  Proposition 1.1.

For any subsystem dimensions $m\leq n$, the exact second moment and third
moment of quantum purity (1.12) over the Bures-Hall ensemble (1.6) are given
respectively by

$\displaystyle\mathbb{E}_{f}\\!\left[\text{P}^{2}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{4\left(n^{2}-1\right)\prod_{i=1}^{3}\left(m^{2}-2mn-2i\right)}\big{(}m^{6}-6m^{5}n+4m^{4}n^{2}+6m^{4}+24m^{3}n^{3}-$
(1.14) $\displaystyle
24m^{3}n-16m^{2}n^{4}+112m^{2}n^{2}-107m^{2}-32mn^{5}-176mn^{3}+214mn-$
$\displaystyle 128n^{4}-32n^{2}+160\big{)}$

and

$\displaystyle\mathbb{E}_{f}\\!\left[\text{P}^{3}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{8n\left(n^{2}-1\right)\left(n^{2}-4\right)\prod_{i=1}^{5}\left(m^{2}-2mn-2i\right)}\big{(}m^{10}n^{2}-2m^{10}-10m^{9}n^{3}+20m^{9}n+$
(1.15) $\displaystyle
28m^{8}n^{4}-35m^{8}n^{2}+22m^{8}+16m^{7}n^{5}-200m^{7}n^{3}-176m^{7}n-160m^{6}n^{6}+$
$\displaystyle
992m^{6}n^{4}-697m^{6}n^{2}+642m^{6}+64m^{5}n^{7}-1808m^{5}n^{5}+6646m^{5}n^{3}-3852m^{5}n+$
$\displaystyle
320m^{4}n^{8}+320m^{4}n^{6}-9208m^{4}n^{4}+13715m^{4}n^{2}-9262m^{4}-128m^{3}n^{9}+$
$\displaystyle
4480m^{3}n^{7}-10760m^{3}n^{5}-29180m^{3}n^{3}+37048m^{3}n-256m^{2}n^{10}-3520m^{2}n^{8}+$
$\displaystyle
30032m^{2}n^{6}-48052m^{2}n^{4}-6520m^{2}n^{2}+31640m^{2}-3072mn^{9}-18944mn^{7}+$
$\displaystyle
144192mn^{5}-61056mn^{3}-63280mn-10240n^{8}+15360n^{6}+132480n^{4}-$
$\displaystyle 114560n^{2}-23040\big{)}.$

###### Proof.

Detailed derivations of the results (1.14) and (1.15) are provided in the next
section. ∎

In the special case $m=n$, the results in proposition 1.1 are reduced
respectively to

$\mathbb{E}_{f}\\!\left[\text{P}^{2}\right]=\frac{5\left(5m^{4}+47m^{2}+32\right)}{4\left(m^{2}+2\right)\left(m^{2}+4\right)\left(m^{2}+6\right)}$
(1.16)

and

$\mathbb{E}_{f}\\!\left[\text{P}^{3}\right]=\frac{5\left(25m^{8}+690m^{6}+6015m^{4}+8750m^{2}+1152\right)}{8m\left(m^{2}+2\right)\left(m^{2}+4\right)\left(m^{2}+6\right)\left(m^{2}+8\right)\left(m^{2}+10\right)}$
(1.17)

as reported222Notice the typo in the second moment expression in [7, equation
(60)], where the correct constant factor in the denominator is $4$ instead of
$2$ as also verified by numerical simulations. in [7, equation (60)]. Besides
the finite-size results, we also point out that the asymptotic distribution of
purity valid for large subsystem dimensions has been obtained in [14] using a
Coulomb gas method.

To illustrate the obtained results, we plot in figure 1 the exact moments as
compared with simulations, where the value of each scatter is computed by
averaging over $10^{6}$ realizations of the Bures-Hall density matrices. We
consider the scenarios of equal subsystem dimensions $n=m$ as well as unequal
subsystem dimensions $n=2m$, where the claimed results match well with the
simulations. By comparing the two subplots, we also observe that the values of
the moments decrease as the dimension of the larger subsystem $n$ increases
from $m$ to $2m$. The observation indicates that this additional subsystem
size ($m$ to $2m$) leads to a system that tends to concentrate on more
entangled states (i.e., with a smaller purity value), as expected.

The rest of the paper is organized as follows. In section 2, we provide
detailed derivations of the main results on the exact second and third moments
of quantum purity. Specifically, in section 2.1 we relate the moment
computation of any order to that over an unconstraint ensemble, whose
correlation functions are given explicitly. In section 2.2, we derive
necessary results on the recurrence relations so as to provide a new proof to
the recently obtained first moment formula. The second and third moments are
derived respectively in section 2.3 and section 2.4, where the recurrence
relation framework is further developed. Potential future work is outlined in
section 3 after summarizing the main findings of the paper.

##  2\. Moment computation and orthogonal polynomial recurrence

###  2.1. Moment relation and correlation functions

The first step is a rather standard procedure, see, e.g., [4, 5, 7, 8, 9, 15,
16], of relating the moment computation to that over an ensemble without the
constraint $\delta\left(1-\sum_{i=1}^{m}\lambda_{i}\right)$ in (1.6). It turns
out that the corresponding unconstraint ensemble is given by [8, 5, 9]

$h\left(\bm{x}\right)=\frac{1}{c^{\prime}}\prod_{1\leq i<j\leq
m}\frac{\left(x_{i}-x_{j}\right)^{2}}{x_{i}+x_{j}}\prod_{i=1}^{m}x_{i}^{\alpha}{\rm
e}^{-x_{i}},$ (2.1)

where $x_{i}\in[0,\infty)$, $i=1,\dots,m$, and the constant $c^{\prime}$
depends on the constant (1.9) through

$c^{\prime}=c~{}\Gamma\left(d\right)$ (2.2)

with $d$ denoting

$d=\frac{1}{2}m\left(m+2\alpha+1\right).$ (2.3)

Beyond the physically relevant $\alpha$ values in (1.8), the results hereafter
are in fact valid for any $\alpha>-1$. It is known that the density
$h\left(\bm{x}\right)$ in (2.1) admits the factorization [5]

$h(\bm{x})\prod_{i=1}^{m}\,\mathrm{d}x_{i}=f(\bm{\lambda})g_{d}(\theta)\,\mathrm{d}\theta\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i},$
(2.4)

where

$g_{d}(\theta)=\frac{1}{\Gamma\left(d\right)}{\rm e}^{-\theta}\theta^{d-1}$
(2.5)

is the density of the trace of the unconstraint ensemble

$\theta=\sum_{i=1}^{m}x_{i},~{}~{}~{}~{}~{}~{}\theta\in[0,\infty).$ (2.6)

The factorization (2.4) implies that the random variable $\theta$ is
independent of each of $\bm{\lambda}$.

Computing $k$-th moment of the purity P can now be converted to computing
$k$-th moment of the induced purity

$\text{T}=\sum_{i=1}^{m}x_{i}^{2}$ (2.7)

over the unconstraint ensemble (2.1) as

$\displaystyle\mathbb{E}_{f}\\!\left[\text{P}^{k}\right]$ $\displaystyle=$
$\displaystyle\int_{\bm{\lambda}}\frac{\text{T}^{k}}{\theta^{2k}}f(\bm{\lambda})\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}$
(2.8) $\displaystyle=$
$\displaystyle\int_{\bm{\lambda}}\frac{T^{k}}{\theta^{2k}}f(\bm{\lambda})\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}\int_{\theta}g_{d+2k}(\theta)\,\mathrm{d}\theta$
(2.9) $\displaystyle=$
$\displaystyle\frac{\Gamma(d)}{\Gamma(d+2k)}\int_{\bm{\lambda}}\text{T}^{k}f(\bm{\lambda})\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}\int_{\theta}g_{d}(\theta)\,\mathrm{d}\theta$
(2.10) $\displaystyle=$
$\displaystyle\frac{\Gamma(d)}{\Gamma(d+2k)}\mathbb{E}_{h}\\!\left[\text{T}^{k}\right],$
(2.11)

where we have used (2.4).

Now the task is to compute the induced moments
$\mathbb{E}_{h}\\!\left[\text{T}^{k}\right]$, the first three of which in
terms of the first three point densities $h_{1}(x)$, $h_{2}(x,y)$, and
$h_{3}(x,y,z)$ are written by definition as

$\displaystyle\mathbb{E}_{h}\\!\left[\text{T}\right]$ $\displaystyle=$
$\displaystyle m\int_{0}^{\infty}x^{2}h_{1}(x)\,\mathrm{d}x$ (2.12a)
$\displaystyle\mathbb{E}_{h}\\!\left[\text{T}^{2}\right]$ $\displaystyle=$
$\displaystyle
m\int_{0}^{\infty}\\!\\!x^{4}h_{1}(x)\,\mathrm{d}x+m(m-1)\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}~{}h_{2}\left(x,y\right)\,\mathrm{d}x\,\mathrm{d}y$
(2.12b) $\displaystyle\mathbb{E}_{h}\\!\left[\text{T}^{3}\right]$
$\displaystyle=$ $\displaystyle
m\int_{0}^{\infty}\\!\\!x^{6}h_{1}(x)\,\mathrm{d}x+3m(m-1)\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{4}y^{2}h_{2}(x,y)\,\mathrm{d}x\,\mathrm{d}y+$
(2.12c) $\displaystyle
m(m-1)(m-2)\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}h_{3}(x,y,z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z.$

In general, any $k$-point ($k\leq m$) density of the unconstraint ensemble is
known to admit a Pfaffian form of a $2k\times 2k$ antisymmetric matrix [17].
In particular, the needed densities can be read off as functions of the four
correlation kernels $K_{00}(x,y)$, $K_{01}(x,y)$, $K_{10}(x,y)$, and
$K_{11}(x,y)$ as [17]

$\displaystyle h_{1}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2m}\left(K_{01}(x,x)+K_{10}(x,x)\right)$ (2.13)
$\displaystyle h_{2}(x,y)$ $\displaystyle=$
$\displaystyle\frac{1}{4m(m-1)}(\left(K_{01}(x,x)+K_{10}(x,x)\right)\left(K_{01}(y,y)+K_{10}(y,y)\right)-2K_{01}(x,y)K_{01}(y,x)$
(2.14)
$\displaystyle-2K_{10}(x,y)K_{10}(y,x)-2K_{00}(x,y)K_{11}(x,y)-2K_{00}(y,x)K_{11}(y,x))$
$\displaystyle h_{3}(x,y,z)$ $\displaystyle=$
$\displaystyle\frac{1}{8m(m-1)(m-2)}\left(h_{\text{A}}+h_{\text{B}}+h_{\text{C}}+h_{\text{D}}\right),$
(2.15)

where

$\displaystyle h_{\text{A}}$ $\displaystyle=$
$\displaystyle(K_{01}(x,x)+K_{10}(x,x))(K_{01}(y,y)+K_{10}(y,y))(K_{01}(z,z)+K_{10}(z,z))$
(2.16) $\displaystyle h_{\text{B}}$ $\displaystyle=$
$\displaystyle-2(K_{01}(x,x)+K_{10}(x,x))(K_{01}(y,z)K_{01}(z,y)+K_{10}(y,z)K_{10}(z,y)+$
(2.17) $\displaystyle
K_{00}(y,z)K_{11}(y,z)+K_{00}(z,y)K_{11}(z,y))-2(K_{01}(y,y)+K_{10}(y,y))\times$
$\displaystyle(K_{01}(x,z)K_{01}(z,x)+K_{10}(x,z)K_{10}(z,x)+K_{00}(x,z)K_{11}(x,z)+$
$\displaystyle
K_{00}(z,x)K_{11}(z,x))-2(K_{01}(z,z)+K_{10}(z,z))(K_{01}(x,y)K_{01}(y,x)+$
$\displaystyle
K_{10}(x,y)K_{10}(y,x)+K_{00}(x,y)K_{11}(x,y)+K_{00}(y,x)K_{11}(y,x))$
$\displaystyle h_{\text{C}}$ $\displaystyle=$ $\displaystyle
2(K_{00}(x,y)K_{01}(y,z)K_{11}(x,z)+K_{00}(x,y)K_{10}(z,x)K_{11}(z,y)+$ (2.18)
$\displaystyle
K_{00}(y,x)K_{01}(x,z)K_{11}(y,z)+K_{00}(y,x)K_{10}(z,y)K_{11}(z,x)+$
$\displaystyle
K_{00}(x,z)K_{01}(z,y)K_{11}(x,y)+K_{00}(x,z)K_{10}(y,x)K_{11}(y,z)+$
$\displaystyle
K_{00}(z,x)K_{01}(x,y)K_{11}(z,y)+K_{00}(z,x)K_{10}(y,z)K_{11}(y,x)+$
$\displaystyle
K_{00}(y,z)K_{01}(y,x)K_{11}(x,z)+K_{00}(y,z)K_{10}(x,z)K_{11}(y,x)+$
$\displaystyle
K_{00}(z,y)K_{01}(z,x)K_{11}(x,y)+K_{00}(z,y)K_{10}(x,y)K_{11}(z,x)-$
$\displaystyle
K_{00}(x,y)K_{01}(x,z)K_{11}(y,z)-K_{00}(x,y)K_{10}(z,y)K_{11}(z,x)-$
$\displaystyle
K_{00}(y,x)K_{01}(y,z)K_{11}(x,z)-K_{00}(y,x)K_{10}(z,x)K_{11}(z,y)-$
$\displaystyle
K_{00}(x,z)K_{01}(x,y)K_{11}(z,y)-K_{00}(x,z)K_{10}(y,z)K_{11}(y,x)-$
$\displaystyle
K_{00}(z,x)K_{01}(z,y)K_{11}(x,y)-K_{00}(z,x)K_{10}(y,x)K_{11}(y,z)-$
$\displaystyle
K_{00}(y,z)K_{01}(z,x)K_{11}(x,y)-K_{00}(y,z)K_{10}(x,y)K_{11}(z,x)-$
$\displaystyle
K_{00}(z,y)K_{01}(y,x)K_{11}(x,z)-K_{00}(z,y)K_{10}(x,z)K_{11}(y,x))$
$\displaystyle h_{\text{D}}$ $\displaystyle=$ $\displaystyle
2(K_{01}(x,y)K_{01}(y,z)K_{01}(z,x)+K_{01}(x,z)K_{01}(z,y)K_{10}(x,y)+$ (2.19)
$\displaystyle
K_{01}(x,z)K_{01}(y,x)K_{10}(y,z)+K_{01}(y,x)K_{01}(z,y)K_{10}(z,x)+$
$\displaystyle
K_{01}(x,y)K_{10}(x,z)K_{10}(z,y)+K_{01}(y,z)K_{10}(x,z)K_{10}(y,x)+$
$\displaystyle
K_{01}(z,x)K_{10}(y,x)K_{10}(z,y)+K_{10}(x,y)K_{10}(y,z)K_{10}(z,x)).$

To make the discussion self-contained, we provide a direct verification of the
normalization constants of the above densities in the appendix A. Note that
the densities (2.14) and (2.15) have been simplified to the current form by
using the following factorization properties333A proof of the factorization
properties of the correlation kernels (2.20) can be found in [16]. [17]

$\displaystyle K_{00}(x,y)+K_{00}(y,x)$ $\displaystyle=$
$\displaystyle\ell_{1}(x)\ell_{1}(y)$ (2.20a) $\displaystyle
K_{01}(x,y)-K_{10}(y,x)$ $\displaystyle=$
$\displaystyle\ell_{2}(x)\ell_{1}(y)$ (2.20b) $\displaystyle
K_{11}(x,y)+K_{11}(y,x)$ $\displaystyle=$
$\displaystyle-\ell_{2}(x)\ell_{2}(y),$ (2.20c)

where

$\displaystyle\ell_{1}(x)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{m-1}\frac{\Gamma(i+m+2\alpha+2)(-x)^{i}}{\Gamma(i+2\alpha+2)\Gamma(i+\alpha+2)\Gamma(m-i)\Gamma(i+1)}$
(2.21a) $\displaystyle\ell_{2}(x)$ $\displaystyle=$
$\displaystyle-x^{2\alpha+1}\sum_{i=0}^{m-1}\frac{\Gamma(i+m+2\alpha+2)\Gamma(-i-\alpha,x)(-x)^{i}}{(i+\alpha+1)\Gamma(i+2\alpha+2)\Gamma(m-i)\Gamma(i+1)}+{\rm
e}^{-x}x^{\alpha}$ (2.21b)

with $\Gamma(a,x)=\int_{x}^{\infty}t^{a-1}{\rm e}^{-t}\,\mathrm{d}t$ denoting
the incomplete Gamma function.

The correlation kernels admit two sets of representations useful in the
subsequent calculations. The first one is the summation representation

$\displaystyle K_{00}(x,y)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{1}{h_{k}}p_{k}(x)q_{k}(y)$ (2.22a)
$\displaystyle K_{01}(x,y)$ $\displaystyle=$
$\displaystyle-x^{\alpha}e^{-x}\sum_{k=0}^{m-1}\frac{1}{h_{k}}p_{k}(y)Q_{k}(-x)$
(2.22b) $\displaystyle K_{10}(x,y)$ $\displaystyle=$
$\displaystyle-y^{\alpha+1}e^{-y}\sum_{k=0}^{m-1}\frac{1}{h_{k}}P_{k}(-y)q_{k}(x)$
(2.22c) $\displaystyle K_{11}(x,y)$ $\displaystyle=$ $\displaystyle
x^{\alpha}y^{\alpha+1}e^{-x-y}\sum_{k=0}^{m-1}\frac{1}{h_{k}}P_{k}(-y)Q_{k}(-x)-W(x,y),$
(2.22d)

where the weight function $W(x,y)$ of the (monic) Cauchy-Laguerre biorthogonal
polynomials $p_{k}(x)$ and $q_{l}(y)$,

$\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}p_{k}(x)q_{l}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=h_{k}\delta_{kl}$
(2.23)

is given by

$W(x,y)=\frac{x^{\alpha}y^{\alpha+1}e^{-x-y}}{x+y}.$ (2.24)

The functions in (2.22) are further related by [19, 17]

$\displaystyle P_{k}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{v^{\alpha}{\rm
e}^{-v}}{x-v}p_{k}(v)\,\mathrm{d}v$ (2.25a) $\displaystyle Q_{k}(y)$
$\displaystyle=$ $\displaystyle\int_{0}^{\infty}\frac{w^{\alpha+1}{\rm
e}^{-w}}{y-w}q_{k}(w)\,\mathrm{d}w$ (2.25b)

known as the Cauchy transforms of $p_{k}(x)$ and $q_{k}(y)$, respectively.
Moreover, they are expressed explicitly via Meijer G-functions as [18, 19, 17]

$\displaystyle p_{k}(x)$
$\displaystyle=(-1)^{k}\frac{\Gamma(k+1)\Gamma(2\alpha+k+2)\Gamma(\alpha+k+1)}{\Gamma(2\alpha+2k+2)}G_{2,3}^{1,1}\left(\left.\begin{array}[]{c}-2\alpha-k-1;~{}k+1\\\
0;~{}-\alpha,\,-2\alpha-1\end{array}\right|x\right)$ (2.26c) $\displaystyle
q_{k}(y)$
$\displaystyle=(-1)^{k}\frac{\Gamma(k+1)\Gamma(2\alpha+k+2)\Gamma(\alpha+k+2)}{\Gamma(2\alpha+2k+2)}G_{2,3}^{1,1}\left(\left.\begin{array}[]{c}-2\alpha-k-1;~{}k+1\\\
0;-\alpha-1,\,-2\alpha-1\end{array}\right|y\right)$ (2.26f) $\displaystyle
P_{k}(x)$
$\displaystyle=(-1)^{k+1}\frac{2\alpha+2k}{\Gamma(k)\Gamma(\alpha+k)}G_{2,3}^{3,1}\left(\left.\begin{array}[]{c}-k;\,k+2\alpha\\\
-1,\alpha-1,2\alpha;\end{array}\right|-x\right)$ (2.26i) $\displaystyle
Q_{k}(y)$
$\displaystyle=(-1)^{k+1}\frac{2\alpha+2k}{\Gamma(k)\Gamma(\alpha+k+1)}G_{2,3}^{3,1}\left(\left.\begin{array}[]{c}-k;\,k+2\alpha\\\
-1,\alpha,2\alpha;\end{array}\right|-y\right),$ (2.26l)

where the Meijer G-function is defined by the contour integral [20]

$\displaystyle
G_{p,q}^{m,n}\left(\left.\begin{array}[]{c}a_{1},\ldots,a_{n};a_{n+1},\ldots,a_{p}\\\
b_{1},\ldots,b_{m};b_{m+1},\ldots,b_{q}\end{array}\right|x.\right)$ (2.29)
$\displaystyle=$
$\displaystyle\frac{1}{2\pi\imath}\int_{\mathcal{L}}{\frac{\prod_{j=1}^{m}\Gamma\left(b_{j}+s\right)\prod_{j=1}^{n}\Gamma\left(1-a_{j}-s\right)x^{-s}}{\prod_{j=n+1}^{p}\Gamma\left(a_{j}+s\right)\prod_{j=m+1}^{q}\Gamma\left(1-b_{j}-s\right)}}\,\mathrm{d}s$
(2.30)

with the contour $\mathcal{L}$ separating the poles of
$\Gamma\left(1-a_{j}-s\right)$ from the poles of $\Gamma\left(b_{j}+s\right)$.
In addition to the summation representation (2.22), the correlation kernels
also admit the integral form [19]

$\displaystyle K_{00}(x,y)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}t^{2\alpha+1}H_{\alpha}(tx)H_{\alpha+1}(ty)\,\mathrm{d}t$
(2.31a) $\displaystyle K_{01}(x,y)$ $\displaystyle=$ $\displaystyle
x^{2\alpha+1}\int_{0}^{1}t^{2\alpha+1}H_{\alpha}(ty)G_{\alpha+1}(tx)\,\mathrm{d}t$
(2.31b) $\displaystyle K_{10}(x,y)$ $\displaystyle=$ $\displaystyle
y^{2\alpha+1}\int_{0}^{1}t^{2\alpha+1}H_{\alpha+1}(tx)G_{\alpha}(ty)\,\mathrm{d}t$
(2.31c) $\displaystyle K_{11}(x,y)$ $\displaystyle=$
$\displaystyle(xy)^{2\alpha+1}\int_{0}^{1}t^{2\alpha+1}G_{\alpha+1}(tx)G_{\alpha}(ty)\,\mathrm{d}t-\frac{x^{\alpha}y^{\alpha+1}}{x+y},$
(2.31d)

where we denote

$\displaystyle H_{q}(x)$ $\displaystyle=$ $\displaystyle
G_{2,3}^{1,1}\left(\left.\begin{array}[]{c}-m-2\alpha-1;~{}m\\\
0;-q,-2\alpha-1\end{array}\right|x\Big{.}\right)$ (2.32c) $\displaystyle
G_{q}(x)$ $\displaystyle=$ $\displaystyle
G_{2,3}^{2,1}\left(\left.\begin{array}[]{c}-m-2\alpha-1;~{}m\\\
0,-q;-2\alpha-1\end{array}\right|x\Big{.}\right).$ (2.32f)

With the results introduced above, the moments (2.12) can now be computed by
integrating over the appropriate form of the correlation kernels. However, the
direct integration may not be an efficient method for the moment computation.
We now discuss how the recurrence coefficients of the Cauchy-Laguerre
biorthogonal polynomials are useful in obtaining the moments of purity. For
positive integers $\beta$ and $\gamma$, the idea is to successively reduce the
power of the average of $x^{\beta}y^{\gamma}$ over the ensemble to averages of
lower order powers. The reduction is substantially simplified by the
recurrence relations of the biorthogonal polynomials. Instead of performing
integrations, the only task now is to collect non-zero contributions as
identified by the biorthogonality condition. As a result, the desired
integrals boil down to summations involving the recurrence coefficients of the
biorthogonal polynomials, which are known explicitly. Besides the considered
purity, the approach in fact works for any polynomial linear spectral
statistics. To some extend, the idea is similar to that in Aomoto’s proof of
the Selberg integral [21, chapter 17.3-17.4], which concerns the computation
of averages over classical random matrix ensembles [21, 22]. The proposed
method also has independent interest in the field of combinatorics [23] as
briefly outlined here. Consider an oriented graph $G=(V,E)$ with vertices
$V=\mathbb{N}^{2}$ and edges
$E=\left\\{(n,k)\to(n+1,m),\,n,k\in\mathbb{N},\,0\leq m\leq k+1\right\\}$,
where the weight on the path is given by $\omega((n,k)\to(n+1,m))=\langle
xp_{k},q_{m}\rangle$. We then have

$\displaystyle xp_{k}(x)=\sum_{\ell=0}^{k+1}\langle xp_{k},q_{\ell}\rangle
p_{\ell}(x),$

and this leads, by induction, to

$\displaystyle
x^{n}p_{k}(x)=\sum_{\ell=0}^{k+n}\left(\sum_{\gamma:(0,k)\to(n,\ell)}\prod_{e\in\gamma}\omega(e)\right)p_{\ell}(x).$

The summation over the oriented paths $\gamma$ on $G$ starts from $(0,k)$ and
ends at $(n,\ell)$, where each path picks the product of the weights along the
edges it crosses.

###  2.2. Computation of the first moment

To illustrate the proposed method, here we rederive the first moment of purity
recently reported in [8, 9]. By the definition (2.12a), one has

$\displaystyle\mathbb{E}_{h}\\!\left[\text{T}\right]=m\int_{0}^{\infty}\\!\\!x^{2}h_{1}(x)\,\mathrm{d}x=\frac{1}{2}\int_{0}^{\infty}\\!\\!x^{2}K_{01}(x,x)\,\mathrm{d}x+\frac{1}{2}\int_{0}^{\infty}\\!\\!x^{2}K_{10}(x,x)\,\mathrm{d}x.$
(2.33)

We now focus on the first integral

$\displaystyle\int_{0}^{\infty}\\!\\!x^{2}K_{01}(x,x)\,\mathrm{d}x,$ (2.34)

where the second one can be similarly computed. By the summation form of the
kernels (2.22), we have

$\displaystyle\int_{0}^{\infty}\\!\\!x^{2}K_{01}(x,x)\,\mathrm{d}x=\sum_{k=0}^{m-1}\frac{1}{h_{k}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}p_{k}(x)q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y.$
(2.35)

To compute the integral (2.35), it is more convenient to consider a general
form

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{\beta}p_{k}(x)q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y$
(2.36)

with $\beta\in\mathbb{Z}_{+}$. To this end, we first state the following
proposition.

######  Proposition 2.1.

Let $\\{a_{k,i}\\}_{i\leq k}$ denote the coefficients of polynomials
$p_{k}(x)$ in (2.22a) for an arbitrary $k\in\mathbb{N}$, i.e.,

$\displaystyle p_{k}(x)=x^{k}+a_{k,k-1}x^{k-1}+\cdots+a_{k,0}x^{0}$

with the convention $a_{k,k}=1$. For the Cauchy-Laguerre biorthogonal
polynomials defined by (2.23)-(2.24), one has

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{\beta}q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=b_{\beta,\beta-k}h_{k},\quad\beta\geq
k,$ (2.37)

where $h_{k}$ is the normalization constant in (2.23) and $b_{\beta,k}$ is
recursively given by

$\displaystyle
b_{\beta,i}=-\sum_{j=0}^{i-1}b_{\beta,j}a_{\beta-j,\beta-i},\quad
b_{\beta,0}=1.$ (2.38)

###### Proof.

The idea is to write $x^{\beta}$ in terms of the orthogonal basis
$\\{p_{k}(x)\\}_{k\leq\beta}$ before applying the orthogonality condition.
When $\beta=k$, it is obvious that

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{k}q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=h_{k}$

by orthogonality. When $\beta=k+1$, since

$\displaystyle x^{k+1}=p_{k+1}(x)-a_{k+1,k}p_{k}(x)+\text{l.o.t.},$

the lower order terms (l.o.t.) do not contribute to the result when taking the
inner product with $q_{k}(y)$ on both sides, one finds that

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{k+1}q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=-a_{k+1,k}h_{k}.$

Continuing this procedure establishes the result of this proposition. ∎

Moreover, if we expand

$\displaystyle
x^{\beta}p_{k}(x)=x^{\beta+k}+a_{k,k-1}x^{\beta+k-1}+\cdots+a_{k,0}x^{\beta}$

and employ the result of above proposition, we then have

$\displaystyle\frac{1}{h_{k}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{\beta}p_{k}(x)q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=\sum_{j=0}^{\beta}a_{k,k-j}b_{\beta+k-j,\beta-j}.$
(2.39)

Similarly, if we expand $q_{k}(y)$ as

$\displaystyle
q_{k}(y)=y^{k}+\hat{a}_{k,k-1}y^{k-1}+\cdots+\hat{a}_{k,0}y^{0},$

then by the recursion

$\displaystyle\hat{b}_{\beta,i}=-\sum_{j=0}^{i-1}\hat{b}_{\beta,j}\hat{a}_{\beta-j,\beta-i},\quad\hat{b}_{\beta,0}=1,$
(2.40)

one has

$\displaystyle\frac{1}{h_{k}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}y^{\beta}p_{k}(x)q_{k}(y)W(x,y)dxdy=\sum_{j=0}^{\beta}\hat{a}_{k,k-j}\hat{b}_{\beta+k-j,\beta-j}.$

Putting the results together, one has

$\displaystyle\mathbb{E}_{h}[\text{T}]=\frac{1}{2}\sum_{k=0}^{m-1}\sum_{j=0}^{2}\left(a_{k,k-j}b_{k+2-j,2-j}+\hat{a}_{k,k-j}\hat{b}_{k+2-j,2-j}\right),$
(2.41)

where the explicit formulas for $\\{a_{k,j},\,\hat{a}_{k,j}\\}_{j\leq k}$ are
given by [17, equation (2.6)]

$\displaystyle\begin{aligned} a_{k,j}&=(-1)^{k-j}{k\choose
j}\frac{\Gamma(2a+k+j+2)\Gamma(2a+k+2)\Gamma(a+k+1)}{\Gamma(2a+2k+2)\Gamma(2a+j+2)\Gamma(a+j+1)}\\\
\hat{a}_{k,j}&=(-1)^{k-j}{k\choose
j}\frac{\Gamma(2a+k+j+2)\Gamma(2a+k+2)\Gamma(a+k+2)}{\Gamma(2a+2k+2)\Gamma(2a+j+2)\Gamma(a+j+2)}.\end{aligned}$
(2.42)

Inserting the above $a_{k,j}$ and $\hat{a}_{k,j}$ into (2.38) and (2.40),
respectively, we claim that

$\displaystyle\begin{aligned} b_{k,j}&={k\choose
j}\frac{\Gamma(2a+k+2)\Gamma(a+k+1)\Gamma(2a+2k+3-2j)}{\Gamma(2a+k+2-j)\Gamma(a+k+1-j)\Gamma(2a+2k+3-j)}\\\
\hat{b}_{k,j}&={k\choose
j}\frac{\Gamma(2a+k+2)\Gamma(a+k+2)\Gamma(2a+2k+3-2j)}{\Gamma(2a+k+2-j)\Gamma(a+k+2-j)\Gamma(2a+2k+3-j)},\end{aligned}$
(2.43)

which is equivalent of stating that

$\displaystyle\sum_{j=0}^{i-1}(-1)^{j}{i\choose
j}\frac{(2{\alpha}+2k+2-2j)\Gamma(2{\alpha}+2k-i-j+2)}{\Gamma(2{\alpha}+2k+3-j)}=(-1)^{i-1}\frac{\Gamma(2{\alpha}+2k+3-2i)}{\Gamma(2{\alpha}+2k+3-i)}.$
(2.44)

The left hand side of the above identity is equal to

$\displaystyle\sum_{j=0}^{i-1}\frac{(-1)^{j}\Gamma(i+1)}{\Gamma(j+1)\Gamma(i+1-j)}\frac{(2{\alpha}+2k+2-2j)\Gamma(2{\alpha}+2k-i-j+2)}{\Gamma(2{\alpha}+2k+3-j)}$
(2.45) $\displaystyle=$
$\displaystyle(-1)^{i-1}\sum_{\ell=0}^{i-1}\frac{(-1)^{\ell}\Gamma(i+1)}{\Gamma(\ell+2)\Gamma(i-\ell)}\frac{(2{\alpha}+2k+4-2i+2\ell)\Gamma(2{\alpha}+2k+3-2i+\ell)}{\Gamma(2{\alpha}+2k+4-i+\ell)}$
$\displaystyle=$
$\displaystyle\frac{(-1)^{i}}{\Gamma(-i)}\sum_{\ell=0}^{i-1}\frac{(2{\alpha}+2k+4-2i+2\ell)\Gamma(2{\alpha}+2k+3-2i+\ell)\Gamma(1-i+\ell)}{\Gamma(2{\alpha}+2k+4-i+\ell)\Gamma(2+\ell)}$
$\displaystyle=$
$\displaystyle\frac{(-1)^{i}}{\Gamma(-i)}\frac{1}{-i}\left(\frac{\Gamma(2{\alpha}+2k+3-i)}{\Gamma(2{\alpha}+2k+3)\Gamma(i+1)}-\frac{\Gamma(1-i)\Gamma(2{\alpha}+2k-2i+3)}{\Gamma(2{\alpha}+2k-i+3)}\right)$
$\displaystyle=$
$\displaystyle(-1)^{i-1}\frac{\Gamma(2{\alpha}+2k-2i+3)}{\Gamma(2{\alpha}+2k-i+3)},$

where the second to last equality is a direct application of [19, lemma 4.1].

Finally, inserting (2.42) and (2.43) into (2.41), one obtains

$\mathbb{E}_{h}\\!\left[\text{T}\right]=\frac{m(2\alpha+m+1)\left(4\alpha^{2}+4\alpha+5m^{2}+10\alpha
m+5m+2\right)}{4(2\alpha+2m+1)},$ (2.46)

which upon applying the moment relation (2.11) for $k=1$ gives

$\displaystyle\mathbb{E}_{f}\\!\left[\text{P}\right]$ $\displaystyle=$
$\displaystyle\frac{4}{m(m+2\alpha+1)\left(m^{2}+2\alpha
m+m+2\right)}\mathbb{E}_{h}\\!\left[\text{T}\right]$ (2.47) $\displaystyle=$
$\displaystyle\frac{5m^{2}+10\alpha
m+5m+4\alpha^{2}+4\alpha+2}{(2m+2\alpha+1)\left(m^{2}+2\alpha m+m+2\right)}.$

This verifies the mean purity formula (1.13) obtained in [8, 9] when replacing
$\alpha$ by (1.8).

###  2.3. Computation of the second moment

We now proceed to the computation of the second moment
$\mathbb{E}_{h}\\!\left[\text{T}^{2}\right]$. Recall that

$\displaystyle\mathbb{E}_{h}\\!\left[\text{T}^{2}\right]=m\int_{0}^{\infty}\\!\\!x^{4}h_{1}(x)\,\mathrm{d}x+m(m-1)\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}~{}h_{2}\left(x,y\right)\,\mathrm{d}x\,\mathrm{d}y,$

where the first term is equal to

$\displaystyle\frac{1}{2}\int_{0}^{\infty}x^{4}(K_{01}(x,x)+K_{10}(x,x))\,\mathrm{d}x=\frac{1}{2}\sum_{k=0}^{m-1}\sum_{j=0}^{4}\left(a_{k,k-j}b_{k+4-j,4-j}+\hat{a}_{k,k-j}\hat{b}_{k+4-j,4-j}\right).$

The second term becomes
$\mathbb{E}_{h}\\!\left[\text{T}^{2}\right]-\mathcal{J}$ if we denote

$\displaystyle\mathcal{J}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}\left(K_{01}(x,y)K_{01}(y,x)+K_{10}(x,y)K_{10}(y,x)\right)\,\mathrm{d}x\,\mathrm{d}y-$
(2.48)
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}K_{00}(x,y)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y.$

The identity

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}K_{01}(x,y)K_{01}(y,x)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\sum_{j,k=0}^{m-1}\left(\frac{1}{h_{j}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}p_{k}(x)q_{j}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y\right)\left(\frac{1}{h_{k}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}p_{j}(x)q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y\right)$

inspires us to consider a more general case

$\displaystyle\frac{1}{h_{j}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{\beta}p_{k}(x)q_{j}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y,\quad\beta\in\mathbb{N}$

for later use.

######  Corollary 2.2.

For arbitrary $\beta,\,k,\,j\in\mathbb{N}$, we have

$\displaystyle\frac{1}{h_{j}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{\beta}p_{k}(x)q_{j}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=\sum_{\ell=0}^{k}a_{k,\ell}b_{\beta+\ell,\beta+\ell-j}$
$\displaystyle\frac{1}{h_{j}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!y^{\beta}p_{j}(x)q_{k}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=\sum_{\ell=0}^{k}\hat{a}_{k,\ell}\hat{b}_{\beta+\ell,\beta+\ell-j}.$

This corollary is a generalization of the result (2.39) and it should be
remarked that when $\beta+\ell<j$, the term $b_{\beta+\ell,\beta+\ell-j}$ is
zero. Therefore, noting that $\\{p_{k}(x)\\}_{k\in\mathbb{N}}$ and
$\\{q_{k}(y)\\}_{k\in\mathbb{N}}$ are dual to each other, we have

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}K_{01}(x,y)K_{01}(y,x)\,\mathrm{d}x\,\mathrm{d}y$
(2.49) $\displaystyle=$
$\displaystyle\sum_{j,k=0}^{m-1}\left(\sum_{\ell=0}^{k}a_{k,\ell}b_{\ell+2,\ell+2-j}\right)\left(\sum_{\ell=0}^{j}a_{j,\ell}b_{\ell+2,\ell+2-k}\right)$

and

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}K_{10}(x,y)K_{10}(y,x)\,\mathrm{d}x\,\mathrm{d}y$
(2.50) $\displaystyle=$
$\displaystyle\sum_{j,k=0}^{m-1}\left(\sum_{\ell=0}^{k}\hat{a}_{k,\ell}\hat{b}_{\ell+2,\ell+2-j}\right)\left(\sum_{\ell=0}^{j}\hat{a}_{j,\ell}\hat{b}_{\ell+2,\ell+2-k}\right).$

We proceed to evaluate (2.49). By changing the order of summation, one obtains

$\displaystyle\sum_{\ell_{1}=0}^{m-1}\sum_{\ell_{2}=0}^{m-1}\left(\sum_{j=\ell_{1}}^{m-1}a_{j,\ell_{2}}b_{\ell_{1}+2,\ell_{1}+2-j}\right)\left(\sum_{k=\ell_{2}}^{m-1}a_{k,\ell_{1}}b_{\ell_{2}+2,\ell_{2}+2-k}\right),$

and by substituting (2.42) and (2.43) into above equation, it leads to

$\displaystyle\sum_{j=\ell_{1}}^{m-1}a_{j,\ell_{2}}b_{\ell_{1}+2,\ell+2-j}$
$\displaystyle=$
$\displaystyle(-1)^{m+\ell_{1}}\frac{\Gamma(\ell_{2}+3)\Gamma(2{\alpha}+\ell_{2}+4)\Gamma({\alpha}+\ell_{2}+3)}{\Gamma(\ell_{1}+1)\Gamma(2{\alpha}+\ell_{1}+2)\Gamma({\alpha}+\ell_{1}+1)}\times$
$\displaystyle\frac{(m-\ell_{1})(4+2{\alpha}+m+\ell_{2})\Gamma(2a+\ell_{1}+m+2)}{(\ell_{1}-\ell_{2}-2)\Gamma(m-\ell_{1}+1)\Gamma(\ell_{2}-m+3)\Gamma(2{\alpha}+m+\ell_{2}+5)}.$

Therefore, the summation (2.49) is evaluated to

$\displaystyle-\frac{m^{2}(2{\alpha}+m+1)^{2}}{48(2{\alpha}+2m+1)^{2}}\big{(}8{\alpha}^{2}-56{\alpha}^{4}+12{\alpha}m-28{\alpha}^{2}m-264{\alpha}^{3}m+m^{2}-66{\alpha}m^{2}-454{\alpha}^{2}m^{2}+$
$\displaystyle
8{\alpha}^{4}m^{2}-32m^{3}-324{\alpha}m^{3}+4{\alpha}^{2}m^{3}+24{\alpha}^{3}m^{3}-81m^{4}+6{\alpha}m^{4}+26{\alpha}^{2}m^{4}+2m^{5}+12{\alpha}m^{5}+2m^{6}\big{)},$

where the one (2.50) can be similarly obtained. For the last integral in
(2.48),

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}K_{00}(x,y)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y,$
(2.51)

we proceed to study a generalized one

$\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{\beta_{1}}y^{\beta_{2}}K_{00}(x,y)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y$
(2.52)

with $\beta_{1}$ and $\beta_{2}$ being non-negative integers. By using the
summation form of the kernels (2.22) and the explicit polynomial expressions
of $p_{k}(x)$ and $q_{k}(y)$, the resulting integrals over $x$ and $y$ can be
separately evaluated by the Mellin transform of Meijer G-function [20]

$\displaystyle\int_{0}^{\infty}\\!\\!x^{s-1}G_{p,q}^{m,n}\left(\left.\begin{array}[]{c}a_{1},\ldots,a_{n};a_{n+1},\ldots,a_{p}\\\
b_{1},\ldots,b_{m};b_{m+1},\ldots,b_{q}\end{array}\right|\eta
x\Big{.}\right)\,\mathrm{d}x$ (2.55) $\displaystyle=$
$\displaystyle\frac{\eta^{-s}\prod_{j=1}^{m}\Gamma\left(b_{j}+s\right)\prod_{j=1}^{n}\Gamma\left(1-a_{j}-s\right)}{\prod_{j=n+1}^{p}\Gamma\left(a_{j}+s\right)\prod_{j=m+1}^{q}\Gamma\left(1-b_{j}-s\right)}.$
(2.56)

This leads to

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{\beta_{1}}y^{\beta_{2}}K_{00}(x,y)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{m-1}\sum_{k=0}^{m-1}\sum_{i=0}^{j}\sum_{s=0}^{j}\frac{4(-1)^{i+s}(\alpha+j+1)(\alpha+k+1)\Gamma(2\alpha+i+j+2)\Gamma(2\alpha+j+s+2)}{i!s!\Gamma(\alpha+i+1)\Gamma(2\alpha+i+2)\Gamma(\alpha+s+2)\Gamma(2\alpha+s+2)}\times$
$\displaystyle\frac{\Gamma\left(\beta_{1}+i+1\right)\Gamma\left(\alpha+\beta_{1}+i+1\right)\Gamma\left(2\alpha+\beta_{1}+i+2\right)}{\Gamma\left(\beta_{1}+i-k+1\right)\Gamma\left(2\alpha+\beta_{1}+i+k+3\right)\Gamma(j-i+1)}\times$
$\displaystyle\frac{\Gamma\left(\beta_{2}+s+1\right)\Gamma\left(\alpha+\beta_{2}+s+2\right)\Gamma\left(2\alpha+\beta_{2}+s+2\right)}{\Gamma\left(\beta_{2}-k+s+1\right)\Gamma\left(2\alpha+\beta_{2}+k+s+3\right)\Gamma(j-s+1)}-\sum_{j=0}^{m-1}\sum_{i=0}^{j}\sum_{k=0}^{j}\frac{2(-1)^{i+k}}{i!k!}\times$
$\displaystyle\frac{(\alpha+j+1)\Gamma(2\alpha+i+j+2)\Gamma(2\alpha+j+k+2)}{\Gamma(j-i+1)\Gamma(j-k+1)\Gamma(\alpha+i+1)\Gamma(2\alpha+i+2)\Gamma(\alpha+k+2)\Gamma(2\alpha+k+2)}\times$
$\displaystyle\frac{\Gamma\left(\alpha+\beta_{1}+i+1\right)\Gamma\left(\alpha+\beta_{2}+k+2\right)}{2\alpha+\beta_{1}+\beta_{2}+i+k+2}$
$\displaystyle=$
$\displaystyle\sum_{i=m-\beta_{1}}^{m-1}\sum_{k=m-\beta_{2}}^{m-1}\frac{(-1)^{i+k-1}\Gamma(2\alpha+i+m+2)\Gamma(2\alpha+k+m+2)}{(2\alpha+i+k+2)\left(2\alpha+\beta_{1}+\beta_{2}+i+k+2\right)i!k!\Gamma(m-i)\Gamma(m-k)}\times$
$\displaystyle\frac{\Gamma\left(\beta_{1}+i+1\right)\Gamma\left(\beta_{2}+k+1\right)\Gamma\left(\alpha+\beta_{1}+i+1\right)\Gamma\left(\alpha+\beta_{2}+k+2\right)}{\Gamma(\alpha+i+1)\Gamma(2\alpha+i+2)\Gamma(\alpha+k+2)\Gamma(2\alpha+k+2)}\times$
$\displaystyle\frac{\Gamma\left(2\alpha+\beta_{1}+i+2\right)\Gamma\left(2\alpha+\beta_{2}+k+2\right)}{\Gamma\left(\beta_{1}+i-m+1\right)\Gamma\left(\beta_{2}+k-m+1\right)\Gamma\left(2\alpha+\beta_{1}+i+m+2\right)\Gamma\left(2\alpha+\beta_{2}+k+m+2\right)},$

where the summations over $j$ and $k$ in (2.3) have been evaluated by the
identities, cf. [5, equations (71)-(72)],

$\sum_{j=i}^{m-1}\frac{(\alpha+j+1)\Gamma(2\alpha+i+j+2)\Gamma(2\alpha+j+s+2)}{\Gamma(j-i+1)\Gamma(j-s+1)}=\frac{\Gamma(i+m+2\alpha+2)\Gamma(s+m+2\alpha+2)}{2(2\alpha+i+s+2)\Gamma(m-i)\Gamma(m-s)},\\\
$ (2.59)
$\displaystyle\\!\\!\\!\\!\\!\\!\sum_{k=0}^{m-1}\frac{(\alpha+k+1)\left(\Gamma\left(\beta_{1}+i-k+1\right)\Gamma\left(\beta_{2}-k+s+1\right)\right)^{-1}}{\Gamma\left(2\alpha+\beta_{1}+i+k+3\right)\Gamma\left(2\alpha+\beta_{2}+k+s+3\right)}=\frac{1}{2\left(2\alpha+\beta_{1}+\beta_{2}+i+s+2\right)}\times$
$\displaystyle\\!\\!\\!\\!\\!\\!\left(\frac{\left(\Gamma\left(i+\beta_{1}+1\right)\Gamma\left(s+\beta_{2}+1\right)\right)^{-1}}{\Gamma\left(i+2\alpha+\beta_{1}+2\right)\Gamma\left(s+2\alpha+\beta_{2}+2\right)}-\frac{\left(\Gamma\left(i-m+\beta_{1}+1\right)\Gamma\left(s-m+\beta_{2}+1\right)\right)^{-1}}{\Gamma\left(i+m+2\alpha+\beta_{1}+2\right)\Gamma\left(s+m+2\alpha+\beta_{2}+2\right)}\right).$

In the interested case $\beta_{1}=\beta_{2}=2$, one readily obtains

$\displaystyle\int_{(0,\infty)^{2}}x^{2}y^{2}K_{00}(x,y)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle-\frac{m^{2}(2\alpha+m+1)^{2}}{8(2\alpha+2m-1)(2\alpha+2m+1)^{2}(2\alpha+2m+3)}\big{(}32\alpha^{4}+64\alpha^{3}+12\alpha^{2}-20\alpha+33m^{4}+$
$\displaystyle 132\alpha m^{3}+66m^{3}+196\alpha^{2}m^{2}+196\alpha
m^{2}+15m^{2}+128\alpha^{3}m+192\alpha^{2}m+28\alpha m-18m-6\big{)}.$

Combining results of (2.49), (2.50), and (2.3), we finally arrive at

$\displaystyle\mathbb{E}_{h}\\!\left[\text{T}^{2}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{k=0}^{m-1}\sum_{j=0}^{4}\left(a_{k,k-j}b_{k+4-j,4-j}+\hat{a}_{k,k-j}\hat{b}_{k+4-j,4-j}\right)+\mathbb{E}_{h}^{2}\\!\left[\text{T}\right]-\mathcal{J}$
(2.62) $\displaystyle=$
$\displaystyle\frac{m(2\alpha+m+1)}{16(2\alpha+2m-1)(2\alpha+2m+3)}\big{(}128\alpha^{4}+256\alpha^{3}+224\alpha^{2}+96\alpha+25m^{6}+$
$\displaystyle 150\alpha m^{5}+75m^{5}+340\alpha^{2}m^{4}+340\alpha
m^{4}+295m^{4}+360\alpha^{3}m^{3}+540\alpha^{2}m^{3}+$ $\displaystyle
1110\alpha
m^{3}+465m^{3}+176\alpha^{4}m^{2}+352\alpha^{3}m^{2}+1448\alpha^{2}m^{2}+1272\alpha
m^{2}+232m^{2}+$ $\displaystyle
32\alpha^{5}m+80\alpha^{4}m+768\alpha^{3}m+1072\alpha^{2}m+376\alpha
m+12m-144\big{)}.$

Employing the moment relation (2.11) for $k=2$ gives the second moment of
quantum purity

$\mathbb{E}_{f}\\!\left[\text{P}^{2}\right]=\frac{\Gamma(d)}{\Gamma(d+4)}\mathbb{E}_{h}\\!\left[\text{T}^{2}\right],$
(2.63)

which upon inserting $\alpha=n-m-1/2$ leads to the claimed result (1.14).

###  2.4. Computation of the third moment

As shown in (2.12c), computing the third moment boils down to computing the
following three integrals

$\displaystyle\text{I}_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!x^{6}h_{1}(x)\,\mathrm{d}x$ (2.64)
$\displaystyle\text{I}_{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{4}y^{2}h_{2}(x,y)\,\mathrm{d}x\,\mathrm{d}y$
(2.65) $\displaystyle\text{I}_{3}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}h_{3}(x,y,z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z,$
(2.66)

where the first two integrals can be obtained in the same manner as in the
first two moments calculation. They are given by

$\displaystyle\text{I}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{k=0}^{m-1}\sum_{j=0}^{6}\left(a_{k,k-j}b_{k+6-j,6-j}+\hat{a}_{k,k-j}\hat{b}_{k+6-j,6-j}\right)$
(2.68) $\displaystyle=$
$\displaystyle\frac{1}{32\prod_{i=1}^{5}(2\alpha+2m+2i-5)}\big{(}1024\alpha^{11}+5632\alpha^{10}+48640\alpha^{9}+176640\alpha^{8}+278592\alpha^{7}+$
$\displaystyle
190176\alpha^{6}-522880\alpha^{5}-1398640\alpha^{4}-930016\alpha^{3}-10608\alpha^{2}+260640\alpha+7293m^{11}+$
$\displaystyle 87516\alpha
m^{10}+43758m^{10}+463320\alpha^{2}m^{9}+463320\alpha
m^{9}+139425m^{9}+1424280\alpha^{3}m^{8}+$ $\displaystyle
2136420\alpha^{2}m^{8}+1304160\alpha
m^{8}+296010m^{8}+2814240\alpha^{4}m^{7}+5628480\alpha^{3}m^{7}+$
$\displaystyle 5268120\alpha^{2}m^{7}+2453880\alpha
m^{7}+253539m^{7}+3734016\alpha^{5}m^{6}+9335040\alpha^{4}m^{6}+$
$\displaystyle 12046320\alpha^{3}m^{6}+8734440\alpha^{2}m^{6}+1728012\alpha
m^{6}-280566m^{6}+3370752\alpha^{6}m^{5}+$ $\displaystyle
10112256\alpha^{5}m^{5}+17149440\alpha^{4}m^{5}+17445120\alpha^{3}m^{5}+5120808\alpha^{2}m^{5}-1916376\alpha
m^{5}-$ $\displaystyle
1049565m^{5}+2048640\alpha^{7}m^{4}+7170240\alpha^{6}m^{4}+15708000\alpha^{5}m^{4}+21344400\alpha^{4}m^{4}+$
$\displaystyle 8657880\alpha^{3}m^{4}-4772460\alpha^{2}m^{4}-5481960\alpha
m^{4}-1404810m^{4}+808320\alpha^{8}m^{3}+$ $\displaystyle
3233280\alpha^{7}m^{3}+9199680\alpha^{6}m^{3}+16282560\alpha^{5}m^{3}+9014040\alpha^{4}m^{3}-5337360\alpha^{3}m^{3}-$
$\displaystyle 10695600\alpha^{2}m^{3}-5852040\alpha
m^{3}-806532m^{3}+192000\alpha^{9}m^{2}+864000\alpha^{8}m^{2}+$ $\displaystyle
3285120\alpha^{7}m^{2}+7465920\alpha^{6}m^{2}+5693184\alpha^{5}m^{2}-2415840\alpha^{4}m^{2}-9654240\alpha^{3}m^{2}-$
$\displaystyle 8908560\alpha^{2}m^{2}-2411328\alpha
m^{2}+49608m^{2}+23552\alpha^{10}m+117760\alpha^{9}m+635520\alpha^{8}m+$
$\displaystyle
1835520\alpha^{7}m+1970304\alpha^{6}m-18816\alpha^{5}m-3921560\alpha^{4}m-5905840\alpha^{3}m-$
$\displaystyle 2513736\alpha^{2}m+165456\alpha m+246240m+86400\big{)}$

and

$\displaystyle\text{I}_{2}$ $\displaystyle=$
$\displaystyle\frac{(2\alpha+m)(2\alpha+m+1)}{64(2\alpha+2m-3)(2\alpha+2m-1)(2\alpha+2m+1)(2\alpha+2m+5)}\big{(}256\alpha^{8}-1152\alpha^{7}+$
(2.69) $\displaystyle
4032\alpha^{6}-1440\alpha^{5}-22656\alpha^{4}+55152\alpha^{3}-60832\alpha^{2}+33840\alpha+1155m^{8}+9240\alpha
m^{7}-$ $\displaystyle 1518m^{7}+31164\alpha^{2}m^{6}-11802\alpha
m^{6}+1506m^{6}+57624\alpha^{3}m^{5}-37896\alpha^{2}m^{5}+$ $\displaystyle
10074\alpha
m^{5}+6024m^{5}+63504\alpha^{4}m^{4}-64992\alpha^{3}m^{4}+30456\alpha^{2}m^{4}+32502\alpha
m^{4}-$ $\displaystyle
32145m^{4}+42336\alpha^{5}m^{3}-63936\alpha^{4}m^{3}+51168\alpha^{3}m^{3}+60576\alpha^{2}m^{3}-137754\alpha
m^{3}+$ $\displaystyle
52158m^{3}+16448\alpha^{6}m^{2}-35712\alpha^{5}m^{2}+47744\alpha^{4}m^{2}+44688\alpha^{3}m^{2}-203716\alpha^{2}m^{2}+$
$\displaystyle 174276\alpha
m^{2}-51156m^{2}+3328\alpha^{7}m-10304\alpha^{6}m+22432\alpha^{5}m+8944\alpha^{4}m-$
$\displaystyle 119864\alpha^{3}m+179416\alpha^{2}m-116424\alpha
m+28296m-7200\big{)}.$

The $\text{I}_{3}$ integral is divided into four parts

$\text{I}_{3}=A+B+C+D$ (2.70)

that correspond to integrals over the blocks of terms of (2.16), (2.17),
(2.18), and (2.19), respectively. The integrals $A$ and $B$ are readily
available based on the results of the first two moments and are given by

$\displaystyle A$ $\displaystyle=$
$\displaystyle\frac{m^{3}(2\alpha+m+1)^{3}\left(4\alpha^{2}+4\alpha+5m^{2}+10\alpha
m+5m+2\right)^{3}}{64m(m-1)(m-2)(2\alpha+2m+1)^{3}}$ (2.71) $\displaystyle B$
$\displaystyle=$
$\displaystyle-\frac{3m(2\alpha+m+1)^{2}\left(4\alpha^{2}+4\alpha+5m^{2}+10\alpha
m+5m+2\right)}{64(m-2)(m-1)(2\alpha+2m-1)(2\alpha+2m+1)^{3}(2\alpha+2m+3)}\big{(}128\alpha^{7}-64\alpha^{6}-$
(2.72) $\displaystyle
160\alpha^{5}+272\alpha^{4}-160\alpha^{3}-64\alpha^{2}+48\alpha+462m^{7}+3234\alpha
m^{6}+717m^{6}+9324\alpha^{2}m^{5}+$ $\displaystyle 3924\alpha
m^{5}+219m^{5}+14280\alpha^{3}m^{4}+8316\alpha^{2}m^{4}+546\alpha
m^{4}+243m^{4}+12432\alpha^{4}m^{3}+$ $\displaystyle
8448\alpha^{3}m^{3}+72\alpha^{2}m^{3}+1008\alpha
m^{3}-117m^{3}+6048\alpha^{5}m^{2}+3984\alpha^{4}m^{2}-768\alpha^{3}m^{2}+$
$\displaystyle 1632\alpha^{2}m^{2}-336\alpha
m^{2}-96m^{2}+1472\alpha^{6}m+576\alpha^{5}m-688\alpha^{4}m+1152\alpha^{3}m-$
$\displaystyle 328\alpha^{2}m-168\alpha m+12m\big{)}.$

By the symmetry, the integral $C$ is written as the following eight distinct
integrals

$C=\frac{1}{2m(m-1)(m-2)}\left(2C_{1}+C_{2}-2C_{3}-C_{4}+2C_{5}+C_{6}-2C_{7}-C_{8}\right),$
(2.73)

where

$\displaystyle C_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(x,y)K_{01}(y,z)K_{11}(x,z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.74) $\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(y,z)K_{01}(y,x)K_{11}(x,z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.75) $\displaystyle C_{3}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(x,y)K_{01}(x,z)K_{11}(y,z)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.76) $\displaystyle C_{4}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(y,z)K_{01}(z,x)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.77) $\displaystyle C_{5}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(x,y)K_{10}(z,x)K_{11}(z,y)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.78) $\displaystyle C_{6}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(y,z)K_{10}(x,z)K_{11}(y,x)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.79) $\displaystyle C_{7}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(x,y)K_{10}(z,y)K_{11}(z,x)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.80) $\displaystyle C_{8}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{00}(y,z)K_{10}(x,y)K_{11}(z,x)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z.$
(2.81)

Similarly, the integral $D$ is divided into

$D=\frac{1}{4m(m-1)(m-2)}\left(D_{1}+3D_{2}+3D_{3}+D_{4}\right),$ (2.82)

where

$\displaystyle D_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{01}(x,y)K_{01}(y,z)K_{01}(z,x)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.83) $\displaystyle D_{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{01}(x,z)K_{01}(z,y)K_{10}(x,y)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.84) $\displaystyle D_{3}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{01}(x,y)K_{10}(x,z)K_{10}(z,y)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z$
(2.85) $\displaystyle D_{4}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!x^{2}y^{2}z^{2}K_{10}(x,y)K_{10}(y,z)K_{10}(z,x)\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z.$
(2.86)

The integrals in $C$ and $D$ can be written in terms of biorthogonal
polynomials and the obtained recurrence relations play an important role in
the evaluation of these integrals. On the other hand, they can also be
directly calculated by integrations over the correlation kernels as
illustrated in the calculation of (2.3). In doing so, it will be efficient to
utilize the integral representation (2.31b) and (2.31c) for the kernels
$K_{01}(x,y)$ and $K_{10}(x,y)$ and the summation form (2.22a) and (2.22d) for
the kernels $K_{00}(x,y)$ and $K_{11}(x,y)$, respectively. The resulting
integrals are evaluated by using the identity [20]

$\displaystyle\int_{0}^{1}\\!x^{a-1}G_{p,q}^{m,n}\left(\left.\begin{array}[]{c}a_{1},\ldots,a_{n};a_{n+1},\ldots,a_{p}\\\
b_{1},\ldots,b_{m};b_{m+1},\ldots,b_{q}\end{array}\right|\eta
x\Big{.}\right)\,\mathrm{d}x$ (2.89) $\displaystyle=$ $\displaystyle
G_{p+1,q+1}^{m,n+1}\left(\left.\begin{array}[]{c}1-a,a_{1},\ldots,a_{n};a_{n+1},\ldots,a_{p}\\\
b_{1},\ldots,b_{m};b_{m+1},\ldots,b_{q},-a\end{array}\right|\eta\Big{.}\right)$
(2.92)

as well as the the Mellin transform of Meijer G-function (2.55). The results
are given by

$C=-\frac{3m(m+1)(2\alpha+m+1)^{2}\gamma}{8(m-2)(2\alpha+2m-3)(2\alpha+2m-1)(2\alpha+2m+1)^{3}(2\alpha+2m+3)}$
(2.93)

and

$D=\frac{\delta}{16(m-2)(m-1)(2\alpha+2m-3)(2\alpha+2m-1)(2\alpha+2m+1)^{3}},$
(2.94)

where $\gamma$ and $\delta$ are respectively

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle
32\alpha^{6}-96\alpha^{5}-40\alpha^{4}+120\alpha^{3}+8\alpha^{2}-24\alpha+15m^{6}+120\alpha
m^{5}+384\alpha^{2}m^{4}-36\alpha m^{4}-$ (2.95) $\displaystyle
48m^{4}+624\alpha^{3}m^{3}-192\alpha^{2}m^{3}-234\alpha
m^{3}+27m^{3}+536\alpha^{4}m^{2}-368\alpha^{3}m^{2}-386\alpha^{2}m^{2}+$
$\displaystyle 128\alpha
m^{2}+224\alpha^{5}m-304\alpha^{4}m-240\alpha^{3}m+208\alpha^{2}m+4\alpha
m-12m$

and

$\displaystyle\delta$ $\displaystyle=$ $\displaystyle
1024\alpha^{11}-9728\alpha^{10}+36352\alpha^{9}-64512\alpha^{8}+41280\alpha^{7}+30624\alpha^{6}-51136\alpha^{5}+1712\alpha^{4}+$
(2.96) $\displaystyle
17696\alpha^{3}-1296\alpha^{2}-2016\alpha+7293m^{11}+87516\alpha
m^{10}-12726m^{10}+463320\alpha^{2}m^{9}-$ $\displaystyle 158004\alpha
m^{9}+10539m^{9}+1424280\alpha^{3}m^{8}-844524\alpha^{2}m^{8}+140976\alpha
m^{8}-5268m^{8}+$ $\displaystyle
2814240\alpha^{4}m^{7}-2560296\alpha^{3}m^{7}+778452\alpha^{2}m^{7}-78504\alpha
m^{7}-20997m^{7}+3734016\alpha^{5}m^{6}-$ $\displaystyle
4871040\alpha^{4}m^{6}+2355888\alpha^{3}m^{6}-441336\alpha^{2}m^{6}-140292\alpha
m^{6}+22338m^{6}+3370752\alpha^{6}m^{5}-$ $\displaystyle
6061440\alpha^{5}m^{5}+4330752\alpha^{4}m^{5}-1295856\alpha^{3}m^{5}-365136\alpha^{2}m^{5}+164100\alpha
m^{5}-5331m^{5}+$ $\displaystyle
2048640\alpha^{7}m^{4}-4972224\alpha^{6}m^{4}+5023200\alpha^{5}m^{4}-2226960\alpha^{4}m^{4}-438984\alpha^{3}m^{4}+$
$\displaystyle 485940\alpha^{2}m^{4}-47952\alpha
m^{4}-8424m^{4}+808320\alpha^{8}m^{3}-2632704\alpha^{7}m^{3}+3668544\alpha^{6}m^{3}-$
$\displaystyle
2317824\alpha^{5}m^{3}-169896\alpha^{4}m^{3}+737688\alpha^{3}m^{3}-159732\alpha^{2}m^{3}-37032\alpha
m^{3}+6048m^{3}+$ $\displaystyle
192000\alpha^{9}m^{2}-848640\alpha^{8}m^{2}+1615488\alpha^{7}m^{2}-1430208\alpha^{6}m^{2}+134400\alpha^{5}m^{2}+$
$\displaystyle
594912\alpha^{4}m^{2}-248640\alpha^{3}m^{2}-52704\alpha^{2}m^{2}+25704\alpha
m^{2}+624m^{2}+23552\alpha^{10}m-$ $\displaystyle
146432\alpha^{9}m+383616\alpha^{8}m-476160\alpha^{7}m+152448\alpha^{6}m+232512\alpha^{5}m-182552\alpha^{4}m-$
$\displaystyle 24064\alpha^{3}m+36216\alpha^{2}m+576\alpha m-1008m.$

Putting everything together, one obtains the third moment formula as

$\displaystyle\mathbb{E}_{f}\\!\left[\text{P}^{3}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{(2\alpha+m+1)\prod_{i=1}^{5}(2\alpha+2m+2i-5)(m^{2}+2\alpha
m+m+2i)}\big{(}81920\alpha^{9}+$ (2.97) $\displaystyle
368640\alpha^{8}+614400\alpha^{7}+430080\alpha^{6}-1059840\alpha^{5}-2864640\alpha^{4}-1894400\alpha^{3}-$
$\displaystyle 7680\alpha^{2}+529920\alpha+500m^{13}+7000\alpha
m^{12}+3500m^{12}+43700\alpha^{2}m^{11}+$ $\displaystyle 43700\alpha
m^{11}+22225m^{11}+160400\alpha^{3}m^{10}+240600\alpha^{2}m^{10}+255900\alpha
m^{10}+$ $\displaystyle
87850m^{10}+384160\alpha^{4}m^{9}+768320\alpha^{3}m^{9}+1294380\alpha^{2}m^{9}+910220\alpha
m^{9}+$ $\displaystyle
256845m^{9}+629600\alpha^{5}m^{8}+1574000\alpha^{4}m^{8}+3783400\alpha^{3}m^{8}+4101100\alpha^{2}m^{8}+$
$\displaystyle 2386800\alpha
m^{8}+562350m^{8}+719616\alpha^{6}m^{7}+2158848\alpha^{5}m^{7}+7068480\alpha^{4}m^{7}+$
$\displaystyle 10538880\alpha^{3}m^{7}+9593004\alpha^{2}m^{7}+4683372\alpha
m^{7}+503935m^{7}+572928\alpha^{7}m^{6}+$ $\displaystyle
2005248\alpha^{6}m^{6}+8807232\alpha^{5}m^{6}+17004960\alpha^{4}m^{6}+21867552\alpha^{3}m^{6}+$
$\displaystyle 16798992\alpha^{2}m^{6}+3491468\alpha
m^{6}-535010m^{6}+311040\alpha^{8}m^{5}+1244160\alpha^{7}m^{5}+$
$\displaystyle
7396992\alpha^{6}m^{5}+17836416\alpha^{5}m^{5}+31093200\alpha^{4}m^{5}+33910560\alpha^{3}m^{5}+$
$\displaystyle 10570948\alpha^{2}m^{5}-3618956\alpha
m^{5}-2077165m^{5}+109568\alpha^{9}m^{4}+493056\alpha^{8}m^{4}+$
$\displaystyle
4122624\alpha^{7}m^{4}+12128256\alpha^{6}m^{4}+28475808\alpha^{5}m^{4}+42019344\alpha^{4}m^{4}+$
$\displaystyle 18278216\alpha^{3}m^{4}-8866596\alpha^{2}m^{4}-10838736\alpha
m^{4}-2813290m^{4}+22528\alpha^{10}m^{3}+$ $\displaystyle
112640\alpha^{9}m^{3}+1452288\alpha^{8}m^{3}+5133312\alpha^{7}m^{3}+16659456\alpha^{6}m^{3}+32484864\alpha^{5}m^{3}+$
$\displaystyle
19378192\alpha^{4}m^{3}-9621472\alpha^{3}m^{3}-21137424\alpha^{2}m^{3}-11734304\alpha
m^{3}-1622340m^{3}+$ $\displaystyle
2048\alpha^{11}m^{2}+11264\alpha^{10}m^{2}+289792\alpha^{9}m^{2}+1219584\alpha^{8}m^{2}+5914752\alpha^{7}m^{2}+$
$\displaystyle
15089088\alpha^{6}m^{2}+12371840\alpha^{5}m^{2}-4003744\alpha^{4}m^{2}-19106720\alpha^{3}m^{2}-$
$\displaystyle 17907952\alpha^{2}m^{2}-4858272\alpha
m^{2}+102600m^{2}+24576\alpha^{10}m+122880\alpha^{9}m+$ $\displaystyle
1124352\alpha^{8}m+3760128\alpha^{7}m+4305408\alpha^{6}m+271872\alpha^{5}m-7811712\alpha^{4}m-$
$\displaystyle 11935488\alpha^{3}m-5078784\alpha^{2}m+344448\alpha
m+496800m+172800\big{)},$

which upon inserting $\alpha=n-m-1/2$ leads to the claimed result (1.15).
Finally, we note that the higher order moments beyond the first three can also
be obtained in a similar manner but with an increased effort. Another
potential approach to obtain the exact moments is via the method of integrable
systems, where, for example, it was established in [24] that the partition
function of Bures-Hall ensemble is the $\tau$-function of certain integrable
hierarchies.

##  3\. Conclusion and outlook

In this work, we derived the exact second and third moments of quantum purity
of arbitrary subsystem dimensions over the Bures-Hall ensemble, generalizing
the results of equal subsystem dimensions in the literature. The derivations
are obtained based on the direct integrations over the correlation kernels and
the established recurrence relations of integrals involving biorthogonal
polynomials. Future work includes devising an integrable system approach to
systematically produce moments of any order as well as the study of asymptotic
distribution of purity of large subsystem dimensions.

## Acknowledgments

We thank Santosh Kumar for providing the simulation codes and Jiyuan Zhang for
correspondence. The work of Lu Wei is supported in part by the U.S. National
Science Foundation (CNS-2006612).

##  Appendix A Normalizations of Bures-Hall $k$-point densities

In this appendix, we verify the normalization constants of the first three
point densities (2.13), (2.14), and (2.15) of the Bures-Hall ensemble, which
in particular makes use of the fact that the kernels $K_{01}(x,y)$,
$K_{10}(x,y)$, and $K_{11}(x,y)$ can be written in terms of the kernel
$K_{00}(x,y)$ as [17]

$\displaystyle K_{01}(x,y)$ $\displaystyle=$ $\displaystyle
x^{\alpha}e^{-x}\int_{0}^{\infty}\frac{v^{\alpha+1}{\rm
e}^{-v}}{x+v}K_{00}(y,v)\,\mathrm{d}v$ (A.1a) $\displaystyle K_{10}(x,y)$
$\displaystyle=$ $\displaystyle
y^{\alpha+1}e^{-y}\int_{0}^{\infty}\frac{w^{\alpha}{\rm
e}^{-w}}{y+w}K_{00}(w,x)\,\mathrm{d}w$ (A.1b) $\displaystyle K_{11}(x,y)$
$\displaystyle=$ $\displaystyle
x^{\alpha}y^{\alpha+1}e^{-x-y}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{v^{\alpha}{\rm
e}^{-v}}{y+v}\frac{w^{\alpha+1}{\rm
e}^{-w}}{x+w}K_{00}(v,w)\,\mathrm{d}v\,\mathrm{d}w-W(x,y).$ (A.1c)

We first note the normalization of $K_{00}(x,y)$,

$\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!K_{00}(x,y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=\sum_{k=0}^{m-1}1=m$
(A.2)

and the reproducing properties of the correlation kernels

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!K_{00}(x,z)K_{00}(w,y)W(w,z)\,\mathrm{d}w\,\mathrm{d}z$
$\displaystyle=$ $\displaystyle K_{00}(x,y)$ (A.3)
$\displaystyle\int_{0}^{\infty}K_{01}(x,y)K_{01}(y,z)\,\mathrm{d}y$
$\displaystyle=$ $\displaystyle K_{01}(x,z)$ (A.4)
$\displaystyle\int_{0}^{\infty}K_{10}(x,y)K_{10}(y,z)\,\mathrm{d}y$
$\displaystyle=$ $\displaystyle K_{10}(x,z).$ (A.5)

For the one-point density (2.13), the normalization constant $2m$ is obtained
as

$\displaystyle\int_{0}^{\infty}K_{01}(x,x)\,\mathrm{d}x$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{\alpha}e^{-x}\frac{v^{\alpha+1}{\rm
e}^{-v}}{x+v}K_{00}(x,v)\,\mathrm{d}v\,\mathrm{d}x$ (A.6) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(x,v)W(x,v)\,\mathrm{d}v\,\mathrm{d}x=m$
(A.7) $\displaystyle\int_{0}^{\infty}K_{10}(x,x)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{\alpha+1}e^{-x}\frac{w^{\alpha}{\rm
e}^{-w}}{x+w}K_{00}(w,x)\,\mathrm{d}w\,\mathrm{d}x$ (A.8) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(w,x)W(w,x)\,\mathrm{d}w\,\mathrm{d}x=m.$
(A.9)

To verify the normalization constant $4m(m-1)$ of the two-point density
(2.14), we first note that

$\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}(\left(K_{01}(x,x)+K_{10}(x,x)\right)\left(K_{01}(y,y)+K_{10}(y,y)\right)\,\mathrm{d}x\,\mathrm{d}y=4m^{2}.$
(A.10)

We now compute

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{01}(x,y)K_{01}(y,x)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{\alpha}e^{-x}\frac{v^{\alpha+1}{\rm
e}^{-v}}{x+v}K_{00}(y,v)y^{\alpha}e^{-y}\frac{w^{\alpha+1}{\rm
e}^{-w}}{y+w}K_{00}(x,w)\,\mathrm{d}v\,\mathrm{d}w\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(y,w)y^{\alpha}e^{-y}\frac{w^{\alpha+1}{\rm
e}^{-w}}{y+w}\,\mathrm{d}w\,\mathrm{d}y=m.$ (A.12)

Similarly, we have

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{10}(x,y)K_{10}(y,x)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}y^{\alpha+1}e^{-y}\frac{w^{\alpha}{\rm
e}^{-w}}{y+w}K_{00}(w,x)x^{\alpha+1}e^{-x}\frac{v^{\alpha}{\rm
e}^{-v}}{x+v}K_{00}(v,y)\,\mathrm{d}w\,\mathrm{d}v\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}y^{\alpha+1}e^{-y}\frac{w^{\alpha}{\rm
e}^{-w}}{y+w}K_{00}(w,y)\,\mathrm{d}w\,\mathrm{d}y=m.$ (A.14)

Finally, one computes

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(x,y)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y\left(=\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(y,x)K_{11}(y,x)\,\mathrm{d}x\,\mathrm{d}y\right)$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(x,y)\left(x^{\alpha}y^{\alpha+1}e^{-x-y}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{v^{\alpha}{\rm
e}^{-v}}{y+v}\frac{w^{\alpha+1}{\rm
e}^{-w}}{x+w}K_{00}(v,w)\,\mathrm{d}v\,\mathrm{d}w-W(x,y)\right)\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\sum_{i,j=0}^{m-1}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\widetilde{p}_{i}(x)\widetilde{q}_{j}(w)\frac{x^{\alpha}w^{\alpha+1}{\rm
e}^{-x-w}}{x+w}\,\mathrm{d}x\,\mathrm{d}w\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\widetilde{p}_{j}(v)\widetilde{q}_{i}(y)\frac{v^{\alpha}y^{\alpha+1}{\rm
e}^{-v-y}}{v+y}\,\mathrm{d}v\,\mathrm{d}y-m$ $\displaystyle=$
$\displaystyle\sum_{i=j=0}^{m-1}1-m=0,$

where for convenience we used the non-monic version [19] of the biorthogonal
polynomials $\widetilde{p}_{k}(x)$ and $\widetilde{q}_{l}(y)$, i.e.,

$K_{00}(x,y)=\sum_{k=0}^{m-1}\widetilde{p}_{k}(x)\widetilde{q}_{k}(y)$ (A.15)

and

$\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\widetilde{p}_{k}(x)\widetilde{q}_{l}(y)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=\delta_{kl}.$
(A.16)

We now arrive at the normalization constant

$4m^{2}-2m-2m-0=4m(m-1).$ (A.17)

Before discussing the $3$-point density normalization, we compute two
additional integral identities. The first one is

$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{01}(x,y)K_{10}(x,y)\,\mathrm{d}x\,\mathrm{d}y\left(=\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{01}(y,x)K_{10}(y,x)\,\mathrm{d}x\,\mathrm{d}y\right)$
(A.18) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{\alpha}{\rm
e}^{-x}\frac{v^{\alpha+1}{\rm e}^{-v}}{x+v}K_{00}(y,v)y^{\alpha+1}{\rm
e}^{-y}\frac{w^{\alpha}{\rm
e}^{-w}}{y+w}K_{00}(w,x)\,\mathrm{d}v\,\mathrm{d}w\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\sum_{i,j=0}^{m-1}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}x^{\alpha}{\rm
e}^{-x}\frac{v^{\alpha+1}{\rm
e}^{-v}}{x+v}\widetilde{p}_{i}(y)\widetilde{q}_{i}(v)y^{\alpha+1}{\rm
e}^{-y}\frac{w^{\alpha}{\rm
e}^{-w}}{y+w}\widetilde{p}_{j}(w)\widetilde{q}_{j}(x)\,\mathrm{d}v\,\mathrm{d}w\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\sum_{i,j=0}^{m-1}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{x^{\alpha}v^{\alpha+1}{\rm
e}^{-x-v}}{x+v}\widetilde{q}_{j}(x)\widetilde{q}_{i}(v)\,\mathrm{d}x\,\mathrm{d}v\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{w^{\alpha}y^{\alpha+1}{\rm
e}^{-w-y}}{w+y}\widetilde{p}_{j}(w)\widetilde{p}_{i}(y)\,\mathrm{d}w\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\sum_{i,j=0}^{m-1}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{x^{\alpha}v^{\alpha+1}{\rm
e}^{-x-v}}{x+v}(\widetilde{p}_{j}(x)+\text{l.o.t.})\widetilde{q}_{i}(v)\,\mathrm{d}x\,\mathrm{d}v$
$\displaystyle\times\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{w^{\alpha}y^{\alpha+1}{\rm
e}^{-w-y}}{w+y}\widetilde{p}_{j}(w)(\widetilde{q}_{i}(y)+\text{l.o.t.})\,\mathrm{d}w\,\mathrm{d}y$
$\displaystyle=$ $\displaystyle\sum_{i<j=0}^{m-1}\text{(non-zero
contribution)}\times 0+\sum_{i>j=0}^{m-1}0\times\text{(non-zero
contribution)}+\sum_{i=j=0}^{m-1}1$ $\displaystyle=$ $\displaystyle m,$

where we have used the fact that the coefficients of the highest degree terms
of the polynomials $\widetilde{p}_{j}(x)$ and $\widetilde{q}_{j}(y)$ are the
same [19]. By the relations (2.20), the correlation kernels satisfy [19]

$\displaystyle
K_{01}(x,y)K_{01}(y,x)+K_{10}(x,y)K_{10}(y,x)+K_{00}(x,y)K_{11}(x,y)+K_{00}(y,x)K_{11}(y,x)$
$\displaystyle=$ $\displaystyle
K_{01}(x,y)K_{10}(x,y)+K_{01}(y,x)K_{10}(y,x)-K_{00}(x,y)K_{11}(y,x)-K_{00}(y,x)K_{11}(x,y),$

which, with the previously obtained results, leads to the identities

$\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(x,y)K_{11}(y,x)\,\mathrm{d}x\,\mathrm{d}y=\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(y,x)K_{11}(x,y)\,\mathrm{d}x\,\mathrm{d}y=0.$
(A.19)

To verify the normalization constant $8m(m-1)(m-2)$ of the three-point density
(2.15), we start from the integrals over the terms (2.16) and (2.17), which
are given by $(2m)^{3}=8m^{3}$ and $3\times(-2\times 2m\times 2m)=-24m^{2}$,
respectively. To compute integrals of the terms (2.18), we first compute some
intermediate results. One has

$\displaystyle\int_{0}^{\infty}K_{00}(x,y)K_{11}(z,y)\,\mathrm{d}y$ (A.20)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}K_{00}(x,y)\left(z^{\alpha}y^{\alpha+1}e^{-z-y}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{v^{\alpha}{\rm
e}^{-v}}{y+v}\frac{w^{\alpha+1}{\rm
e}^{-w}}{z+w}K_{00}(v,w)\,\mathrm{d}v\,\mathrm{d}w-W(z,y)\right)\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{z^{\alpha}w^{\alpha+1}{\rm
e}^{-z-w}}{z+w}K_{01}(v,x)K_{00}(v,w)\,\mathrm{d}v\,\mathrm{d}w-\int_{0}^{\infty}K_{00}(x,y)W(z,y)\,\mathrm{d}y$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}K_{01}(z,v)K_{01}(v,x)\,\mathrm{d}v-K_{01}(z,x)=K_{01}(z,x)-K_{01}(z,x)=0$

and

$\displaystyle\int_{0}^{\infty}K_{00}(x,y)K_{11}(x,z)\,\mathrm{d}x$ (A.21)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}K_{00}(x,y)\left(x^{\alpha}z^{\alpha+1}e^{-x-z}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{v^{\alpha}{\rm
e}^{-v}}{z+v}\frac{w^{\alpha+1}{\rm
e}^{-w}}{x+w}K_{00}(v,w)\,\mathrm{d}v\,\mathrm{d}w-W(x,z)\right)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(v,w)\frac{v^{\alpha}z^{\alpha+1}{\rm
e}^{-v-z}}{v+z}\,\mathrm{d}v\int_{0}^{\infty}K_{00}(x,y)\frac{x^{\alpha}w^{\alpha+1}{\rm
e}^{-x-w}}{x+w}\,\mathrm{d}x\,\mathrm{d}w-$
$\displaystyle\int_{0}^{\infty}K_{00}(x,y)W(x,z)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}K_{10}(w,z)K_{10}(y,w)\,\mathrm{d}w-K_{10}(y,z)=K_{10}(y,z)-K_{10}(y,z)=0,$

which respectively lead to

$0=\int_{0}^{\infty}K_{00}(x,y)K_{11}(z,y)\,\mathrm{d}y\int_{0}^{\infty}\frac{x^{\alpha}v^{\alpha+1}{\rm
e}^{-x-v}}{x+v}\,\mathrm{d}x=\int_{0}^{\infty}K_{10}(y,v)K_{11}(z,y)\,\mathrm{d}y$
(A.22)

and

$0=\int_{0}^{\infty}K_{00}(x,y)K_{11}(x,z)\,\mathrm{d}x\int_{0}^{\infty}\frac{v^{\alpha}y^{\alpha+1}{\rm
e}^{-v-y}}{v+y}\,\mathrm{d}y=\int_{0}^{\infty}K_{01}(v,x)K_{11}(x,z)\,\mathrm{d}x.$
(A.23)

One also has

$\displaystyle\int_{0}^{\infty}K_{00}(x,y)K_{01}(x,z)\,\mathrm{d}x$ (A.24)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(x,y)\frac{x^{\alpha}v^{\alpha+1}{\rm
e}^{-x-v}}{x+v}K_{00}(z,v)\,\mathrm{d}x\,\mathrm{d}v$ $\displaystyle=$
$\displaystyle\sum_{i,j=0}^{m-1}\widetilde{p}_{j}(z)\widetilde{q}_{i}(y)\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\widetilde{p}_{i}(x)\widetilde{q}_{j}(v)\frac{x^{\alpha}v^{\alpha+1}{\rm
e}^{-x-v}}{x+v}\,\mathrm{d}x\,\mathrm{d}v$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{m-1}\widetilde{p}_{i}(z)\widetilde{q}_{i}(y)=K_{00}(z,y)$

and

$\displaystyle\int_{0}^{\infty}K_{00}(x,y)K_{10}(z,y)\,\mathrm{d}y$ (A.25)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}K_{00}(x,y)\frac{w^{\alpha}y^{\alpha+1}{\rm
e}^{-w-y}}{w+y}K_{00}(w,z)\,\mathrm{d}w\,\mathrm{d}y$ $\displaystyle=$
$\displaystyle\sum_{i,j=0}^{m-1}\widetilde{p}_{i}(x)\widetilde{q}_{j}(z)\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\widetilde{p}_{j}(w)\widetilde{q}_{i}(y)\frac{w^{\alpha}y^{\alpha+1}{\rm
e}^{-w-y}}{w+y}\,\mathrm{d}w\,\mathrm{d}y$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{m-1}\widetilde{p}_{i}(x)\widetilde{q}_{i}(z)=K_{00}(x,z).$

With the above integral identities, it can be verified that each of the $24$
terms in (2.18) has a zero contribution to the normalization. At the same
time, each of the $8$ terms in (2.19) integrates to $m$ by using the
reproducing properties (A.4) and (A.5) as well as the identity (A.18). The
normalization constant is therefore obtained as

$8m^{3}-24m^{2}+2\times 8m=8m(m-1)(m-2).$ (A.26)

## References

  * [1] Hall M J W 1998 Random quantum correlations and density operator distributions Phys. Lett. A 242 123-9
  * [2] Życzkowski K and Sommers H J 2001 Induced measures in the space of mixed quantum states J. Phys. A: Math. Gen. 34 7111-25
  * [3] Sommers H-J and Życzkowski K 2003 Random Bures miexed states and the distribution of their purity J. Phys. A: Math. Gen. 36 10083-100
  * [4] Wei L 2020 Skewness of von Neumann entanglement entropy J. Phys. A: Math. Theor. 53 075302
  * [5] Wei L 2020 Exact variance of von Neumann entanglement entropy over the Bures-Hall measure Phys. Rev. E 102 062128
  * [6] Sommers H J and Życzkowski K 2004 Statistical properties of random density matrices J. Phys. A: Math. Gen. 37 8457-66
  * [7] Osipov V A, Sommers H J and Życzkowski K 2010 Random Bures mixed states and the distribution of their purity J. Phys. A: Math. Theor. 43 055302
  * [8] Sarkar A and Kumar S 2019 Bures-Hall ensemble: Spectral densities and average entropies J. Phys. A: Math. Theor. 52 295203
  * [9] Wei L 2020 Proof of Sarkar-Kumar’s conjectures on average entanglement entropies over the Bures-Hall ensemble J. Phys. A: Math. Theor. 53 235203
  * [10] von Neumann J 1927 Wahrscheinlichkeitstheoretischer aufbau der quantenmechanik Nachr. Ges. Wiss. Göttingen 1 245
  * [11] Bengtsson I and Życzkowski K 2017 Geometry of Quantum States: An Introduction to Quantum Entanglement 2nd ed (Cambridge: Cambridge University Press)
  * [12] Życzkowski K and Sommers H-J 2005 Average fidelity between random quantum states Phys. Rev. A 71 032313
  * [13] Islam R, Ma R, Preiss P M, Tai M E, Lukin A, Rispoli M and Greiner M 2015 Measuring entanglement entropy in a quantum many-body system Nature 528 77
  * [14] Borot G and Nadal C 2012 Purity distribution for generalized random Bures mixed states J. Phys. A: Math. Theor. 45 075209
  * [15] Wei L 2017 Proof of Vivo-Pato-Oshanin’s conjecture on the fluctuation of von Neumann entropy Phys. Rev. E 96 022106
  * [16] Wei L and Witte N 2021 Quantum interpolating ensemble: Average entropies and orthogonal polynomials (arXiv:2103.04231)
  * [17] Forrester P J and Kieburg M 2016 Relating the Bures measure to the Cauchy two-matrix model Commun. Math. Phys. 342 151-87
  * [18] Bertola M, Gekhtman M and Szmigielski J 2010 Cauchy biorthogonal polynomials J. Approx. Theory 162 832-67
  * [19] Bertola M, Gekhtman M and Szmigielski J 2014 Cauchy-Laguerre two-matrix model and the Meijer-G random point field Commun. Math. Phys. 326 111-44
  * [20] Prudnikov A P, Brychkov Y A and Marichev O I 1986 Integrals and Series: More special functions, Integrals and Series vol 3 (London: Gordon and Breach Science Publishers)
  * [21] Mehta M L 2004 Random Matrices 3rd ed (Amsterdam: Elsevier)
  * [22] Forrester P J 2010 Log-gases and Random Matrices (Princeton: Princeton University Press)
  * [23] Hardy A 2018 Polynomial ensembles and recurrence coefficients Constr. Approx. 48 137-62
  * [24] Hu X-B and Li S-H 2017 The partition function of the Bures ensemble as the $\tau$-function of BKP and DKP hierarchies: continuous and discrete J. Phys. A: Math. Theor. 50 285201