# Diffuse $\gamma$-ray emission in Cygnus X: Comments to Yan&Pavaskar
https://arxiv.org/abs/2211.17057

Ottavio Fornieri  Gran Sasso Science Institute, Viale F. Crispi 7, 67100
L’Aquila, Italy  INFN - Laboratori Nazionali del Gran Sasso (LNGS), Via G.
Acitelli 22, 67100 Assergi (AQ), Italy  Heshou Zhang  INAF - Osservatorio
Astronomico di Brera, Via E. Bianchi 46, 23807 Merate (LC), Italy.

###### Abstract

Abstract

It has been pointed out in https://arxiv.org/abs/2211.17057 [1] that our
recent (published) paper [2] might be revised, due to an incorrect evaluation
of the diffusion coefficients, $D(E)$, employed in the calculations.
Unfortunately, there is no indication as to where the incorrectness might be.
Here we offer the opportunity to be more specific, by providing the community
with the whole description of the equations involved in the calculation of
$D(E)$, which is missing in the arXiv note [1]. In this context, we mention
that no ad hoc parameterisation has been used in our paper. Furthermore, any
assumption on the injection mechanisms is explicitly described in the paper as
an input factor and is obviously part of the modelisation procedure, hence the
final outcome is subject to it. Finally, we discuss in a more broad context
what, in this calculation of the diffusion coefficient, we believe is the key
message.

MHD — Magnetosonic modes — Cosmic-ray diffusion

##  I Introduction

To compute the diffusion coefficient that cosmic rays (CRs) feel, $D(E)$, as a
function of their energy, it is necessary to estimate the interaction between
the particle — represented by its Larmor radius $r_{L}$ — and the turbulent
waves, which act as scattering centers for the CRs. These waves, depending on
the environment where they are produced, can propagate undisturbed — therefore
contributing to the diffusion coefficient of the particles — or can be damped.
This is formalized in the following equation, to compute the pitch-angle
scattering of the CRs, $D_{\mu\mu}$, where the integral bounds up to a maximum
wavenumber, $k_{\max}$, account for the damping of the fluctuations at a
certain scale [3, 4]:

$\displaystyle D_{\mu\mu}$ $\displaystyle=\Omega^{2}(1-\mu^{2})\int
d^{3}\bm{k}\sum^{+\infty}_{n=-\infty}R_{n}(k_{\parallel}v_{\parallel}-\omega+n\Omega)\times$
(1) $\displaystyle\times\left[\frac{n^{2}J_{n}^{2}(z)}{z^{2}}I^{\rm
A}(\bm{k})+\frac{k^{2}_{\parallel}}{k^{2}}J^{\prime 2}_{n}(z)I^{\rm
M}(\bm{k})\right],$

where $\Omega=qB_{0}/m\gamma c$ is the particle’s Larmor frequency; $\bm{k}$
is the wave-vector of the turbulent fluctuations;
$k_{\parallel}\equiv|\bm{k}|\cos{\alpha_{w}}$ is its field-aligned component
($\alpha_{w}$ being the angle between the wave vector and the direction of the
background magnetic field, $B_{0}$); $\omega=\omega(\bm{k})$ the associated
fluctuations’ frequency. The details of the calculation are described in the
Appendix of Fornieri  _et al._ [5].

Several damping mechanisms have been considered in the literature. Given the
typical environmental parameters considered in our Galaxy, only the
collisionless [6] and the collisional (viscous) [7] dampings have been
isolated and are at play [8, 5].

The parameters of the plasma involved in the calculation of the scale at which
turbulent fluctuations are damped (in wavenumber $k_{\mathrm{max}}$,
normalized to the turbulence injection scale,
$k_{\mathrm{max}}L_{\mathrm{inj}}$) — and therefore cannot act as CR
scattering centers — depend on physical quantities characteristic of the
Cygnus-X environment. In particular, we use the following values, as listed in
our paper:

$\displaystyle T_{\mathrm{Cyg}}=5\cdot 10^{3}\,\mathrm{K}$ (2) $\displaystyle
n_{\mathrm{Cyg}}=10^{-1}\,\mathrm{cm^{-3}}$ $\displaystyle
B_{\mathrm{Cyg}}=10\,\mu\mathrm{G}$ $\displaystyle
L_{\mathrm{inj}}=10\,\mathrm{pc}.$

To compute the truncation scale of the fluctuations, we equate the cascading
timescale, for the fast magnetosonic modes, with the (specific) damping rate
[9]. In particular, provided that for the typical Interstellar Medium (ISM)
conditions $\beta\ll 1$ everywhere, we have:

$\displaystyle
k_{\mathrm{max}}L_{\mathrm{inj}}=x_{c}\,(1-\xi^{2})^{-2/3},\qquad\beta\ll 1$
(3) $\displaystyle
k_{\mathrm{max}}L_{\mathrm{inj}}=\frac{4\,M_{\mathrm{A}}^{4}\,\frac{m_{p}}{m_{e}}\,\xi^{2}}{\pi\beta\,(1-\xi^{2})^{2}}\cdot\exp\left(\frac{2}{\beta\,\frac{m_{p}}{m_{e}}\,\xi^{2}}\right),$

where $\xi\equiv\cos\alpha_{w}$. The equations show the truncation scales for
the viscous and collisionless dampings, respectively, while the other
quantities are discussed below. From here on, $L=L_{\mathrm{inj}}$
interchangeably, to avoid heavy notation.

The plasma parameters for the viscous case can be rearranged as follows:

$\begin{aligned} &x_{c}=\left(\frac{6\rho\delta
u^{2}L_{\mathrm{inj}}}{\eta_{0}v_{A}}\right)^{2/3}=\left(\frac{6\rho
M_{A}^{2}L_{\mathrm{inj}}\,v_{A}}{\eta_{0}}\right)^{2/3},\qquad
M_{A}\equiv\frac{\delta u}{v_{A}}\bigg{|}_{\mathrm{inj}}\\\\[10.0pt] &-\hskip
5.69046pt\rho=n_{\mathrm{Cyg}}\cdot m_{p}\\\ &-\hskip 5.69046pt\eta_{0}=6\cdot
10^{3}\left(\frac{37}{\ln\Lambda}\right)\left(\frac{T_{\mathrm{Cyg}}}{10^{8}\,\mathrm{K}}\right)^{5/2}\,\mathrm{\frac{g}{cm\cdot
s}}\\\ &\hskip
15.6491pt\ln\Lambda\equiv\ln\left(\frac{b_{\mathrm{max}}}{b_{\mathrm{min}}}\right)\simeq\ln\left(\frac{\lambda_{D}}{b_{\mathrm{min}}}\right),\;\;b_{\mathrm{min}}\approx
8.63\cdot
10^{-8}\left(\frac{10^{4}\,\mathrm{K}}{T_{\mathrm{Cyg}}}\right)\,\mathrm{cm}\\\
&-\hskip 5.69046pt\lambda_{D}\approx 0.95\cdot
10^{5}\left(\frac{T_{\mathrm{Cyg}}}{10^{6}\,\mathrm{K}}\right)^{1/2}\left(\frac{10^{-3}\,\mathrm{cm^{-3}}}{n_{\mathrm{Cyg}}}\right)^{1/2}\,\mathrm{cm}\end{aligned}\qquad\qquad\qquad\begin{aligned}
&\\\\[25.0pt] &\textrm{density of the medium}\\\\[10.0pt] &\textrm{Spitzer
1962~{}\cite[citep]{[\@@bibref{Number}{1962pfig.book.....S,
1965RvPP....1..205B}{}{}]}}\\\\[16.0pt] &\textrm{Coulomb
logarithm}\\\\[18.0pt] &\textrm{Debye lenght}.\end{aligned}$ (4)

The details of these equations can be found in Spitzer [10], Braginskii [7],
Ginzburg  _et al._ [6].

The plasma parameters for the collisionless case are:

$\displaystyle\beta\equiv\frac{P_{\mathrm{th}}}{P_{B}}=3.3\cdot\left(\frac{3\,\mu\mathrm{G}}{B_{\mathrm{Cyg}}}\right)^{2}\left(\frac{n_{\mathrm{Cyg}}}{1\,\mathrm{cm^{-3}}}\right)\left(\frac{T_{\mathrm{Cyg}}}{10^{4}\,\mathrm{K}}\right)$
(5) $\displaystyle v_{A}=6.27\cdot
10^{5}\left(\frac{B_{\mathrm{Cyg}}}{3\,\mu\mathrm{G}}\right)\left(\frac{1\,\mathrm{cm^{-3}}}{n_{\mathrm{Cyg}}}\right)^{1/2}\,\mathrm{cm\cdot
s^{-1}}$

where $\beta$ is the ratio between the thermal pressure and the magnetic
pressure in the ISM, and $v_{A}$ is the usual Alfvén speed.

These parameterisations have been used to compute the diffusion coefficients
in the paper in question [2] and in Fornieri  _et al._ [5], and for both of
them have been checked multiple times by multiple persons, other than the
reviewers, obviously. In any case, this is what should be addressed to find
any (certainly possible) problems in the calculation, with an independent
estimation. Again, we would like to stress that no ad hoc parameterisation has
been involved, unlike what claimed in the arXiv note [1].

One thing that is worth noticing is the truncation scales resulting from the
two damping mechanisms computed with the listed values, that can be seen in
Figure 1, where we show them (a) for the case $M_{A}=0.5$ (Alfvén-dominated
region in Fornieri and Zhang [2]) and (b) for the case $M_{A}=1$ (fast-
dominated region in Fornieri and Zhang [2]). In the figure we also show
(dashed lines) the $k_{\mathrm{max}}L$ scales (of the cascade) corresponding
to fluctuations resonating with $100\,\mathrm{GeV}$ (red dashed) and
$10,\mathrm{TeV}$ (green dashed) CRs, based on their Larmor radii $r_{L}$,
being $E_{\mathrm{CR}}\in[100\,\mathrm{GeV},\,10\,\mathrm{TeV}]$ the energy
range of the particles considered in our paper. As it can be seen, a little
change in the properties of the cascade (its $M_{A}$) makes, e.g.,
$10\,\mathrm{TeV}$ CRs scatter off all the fluctuations of the fast cascade in
the $M_{A}=1$ case, while less waves are available to scatter off, for the
same CRs, in the $M_{A}=0.5$ case. This difference makes our two diffusion
coefficients shown in Figure 1a of Fornieri and Zhang [2] to be different by
quite large amounts ($\sim 5$ at $E_{\mathrm{CR}}\simeq 10\,\mathrm{TeV}$ and
$\sim 7.5$ at $E_{\mathrm{CR}}\simeq 100\,\mathrm{GeV}$).

(a)

(b)

Figure 1: The truncation scales of the fluctuations as a function of the wave
pitch-angle, due to viscous (solid orange lines) and collisionless (solid blue
line) dampings are shown, (a) for the case $M_{A}=0.5$ and (b) for the case
$M_{A}=1$. The dashed horizontal lines correspond to the values of
$k_{\mathrm{max}}L$ that resonate with Larmor radii of $100\,\mathrm{GeV}$
particles (dashed red) and of $10\,\mathrm{TeV}$ particles (dashed green). As
an example, it can be seen that, in the case $M_{A}=1$, $10\,\mathrm{TeV}$ CRs
can scatter against all the waves, since they have not been damped by the
medium, while this is not true anymore where $M_{A}=0.5$, since they can only
scatter off waves with a pitch angle up to $\sim 50^{\circ}$. More scattering
implies obviously more confinement, namely a larger $D_{\mu\mu}$ and a smaller
(spatial) $D(E)$.

##  II Effect of a different normalisation on the resulting $\gamma$-ray
emission

As a matter of fact, it has been pointed out in our paper [2] (Section IIb,
page 5, right column) that the final result is, in terms of the diffusion
properties of the CRs, nearly exclusively due to the relative normalisations
of the diffusion coefficients in the two regions, up to the point that even a
different rigidity scaling does not affect it significantly. This is due to
the significance maps that we are plotting to show the $\gamma$-ray
distribution, that involve relative differences. The referee raised this point
and was satisfied with this answer.

A different matter would be if the energetics of the sources were analyzed in
the paper, to study the broken power-law spectrum shown by Fermi-LAT [11] \+
HAWC [12]: in that case, a careful evaluation of the age of the sources and
the CR injection spectrum would have been required. However, we found it to be
pointless, due to the high level of uncertainties involved in those
estimations, which would allow us to properly readjust the environmental
parameters, still within the reasonable uncertainty band.

##  III Effect of a different injection mechanism on the resulting
$\gamma$-ray emission

For what concerns the injection mechanism, we explicitly describe in the paper
that $\sim 90\%$ of the particle emission occurs as an initial burst (Section
IIb, page 4, left column). We believe that this is physically motivated by two
main reasons:

  * •

the acceleration of CRs occurs inside the bubble in a reverse shock, so that
the region upstream of the shock can trap the particles inside the bubble
until dissipation of the shock itself, a process that can take as long as
$\sim\mathcal{O}(\mathrm{Myr})$ [13], after which the bulk of the accelerated
CRs can be released into the ISM,

  * •

there is evidence for a $\sim 1\big{/}r$ radial profile of the surface
brightness detected from the OB cluster [14], which should indicate a
continuous injection of CR particles from the region during (at least) its
most recent stage — this is a standard calculation: see, e.g. De La Torre
Luque  _et al._ [15], Appendix B.

We implemented the above considerations as described in the paper.
Furthermore, we remark that this is part of the modeling procedure, and as
such is considered to obtain the expected outcome.

##  IV General discussion

It is probably worth at this point to clarify what are — in our opinion — the
key points to learn from any calculation of the diffusion coefficient in these
contexts. As shown in Equation (4), the calculation of the damping-parameter
scales, which is of primary importance in this framework (see e.g. Figure 1
and its discussion above), is based on the results from previous research.
These quantities could be updated when we have a better knowledge of the
composition of the background plasma in the ISM. Furthermore, to compute the
truncation scales in Equation (3), we had to estimate the developement
timescale of the fast magnetosonic cascade, which involves the estimation of
the fluid velocity of the turbulence, at each scale — see, e.g., Schekochihin
[16] for a recent review on the topic, where the scalings of the magneto-
hydro-dynamic (MHD) cascades are discussed. Besides, to the best of our
knowledge, there is no analytic theory on the MHD compressible cascade, but
only on the Alfvénic one [16]. This quantities, therefore, should be taken
with a grain of salt.

With this in mind, what we believe is really relevant in the context of our
paper — and explicitly stated in it — is the difference in the propagation
properties of the CRs between two regions where different modes dominate the
diffusion process: (i) where Alfvén modes dominate in the background plasma,
particles are less trapped and propagate faster, (ii) where fast modes
dominate, particles are confined for longer. This has the implications in the
$\gamma$-ray emission of the Cygnus-X region that are described in our paper
[2]. In particular, the diffuse $\gamma$-ray emission in the region is
correlated with such difference, rather than with the baryonic-gas
distribution. Within reasonable physical uncertainties, all the parameters can
be, instead, appropriately tuned. This is the reason why we believe an
accurate analysis of the energetics of the sources might not be based on solid
grounds.

We remain open to possible quantitative criticism on the results of the paper
Fornieri and Zhang [2], which we certainly have already received during the
review process from the Journal. Until then, we sincerely hope that everything
is clearer, now.

## References

  * Yan and Pavaskar [2022] H. Yan and P. Pavaskar, (2022), arXiv:2211.17057 [astro-ph.HE] .
  * Fornieri and Zhang [2022] O. Fornieri and H. Zhang, Phys. Rev. D 106, 103015 (2022).
  * Kulsrud and Pearce [1969] R. Kulsrud and W. P. Pearce, ApJ 156, 445 (1969).
  * Voelk [1975] H. J. Voelk, Reviews of Geophysics and Space Physics 13, 547 (1975).
  * Fornieri  _et al._ [2021] O. Fornieri, D. Gaggero, S. S. Cerri, P. De la Torre Luque, and S. Gabici, Mon. Not. Roy. Astron. Soc. 502, 5821 (2021), arXiv:2011.09197 [astro-ph.HE] .
  * Ginzburg  _et al._ [1970] V. Ginzburg, V. Ginzburg, J. Sykes, and R. Tayler,  _The Propagation of Electromagnetic Waves in Plasmas_, Commonwealth and International Library (Pergamon Press, 1970).
  * Braginskii [1965] S. I. Braginskii, Reviews of Plasma Physics 1, 205 (1965).
  * Petrosian  _et al._ [2006] V. Petrosian, H. Yan, and A. Lazarian, ApJ 644, 603 (2006), arXiv:astro-ph/0508567 [astro-ph] .
  * Yan and Lazarian [2008] H. Yan and A. Lazarian, Astrophys. J. 673, 942 (2008), arXiv:0710.2617 [astro-ph] .
  * Spitzer [1962] L. Spitzer,  _Physics of Fully Ionized Gases_ (1962).
  * Ackermann  _et al._ [2011] M. Ackermann _et al._ , Science 334, 1103 (2011).
  * Abeysekara  _et al._ [2021] A. U. Abeysekara _et al._ , Nature Astron. 5, 465 (2021), arXiv:2103.06820 [astro-ph.HE] .
  * 1  _et al._ [2021] G. M. 1, P. Blasi, E. Peretti, and P. Cristofari,  (2021), 10.1093/mnras/stab690, arXiv:2102.09217 [astro-ph.HE] .
  * Aharonian  _et al._ [2019] F. Aharonian, R. Yang, and E. de Oña Wilhelmi, Nature Astron. 3, 561 (2019), arXiv:1804.02331 [astro-ph.HE] .
  * De La Torre Luque  _et al._ [2022] P. De La Torre Luque, O. Fornieri, and T. Linden, (2022), arXiv:2205.08544 [astro-ph.HE] .
  * Schekochihin [2020] A. A. Schekochihin, (2020), arXiv:2010.00699 [physics.plasm-ph] .


  * March 2021

# Quantum interpolating ensemble: Average entropies and orthogonal polynomials

Lu Wei1 and Nicholas Witte2 1 Department of Electrical and Computer
Engineering, University of Michigan, Dearborn, MI 48128, USA  2 School of
Mathematics and Statistics, Victoria University of Wellington, Wellington
6012, New Zealand<EMAIL_ADDRESS>
###### Abstract

The density matrix formalism is a fundamental tool in studying various
problems in quantum information processing. In the space of density matrices,
the most well-known and physically relevant measures are the Hilbert-Schmidt
ensemble and the Bures-Hall ensemble. In this work, we propose a generalized
ensemble of density matrices, termed quantum interpolating ensemble, which is
able to interpolate between these two seemingly unrelated ensembles. As a
first step to understand the proposed ensemble, we derive the exact mean
formulas of entanglement entropies over such an ensemble generalizing several
recent results in the literature. We also derive some key properties of the
corresponding orthogonal polynomials relevant to obtaining other statistical
information of the entropies. Numerical results demonstrate the usefulness of
the proposed ensemble in estimating the degree of entanglement of quantum
states.

Keywords: quantum entanglement, interpolating ensemble, entanglement entropy,
orthogonal polynomials, special functions

##  1 Introduction and quantum interpolating ensemble

Quantum information theory aims at understanding the theoretical underpinnings
of quantum technologies including quantum computing and quantum
communications. It is based on probabilistic interpretations of quantum states
to explain various quantum effects. The density matrix formalism introduced by
von Neumann [1] provides a natural framework to describe density matrices of
quantum states. The density matrix is a fundamental object that encodes all
the information of a quantum state. Among the different measures of density
matrices, the most well-known and physically relevant ones [2] are the
Hilbert-Schmidt measure and the Bures-Hall measure.

The Hilbert-Schmidt measure is formulated as follows. Consider a bipartite
quantum system consisting of two subsystems $A$ and $B$ in the Hilbert space
$\mathcal{H}_{m}$ and $\mathcal{H}_{n}$ (with $m\leq n$), respectively. A
random pure state $\Ket{\psi}$, defined as a linear combination of the
complete basis of the subsystems, belongs to the composite Hilbert space
$\Ket{\psi}\in\mathcal{H}_{m}\otimes\mathcal{H}_{n}$. The reduced density
matrix is obtained by partial tracing over the larger system of the full
density matrix $\rho=\Ket{\psi}\Bra{\psi}$ as $\rho_{A}=\tr_{B}\rho$. The
resulting density of eigenvalues of $\rho_{A}$ is the Hilbert-Schmidt measure
[2]

$f_{\rm{HS}}\left(\bm{\lambda}\right)\propto\delta\left(1-\sum_{i=1}^{m}\lambda_{i}\right)\prod_{1\leq
i<j\leq
m}\left(\lambda_{i}-\lambda_{j}\right)^{2}\prod_{i=1}^{m}\lambda_{i}^{n-m},$
(1)

where $\delta(\cdot)$ is the Dirac delta function. The joint density (1) is
also referred to as the Hilbert-Schmidt (random matrix) ensemble, which is the
eigenvalue density of a normalized Wishart matrix

$\frac{\mathbf{GG}^{{\dagger}}}{\tr\left(\mathbf{GG}^{{\dagger}}\right)}$ (2)

with $\mathbf{G}$ being an $m\times n$ complex Gaussian matrix. For the Bures-
Hall measure, the random pure state is given by a superposition of that of the
Hilbert-Schmidt measure as $\Ket{\varphi}=\Ket{\psi}+\left(\mathbf{U}\otimes
I_{m}\right)\Ket{\psi}$, where $\mathbf{U}$ is an $m\times m$ unitary matrix
with the measure proportional to $\det\left(I_{m}+\mathbf{U}\right)^{2(n-m)}$.
The resulting density of eigenvalues of the reduced density matrix
$\rho_{A}=\tr_{B}\Ket{\varphi}\Bra{\varphi}$ is the (generalized) Bures-Hall
measure [3, 4]

$f_{\rm{BH}}\left(\bm{\lambda}\right)\propto\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}^{n-m-\frac{1}{2}}.$
(3)

The joint density (3) is also known as the Bures-Hall ensemble in random
matrix theory, which is the eigenvalue density of the normalized product of
the matrix $I_{m}+\mathbf{U}$ with a complex Gaussian matrix $\mathbf{G}$,

$\frac{(I_{m}+\mathbf{U})\mathbf{GG}^{{\dagger}}(I_{m}+\mathbf{U}^{{\dagger}})}{\tr\left((I_{m}+\mathbf{U})\mathbf{GG}^{{\dagger}}(I_{m}+\mathbf{U}^{{\dagger}})\right)}.$
(4)

The Hilbert-Schmidt measure (1) and the Bures-Hall measure (3) are supported
in the probability simplex

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

which reflects the constraint $\tr\rho_{A}=1$ of density matrices. Note also
that the normalization constants in the densities (1) and (3) are omitted.

The study of the Hilbert-Schmidt measure has received substantial attention,
see, for example, the results in [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23]. These results include information-theoretic
studies of different entanglement entropies [5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19, 20] as well as applications to quantum information
processing [21, 22, 23]. The relatively less-studied Bures-Hall ensemble [3,
4, 24, 25, 26, 27, 28, 29] gains renewed interest very recently [30, 31, 32].
This is owing to the recent breakthrough in probability theory in
understanding various aspects of the Bures-Hall ensemble [33, 34, 35, 36, 37].
Despite the distinct behavior of Hilbert-Schmidt measure and Bures-Hall
measure, an interesting question is whether one could propose a measure that
interpolates between the two measures. This question has also been motivated
by the observation in [30] that the Bures-Hall measure tends to be more
conservative than the Hilbert-Schmidt measure in estimating entanglement
entropies. Namely, the Bures-Hall measure leads towards an estimate of less
entangled states than the Hilbert-Schmidt measure does. In this context, one
tries to control the appropriate amount of entanglement as a resource for
quantum information processing by constructing new measures that interpolates
between the two major measures.

In this work, we consider the following measure, supported in (5),

$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{\lambda_{i}-\lambda_{j}}{\lambda_{i}+\lambda_{j}}\left(\lambda_{i}^{\theta}-\lambda_{j}^{\theta}\right)\prod_{i=1}^{m}\lambda_{i}^{a}$
(6)

termed the quantum interpolating ensemble, where $\theta$ is assumed to be a
positive real parameter and $a>-1$. Clearly, the proposed ensemble (6) reduces
to the Hilbert-Schmidt ensemble (1) and the Bures-Hall ensemble (3) as special
cases,

$f\left(\bm{\lambda}\right)=\cases{f_{\rm{BH}}\left(\bm{\lambda}\right)&for
$\theta=1,~{}a=n-m-\frac{1}{2}$\\\ f_{\rm{HS}}\left(\bm{\lambda}\right)&for
$\theta=2,~{}a=n-m$.\\\ }$ (7)

Namely, as $\theta$ varies from $\theta=1$ to $\theta=2$, the quantum
interpolating ensemble interpolates between the Bures-Hall ensemble and the
Hilbert-Schmidt ensemble. Due to Schur’s Pfaffian identity [38, 39]

$\prod_{1\leq i<j\leq
2m}\frac{\lambda_{i}-\lambda_{j}}{\lambda_{i}+\lambda_{j}}={\rm{Pf}}\left(\frac{\lambda_{i}-\lambda_{j}}{\lambda_{i}+\lambda_{j}}\right)_{1\leq
i,j\leq 2m},$ (8)

the proposed ensemble (6) is described by a Pfaffian point process for any
$\theta>0$ except for the special value $\theta=2$ when the ensemble becomes a
determinantal point process. Therefore, in the interested interval
$\theta\in[1,2]$ the new ensemble (6) corresponds to the transition from a
Pfaffian point process to a determinantal point process. It is also worth
mentioning that besides the half integers values of $a$ for the Bures-Hall
ensemble and the integers values of $a$ for the Hilbert-Schmidt ensemble, the
proposed ensemble is valid for any $a>-1$. Therefore, in addition to $\theta$,
the parameter $a$ is also considered as a deformation parameter that defines
the interpolating ensemble.

With an interpolating ensemble being identified, a natural question is what
will be the statistical behavior of entanglement entropies over such an
ensemble? Will the values of entropies also interpolate between those of the
Hilbert-Schmidt ensemble and the Bures-Hall ensemble? In addition, what are
the properties of the underlying orthogonal polynomial systems? We address
these questions in this work. The rest of the paper is organized as follows.
In section 2, we present the mains results of the paper on the mean values of
quantum purity and von Neumann entanglement entropy of the interpolating
ensemble (section 2.1) and some basic properties of the corresponding
orthogonal polynomials relevant to higher moment calculation (section 2.2). In
section 3, we perform extensive numerical studies on the impact of the
deformation parameters of the ensemble on the behavior of the entanglement.
Some proofs of the results can be found in the two appendices.

##  2 Entropies and orthogonal polynomials

###  2.1 Average entropies over interpolating ensemble

The degree of entanglement of the subsystems $A$ and $B$ is estimated by
entanglement entropies, which are functions of the eigenvalues (entanglement
spectrum) of a given ensemble. Any function that satisfies a list of axioms
can be considered as an entanglement entropy. In particular, an entropy should
monotonically change from the separable state

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

to the maximally-entangled state

$\lambda_{1}=\lambda_{2}=\dots\lambda_{m}=1/m.$ (10)

A standard one we consider here is quantum purity [2]

$S_{\rm{P}}=\sum_{i=1}^{m}\lambda_{i}^{2},$ (11)

supported in $S_{\rm{P}}\in[1/m,1]$, which attains the separable state and
maximally-entangled state when $S_{\rm{P}}=1$ and when $S_{\rm{P}}=1/m$,
respectively. Quantum purity (11) is an example of polynomial entropies,
whereas a well-known non-polynomial entropy that we also consider here is von
Neumann entropy [2]

$S_{\rm{vN}}=-\sum_{i=1}^{m}\lambda_{i}\ln\lambda_{i}.$ (12)

The von Neumann entropy (12) is supported in $S_{\rm{vN}}\in[0,\ln{m}]$ that
achieves the separable state and maximally-entangled state when
$S_{\rm{vN}}=0$ and when $S_{\rm{vN}}=\ln{m}$, respectively.

Statistical information of entanglement entropies is encoded through their
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 distributions. In this work, we focus on
the first moments of purity (11) and von Neumann entropy (12) over the quantum
interpolating ensemble (6). Moment computation over an ensemble with the
probability constraint $\delta\left(1-\sum_{i=1}^{m}\lambda_{i}\right)$ is
typically performed over an ensemble without the constraint [7, 9, 11, 12, 30,
31, 32]. As will be seen, the unconstrained ensemble that corresponds to (6)
is given by

$h\left(\bm{x}\right)=\frac{1}{C^{\prime}}\prod_{1\leq i<j\leq
m}\frac{x_{i}-x_{j}}{x_{i}+x_{j}}\left(x_{i}^{\theta}-x_{j}^{\theta}\right)\prod_{i=1}^{m}x_{i}^{a}\rme^{-x_{i}},$
(13)

where $x_{i}\in[0,\infty)$, $i=1,\dots,m$. This ensemble has been recently
studied in [40] in connection to a $\theta$-deformed Cauchy-Laguerre two-
matrix model. In the case when $\theta=1$, the corresponding ensembles have
been studied in [33, 34, 35, 36, 37]. Among other results in [40], the
$k$-point ($1\leq k\leq m$) correlation function of the unconstrained ensemble
is shown to follow a Pfaffian point process

$\rho_{k}(x_{1},\dots,x_{k})\propto{\rm{Pf}}\left(\begin{array}[]{cc}\Delta
K_{11}(x_{i},x_{j})&\Sigma K_{01}(x_{i},x_{j})\\\ -\Sigma
K_{01}(x_{j},x_{i})&\Delta K_{00}(x_{i},x_{j})\\\ \end{array}\right)_{1\leq
i,j\leq k},$ (14)

where we denote

$\displaystyle\Delta K_{00}(x,y)=K_{00}(x,y)-K_{00}(y,x)$ (15a)
$\displaystyle\Sigma K_{01}(x,y)=K_{01}(x,y)+K_{10}(y,x)$ (15b)
$\displaystyle\Delta K_{11}(x,y)=K_{11}(x,y)-K_{11}(y,x).$ (15c)

The above kernel functions can be expressed via the following Fox H-functions
[40]

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

as

$\displaystyle
K_{00}(x,y)=\theta\int_{0}^{1}t^{\alpha}H_{a}(x)H_{a+1}(y)\,\mathrm{d}t$
(15pqa) $\displaystyle K_{01}(x,y)=\theta
x^{2a+1}\int_{0}^{1}t^{\alpha}H_{a}(y)G_{a+1}(x)\,\mathrm{d}t$ (15pqb)
$\displaystyle K_{10}(x,y)=\theta
y^{2a+1}\int_{0}^{1}t^{\alpha}H_{a+1}(x)G_{a}(y)\,\mathrm{d}t$ (15pqc)
$\displaystyle
K_{11}(x,y)=\theta(xy)^{2a+1}\int_{0}^{1}t^{\alpha}G_{a+1}(x)G_{a}(y)\,\mathrm{d}t-\frac{x^{a}y^{a+1}}{x+y},$
(15pqd)

where we denote

$\alpha=\frac{2(a+1)}{\theta}-1.$ (15pqr)

In general, the Fox H-function is defined through the following contour
integral [41]

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

where the contour $\mathcal{L}$ separates the poles of
$\Gamma\left(b_{j}+B_{j}s\right)$ from the poles of
$\Gamma\left(1-a_{j}-A_{j}s\right)$. In the special case $A_{1}=\dots=A_{p}=1$
and $B_{1}=\dots=B_{q}=1$, the Fox H-function reduces to the Meijer G-function
[41]. The integral forms of the kernel functions (15pqa)–(15pqd) are useful
for the computation of mean entropies, whereas several other distinct
representations of the kernel functions exist [35, 36, 37, 40, 32]. In
particular, we present the following biorthogonal polynomial forms [40] useful
for a later discussion

$\displaystyle
K_{00}(x,y)=\sum_{k=0}^{m-1}p_{k}\left(x^{\theta}\right)q_{k}\left(y^{\theta}\right)$
(15pqwa) $\displaystyle
K_{01}(x,y)=x^{a}\rme^{-x}\int_{0}^{\infty}\frac{v^{a+1}\rme^{-v}}{x+v}K_{00}(y,v)\,\mathrm{d}v$
(15pqwb) $\displaystyle
K_{10}(x,y)=y^{a+1}\rme^{-y}\int_{0}^{\infty}\frac{w^{a}\rme^{-w}}{y+w}K_{00}(w,x)\,\mathrm{d}w$
(15pqwc) $\displaystyle
K_{11}(x,y)=x^{a}y^{a+1}\rme^{-x-y}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{v^{a}\rme^{-v}}{y+v}\frac{w^{a+1}\rme^{-w}}{x+w}K_{00}(v,w)\,\mathrm{d}v\,\mathrm{d}w-W(x,y),$
(15pqwd)

where the biorthogonal polynomials

$\displaystyle p_{j}\left(x^{\theta}\right)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{j}\frac{\sqrt{2}(-1)^{k+j}\Gamma(k+j+\alpha+1)x^{\theta
k}}{\Gamma(\theta k+a+1)\Gamma(k+\alpha+1)(j-k)!k!}$ (15pqwxa) $\displaystyle
q_{j}\left(y^{\theta}\right)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{j}\frac{\sqrt{2}(-1)^{k+j}(\theta
j+a+1)\Gamma(k+j+\alpha+1)y^{\theta k}}{\Gamma(\theta
k+a+2)\Gamma(k+\alpha+1)(j-k)!k!}$ (15pqwxb)

are orthogonal with respect to the weight function

$W(x,y)=\frac{x^{a}y^{a+1}\rme^{-x-y}}{x+y}$ (15pqwxy)

as

$\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}p_{k}\left(x^{\theta}\right)q_{l}\left(y^{\theta}\right)W(x,y)\,\mathrm{d}x\,\mathrm{d}y=\delta_{kl}.$
(15pqwxz)

We now study the moment relations of entanglement entropies between the
proposed ensemble (6) and its unconstrained version (13). Firstly, the density
$g_{d}(r)$ of the trace

$r=\sum_{i=1}^{m}x_{i},~{}~{}~{}~{}r\in[0,\infty)$ (15pqwxaa)

of the unconstrained ensemble (13) is obtained as

$\displaystyle g_{d}(r)$ $\displaystyle=$
$\displaystyle\int_{\bm{x}}h(\bm{x})\delta\left(r-\sum_{i=1}^{m}x_{i}\right)\prod_{i=1}^{m}\,\mathrm{d}x_{i}$
(15pqwxab) $\displaystyle=$
$\displaystyle\frac{C}{C^{\prime}}\rme^{-r}r^{d-1}\int_{\bm{\lambda}}f(\bm{\lambda})\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}$
(15pqwxac) $\displaystyle=$
$\displaystyle\frac{1}{\Gamma\left(d\right)}\rme^{-r}r^{d-1},$ (15pqwxad)

where we have used the change of variables

$x_{i}=r\lambda_{i},~{}~{}~{}~{}i=1,\ldots,m$ (15pqwxae)

and the resulting Jacobian calculation leads to the normalization $\Gamma(d)$
with

$d=\frac{m}{2}(m\theta-\theta+2a+2).$ (15pqwxaf)

The above calculation implies that the density $h(\bm{x})$ can be factored as

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

i.e., the random variable $r$ is independent of each $\lambda_{i}$ (hence
independent of $S_{\rm{P}}$ and $S_{\rm{vN}}$). Similar factorizations also
exist for the Hilbert-Schmidt ensemble [7] and the Bures-Hall ensemble [30,
32]. Introducing the corresponding quantum purity of the unconstrained
ensemble

$T_{\rm{P}}=\sum_{i=1}^{m}x_{i}^{2},$ (15pqwxah)

the $k$-th moment of quantum purity $S_{\rm{P}}$ is represented as

$\displaystyle\mathbb{E}_{f}\\!\left[S_{\rm{P}}^{k}\right]$ $\displaystyle=$
$\displaystyle\int_{\bm{\lambda}}S_{\rm{P}}^{k}~{}f(\bm{\lambda})\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}$
(15pqwxai) $\displaystyle=$
$\displaystyle\int_{\bm{\lambda}}\frac{T_{\rm{P}}^{k}}{r^{2k}}f(\bm{\lambda})\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}\int_{r}g_{d+2k}(r)\,\mathrm{d}r$
(15pqwxaj) $\displaystyle=$
$\displaystyle\frac{\Gamma(d)}{\Gamma(d+2k)}\int_{\bm{\lambda}}\int_{r}T_{\rm{P}}^{k}~{}f(\bm{\lambda})g_{d}(r)\,\mathrm{d}r\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}$
(15pqwxak) $\displaystyle=$
$\displaystyle\frac{\Gamma(d)}{\Gamma(d+2k)}\mathbb{E}_{h}\\!\left[T_{\rm{P}}^{k}\right],$
(15pqwxal)

where we have used the change of variables (15pqwxae) and the independence
property (15pqwxag). Therefore, computing the $k$-th moment of $S_{\rm{P}}$
can be converted to computing the $k$-th moment of $T_{\rm{P}}$. In
particular, the first moments are related by

$\mathbb{E}_{f}\\!\left[S_{\rm{P}}\right]=\frac{1}{d(d+1)}\mathbb{E}_{h}\\!\left[T_{\rm{P}}\right].$
(15pqwxam)

We now introduce von Neumann entropy of the unconstrained ensemble

$T_{\rm{vN}}=\sum_{i=1}^{m}x_{i}\ln x_{i}$ (15pqwxan)

that leads to the identity

$S_{\rm{vN}}=\ln r-r^{-1}T_{\rm{vN}},$ (15pqwxao)

then the first moment relation is similarly obtained as

$\displaystyle\mathbb{E}_{f}\\!\left[S_{\rm{vN}}\right]$ $\displaystyle=$
$\displaystyle\int_{\bm{\lambda}}S_{\rm{vN}}~{}f(\bm{\lambda})\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}\int_{r}g_{d+1}(r)\,\mathrm{d}r$
$\displaystyle=$
$\displaystyle\frac{\Gamma(d)}{\Gamma(d+1)}\left(\int_{r}g_{d}(r)r\ln
r\,\mathrm{d}r-\int_{\bm{\lambda}}\int_{r}T_{\rm{vN}}~{}f(\bm{\lambda})g_{d}(r)\,\mathrm{d}r\prod_{i=1}^{m}\,\mathrm{d}\lambda_{i}\right)$
$\displaystyle=$
$\displaystyle\psi_{0}(d+1)-\frac{1}{d}\mathbb{E}_{h}\\!\left[T_{\rm{vN}}\right],$
(15pqwxaq)

where we have also used

$\int_{0}^{\infty}\rme^{-r}r^{a-1}\ln{r}\,\mathrm{d}r=\Gamma(a)\psi_{0}(a),~{}~{}~{}~{}\Re(a)>0$
(15pqwxar)

with $\psi_{0}(x)=\,\mathrm{d}\ln\Gamma(x)/\,\mathrm{d}x$ denoting the digamma
function [42]. For a positive integer $l$, the digamma function admits the
following useful identities

$\displaystyle\psi_{0}(l)=-\gamma+\sum_{k=1}^{l-1}\frac{1}{k}$ (15pqwxasa)
$\displaystyle\psi_{0}\left(l+\frac{1}{2}\right)=-\gamma-2\ln
2+2\sum_{k=0}^{l-1}\frac{1}{2k+1},$ (15pqwxasb)

where $\gamma\approx 0.5772$ is the Euler’s constant.

With the above preparation, we now state the results of average entropies over
the interpolating ensemble.

######  Proposition 1.

The average value of quantum purity (11) under the quantum interpolating
ensemble (6), for any $\theta>0$ and $a>-1$, is given by

$\mathbb{E}_{f}\\!\left[S_{\rm{P}}\right]=\sum_{k=0}^{m-1}\frac{(-1)^{k+m-1}\theta(\theta
k+a+2)^{2}}{2d(d+1)(m-1-k)!k!}A_{\theta}(k),$ (15pqwxasat)

where

$A_{\theta}(k)=\frac{\Gamma\left(k+\frac{2(a+2)}{\theta}\right)\Gamma\left(k+m+\frac{2(a+1)}{\theta}\right)\Gamma\left(k+1+\frac{2}{\theta}\right)}{\Gamma\left(k+\frac{2(a+1)}{\theta}\right)\Gamma\left(k+m+\frac{2(a+2)}{\theta}\right)\Gamma\left(k+1+\frac{2}{\theta}-m\right)}.$
(15pqwxasau)

Before proving proposition 1, two remarks on its special cases are in order.

######  Remark 1.

In the special case $\theta=1$, $a=n-m-\frac{1}{2}$ that corresponds to the
Bures-Hall ensemble (3), by recognizing $A_{1}(m-2)$ and $A_{1}(m-1)$ as the
only non-vanishing terms in (15pqwxasau) with

$A_{1}(k)=\frac{(k+2n-2m+2)(k+2n-2m+1)\Gamma(k+3)}{(k+2n-m+2)(k+2n-m+1)\Gamma(k+3-m)},$
(15pqwxasav)

the expression (15pqwxasat) simplifies to

$\displaystyle\mathbb{E}_{f}\\!\left[S_{\rm{P}}\right]$ $\displaystyle=$
$\displaystyle\sum_{k=m-2}^{m-1}\frac{2(-1)^{k+m-1}(k+n-m+3/2)^{2}}{m(2n-m)(2mn-m^{2}+2)(m-1-k)!k!}A_{1}(k)$
(15pqwxasaw) $\displaystyle=$
$\displaystyle\frac{2n(2n+m)-m^{2}+1}{2n\left(2mn-m^{2}+2\right)}.$
(15pqwxasax)

This recovers the mean purity formula of the Bures-Hall ensemble recently
reported in [30, 31].

######  Remark 2.

In the special case $\theta=2$, $a=n-m$ that corresponds to the Hilbert-
Schmidt ensemble (1), by recognizing $A_{2}(m-1)$ as the only non-vanishing
term in (15pqwxasau) with

$A_{2}(k)=\frac{(k+n-m+1)\Gamma(k+2)}{(k+n+1)\Gamma(k+2-m)},$ (15pqwxasay)

the expression (15pqwxasat) simplifies to

$\displaystyle\mathbb{E}_{f}\\!\left[S_{\rm{P}}\right]$ $\displaystyle=$
$\displaystyle\sum_{k=m-1}^{m-1}\frac{(-1)^{k+m-1}(2k+n-m+2)^{2}}{mn(mn+1)(m-1-k)!k!}A_{2}(k)$
(15pqwxasaz) $\displaystyle=$ $\displaystyle\frac{m+n}{mn+1}.$ (15pqwxasba)

We recover the mean purity formula of the Hilbert-Schmidt ensemble obtained in
[5].

We now prove the proposition 1.

######  Proof.

The essential task is to compute $\mathbb{E}_{h}\\!\left[T_{\rm{P}}\right]$,
which, after inserting into the moment relation (15pqwxam), will establish the
proposition 1. The required single eigenvalue density $h_{1}(x)$ of the
unconstrained ensemble (13) can be read off from the correlation function (14)
that corresponds to a Pfaffian of a $2\times 2$ matrix as

$h_{1}(x)=\frac{1}{m}\rho_{1}(x)=\frac{1}{2m}\left(K_{01}(x,x)+K_{10}(x,x)\right).$
(15pqwxasbb)

The computation now boils down to computing two integrals

$\displaystyle\mathbb{E}_{h}\\!\left[T_{\rm{P}}\right]$ $\displaystyle=$
$\displaystyle m\int_{0}^{\infty}x^{2}h_{1}(x)\,\mathrm{d}x$ (15pqwxasbc)
$\displaystyle=$
$\displaystyle\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.$
(15pqwxasbd)

The starting point to calculate the above integrals is the fact that the
contour form (15pqu) of the Fox H-function (15pc) admits a finite number of
single poles, which by residue calculation gives a finite sum

$H_{q}(x)=\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\left(tx^{\theta}\right)^{k}}{\Gamma(k+\alpha+1)\Gamma(\theta
k+q+1)(m-1-k)!k!}.$ (15pqwxasbe)

Therefore, we have

$\displaystyle\int_{0}^{\infty}x^{2}K_{01}(x,x)\,\mathrm{d}x$ (15pqwxasbf)
$\displaystyle=$ $\displaystyle\int_{0}^{\infty}\theta
x^{2a+3}\int_{0}^{1}t^{\alpha}H_{a}(y)G_{a+1}(x)\,\mathrm{d}t\,\mathrm{d}x$
(15pqwxasbg) $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{\theta(-1)^{k}\Gamma(k+\alpha+m+1)}{\Gamma(k+\alpha+1)\Gamma(\theta
k+a+1)(m-1-k)!k!}\int_{0}^{\infty}\\!\\!x^{\theta
k+2a+3}\int_{0}^{1}\\!t^{\alpha+k}G_{a+1}(x)\,\mathrm{d}t\,\mathrm{d}x$
(15pqwxasbh) $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{(-1)^{k+m-1}\theta(\theta k+a+1)(\theta
k+a+2)}{2(m-1-k)!k!}A_{\theta}(k),$ (15pqwxasbi)

where the integrals over $t$ and $x$ in (15pqwxasbh) are evaluated
respectively by the identity [41]

$\displaystyle\int_{0}^{1}x^{\rho-1}H_{p,q}^{m,n}\left(\begin{array}[]{c}(a_{1},A_{1}),\ldots,(a_{n},A_{n});(a_{n+1},A_{n+1}),\ldots,(a_{p},A_{p})\\\
(b_{1},B_{1}),\ldots,(b_{m},B_{m});(b_{m+1},B_{m+1}),\ldots,(b_{q},B_{q})\end{array}\Big{|}\eta
x\Big{.}\right)\,\mathrm{d}x$ (15pqwxasbl) $\displaystyle=$ $\displaystyle
H_{p+1,q+1}^{m,n+1}\left(\begin{array}[]{c}(1-\rho,1),(a_{1},A_{1}),\ldots,(a_{n},A_{n});(a_{n+1},A_{n+1}),\ldots,(a_{p},A_{p})\\\
(b_{1},B_{1}),\ldots,(b_{m},B_{m});(b_{m+1},B_{m+1}),\ldots,(b_{q},B_{q}),(-\rho,1)\end{array}\Big{|}\eta\Big{.}\right)$
(15pqwxasbo)

and the Mellin transform of the Fox H-function [41], cf. (15pqu),

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

In the same manner, the second integral in (15pqwxasbd) is evaluated to

$\int_{0}^{\infty}x^{2}K_{10}(x,x)\,\mathrm{d}x=\sum_{k=0}^{m-1}\frac{(-1)^{k+m-1}\theta(\theta
k+a+2)(\theta k+a+3)}{2(m-1-k)!k!}A_{\theta}(k).$ (15pqwxasbt)

Putting together the results (15pqwxam), (15pqwxasbd), (15pqwxasbi), and
(15pqwxasbt), we complete the proof of proposition 1. $\square$ $\blacksquare$

The next result on the average von Neumann entropy is summarized in the
following proposition.

######  Proposition 2.

The average value of von Neumann entropy (12) under the quantum interpolating
ensemble (6), for any $\theta>0$ and $a>-1$, is given by

$\mathbb{E}_{f}\\!\left[S_{\rm{vN}}\right]=\psi_{0}(d+1)-\sum_{k=0}^{m-1}\frac{(-1)^{k+m-1}\left(\theta
k+a+3/2\right)}{d(m-1-k)!k!}B_{\theta}(k),$ (15pqwxasbu)

where

$\displaystyle B_{\theta}(k)$ $\displaystyle=$
$\displaystyle\frac{\Gamma\left(k+\frac{2a+3}{\theta}\right)\Gamma\left(k+m+\frac{2a+2}{\theta}\right)\Gamma\left(k+1+\frac{1}{\theta}\right)}{\Gamma\left(k+\frac{2a+2}{\theta}\right)\Gamma\left(k+m+\frac{2a+3}{\theta}\right)\Gamma\left(k+1+\frac{1}{\theta}-m\right)}\Bigg{(}\psi_{0}\left(k+\frac{2a+3}{\theta}\right)$
(15pqwxasbv)
$\displaystyle+\psi_{0}\left(k+1+\frac{1}{\theta}\right)-\psi_{0}\left(k+m+\frac{2a+3}{\theta}\right)-\psi_{0}\left(k+1+\frac{1}{\theta}-m\right)$
$\displaystyle+\theta\left(\psi_{0}(\theta k+a+2)-\frac{\theta k+a+1}{\theta
k+a+3/2}\right)\\!\Bigg{)}.$

######  Proof.

The main task is to compute the average
$\mathbb{E}_{h}\\!\left[T_{\rm{vN}}\right]$, which, after inserting into the
moment relation (15pqwxaq), establishes the proposition 2. By employing the
one arbitrary eigenvalue density (15pqwxasbb), the task boils down to
computing two integrals

$\mathbb{E}_{h}\\!\left[T_{\rm{vN}}\right]=\frac{1}{2}\int_{0}^{\infty}x\ln
xK_{01}(x,x)\,\mathrm{d}x+\frac{1}{2}\int_{0}^{\infty}x\ln
xK_{10}(x,x)\,\mathrm{d}x.$ (15pqwxasbw)

We first compute the integral

$\int_{0}^{\infty}x^{\beta}K_{01}(x,x)\,\mathrm{d}x,~{}~{}~{}~{}\beta>0$
(15pqwxasbx)

by using the results (15pqwxasbe), (15pqwxasbl), and (15pqwxasbr) as

$\displaystyle\int_{0}^{\infty}x^{\beta}K_{01}(x,x)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{(-1)^{k+m}\Gamma\left(k+\alpha+m+1\right)}{2\theta(m-1-k)!k!\Gamma\left(k+\alpha+1\right)\Gamma(\theta
k+a+1)}$ (15pqwxasby)
$\displaystyle\times\frac{\Gamma(s)\Gamma(s-\alpha)\Gamma(k+\alpha-s+1)\Gamma(\theta
s-a-1)}{\Gamma(s+m)\Gamma(s-\alpha-m)\Gamma(k+\alpha-s+2)},$

where $\alpha$ is given by (15pqr) and we denote

$s=\beta-1+k+\frac{2a+3}{\theta}.$ (15pqwxasbz)

Similarly, we also obtain

$\displaystyle\int_{0}^{\infty}x^{\beta}K_{10}(x,x)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{(-1)^{k+m}\Gamma\left(k+\alpha+m+1\right)}{2\theta(m-1-k)!k!\Gamma\left(k+\alpha+1\right)\Gamma(\theta
k+a+2)}$ (15pqwxasca)
$\displaystyle\times\frac{\Gamma(s)\Gamma(s-\alpha)\Gamma(k+\alpha-s+1)\Gamma(\theta
s-a)}{\Gamma(s+m)\Gamma(s-\alpha-m)\Gamma(k+\alpha-s+2)}.$

Taking the derivative of (15pqwxasby) and (15pqwxasca) with respect to $\beta$
before setting $\beta\to 1$ leads to the desired expression of (15pqwxasbw),
which upon inserting into the moment relation (15pqwxaq) completes the proof
of proposition 2. $\square$ $\blacksquare$

######  Remark 3.

In the special case $\theta=1$, $a=n-m-\frac{1}{2}$ that corresponds to the
Bures-Hall ensemble (3), the result (15pqwxasbu) in proposition 2 reduces to

$\displaystyle\mathbb{E}_{f}\\!\left[S_{\rm{vN}}\right]$ $\displaystyle=$
$\displaystyle\psi_{0}\left(mn-\frac{m^{2}}{2}+1\right)-\psi_{0}\left(n+\frac{1}{2}\right)+\psi_{0}(2n+1)-\psi_{0}(2n-m+1)+\psi_{0}(1)$
(15pqwxascb)
$\displaystyle-\psi_{0}(m+1)+\frac{2n-1}{2n}+\sum_{k=0}^{m-2}\frac{2(k+1)(k+n-m+1)(k+2n-2m+1)}{m(2n-m)(m-1-k)(k+2n-m+1)}$
$\displaystyle=$
$\displaystyle\psi_{0}\left(mn-\frac{m^{2}}{2}+1\right)-\psi_{0}\left(n+\frac{1}{2}\right).$

This recovers the mean formula of von Neumann entropy under the Bures-Hall
ensemble recently obtained in [30, 31].

######  Remark 4.

In the special case $\theta=2$, $a=n-m$ that corresponds to the Hilbert-
Schmidt ensemble (1), the mean formula of von Neumann entropy is well-known

$\mathbb{E}_{f}\\!\left[S_{\rm{vN}}\right]=\psi_{0}(mn+1)-\psi_{0}(n)-\frac{m+1}{2n},$
(15pqwxascc)

which was conjectured by Page [7] and later proved in [8, 9]. For this special
case, to show that the proposition 2 can be reduced to the above result
requires showing the identity

$\displaystyle\frac{1}{mn}\sum_{k=0}^{m-1}\frac{(-1)^{k+m-1}\left(2k+n-m+\frac{3}{2}\right)\Gamma\left(k+\frac{3}{2}\right)\Gamma(k+n+1)\Gamma\left(k+n-m+\frac{3}{2}\right)}{(m-1-k)!k!\Gamma\left(k-m+\frac{3}{2}\right)\Gamma\left(k+n+\frac{3}{2}\right)\Gamma(k+n-m+1)}$
$\displaystyle\times\Bigg{(}\psi_{0}\left(k+\frac{3}{2}\right)+\psi_{0}\left(k+n-m+\frac{3}{2}\right)-\psi_{0}\left(k+n+\frac{3}{2}\right)-\psi_{0}\left(k-m+\frac{3}{2}\right)$
$\displaystyle+2\psi_{0}(2k+n-m+2)-\frac{2(2k+n-m+1)}{2k+n-m+3/2}\Bigg{)}=\psi_{0}(n)+\frac{m+1}{2n},$
(15pqwxascd)

which, despite strong numerical evidence, remains unproven.

###  2.2 Some properties of orthogonal polynomials

Before illustrating the behavior of average entropies via numerical
simulations, we present some results on the $\theta$-deformed orthogonal
polynomials (15pqwxa), (15pqwxb). Our first result is the recurrence relation
of $p_{j}\left(x^{\theta}\right)$ for $\theta=2$, which will be useful in
computing the higher order moments of entanglement entropies.

######  Proposition 3.

The $\theta$-deformed orthogonal polynomials $p_{j}\left(x^{\theta}\right)$ in
(15pqwxa) for $\theta=2$ satisfy a five-term recurrence relation as

$\displaystyle
x^{2}\left(a_{2}p_{j+2}\left(x^{2}\right)+a_{1}p_{j+1}\left(x^{2}\right)+a_{0}p_{j}\left(x^{2}\right)\right)=$
$\displaystyle
r_{3}p_{j+3}\left(x^{2}\right)+r_{2}p_{j+2}\left(x^{2}\right)+r_{1}p_{j+1}\left(x^{2}\right)+r_{0}p_{j}\left(x^{2}\right)+r_{-1}p_{j-1}\left(x^{2}\right),$
(15pqwxasce)

where the coefficients are explicitly given by

$\displaystyle a_{2}=\frac{2aj+3a+2j^{2}+6j+5}{a+2j+4}$ (15pqwxascfa)
$\displaystyle
a_{1}=-\frac{2(a+2j+3)\left(2aj+3a+2j^{2}+6j+5\right)}{(a+2j+2)(a+2j+4)}$
(15pqwxascfb) $\displaystyle a_{0}=\frac{2aj+3a+2j^{2}+6j+5}{a+2j+2}$
(15pqwxascfc) $\displaystyle
r_{3}=\frac{(j+3)(a+j+3)\left(2aj+3a+2j^{2}+6j+5\right)}{a+2j+4}$
(15pqwxascfd) $\displaystyle
r_{2}=\frac{a^{3}+6a^{2}j+12a^{2}+12aj^{2}+46aj+41a+8j^{3}+46j^{2}+82j+42}{\left(2aj+3a+2j^{2}+6j+5\right)^{-1}(a+2j+2)(a+2j+4)}$
(15pqwxascfe) $\displaystyle
r_{1}=\frac{(a+2j+3)\left(2aj+3a+2j^{2}+6j+5\right)\left(2a^{2}+6aj+9a+6j^{2}+18j+10\right)}{(a+2j+2)(a+2j+4)}$
(15pqwxascff) $\displaystyle
r_{0}=\frac{a^{3}+6a^{2}j+6a^{2}+12aj^{2}+26aj+11a+8j^{3}+26j^{2}+22j+6}{\left(2aj+3a+2j^{2}+6j+5\right)^{-1}(a+2j+2)(a+2j+4)}$
(15pqwxascfg) $\displaystyle
r_{-1}=\frac{j(a+j)\left(2aj+3a+2j^{2}+6j+5\right)}{a+2j+2}.$ (15pqwxascfh)

The proof of proposition 3 is in A.

Figure 1: Average purity (15pqwxasat) as a function of subsystem dimensions:
the impact of $\theta$-deformation. The two solid curves represent the cases
of the standard ensembles (7) with no deformations and the other curves
represent the cases of the deformed ensemble (6). Figure 2: Average purity
(15pqwxasat) as a function of the parameters $a$ and $\theta$. The data points
marked by diamond and square shapes represent the special cases of the
undeformed ensembles (7). Figure 3: Average von Neumann entropy (15pqwxasbu)
as a function of subsystem dimensions: the impact of $\theta$-deformation. The
two solid curves represent the cases of the standard ensembles (7) with no
deformations and the other curves represent the cases of the deformed ensemble
(6). Figure 4: Average von Neumann entropy (15pqwxasbu) as a function of the
parameters $a$ and $\theta$. The data points marked by diamond and square
shapes represent the special cases of the undeformed ensembles (7).

######  Remark 5.

The corresponding recurrence relation of the dual polynomial
$q_{j}\left(y^{\theta}\right)$ can be also similarly derived. We defer to a
future publication a more complete study of the orthogonal polynomials
including the recurrence relations for an arbitrary integer $\theta$ and the
resulting Christoffel-Darboux formulas. Note that the recurrence relations for
the undeformed case $\theta=1$ has been established in [34, 37] in terms of a
four-term recurrence.

The second result is a factorization property of the kernels in the following
lemma. This property is useful in simplifying the $k$-point densities for the
higher moment calculations as demonstrated in [32] for the case $\theta=1$.

######  Lemma 1.

For any $\theta>0$, the correlation kernels (15pqwa)–(15pqwd) can be
factorized as

$\displaystyle K_{00}(x,y)+K_{00}(y,x)$ $\displaystyle=$ $\displaystyle
w(x)w(y)$ (15pqwxascga) $\displaystyle K_{01}(x,y)-K_{10}(y,x)$
$\displaystyle=$ $\displaystyle v(x)w(y)$ (15pqwxascgb) $\displaystyle
K_{11}(x,y)+K_{11}(y,x)$ $\displaystyle=$ $\displaystyle-v(x)v(y),$
(15pqwxascgc)

where

$\displaystyle w(x)$ $\displaystyle=$
$\displaystyle\theta\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)x^{\theta
k}}{\Gamma(\theta k+a+2)\Gamma(k+\alpha+1)\Gamma(m-k)k!}$ (15pqwxascgcha)
$\displaystyle v(x)$ $\displaystyle=$ $\displaystyle\rme^{-x}x^{a}-\theta
x^{2a+1}\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,x)x^{\theta k}}{(\theta k+a+1)\Gamma(k+\alpha+1)\Gamma(m-k)k!}$
(15pqwxascgchb)

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

The proof of lemma 1 is in B. Note that the corresponding factorization
property of the special case $\theta=1$ has been established in [36].

##  3 Numerical results

In this section, we perform numerical studies of the main results on the
average entanglement entropies of the interpolating ensemble. We first focus
on the result of average purity in proposition 1. In figure 1, we plot the
numerical values of average purity (15pqwxasat) as a function of the subsystem
dimensions. For different values of $\theta$, we consider both the
$\theta$-deformed Hilbert-Schmidt ensemble assuming $a=0$ and and the
$\theta$-deformed Bures-Hall ensemble assuming $a=-1/2$ as shown in the left
subfigure and the right subfigure, respectively. Note that the choices of $a$
imply equal subsystem dimensions $m=n$ in both cases. The solid curves in
figure 1 describe the behavior of the standard ensembles (7) with no
deformations, whereas the other curves represent the corresponding
$\theta$-deformed ones. It is observed that as the deformation parameter
$\theta$ increases, the values of average purity decrease monotonically
resulting in estimations of entanglement towards more entangled states. The
observation suggests that the proposed interpolating ensemble (6) is indeed
able to continuously interpolate among the possible values of purity by
varying the $\theta$ parameter. It is also observed in figure 1 that for a
given $\theta$ the average purity under the Bures-Hall ensemble tends to an
estimate of more separable state (i.e., a larger purity value) than that of
the Hilbert-Schmidt ensemble. On the other hand, the differences are
diminishing as the dimension increases. This behavior has also been recently
observed in [30].

Since the parameter $a$ of the proposed ensemble (6) can be also considered as
a deformation variable, we wish to understand its impact on the quantum
purity. In figure 2, we plot the average purity (15pqwxasat) as a function of
the parameter $a$ for different values of $\theta$, where the dimension of
subsystem is assumed to be $m=8$. The data points marked by diamond shape for
$\theta=1$ and square shape for $\theta=2$ corresponds to the special case of
Bures-Hall ensemble and Hilbert-Schmidt ensemble, respectively. It is observed
in figure 2 that as $a$ increases, the values of average purity decrease
monotonically indicating more entangled states. In particular, for the cases
$\theta=1$ and $\theta=2$, the average purity is seen to interpolate
continuously among the permissible values (7) of the parameter $a$.

We now turn to the numerical study of the von Neumann entropy in proposition
2. In figure 3, we plot the average von Neumann entropy (15pqwxasbu) as a
function of the subsystem dimensions for different values of $\theta$. We
consider both the $\theta$-deformed Hilbert-Schmidt and Bures-Hall ensembles
with the same values of $a$ as in figure 1. It is seen that as the deformation
parameter $\theta$ increases, the average von Neumann entropy increases
monotonically, which also results in estimations of entanglement towards more
entangled states as in figure 1. In particular, the proposed interpolating
ensemble (6) continuously interpolates among the possible values of the von
Neumann entropy. Similar to figure 1, we also observe in figure 3 that the
average von Neumann entropy under the Bures-Hall ensemble tends to an estimate
of more separable state (i.e., a smaller value of von Neumann entropy) than
that of the Hilbert-Schmidt ensemble. The differences, however, diminish as
the dimension increases, which is in line with the recent observation [30]. To
understand the impact of parameter $a$, we plot in figure 4 the average von
Neumann entropy (15pqwxasbu) as a function of the parameter $a$ for different
values of $\theta$. The dimension of subsystem is also assumed to be $m=8$.
Similarly as observed in figure 2, the values of average von Neumann entropy
increase monotonically indicating more entangled states as $a$ increases.
Finally, we point out that various other numerical simulations have been
performed, where the same relative behavior as discussed in above four figures
persists.

##  4 Conclusion and outlook

In this work, we proposed and studied a generalized measure that interpolates
between the two major measures of density matrices - the Hilbert-Schmidt
ensemble and the Bures-Hall ensemble. In particular, we derived the mean
quantum purity and the mean von Neumann entropy over the proposed ensemble as
summarized in proposition $1$ and proposition $2$, respectively. We also
studied the corresponding orthogonal polynomial system including its
recurrence relation in proposition $2$ and a kernel factorization property in
lemma $1$. Extensive numerical simulations show that the proposed ensemble
provides additional power in estimating the degree of entanglement by varying
the deformation parameters. Future work includes further study of the
statistical information of the ensemble such as higher order moments of
entropies, fidelity, and volumes as well as further study of the underlying
orthogonal polynomials.

## Acknowledgments

We wish to thank Shi-Hao Li for correspondence. The work of Lu Wei is
supported in part by the U.S. National Science Foundation (CNS-2006612).

##  Appendix A Proof of proposition 3

######  Proof.

We give a proof of the recurrence relation for the general bi-orthogonal
system $\\{p_{j},q_{j}\\}^{\infty}_{j=0}$ for all positive integer
$\theta\in\mathbb{N}$ with respect to the generalized weight function

$\widetilde{W}(x,y)=\frac{x^{a}y^{b}\rme^{-x-y}}{x+y},$ (15pqwxasci)

and all $\Re{(a,b)}>-1$. This reduces to the weight function (15pqwxy) when
$b=a+1$. The essence of our proof is a generalization of that used in [34]
from the rank one shift condition to a rank-$\theta$ condition. Towards the
end of proof we will specialize to $\theta=2$ only when evaluating the
explicit recurrence coefficients, without loss of generality.

To begin we define the parameter

$\beta=\frac{a+b+1}{\theta},$ (15pqwxascj)

which reduces to $\alpha+1$ in (15pqr) when $b=a+1$. We will see that the
simplest form of the recurrence will be expressed in terms of the bi-
orthogonal polynomials

$\tilde{p}_{j}(x)=\frac{(-1)^{j}}{\sqrt{2}}p_{j}(x)$ (15pqwxasck)

of argument $x$, which we call hybrid polynomials, rather than those defined
in (15pqwxa) of argument $x^{\theta}$. The normalized bi-orthogonal
polynomials are denoted by $\\{P_{j},Q_{j}\\}^{\infty}_{j=0}$ and the monic
system by $\\{p_{j},q_{j}\\}^{\infty}_{j=0}$ with the inter-relations

$\displaystyle P_{j}(x)=\frac{1}{\sqrt{h_{j}}}p_{j}(x)$ (15pqwxascl)
$\displaystyle p_{j}(x)=(-1)^{j}j!\frac{\Gamma(a+1+\theta
j)\Gamma(j+\beta)}{\Gamma(2j+\beta)}\tilde{p}_{j}(x)$ (15pqwxascm)

along with the squared norm

$h_{j}=\theta^{-1}\frac{(j!)^{2}\Gamma^{2}(j+\beta)}{\Gamma(2j+\beta)\Gamma(2j+\beta+1)}\Gamma(a+1+\theta
j)\Gamma(b+1+\theta j).$ (15pqwxascn)

The bi-moments are defined and evaluated as

$I_{k,l}=\int_{\mathbb{R}_{+}^{2}}\frac{e^{-x-y}}{x+y}x^{a+\theta
k}y^{b+\theta l}\,\mathrm{d}x\,\mathrm{d}y=\frac{\Gamma(a+1+\theta
k)\Gamma(b+1+\theta l)}{a+b+1+\theta(k+l)},\quad k,l=0,1,\ldots.$ (15pqwxasco)

Now a key observation is that $x+y$ divides
$x^{\theta}-e^{\pi\imath\theta}y^{\theta}$ without remainder and therefore

$I_{k+1,l}-e^{\pi\imath\theta}I_{k,l+1}=\sum_{s=0}^{\theta-1}(-1)^{s}\Gamma(a+\theta-s+\theta
k)\Gamma(b+1+s+\theta l).$ (15pqwxascp)

Employing semi-infinite matrices $\mathbf{I}=(I_{k,l})_{k,l\geq 0}$ and the
shift matrix $\mathbf{\Lambda}=(\delta_{k+1,l})_{k,l\geq 0}$, this is written
as the rank-$\theta$ decomposition

$\mathbf{\Lambda}\mathbf{I}-e^{\pi
i\theta}\mathbf{I}\mathbf{\Lambda}^{T}=\sum_{s=0}^{\theta-1}(-1)^{s}\bm{\alpha}_{s}\bm{\beta}_{s}^{T}$
(15pqwxascq)

with semi-infinite column vectors $\bm{\alpha}_{s}=(\Gamma(a+\theta-s+\theta
k))_{k\geq 0}$, $\bm{\beta}_{s}=(\Gamma(b+1+s+\theta k))_{k\geq 0}$. Now
denote the vectors of monomials $\bm{x}=(x^{k})_{k\geq 0}$,
$\bm{y}=(y^{k})_{k\geq 0}$ and the vector of normalised bi-orthogonal
polynomials $\bm{P}=(P_{k}(x))_{k\geq 0}$, $\bm{Q}=(Q_{k}(y))_{k\geq 0}$.
These basis vectors are related by lower triangular matrices $\bm{S}_{P}$ and
$\bm{S}_{Q}$ as

$\bm{P}=\bm{S}_{P}\bm{x},\qquad\bm{Q}=\bm{S}_{Q}\bm{y}.$ (15pqwxascr)

Thus, we have the L-U decomposition of the bi-moment matrix

$\mathbf{I}=\bm{S}_{P}^{-1}\left(\bm{S}_{Q}^{-1}\right)^{T}.$ (15pqwxascs)

The monomial bases are the right and left eigenvectors of the shift matrix and
its transpose

$\mathbf{\Lambda}\bm{x}=x\bm{x},\qquad\bm{y}^{T}\mathbf{\Lambda}^{T}=y\bm{y}^{T}.$
(15pqwxasct)

Thinking of the action of a multiplication operator on the normalized
polynomial bases, it is clear that it can be written generally as

$x\bm{P}=\bm{X}\bm{P},\qquad y\bm{Q}^{T}=\bm{Q}^{T}\bm{Y}^{T}$ (15pqwxascu)

for some lower Hessenberg matrices $\bm{X}$ and $\bm{Y}$. These can be
explicitly calculated although they are not needed and so we will not exhibit
them here. These multiplication matrices now have the decomposition

$\bm{X}=\bm{S}_{P}\bm{\Lambda}\bm{S}_{P}^{-1},\qquad\bm{Y}=\bm{S}_{Q}\bm{\Lambda}\bm{S}_{Q}^{-1},$
(15pqwxascv)

and by premultiplying by $\bm{S}_{P}$ and postmultiplying (15pqwxascq) by
$\bm{S}_{Q}^{T}$ we deduce

$\bm{X}-e^{\pi
i\theta}\bm{Y}^{T}=\sum_{s=0}^{\theta-1}(-1)^{s}\bm{\pi}_{s}\bm{\eta}_{s}^{T},$
(15pqwxascw)

where $\bm{\pi}_{s}=\bm{S}_{P}\bm{\alpha}_{s}$ and
$\bm{\eta}_{s}=\bm{S}_{Q}\bm{\beta}_{s}$.

In order to proceed further it is useful to construct rank-$\theta$
annihilators, along the lines that was done in [34] for $\theta=1$, for the
right-hand side of (15pqwxascw). We will do this recursively in $\theta$
steps, but will show the details for $\theta=2$ only. For any vector
$\bm{\pi}$ let us construct the semi-infinite diagonal matrix
$\bm{D}_{\bm{\pi}}$ so that $\bm{\pi}=\bm{D}_{\bm{\pi}}\bm{1}$. From the
knowledge that the unit vector $\bm{1}$ (or any constant vector) is left-
annihilated by $\bm{\Lambda}-Id$ we can left-annihilate the $s=0$ term on the
right-hand side of (15pqwxascw) by premultiplying with
$\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\pi}_{0}}^{-1}$. To annihilate the
remaining $s=1$ term we need to calculate
$\bm{\psi}=\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\pi}_{0}}^{-1}\bm{D}_{\bm{\pi}_{1}}\bm{1}$.
A simple calculation gives the components of $\bm{\psi}$ as

$\psi_{j}=\frac{\pi_{1,j}}{\pi_{0,j}}-\frac{\pi_{1,j+1}}{\pi_{0,j+1}},\qquad
j=0,1,\dots,$ (15pqwxascx)

assuming $\pi_{0,j}\neq 0$. Our required second left-annihilator is therefore
$\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\psi}}^{-1}$ and so the composite
operator is the second order difference operator

$\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\psi}}^{-1}\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\pi}_{0}}^{-1}.$
(15pqwxascy)

This recursion can be repeated up to $\theta$ levels leaving us with a
$\theta$-order difference operator in the general case, modulo the non-
vanishing condition given above. Before we compute this our final step in
deriving the recurrence relation for $P_{j}$ is to put some of these pieces
together. Let us act on the first equation of (15pqwxascu) with this operator
- doing so on the left-hand side gives

$x\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\psi}}^{-1}\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\pi}_{0}}^{-1}\bm{P},$
(15pqwxascz)

whereas acting on the right-hand side gives

$\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\psi}}^{-1}\left(\bm{\Lambda}-Id\right)\bm{D}_{\bm{\pi}_{0}}^{-1}\bm{X}\bm{P}=\bm{A}\bm{P}$
(15pqwxasda)

that defines a banded matrix $\bm{A}$. This banded matrix has non-zero
elements only for $\theta+1$ super-diagonals above the diagonal - the
$\theta$-order difference operator adds $\theta$ super-diagonals to the
initial single one of the lower Hessenberg $\bm{X}$ \- and a single sub-
diagonal - this operator does not add any additional sub-diagonals to the
upper Hessenberg matrix $\bm{Y}^{T}$. For $\theta=2$ the second order
difference operator acting on $\bm{P}$ has $j$-th component

$\displaystyle\frac{\pi_{0,j+1}}{\pi_{0,j+2}\pi_{1,j+1}-\pi_{0,j+1}\pi_{1,j+2}}P_{j+2}+\frac{\pi_{0,j+1}}{\pi_{0,j+1}\pi_{1,j}-\pi_{0,j}\pi_{1,j+1}}P_{j}$
$\displaystyle+\frac{\pi_{0,j+1}\left(\pi_{0,j}\pi_{1,j+2}-\pi_{0,j+2}\pi_{1,j}\right)}{\left(\pi_{0,j+1}\pi_{1,j}-\pi_{0,j}\pi_{1,j+1}\right)\left(\pi_{0,j+2}\pi_{1,j+1}-\pi_{0,j+1}\pi_{1,j+2}\right)}P_{j+1}.$
(15pqwxasdb)

In the case at hand we compute the components of $\bm{\pi}_{0}$ and
$\bm{\pi}_{1}$ to be

$\displaystyle\pi_{0,j}$
$\displaystyle=\frac{1}{\sqrt{h_{j}}}j!\frac{\Gamma(a+1+2j)\Gamma(\beta+j)}{\Gamma(\beta+2j)}\left(a+1+2j(j+\beta)\right)$
(15pqwxasdc) $\displaystyle\pi_{1,j}$
$\displaystyle=\frac{1}{\sqrt{h_{j}}}j!\frac{\Gamma(a+1+2j)\Gamma(\beta+j)}{\Gamma(\beta+2j)}.$
(15pqwxasdd)

Thus, we arrive at a structure for the recurrence for our hybrid system

$\displaystyle
x\left(a_{2}\tilde{p}_{j+2}(x)+a_{1}\tilde{p}_{j+1}(x)+a_{0}\tilde{p}_{j}(x)\right)=$
$\displaystyle
r_{3}\tilde{p}_{j+3}(x)+r_{2}\tilde{p}_{j+2}(x)+r_{1}\tilde{p}_{j+1}(x)+r_{0}\tilde{p}_{j}(x)+r_{-1}\tilde{p}_{j-1}(x),$
(15pqwxasde)

along with

$\displaystyle a_{2}$
$\displaystyle=\frac{\left(a+1+2(j+1)(j+1+\beta)\right)}{2j+3+\beta}$
(15pqwxasdf) $\displaystyle a_{1}$
$\displaystyle=\frac{2\left(a+1+2(j+1)(j+1+\beta)\right)(2j+2+\beta)}{(2j+1+\beta)(2j+3+\beta)}$
(15pqwxasdg) $\displaystyle a_{0}$
$\displaystyle=\frac{\left(a+1+2(j+1)(j+1+\beta)\right)}{2j+1+\beta}.$
(15pqwxasdh)

On the other hand, the $r$ coefficients can be deduced in a number of ways,
such as employing the explicit series form and peeling off the leading terms
from the highest degree ($j+3$) down in four successive iterations. Either way
we find

$\displaystyle
r_{3}=-\frac{\left(a+1+2(j+1)(j+1+\beta)\right)(j+3)(a+2j+5)(a+2j+6)\left(\beta+j+2\right)}{\left(\beta+2j+3\right)\left(\beta+2j+4\right)\left(\beta+2j+5\right)}$
(15pqwxasdi) $\displaystyle r_{2}$ $\displaystyle=$
$\displaystyle\frac{\left(a+1+2(j+1)(j+1+\beta)\right)}{(\beta+2j+1)(\beta+2j+3)}\Bigg{(}\\!\\!-(j+2)(a+2j+3)(a+2j+4)\left(\beta+j+1\right)$
(15pqwxasdj)
$\displaystyle+\frac{(j+3)(a+2j+5)(a+2j+6)(\beta+j+2)(\beta+2j+1)}{\beta+j+5}\Bigg{)}$
$\displaystyle r_{1}$ $\displaystyle=$
$\displaystyle-(\beta+2j+2)\left(a+1+2(j+1)(j+1+\beta)\right)\Bigg{(}\frac{1}{8}\beta(\beta-2)(a-\beta+1)(a-\beta+2)$
(15pqwxasdk)
$\displaystyle\times\left(\frac{1}{\beta+2j+4}-\frac{1}{\beta+2j}\right)+\frac{1}{4}\Big{(}\beta^{4}-(2a+5)\beta^{3}+(a(a+7)+8)\beta^{2}-2a(a+2)\beta$
$\displaystyle-a(a+7)+\beta-9\Big{)}\left(\frac{1}{\beta+2j+1}-\frac{1}{\beta+2j+3}\right)+\frac{3}{2}\Bigg{)}$
$\displaystyle r_{0}$ $\displaystyle=$
$\displaystyle\left(a+1+2(j+1)(j+1+\beta)\right)\Bigg{(}\frac{4\beta-2a-3}{2}+\frac{(\beta-1)^{2}(a-\beta)(a-\beta+1)}{8(\beta+2j-1)}$
(15pqwxasdl)
$\displaystyle-\frac{\beta^{4}-(2a+5)\beta^{3}+(a(a+7)+8)\beta^{2}-2a(a+2)\beta-a(a+7)+\beta-9}{4(\beta+2j+1)}$
$\displaystyle+\frac{(\beta-3)(\beta+1)(a-\beta+2)(a-\beta+3)}{8(\beta+2j+3)}+2j\Bigg{)}$
$\displaystyle r_{-1}$ $\displaystyle=$
$\displaystyle-\frac{j(\beta+j-1)(a-2\beta-2j+1)(a-2(\beta+j-1))}{(\beta+2j-1)(\beta+2j)(\beta+2j+1)}$
(15pqwxasdm) $\displaystyle\times\left(a+1+2(j+1)(j+1+\beta)\right).$

Finally, upon the substitution $x\to x^{2}$ and specialization $b=a+1$,
$\beta=\alpha+1$ while keeping in mind the relation (15pqwxasck) prove the
proposition 3. $\square$ $\blacksquare$

##  Appendix B Proof of lemma 1

######  Proof.

The starting point of the proof is the orthogonal polynomial forms of the
kernel functions (15pqwa)–(15pqwd). To show (15pqwxascga), we first represent
(15pqwa) via (15pqwxa) and (15pqwxb) as

$\displaystyle K_{00}(x,y)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\sum_{i=0}^{m-1}\frac{\theta(-1)^{i+k}x^{\theta
k}y^{\theta i}}{\Gamma(\theta i+a+2)\Gamma(\theta
k+a+1)\Gamma(i+\alpha+1)\Gamma(k+\alpha+1)i!k!}$ (15pqwxasdo)
$\displaystyle\times\sum_{j=i}^{m-1}(2j+\alpha+1)\frac{\Gamma(j+i+\alpha+1)\Gamma(j+k+\alpha+1)}{\Gamma(j-i+1)\Gamma(j-k+1)}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\sum_{i=0}^{m-1}\frac{\theta(-1)^{i+k}x^{\theta
k}y^{\theta i}}{\Gamma(\theta i+a+2)\Gamma(\theta
k+a+1)\Gamma(i+\alpha+1)\Gamma(k+\alpha+1)i!k!}$
$\displaystyle\times\frac{\Gamma(i+\alpha+m+1)\Gamma(k+\alpha+m+1)}{(i+k+\alpha+1)\Gamma(m-i)\Gamma(m-k)},$

where the last step is obtained by lemma 4.1 in [35]. We then have

$\displaystyle K_{00}(x,y)+K_{00}(y,x)$ (15pqwxasdq) $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\sum_{i=0}^{m-1}\frac{\theta^{2}(-1)^{i+k}x^{\theta
k}y^{\theta i}}{\Gamma(\theta i+a+2)\Gamma(\theta
k+a+2)\Gamma(i+\alpha+1)\Gamma(k+\alpha+1)i!k!}$
$\displaystyle\times\frac{\Gamma(i+\alpha+m+1)\Gamma(k+\alpha+m+1)}{(i+k+\alpha+1)\Gamma(m-i)\Gamma(m-k)}\left(i+k+\frac{2(a+1)}{\theta}\right)$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{\theta(-1)^{k}\Gamma(k+\alpha+m+1)x^{\theta
k}}{\Gamma(\theta k+a+2)\Gamma(k+\alpha+1)\Gamma(m-k)k!}$ (15pqwxasds)
$\displaystyle\times\sum_{i=0}^{m-1}\frac{\theta(-1)^{i}\Gamma(i+\alpha+m+1)y^{\theta
i}}{\Gamma(\theta i+a+2)\Gamma(i+\alpha+1)\Gamma(m-i)i!}=w(x)w(y),$

which establishes (15pqwxascga). To show (15pqwxascgb), we insert (15pqwa)
into (15pqwb) and (15pqwc) that gives

$\displaystyle K_{01}(x,y)-K_{10}(y,x)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{x^{a}v^{a}\rme^{-x-v}}{x+v}\left(vK_{00}(y,v)-xK_{00}(v,y)\right)\,\mathrm{d}v$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\sum_{i=0}^{m-1}\frac{\theta(-1)^{i+k}\Gamma(i+\alpha+m+1)\Gamma(k+\alpha+m+1)\left(\Gamma(m-i)\Gamma(m-k)i!k!\right)^{-1}y^{\theta
i}}{(i+k+\alpha+1)\Gamma(\theta i+a+2)\Gamma(\theta
k+a+2)\Gamma(i+\alpha+1)\Gamma(k+\alpha+1)}$
$\displaystyle\times\int_{0}^{\infty}\frac{x^{a}v^{a}\rme^{-x-v}}{x+v}v^{\theta
k}\left((\theta i+a+1)v-x(\theta k+a+1)\right)\,\mathrm{d}v$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{m-1}\frac{\theta(-1)^{i}\Gamma(i+\alpha+m+1)\rme^{-x}x^{a}y^{\theta
i}}{\Gamma(\theta
i+a+1)\Gamma(i+\alpha+1)\Gamma(m-i)i!}\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\left(\Gamma(m-k)k!\right)^{-1}}{(i+k+\alpha+1)(\theta
k+a+1)\Gamma(k+\alpha+1)}$ $\displaystyle-$
$\displaystyle\sum_{i=0}^{m-1}\frac{\theta^{2}(-1)^{i}\Gamma(i+\alpha+m+1)x^{2a+1}y^{\theta
i}}{\Gamma(\theta
i+a+2)\Gamma(i+\alpha+1)\Gamma(m-i)i!}\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,x)x^{\theta k}}{(\theta k+a+1)\Gamma(k+\alpha+1)\Gamma(m-k)k!}$
$\displaystyle=$ $\displaystyle\left(\rme^{-x}x^{a}-\theta
x^{2a+1}\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,x)x^{\theta k}}{(\theta k+a+1)\Gamma(k+\alpha+1)\Gamma(m-k)k!}\right)$
$\displaystyle\times\sum_{i=0}^{m-1}\frac{\theta(-1)^{i}\Gamma(i+\alpha+m+1)y^{\theta
i}}{\Gamma(\theta i+a+2)\Gamma(i+\alpha+1)\Gamma(m-i)i!}=v(x)w(y),$

where the second to last equality is obtained by the identity

$\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)}{(i+k+\alpha+1)(\theta
k+a+1)\Gamma(k+\alpha+1)\Gamma(m-k)k!}=\frac{1}{\theta i+a+1}.$ (15pqwxasdt)

This identity is established by the fact that the sum can be written in terms
of a unit argument terminating hypergeometric function of Saalschützian type
[42] as

${}_{4}F_{3}\left(1-m,\frac{\alpha+1}{2},\alpha+m+1,i+\alpha+1;\frac{\alpha+3}{2},\alpha+1,i+\alpha+2;1\right),$
(15pqwxasdu)

which in general admits [42]

${}_{4}F_{3}(-n,b,c,d;b+1,c-l,d-m;1)=\frac{n!(b-c+1)_{l}(b-d+1)_{m}}{(b+1)_{n}(1-c)_{l}(1-d)_{m}},~{}~{}~{}~{}0\leq
l+m\leq n$

with

$(a)_{n}=\frac{\Gamma(a+n)}{\Gamma(a)}$ (15pqwxasdv)

denoting the Pochhammer’s symbol. This completes the proof of (15pqwxascgb).
To show (15pqwxascgc), we insert (15pqwa) into (15pqwd) that leads to

$\displaystyle K_{11}(x,y)+K_{11}(y,x)$ $\displaystyle=$ $\displaystyle
x^{a}y^{a}\rme^{-x-y}\left(\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{v^{a}w^{a}\rme^{-v-w}}{(x+w)(y+v)}\left(ywK_{00}(v,w)+xvK_{00}(w,v)\,\mathrm{d}w\,\mathrm{d}v\right)-1\right)$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\sum_{i=0}^{m-1}\frac{\theta(-1)^{i+k}\Gamma(i+\alpha+m+1)\Gamma(k+\alpha+m+1)\left(\Gamma(m-i)\Gamma(m-k)\right)^{-1}}{(i+k+\alpha+1)\Gamma(\theta
i+a+2)\Gamma(\theta k+a+1)\Gamma(i+\alpha+1)\Gamma(k+\alpha+1)i!k!}$
$\displaystyle\times
x^{a}y^{a}\rme^{-x-y}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\frac{v^{a}w^{a}\rme^{-v-w}\left(yv^{\theta
k}w^{\theta i+1}+xw^{\theta k}v^{\theta
i+1}\right)}{(x+w)(y+v)}\,\mathrm{d}w\,\mathrm{d}v-x^{a}y^{a}\rme^{-x-y}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{\theta(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,x)x^{\theta k+2a+1}y^{a}e^{-y}}{\Gamma(k+\alpha+1)\Gamma(m-k)k!}$
$\displaystyle\times\sum_{i=0}^{m-1}\frac{(-1)^{i}\Gamma(i+\alpha+m+1)}{(i+k+\alpha+1)(\theta
i+a+1)\Gamma(i+\alpha+1)\Gamma(m-i)i!}$
$\displaystyle+\sum_{k=0}^{m-1}\frac{\theta(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,y)y^{\theta k+2a+1}x^{a}e^{-x}}{\Gamma(k+\alpha+1)\Gamma(m-k)k!}$
$\displaystyle\times\sum_{i=0}^{m-1}\frac{(-1)^{i}\Gamma(i+\alpha+m+1)}{(i+k+\alpha+1)(\theta
i+a+1)\Gamma(i+\alpha+1)\Gamma(m-i)i!}$ $\displaystyle-\theta
x^{2a+1}\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,x)x^{\theta k}}{(\theta k+a+1)\Gamma(k+\alpha+1)\Gamma(m-k)k!}$
$\displaystyle\times\theta
y^{2a+1}\sum_{i=0}^{m-1}\frac{(-1)^{i}\Gamma(i+\alpha+m+1)\Gamma(-\theta
i-a,y)y^{\theta i}}{(\theta
i+a+1)\Gamma(i+\alpha+1)\Gamma(m-i)i!}-x^{a}y^{a}\rme^{-x-y}$ $\displaystyle=$
$\displaystyle-\left(\rme^{-x}x^{a}-\theta
x^{2a+1}\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,x)x^{\theta k}}{(\theta k+a+1)\Gamma(k+\alpha+1)\Gamma(m-k)k!}\right)$
$\displaystyle\times\left(\rme^{-y}y^{a}-\theta
y^{2a+1}\sum_{k=0}^{m-1}\frac{(-1)^{k}\Gamma(k+\alpha+m+1)\Gamma(-\theta
k-a,y)y^{\theta k}}{(\theta
k+a+1)\Gamma(k+\alpha+1)\Gamma(m-k)k!}\right)=-v(x)v(y),$

where the second to last equality is obtained by applying twice the identity
(15pqwxasdt). This completes the proof of lemma 1. $\square$ $\blacksquare$

## References

## References

  * [1] von Neumann J 1927 Wahrscheinlichkeitstheoretischer aufbau der quantenmechanik Nachr. Ges. Wiss. Göttingen 1 245
  * [2] Bengtsson I and Życzkowski K 2017 Geometry of Quantum States: An Introduction to Quantum Entanglement 2nd ed (Cambridge: Cambridge University Press)
  * [3] Hall M J W 1998 Random quantum correlations and density operator distributions Phys. Lett.A 242 123-9
  * [4] Życzkowski K and Sommers H-J 2001 Induced measures in the space of mixed quantum states J. Phys. A: Math. Gen.34 7111-25
  * [5] Lubkin E 1978 Entropy of an n-system from its correlation with a k-reservoir J. Math. Phys.19 1028-31
  * [6] Giraud O 2007 Distribution of bipartite entanglement for random pure states 40 2793-801
  * [7] Page D N 1993 Average entropy of a subsystem Phys. Rev. Lett.71 1291-4
  * [8] Foong S K and Kanno S 1994 Proof of Page’s conjecture on the average entropy of a subsystem Phys. Rev. Lett.72 1148-51
  * [9] Sánchez-Ruiz J 1995 Simple proof of Page’s conjecture on the average entropy of a subsystem Phys. Rev.E 52 5653-5
  * [10] Vivo P, Pato M P and Oshanin G 2016 Random pure states: Quantifying bipartite entanglement beyond the linear statistics Phys. Rev.E 93 052106
  * [11] Wei L 2017 Proof of Vivo-Pato-Oshanin’s conjecture on the fluctuation of von Neumann entropy Phys. Rev.E 96 022106
  * [12] Wei L 2020 Skewness of von Neumann entanglement entropy 53 075302
  * [13] Malacarne L C, Mendes R S and Lenzi E K 2002 Average entropy of a subsystem from its average Tsallis entropy Phys. Rev.E 65 046131
  * [14] Wei L 2019 On the exact variance of Tsallis entanglement entropy in a random pure state Entropy 21 539
  * [15] Życzkowski K, Penson K A, Nechita I and Collins B 2011 Generating random density matrices J. Math. Phys.52 06220
  * [16] Nadal C, Majumdar S N and Vergassola M 2011 Statistical distribution of quantum entanglement for a random bipartite state J. Stat. Phys. 142 403-38
  * [17] Chen Y, Liu D-Z and Zhou D-S 2010 Smallest eigenvalue distribution of the fixed-trace Laguerre beta-ensemble 43 315303
  * [18] Liu Z-W, Lloyd S, Zhu E Y and Zhu H 2018 Generalized entanglement entropies of quantum designs Phys. Rev. Lett.120 130502
  * [19] Hayden P, Leung D W and Winter A 2006 Aspects of generic entanglement Commun. Math. Phys. 265 95-117
  * [20] Aubrun G, Szarek S J and Ye D 2014 Entanglement thresholds for random induced states Comm. Pure Appl. Math. 67 129-71
  * [21] Szarek S J 2005 Volume of separable states is super-doubly-exponentially small in the number of qubits Phys. Rev.A 72 032304
  * [22] Kendon V M, Życzkowski K and Munro W J 2002 Bounds on entanglement in qudit subsystems Phys. Rev.A 66 062310
  * [23] Aubrun G 2012 Partial transposition of random states and non-centered semicircular distributions Random Matrices: Theory Appl. 1 1-29
  * [24] 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
  * [25] Sommers H-J and Życzkowski K 2004 Statistical properties of random density matrices J. Phys. A: Math. Gen.37 8457-66
  * [26] Osipov V A, Sommers H-J and Życzkowski K 2010 Random Bures mixed states and the distribution of their purity 43 055302
  * [27] Borot G and Nadal C 2012 Purity distribution for generalized random Bures mixed states 45 075209
  * [28] Życzkowski K and Sommers H-J 2005 Average fidelity between random quantum states Phys. Rev.A 71 032313
  * [29] Ye D 2009 On the Bures volume of separable quantum states J. Math. Phys.50 083502
  * [30] Sarkar A and Kumar S 2019 Bures-Hall ensemble: Spectral densities and average entropies 52 295203
  * [31] Wei L 2020 Proof of Sarkar-Kumar’s conjectures on average entanglement entropies over the Bures-Hall ensemble 53 235203
  * [32] Wei L 2020 Exact variance of von Neumann entanglement entropy over the Bures-Hall measure Phys. Rev.E 102 062128
  * [33] Bertola M, Gekhtman M and Szmigielski J 2009 The Cauchy two-matrix model Commun. Math. Phys. 287 983-1014
  * [34] Bertola M, Gekhtman M and Szmigielski J 2010 Cauchy biorthogonal polynomials J. Approx. Theory 162 832-67
  * [35] 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
  * [36] Forrester P J and Kieburg M 2016 Relating the Bures measure to the Cauchy two-matrix model Commun. Math. Phys. 342 151-87
  * [37] Li C and Li S-H 2019 The Cauchy two-matrix model, C-Toda lattice and CKP hierarchy J. Nonlinear Sci. 29 3-27
  * [38] Mehta M L 2004 Random Matrices 3rd ed (Amsterdam: Elsevier)
  * [39] Forrester P J 2010 Log-gases and Random Matrices (Princeton: Princeton University Press)
  * [40] Forrester P J and Li S-H 2019 Fox H-kernel and $\theta$-deformation of the Cauchy two-matrix model and Bures ensemble Int. Math. Res. Not. 1-34
  * [41] Mathai A M, Saxena R K and Haubold H J 2010 The H-Function: Theory and Applications (Heidelberg: Springer)
  * [42] 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)