  * November 2021

# Second-order statistics of fermionic Gaussian states

Youyi Huang and Lu Wei  Department of Computer Science, Texas Tech University,
Texas 79409, USA<EMAIL_ADDRESS>
###### Abstract

We study the statistical behavior of entanglement in quantum bipartite systems
over fermionic Gaussian states as measured by von Neumann entropy and
entanglement capacity. The focus is on the variance of von Neumann entropy and
the mean entanglement capacity that belong to the so-defined second-order
statistics. The main results are the exact yet explicit formulas of the two
considered second-order statistics for fixed subsystem dimension differences.
We also conjecture the exact variance of von Neumann entropy valid for
arbitrary subsystem dimensions. Based on the obtained results, we analytically
study the numerically observed phenomena of Gaussianity of von Neumann entropy
and linear growth of average capacity.

Keywords: fermionic Gaussian states, quantum entanglement, von Neumann
entropy, entanglement capacity, random matrix theory

##  1 Introduction

Quantum information theory aims to construct the theoretical groundwork of
quantum technologies such as quantum computing and quantum communications. In
the exploitation of the revolutionary advances of quantum mechanical systems,
the understanding of the non-classical phenomenon of quantum entanglement is
quite crucial. Quantum entanglement is also the resource and medium that
enable quantum technologies.

In this work, we study the statistical behavior of entanglement of quantum
bipartite systems of fermionic Gaussian states. In the literature, the degree
of entanglement as measured by entanglement entropies such as von Neumann
entropy, quantum purity, and Tsallis entropy has been studied over different
ensembles, e.g. the Hilbert-Schmidt ensemble [1, 2, 5, 3, 4, 6, 7, 8, 9, 10,
11] and the Bures-Hall ensemble [12, 13, 14, 15, 16, 17, 18]. Similar
investigations are now being carried out over the fermionic Gaussian ensemble,
where the average von Neumann entropy is obtained in [19]. As an important
step towards further understanding its statistical distribution, we aim to
derive the variance. The variance describes the fluctuation of the entropy
around its mean value. It also provides crucial information such as if the
average entropy is typical. In addition, we study the capacity of entanglement
over the fermionic Gaussian ensemble. Entanglement capacity is another
entanglement measure that possesses other distinctive properties as compared
to entanglement entropies [20].

We now introduce the formulation that leads to the fermionic Gaussian states
as follows. Consider a phase space $V$ in $\mathbb{R}^{2N}$ given by a
$2N$-dimensional real vector space that corresponds to a system of $N$
fermionic degrees of freedom. Its dual space $V^{*}$ in $\mathbb{R}^{2N}$ and
$V$ are equipped with a positive-definite bilinear form $G^{ab}:V^{*}\times
V^{*}\to\mathbb{R}$ and its inverse $G_{ab}^{-1}:V\times V\to\mathbb{R}$.
Observables in quantum theory are linear operators on a Hilbert space. Here,
it is convenient to introduce an operator-valued vector $\hat{\xi}^{a}$, i.e.,
a quantization map [21, 22]

$\\{\hat{\xi}^{a},\hat{\xi}^{b}\\}=G^{ab}\mathbb{I}$ (1)

with
$\\{\hat{\xi}^{a},\hat{\xi}^{b}\\}=\hat{\xi}^{a}\hat{\xi}^{b}+\hat{\xi}^{b}\hat{\xi}^{a}$
denoting the anti-commutation relation and $\mathbb{I}$ being an identity
operator. The operators $\hat{\xi}^{a}$ can be represented via the bases of
the Majorana operators or fermionic creation and annihilation operators as
discussed in details in [22]. A fermionic Gaussian state is labeled by its
covariance matrix

$\Omega^{ab}=\Bra{\Omega}\hat{\xi}^{a}\hat{\xi}^{b}-\hat{\xi}^{b}\hat{\xi}^{a}\Ket{\Omega},$
(2)

which is clearly anti-symmetric.

Given a fermionic Gaussian state, the phase space $V$ of a system of $N$
fermionic degrees of freedom can be decomposed into two orthogonal
complementary subsystems $A$ and $B$ of the dimensions $m$ and $n$,
respectively,

$V=A\oplus B.$ (3)

Here, $m+n=N$ and we assume $m\leq n$ without loss of generality. One can
always choose an orthonormal basis
$\left(q_{1}^{A},p_{1}^{A},\dots,q_{m}^{A},p_{m}^{A},q_{1}^{B},p_{1}^{B},\dots,q_{n}^{B},p_{n}^{B}\right)$
such that the covariance matrix is transformed into the form

$\Omega=\left(\begin{array}[]{cc}\Omega_{A}&\Omega_{AB}\\\
\Omega_{BA}&\Omega_{B}\\\ \end{array}\right).$ (4)

In the above decomposition, the mode $(q_{i}^{A},p_{i}^{A})$ is entangled with
the mode $(q_{i}^{B},p_{i}^{B})$ for each $i=1,\dots,m$ and the $2m\times 2m$
subsystem $\Omega_{A}$ and the $2n\times 2n$ subsystem $\Omega_{B}$ have the
same spectrum except for the $2n-2m$ eigenvalues of the values $\pm\imath$,
$\imath=\sqrt{-1}$. Namely, $\Omega_{A}$ and $\Omega_{B}$ contain the same
information regarding the entanglement between the two subsystems. The matrix
$\Omega_{A}$ takes the block diagonal form [21, 22]

$\Omega_{A}=\left(\begin{array}[]{ccc}\cos(2r_{1})\mathbb{A}&\dots&0\\\
\vdots&\ddots&\vdots\\\ 0&\dots&\cos(2r_{m})\mathbb{A}\\\ \end{array}\right),$
(5)

where

$\mathbb{A}=\left(\begin{array}[]{cc}0&1\\\ -1&0\\\ \end{array}\right)$ (6)

and the amount of entanglement is now encoded in the parameters
$r_{i}\in[0,\pi/4]$. Therefore, by defining the new variables
$x_{i}=\cos(2r_{i})$, $i=1,\dots,m$, the von Neumann entropy of the considered
fermionic system can be represented as [19]

$\displaystyle S$ $\displaystyle=$
$\displaystyle-\tr\left(\frac{\mathrm{I}_{m}+\imath\Omega_{A}}{2}\ln\frac{\mathrm{I}_{m}+\imath\Omega_{A}}{2}\right)$
(7) $\displaystyle=$ $\displaystyle-\sum_{i=1}^{m}v(x_{i}),\qquad
x_{i}\in[0,1],$ (8)

where

$v(x)=\frac{1-x}{2}\ln\frac{1-x}{2}+\frac{1+x}{2}\ln\frac{1+x}{2},$ (9)

and $\mathrm{I}_{m}$ is an identity matrix of dimension $m$. The density of
the eigenvalues $x_{i},i=1,\dots,m$, follows a non-standard Jacobi unitary
ensemble, which is proportional to [19]

$\prod_{1\leq i<j\leq
m}\left(x_{i}^{2}-x_{j}^{2}\right)^{2}\prod_{i=1}^{m}\left(1-x_{i}^{2}\right)^{n-m}.$
(10)

The above random matrix ensemble is referred to as the fermionic Gaussian
ensemble. Its eigenvalues are supported in

$x_{i}\in[0,1]$ (11)

instead of the natural support $x_{i}\in[-1,1]$ of a standard Jacobi ensemble.
This corresponds to a bipartite quantum system in the fermionic Gaussian
states of subsystem dimensions $m$ and $n$ with the assumption $m\leq n$. For
convenience, we denote

$a=n-m$ (12)

as the difference of subsystem dimensions.

We define the $k$-th order statistic over the fermionic Gaussian ensemble as
the average of random variables consisting of sum of all eigenvalues
$x_{i},i=1,\dots,m$ of the form

$\ln^{k}(1\pm x_{i}).$ (13)

In this work, we focus on the study of two major second-order statistics of
fermionic Gaussian states. The first one is the variance of the von Neumann
entanglement entropy (8). The von Neumann entropy is the most fundamental
entanglement measure that satisfies various interesting properties. It is an
important theoretic quantity in quantum information theory and has
applications in different fields, including the many-body systems and black
hole radiation [23]. Recently, the average von Neumann entropy over the
fermionic Gaussian ensemble (10) has been obtained as [19]

$\displaystyle\mathbb{E}\\!\left[S\right]$ $\displaystyle=$
$\displaystyle\left(m+n-\frac{1}{2}\right)\psi_{0}(2m+2n)+\left(\frac{1}{4}-m\right)\psi_{0}(m+n)+$
(14)
$\displaystyle\left(\frac{1}{2}-n\right)\psi_{0}(2n)-\frac{1}{4}\psi_{0}(n)-m,$

where

$\psi_{0}(x)=\frac{\,\mathrm{d}\ln\Gamma(x)}{\,\mathrm{d}x}$ (15)

is the digamma function and for a positive integer [27]

$\psi_{0}(l)=-\gamma+\sum_{k=1}^{l-1}\frac{1}{k}$ (16)

with $\gamma\approx 0.5772$ being the Euler’s constant. As one step further
towards understanding the statistical behaviour of von Neumann entropy, we
study the exact variance $\mathbb{V}\\!\left[S\right]$ over the fermionic
Gaussian states. Specifically, we are able to derive the exact variance for
equal subsystem dimensions $n=m$. We also conjecture the exact variance for
arbitrary subsystem dimensions $n\geq m$.

Another second-order statistic we are considering is the mean value of the
entanglement capacity. The capacity of entanglement over the fermionic
Gaussian ensemble (10) is defined as [24, 20]

$C=\sum_{i=1}^{m}\frac{1-x_{i}^{2}}{4}\ln^{2}\frac{1+x_{i}}{1-x_{i}}.$ (17)

The mean capacity is a second-order statistics because the above definition
can be rewritten as

$\displaystyle C$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\Bigg{(}\frac{1-x_{i}}{2}\ln^{2}\frac{1-x_{i}}{2}+\frac{1+x_{i}}{2}\ln^{2}\frac{1+x_{i}}{2}-$
(18)
$\displaystyle\left(\frac{1-x_{i}}{2}\ln\frac{1-x_{i}}{2}+\frac{1+x_{i}}{2}\ln\frac{1+x_{i}}{2}\right)^{2}\Bigg{)},$

which corresponds to our definition of $k$-th order statistic (13) when $k=2$.
As the quantum information theoretic counterpart of heat capacity, the
capacity of entanglement is known to capture the deviation from the maximal
entanglement [24, 25]. The mean of entanglement capacity under the Hilbert-
Schmidt ensemble has been recently computed in [26]. For the fermionic
Gaussian ensemble, we compute the exact average capacity when the subsystem
dimension differences are $a=0,1,2,3$.

In proposition 19 below, we summarize the exact variance of von Neumann
entanglement entropy of equal subsystem dimensions $n=m$.

######  Proposition 1.

For equal subsystem dimensions $n=m$, i.e., $a=0$, the exact variance of the
von Neumann entropy (8) of fermionic Gaussian states (10) is given by

$\displaystyle\mathbb{V}\\!\left[S\right]$ $\displaystyle=$
$\displaystyle\\!\left(\frac{1}{2}-2n\right)\psi_{1}(4n)+\frac{56n^{2}-36n+5}{8(4n-1)}\psi_{1}(2n)+\frac{\psi_{1}(n)}{8}-\frac{\psi_{0}(4n)}{2}+\frac{\psi_{0}(2n)}{2},$
(19)

where $\psi_{0}(x)$ is the digamma function (15) and
$\psi_{1}(x)=\,\mathrm{d}^{2}\ln\Gamma(x)/\,\mathrm{d}x^{2}$ is the trigamma
function.

The proof of proposition 19 can be found in section 2.2. Here, the trigamma
function also admits the finite sum form as [27]

$\psi_{1}(l)=\frac{\pi^{2}}{6}-\sum_{k=1}^{l-1}\frac{1}{k^{2}}.$ (20)

Our approach in deriving the result of proposition 19 in principle works for
any fixed dimension difference (i.e., $a=n-m=1,2,\dots$). Based on the further
calculations of the variance for the cases of different values of $a$, we come
up with the following conjecture on the exact variance of arbitrary subsystem
dimensions.

######  Conjecture.

For any subsystem dimensions $m\leq n$, the exact variance of the von Neumann
entanglement entropy (8) of fermionic Gaussian states (10) is given by

$\displaystyle\mathbb{V}\\!\left[S\right]$ $\displaystyle=$
$\displaystyle\\!\left(\frac{1}{2}-m-n\right)\psi_{1}(2m+2n)+\left(n-\frac{1}{2}\right)\psi_{1}(2n)+\left(\frac{m(2m+n-1)}{2m+2n-1}-\frac{1}{8}\right)\times$
(21)
$\displaystyle\psi_{1}(m+n)+\frac{\psi_{1}(n)}{8}-\frac{1}{2}(\psi_{0}(2m+2n)-\psi_{0}(2n)).$

Note that in the case of equal subsystem dimensions $n=m$, the above
conjectured result reduces to (19) in proposition 19, as expected. It is also
worth mentioning that the conjectured formula (21) recovers the obtained
asymptotic variance in [19], i.e.,

$\lim_{m,n\rightarrow\infty}\mathbb{V}\\!\left[S\right]=\frac{1}{2}\left(f+f^{2}+\ln(1-f)\right)$
(22)

for a fixed $f=m/(n+m)$. This asymptotic result (22) can be obtained by using
the asymptotic behavior of polygamma functions [27]

$\displaystyle\psi_{0}(x)$ $\displaystyle=$
$\displaystyle\Theta\left(\ln(x)-\frac{1}{2x}\right),\qquad x\to\infty$ (23)
$\displaystyle\psi_{j}(x)$ $\displaystyle=$
$\displaystyle\Theta\left(\frac{1}{x^{j}}\right),\qquad
x\to\infty,~{}~{}~{}~{}~{}~{}j\geq 1,$ (24)

where $\Theta(\cdot)$ is the Big-Theta notation of the Bachmann-Landau
symbols.

We now move on to the results of the entanglement capacity. Here, we also
focus on cases when the differences of the subsystem dimensions are fixed. The
exact formulas of mean capacity of subsystem dimension differences
$a=n-m=0,1,2,3$ are summarized in proposition 2 below, where the corresponding
calculation is provided in section 2.3.

######  Proposition 2.

For subsystem dimension differences $a=n-m=0,1,2,3$, the respective averages
of entanglement capacity (17) of fermionic Gaussian states (10) are

$\displaystyle\mathbb{E}\\!\left[C\right]$ $\displaystyle=$ $\displaystyle
a_{0}\psi_{1}(2n)+a_{1}\psi_{1}(n)+a_{2}$ (25)
$\displaystyle\mathbb{E}\\!\left[C\right]$ $\displaystyle=$ $\displaystyle
b_{0}\psi_{1}(2n)+b_{1}\psi_{1}(n)+b_{2}$ (26)
$\displaystyle\mathbb{E}\\!\left[C\right]$ $\displaystyle=$ $\displaystyle
c_{0}\psi_{1}(2n)+c_{1}\psi_{1}(n)+c_{2}(\psi_{0}(2n)-\psi_{0}(1))+c_{3}$ (27)
$\displaystyle\mathbb{E}\\!\left[C\right]$ $\displaystyle=$ $\displaystyle
d_{0}\psi_{1}(2n)+\\!d_{1}\psi_{1}(n)+\\!d_{2}(\psi_{0}(2n)-\psi_{0}(1))+d_{3},$
(28)

where the coefficients $a$, $b$, $c$, and $d$ are summarized in table 1.

Table 1: Coefficients of $\mathbb{E}\\!\left[C\right]$ of subsystem dimension
differences $a=0,1,2,3$

$a_{0}=$ |  $\displaystyle-\frac{(2n-1)^{2}}{2(4n-1)}$
---|---
$a_{1}=$ |  $\displaystyle-\frac{1}{8}$
$a_{2}=$ |  $\displaystyle\frac{\pi^{2}\left(8n^{2}-4n+1\right)}{16(4n-1)}+\frac{1}{2}(1-2n)$
$a=0$

$b_{0}=$ |  $\displaystyle-\frac{4n^{2}-8n+3}{2(4n-3)}$
---|---
$b_{1}=$ |  $\displaystyle-\frac{1}{8}$
$b_{2}=$ |  $\displaystyle\frac{\pi^{2}\left(8n^{2}-12n+3\right)}{16(4n-3)}-\frac{16n^{3}-36n^{2}+28n-9}{2(2n-1)(4n-3)}$
$a=1$

$c_{0}=$ |  $\displaystyle\frac{-4n^{2}+12n-5}{8n-10}$
---|---
$c_{1}=$ |  $\displaystyle-\frac{1}{8}$
$c_{2}=$ |  $\displaystyle\frac{1}{2n^{2}-5n+3}$
$c_{3}=$ |  $\displaystyle\frac{\pi^{2}\left(8n^{2}-20n+5\right)}{64n-80}-\frac{32n^{4}-152n^{3}+268n^{2}-210n+75}{2(2n-3)(2n-1)(4n-5)}$
$a=2$

$d_{0}=$ |  $\displaystyle-\frac{(2n-7)(2n-1)}{2(4n-7)}$
---|---
$d_{1}=$ |  $\displaystyle-\frac{1}{8}$
$d_{2}=$ |  $\displaystyle\frac{2\left(4n^{2}-14n+11\right)}{(n-2)(n-1)(2n-5)(2n-3)}$
$d_{3}=$ |  $\displaystyle\frac{\pi^{2}\left(8n^{2}-28n+7\right)}{16(4n-7)}-\frac{1}{2(n-2)(n-1)(2n-5)(2n-3)(2n-1)(4n-7)}\times$
|
$\displaystyle(64n^{7}-720n^{6}+3408n^{5}-8736n^{4}+13176n^{3}-11967n^{2}+6258n-1470)$
$a=3$

From the results of proposition 2, we observe that the coefficients of
$\psi_{1}(n)$ remain the same for $a=0,1,2,3$. We also observe that the term
involving digamma functions, i.e., $(\psi_{0}(2n)-\psi_{0}(1))$, starts to
appear when $a\geqslant 2$. Our approach in deriving the above formulas in
principle works for any given subsystem dimension difference $a=n-m$. Note
that in a recent work [20], an exact representation of mean capacity for
arbitrary $a$ has been obtained in terms of a rather complicated triple
summations.

As will be seen in section 2, the computation of the exact second-order
statistics of the fermionic Gaussian states requires one to calculate the
corresponding integrals involving Jacobi polynomials. In principle, one could
obtain different representations of the same integral due to the availability
of different summation representations such as the ones of Jacobi polynomial
(35)–(37). These different choices lead to different summation forms, where
most of the combinations may not permit further simplifications. In practice,
the integral identities utilized in establishing the claimed results have
underwent a sensible choice so as to facilitate the cancellation among the
summations. The cancellation eventually leads to the desired closed-form
results as shown in proposition 19 and proposition 2.

Finally, we point out the fact that our obtained summation representations
(77chcm)–(77chcq) may not seem to permit further simplification that would
lead to the conjectured variance (21) for an arbitrary $a$. This is because
the gamma functions involved in these summation representations, e.g. in
(77chcn) and (77chcp), will not cancel in pairs. In contrast, for a given $a$,
these gamma functions are reduced to summable rational terms leading to the
claimed results.

Figure 1: Variance of von Neumann entropy: analytical results versus
simulations. The solid black line represents the obtained analytical result
(19) of equal subsystem dimensions $n=m$. The solid blue line ($n=2m$) and the
solid red line ($n=3m$) are drawn by the conjectured result (21). The scatters
(diamond, triangular, and square shapes) represent numerical simulations.

To illustrate the derived result (19) and provide numerical evidence of the
conjecture (21), we plot in figure 1111The simulations performed in figure 1–3
make use of the Mathematica codes provided by Santosh Kumar based on the log-
gas approach as discussed in[15, appendix B]. the exact variance of von
Neumann entropy as compared with the simulations. We consider the scenarios of
equal and unequal subsystem dimensions, where the claimed results in solid
lines match well with the corresponding simulations as represented by the
scatters. We also observe that the slope of the curves approaches to zero when
the dimension of the system becomes larger. This phenomenon can be explained
by the asymptotic formula (22) of the variance, i.e., the variance is
dominated by the ratios of subsystem dimensions for large-dimensional systems.

We now consider the standardized von Neumann entropy

$X=\frac{S-\mathbb{E}\\!\left[S\right]}{\sqrt{\mathbb{V}\\!\left[S\right]}},$
(29)

where its average $\mathbb{E}\\!\left[S\right]$ and variance
$\mathbb{V}\\!\left[S\right]$ are given by (14) and (21), respectively.

Figure 2: Probability densities of standardized von Neumann entropy: a
comparison of Gaussian approximation (30) to the simulation results. The
dashed line in blue and the dash-dot line in red refer to the standardized von
Neumann entropy (29) of subsystem dimensions $m=2$, $n=4$, and $m=16$, $n=32$,
respectively. The solid black line represents the Gaussian approximation (30).

In figure 2, we plot the standardized von Neumann entropy $X$ for the cases of
a fixed subsystem dimension ratio $m/n=1/2$ along with the standard Gaussian
density

$\frac{1}{\sqrt{2\pi}}\rme^{-\frac{1}{2}x^{2}},\qquad
x\in\left(-\infty,\infty\right).$ (30)

We observe from this plot that the true distribution of the standardized von
Neumann entropy $X$ is non-symmetric and appears to be left-skewed when $m=2$,
$n=4$. In comparison, when the subsystem dimensions increase to $m=16$, $n=32$
with the same ratio $m/n=1/2$, the distribution of $X$ appears to be closer to
the Gaussian distribution. The observed asymptotic Gaussian behavior is
typical for a wide class of linear spectral statistics (such as the von
Neumann entropy) over different random matrix ensembles as already been
observed in the Hilbert-Schmidt ensemble [10] and Bures-Hall ensemble [17].
Here, we also conjecture that in the limit

$m\to\infty,~{}~{}~{}~{}n\to\infty,~{}~{}~{}~{}\frac{m}{n}\in(0,1],$ (31)

the standardized von Neumann entropy (29) converges in distribution to a
standard Gaussian random variable.

Figure 3: Average entanglement capacity: analytical results versus
simulations. The solid black lines represent the analytical results (25)–(28)
and the diamond scatters represent the corresponding simulations. The
parameter $a$ in the plots denotes the dimension difference $a=n-m$.

For the entanglement capacity, we numerically compare our mean results
(25)–(28) with the simulations in figure 3. As can be seen, the claimed
results match well with the simulations. In addition, we notice the linear
growth of the mean capacity, which can be captured by the following asymptotic
results. For example, when $a=n-m=0$, by taking the appropriate limit of the
result (25) and using the asymptotic behavior (24) of the trigamma function,
one obtains the slope of linear behavior

$\lim_{n\rightarrow\infty}\frac{\mathbb{E}\\!\left[C\right]}{n}=\frac{1}{8}\left(\pi^{2}-8\right).$
(32)

The rest of the paper is organized as follows. In section 2, we perform the
derivation of the considered second-order statistics of fermionic Gaussian
states. Specifically, in section 2.1 we first provide necessary results of
Jacobi polynomials together with two integral identities useful in the
subsequent computations. In section 2.2, we compute the exact variance of von
Neumann entropy based on the computation and simplification of the four
integrals involved. In section 2.3, the exact mean capacity is derived by
calculating one additional integral along with the results already obtained.
The summation forms of the five integrals are listed in appendix A, and the
relevant finite sum identities utilized in simplifying the summations are
listed in appendix B.

##  2 Computation of second-order statistics

In this section, we first introduce some useful results of Jacobi polynomials
in section 2.1. We then present the computation of the variance of von Neumann
entropy $\mathbb{V}\\!\left[S\right]$ in section 2.2 and that of the mean
capacity $\mathbb{E}\\!\left[C\right]$ in section 2.3.

###  2.1 Integrals involving Jacobi polynomials

We list here the necessary random matrix results and integral identities
regarding the Jacobi unitary ensemble. While most of the listed results are
known in the literature [28, 29, 30], some of which seem new.

For the considered Jacobi ensemble (10), its joint density
$g_{l}(x_{1},\dots,x_{l})$ of $l$ (out of $m$) eigenvalues can be written in
terms of an $l\times l$ determinant as [19]

$g_{l}(x_{1},\dots,x_{l})=\frac{(m-l)!}{m!}\det\left(K\left(x_{i},x_{j}\right)\right)_{i,j=1}^{l}.$
(33)

The corresponding correlation kernel is given by

$K\left(x,y\right)=\sqrt{(1-x^{2})^{a}(1-y^{2})^{a}}\sum_{k=0}^{m-1}\frac{p_{k}(x)p_{k}(y)}{h_{k}},$
(34)

where $a$ has been defined in (12). The polynomial $p_{k}(x)$ of degree $2k$
is related to different representations of the standard Jacobi polynomial [19,
28]

$\displaystyle J^{(a,b)}_{k}(x)$ $\displaystyle=$
$\displaystyle\frac{(-1)^{k}(b+1)_{k}}{k!}\sum_{i=0}^{k}\frac{(-k)_{i}(k+a+b+1)_{i}}{(b+1)_{i}\Gamma(i+1)}\left(\frac{1+x}{2}\right)^{i}$
(35) $\displaystyle=$
$\displaystyle\frac{1}{k!}\sum_{i=0}^{k}\frac{(-k)_{i}(k+a+b+1)_{i}(i+a+1)_{k-i}}{\Gamma(i+1)}\left(\frac{1-x}{2}\right)^{i}$
(36) $\displaystyle=$
$\displaystyle\sum_{i=0}^{k}\frac{(-1)^{i}\Gamma(a+k+1)(k+b-i+1)_{i}}{\Gamma(i+1)\Gamma(a+i+1)\Gamma(k-i+1)}\left(\frac{1-x}{2}\right)^{i}\left(\frac{1+x}{2}\right)^{k-i}$
(37)

supported in $x\in[-1,1]$, as

$p_{k}(x)=J^{(n-m,n-m)}_{2k}(x)$ (38)

with

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

denoting the Pochhammer’s symbol. The orthogonality relation of the standard
Jacobi polynomials

$\displaystyle\int_{-1}^{1}(1-x)^{a}(1+x)^{b}J^{(a,b)}_{k}(x)J^{(a,b)}_{l}(x)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\frac{2^{a+b+1}\Gamma(k+a+1)\Gamma(k+b+1)}{k!(2k+a+b+1)\Gamma(k+a+b+1)}\delta_{kl},\quad\Re\\{a,b\\}>-1,$
(40)

leads to the normalization constant $h_{k}$ of the desired polynomials
$p_{k}(x)$,

$\int_{0}^{1}\left(1-x^{2}\right)^{a}p_{k}(x)p_{l}(x)\,\mathrm{d}x=h_{k}\delta_{kl}$
(41)

as

$h_{k}=\frac{2^{2a}\Gamma^{2}(2k+a+1)}{(4k+2a+1)\Gamma(2k+2a+1)\Gamma(2k+1)}.$
(42)

Note that the polynomials (38) are even

$p_{k}(-x)=p_{k}(x),$ (43)

which is a consequence of the parity property of Jacobi polynomials

$J^{(a,b)}_{k}(-x)=(-1)^{k}J^{(b,a)}_{k}(x).$ (44)

We now present two integral identities useful in computing the considered
second-order statistics. The first one is

$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{a}\left(\frac{1+x}{2}\right)^{c}J^{(a,b)}_{k}(x)\,\mathrm{d}x$
(45) $\displaystyle=$
$\displaystyle\frac{2\Gamma(c+1)\Gamma(k+a+1)\Gamma(c-b+1)}{k!\Gamma(k+a+c+2)\Gamma(c-k-b+1)},\quad\Re\\{a,b,c\\}>-1.$
(46)

This identity is obtained by using the definition (35) and the identity

$\int_{-1}^{1}(1-x)^{a}(1+x)^{b}\,\mathrm{d}x=\frac{2^{a+b+1}\Gamma(a+1)\Gamma(b+1)}{\Gamma(a+b+2)},\quad\Re\\{a,b\\}>-1,$
(47)

where the resulting summation can be written as a unit argument hypergeometric
function of Saalschützian type that terminates

${}_{3}F_{2}(-n,a,b;d,a+b-n-d+1;1)=\frac{(d-a)_{n}(d-b)_{n}}{(d)_{n}(-a-b+d)_{n}},\quad
n\in\mathbb{Z^{+}}.$ (48)

The next identity is

$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{d}\left(\frac{1+x}{2}\right)^{c}J^{(a,b)}_{k}(x)\,\mathrm{d}x$
(50) $\displaystyle=$
$\displaystyle\frac{2\Gamma(c-b+1)\Gamma(d-a+1)}{\Gamma(c+d+k+2)}\sum_{i=0}^{k}\frac{(-1)^{i}\Gamma(c+i+1)\Gamma(d-i+k+1)}{\Gamma(i+1)\Gamma(k-i+1)}\times$
$\displaystyle\frac{1}{\Gamma(d-a-i+1)\Gamma(c-b+i-k+1)},\quad\Re\\{a,b,c,d\\}>-1,$

where the integral involves an additional variable $d$ as compared to the
integral (45). Proving the above identity requires the Rodrigues’ formula[30],
i.e., for $\Re\\{a,b\\}>-1$, one has

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!(1-x)^{a}(1+x)^{b}J^{(a,b)}_{k}(x)=\frac{(-1)^{k}}{2^{k}k!}\frac{\,\mathrm{d}^{k}}{\,\mathrm{d}x^{k}}\left((1-x)^{k+a}(1+x)^{k+b}\right).$
(51)

By using the above formula, the integral (50) is evaluated as

$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{d}\left(\frac{1+x}{2}\right)^{c}J^{(a,b)}_{k}(x)\,\mathrm{d}x$
(52) $\displaystyle=$
$\displaystyle\frac{(-1)^{k}}{2^{k+c+d}k!}\int_{-1}^{1}(1-x)^{d-a}(1+x)^{c-b}\frac{\,\mathrm{d}^{k}}{\,\mathrm{d}x^{k}}\left((1-x)^{k+a}(1+x)^{k+b}\right)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\frac{1}{2^{k+c+d}k!}\int_{-1}^{1}(1-x)^{k+a}(1+x)^{k+b}\frac{\,\mathrm{d}^{k}}{\,\mathrm{d}x^{k}}\left((1-x)^{d-a}(1+x)^{c-b}\right)\,\mathrm{d}x$
(55) $\displaystyle=$
$\displaystyle\frac{1}{2^{k+c+d}k!}\sum_{i=0}^{k}{{k}\choose{i}}(-1)^{i}(d-a-i+1)_{i}(c-b+i-k+1)_{k-i}\times$
$\displaystyle\int_{-1}^{1}(1-x)^{d-i+k}(1+x)^{c+i}\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle\frac{2\Gamma(c-b+1)\Gamma(d-a+1)}{\Gamma(c+d+k+2)}\sum_{i=0}^{k}\frac{(-1)^{i}\Gamma(c+i+1)\Gamma(d-i+k+1)}{\Gamma(i+1)\Gamma(k-i+1)}\times$
$\displaystyle\frac{1}{\Gamma(d-a-i+1)\Gamma(c-b+i-k+1)},\quad\Re\\{a,b,c,d\\}>-1,$

where the equality (55) is obtained by applying $k$ times integration by parts
and the last equality (55) is obtained by the identity (47). Finally, we note
that besides the results (46) and (50), the integrals (45) and (50) also have
alternative forms obtained by using the different representations of the
Jacobi polynomials such as (35)–(37). Among the alternative formulas, no other
combinations would lead to the desired results in proposition 19 and
proposition 2.

###  2.2 Computation of von Neumann entropy variance

In this section, we compute the exact variance of von Neumann entropy shown in
proposition 19. The eigenvalue densities required to compute the first two
moments of von Neumann entropy can be read off from the $l$-point correlation
function (33) as

$\displaystyle g_{1}(x_{1})$ $\displaystyle=$
$\displaystyle\frac{1}{m}K\left(x_{1},x_{1}\right)$ (56)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!g_{2}(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{m(m-1)}\left(K(x_{1},x_{1})K(x_{2},x_{2})-K^{2}(x_{1},x_{2})\right).$
(57)

We now first reproduce the mean formula (14) of von Neumann entropy in order
to set up our notations while providing a more detailed derivation.
Calculating the mean value requires the one-point density (56) as

$\displaystyle\mathbb{E}\\!\left[S\right]$ $\displaystyle=$
$\displaystyle-\int_{0}^{1}v(x)K(x,x)\,\mathrm{d}x.$ (58)

Inserting the definitions (9) and (34) into (58), one has

$\mathbb{E}\\!\left[S\right]=-\sum_{k=0}^{m-1}\frac{1}{h_{k}}\int_{-1}^{1}\frac{1+x}{2}\ln\frac{1+x}{2}\left(1-x^{2}\right)^{a}p_{k}^{2}(x)\,\mathrm{d}x,$
(59)

where $a=n-m$ and it becomes clear that the key is to evaluate the integral

$\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{a}\left(\frac{1+x}{2}\right)^{c}p_{k}^{2}(x)\,\mathrm{d}x.$
(60)

By using the identity (46) along with the definitions (35) and (38), the above
integral is evaluated to

$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{a}\left(\frac{1+x}{2}\right)^{c}p_{k}^{2}(x)\,\mathrm{d}x$
(61) $\displaystyle=$
$\displaystyle\frac{2\Gamma^{2}(2k+a+1)}{\Gamma(2k+1)\Gamma(2k+2a+1)}\sum_{i=0}^{2k}\frac{(-1)^{i}\Gamma(i+2k+2a+1)}{i!\Gamma(i+a+1)\Gamma(2k-i+1)}\times$
$\displaystyle\frac{\Gamma(i+c+1)\Gamma(i-a+c+1)}{\Gamma(i-2k-a+c+1)\Gamma(i+2k+a+c+2)}.$

Now taking the derivative of (61) before setting $c=a+1$, the resulting sum
can be evaluated that leads to

$\displaystyle\int_{-1}^{1}\frac{1+x}{2}\ln\frac{1+x}{2}\left(1-x^{2}\right)^{a}p_{k}^{2}(x)\,\mathrm{d}x$
(62) $\displaystyle=$ $\displaystyle
h_{k}\big{(}1+\psi_{0}(2k+a)+\psi_{0}(2k+2a)-2\psi_{0}(4k+2a)\big{)}+$
$\displaystyle\frac{h_{k}}{2}\left(\frac{1}{k+a}-\frac{a}{2k+a}-\frac{a}{2k+a+1}-\frac{2}{4k+2a+1}\right),$

where $h_{k}$ is given by (42). Inserting the above result into (59), the
remaining summation over $k$ is evaluated by repeated use of the identity
(77chcr) as well as the result [31]

$\displaystyle\psi_{0}(mk)=\ln
m+\frac{1}{m}\sum_{i=0}^{m-1}\psi_{0}\left(k+\frac{i}{m}\right),\qquad
m\in\mathbb{Z^{+}}.$ (63)

This completes the proof of the mean entropy formula (14).

We now turn to the variance. By definition, the variance calculation boils
down to computing two integrals, cf. [8, 17, 19],

$\displaystyle\mathbb{V}\\!\left[S\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[S^{2}\right]-\mathbb{E}^{2}\\!\left[S\right]=\mathrm{I_{A}}-\mathrm{I_{B}},$
(64)

where

$\displaystyle\mathrm{I_{A}}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}v^{2}(x)K(x,x)\,\mathrm{d}x$ (65)
$\displaystyle\mathrm{I_{B}}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}v(x)v(y)K^{2}\left(x,y\right)\,\mathrm{d}x\,\mathrm{d}y.$
(66)

According to the definitions (9) and (34), the $\mathrm{I_{A}}$ integral (65)
boils down computing the two parts

$\mathrm{I_{A}}=\mathrm{A_{1}}+\mathrm{A_{2}},$ (67)

where

$\displaystyle\mathrm{A_{1}}$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{1}{h_{k}}\int_{-1}^{1}\left(\frac{1+x}{2}\right)^{2}\ln^{2}\frac{1+x}{2}\left(1-x^{2}\right)^{a}p_{k}^{2}(x)\,\mathrm{d}x$
(68) $\displaystyle\mathrm{A_{2}}$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{1}{h_{k}}\int_{-1}^{1}\frac{1-x}{2}\frac{1+x}{2}\ln\frac{1-x}{2}\ln\frac{1+x}{2}\left(1-x^{2}\right)^{a}p_{k}^{2}(x)\,\mathrm{d}x.$
(69)

Similarly, the $\mathrm{I_{B}}$ integral (66) can be written in terms of the
following two integrals

$\mathrm{I_{B}}=\mathrm{B_{1}}+\mathrm{B_{2}},$ (70)

where

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathrm{B_{1}}$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\frac{1}{h_{k}^{2}}\left(\int_{-1}^{1}\frac{1+x}{2}\ln\frac{1+x}{2}\left(1-x^{2}\right)^{a}p_{k}^{2}(x)\,\mathrm{d}x\right)^{2}$
(71) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathrm{B_{2}}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{m-1}\sum_{k=0}^{m-j-1}\frac{2}{h_{k+j}h_{k}}\Bigg{(}\int_{-1}^{1}\frac{1+x}{2}\ln\frac{1+x}{2}\left(1-x^{2}\right)^{a}p_{k+j}(x)p_{k}(x)\,\mathrm{d}x\Bigg{)}^{2}.$
(72)

The next step is to compute the integrals in $\mathrm{A_{1}}$,
$\mathrm{A_{2}}$, $\mathrm{B_{1}}$, and $\mathrm{B_{2}}$. These integrals can
be calculated into summations by using the corresponding integral identities
of Jacobi polynomials introduced in section 2.1. The strategies in computing
the integrals will be discussed in the following. The integral in
$\mathrm{A_{1}}$ is calculated by taking twice derivatives of (61) with
respect to $c$ before setting $c=a+2$. The calculation of the integral in
$\mathrm{A_{2}}$ requires taking derivatives with respect to both $c$ and $d$
of the identity

$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{d}\left(\frac{1+x}{2}\right)^{c}p_{k}^{2}(x)\,\mathrm{d}x$
(73) $\displaystyle=$
$\displaystyle\frac{2\Gamma(2k+a+1)^{2}}{\Gamma(4k+c+d+2)}\sum_{j=0}^{2k}\sum_{i=0}^{2k}\frac{(-1)^{i+j}\Gamma(2k-j+c-a+1)}{\Gamma(i+1)\Gamma(j+1)\Gamma(2k-j+a+1)}\times$
$\displaystyle\frac{\Gamma(j-a+d+1)\Gamma(2k+i-j+c+1)\Gamma(2k+d-i+j+1)}{\Gamma(j+a+1)\Gamma(2k-i+1)\Gamma(2k-j+1)\Gamma(i-j-a+c+1)}\times$
$\displaystyle\frac{1}{\Gamma(j-i-a+d+1)}$

before setting $c=d=a+1$. The above identity is established by using the
identity (50) along with the definitions (37) and (38). The integral in
$\mathrm{B_{1}}$ is given by (62), which has been obtained in the mean
calculation. The integral in $\mathrm{B_{2}}$ is computed by taking the
derivative of the identity

$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{a}\left(\frac{1+x}{2}\right)^{c}p_{k+j}(x)p_{k}(x)\,\mathrm{d}x$
(74) $\displaystyle=$
$\displaystyle\frac{2\Gamma(2k+a+1)\Gamma(2k+2j+a+1)}{\Gamma(2k+2j+1)\Gamma(2k+2a+1)}\sum_{i=0}^{2k}\frac{(-1)^{i}\Gamma(i+2k+2a+1)}{\Gamma(i+a+1)\Gamma(2k-i+1)}\times$
$\displaystyle\frac{\Gamma(i+c+1)\Gamma(i-a+c+1)}{\Gamma(i+1)\Gamma(i-2k-2j-a+c+1)\Gamma(i+2k+2j+a+c+2)}$

with respect to $c$ and setting $c=a+1$. To derive the above identity, we
first apply the definitions (35) and (36) on $p_{k+j}(x)$ and $p_{k}(x)$,
respectively. Keeping in mind the relation (38), we have

$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{a}\left(\frac{1+x}{2}\right)^{c}p_{k+j}(x)p_{k}(x)\,\mathrm{d}x$
(75) $\displaystyle=$
$\displaystyle\frac{\Gamma(a+2k+1)\Gamma(a+2j+2k+1)}{\Gamma(2a+2k+1)\Gamma(2a+2j+2k+1)}\sum_{s=0}^{2j+2k}\frac{(-1)^{s}}{\Gamma(s+1)\Gamma(a+s+1)}\times$
$\displaystyle\frac{\Gamma(2a+2j+2k+s+1)}{\Gamma(2j+2k-s+1)}\sum_{i=0}^{2k}\frac{(-1)^{i}\Gamma(2a+i+2k+1)}{\Gamma(i+1)\Gamma(a+i+1)\Gamma(2k-i+1)}\times$
$\displaystyle\int_{-1}^{1}\left(\frac{1-x}{2}\right)^{a+s}\left(\frac{1+x}{2}\right)^{c+i}\,\mathrm{d}x.$

The integral in (75) is further evaluated by the identity (47), where one of
the resulting summations can be written as a unit argument ${}_{2}F_{1}$
hypergeometric function

${}_{2}F_{1}(-2j-2k,2a+2j+2k+1;a+c+i+2;1)$ (76) $\displaystyle=$
$\displaystyle\frac{\Gamma(-a+c+i+1)\Gamma(a+c+i+2)}{\Gamma(-a+c+i-2j-2k+1)\Gamma(a+c+i+2j+2k+2)}$

leading to the claimed identity (74). In writing down the summation forms of
$\mathrm{A_{1}}$, $\mathrm{A_{2}}$, $\mathrm{B_{1}}$, and $\mathrm{B_{2}}$,
one will have to solve the indeterminacy by using the following formulas for
$\epsilon\rightarrow 0$,

$\displaystyle\Gamma(-l+\epsilon)$ $\displaystyle\\!=$
$\displaystyle\frac{(-1)^{l}}{l!\epsilon}\left(1+\psi_{0}(l+1)\epsilon+o\left(\epsilon^{2}\right)\right)$
(77a) $\displaystyle\psi_{0}(-l+\epsilon)$ $\displaystyle=$
$\displaystyle-\frac{1}{\epsilon}+\psi_{0}(l+1)+\left(2\psi_{1}(1)-\psi_{1}(l+1)\right)\epsilon+o\left(\epsilon^{2}\right)$
(77b) $\displaystyle\psi_{1}(-l+\epsilon)$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon^{2}}-\psi_{1}(l+1)+\psi_{1}(1)+\zeta(2)+o\left(\epsilon\right)$
(77c)

with $l$ being a non-positive integer. The resulting summations forms of
$\mathrm{A_{1}}$, $\mathrm{A_{2}}$, $\mathrm{B_{1}}$, and $\mathrm{B_{2}}$ are
summarized in appendix A. The rest of the task to show proposition 19 is to
simplify these summations.

The formulas in appendix A involve summations over gamma, digamma, and
trigamma functions. The identities utilized in evaluating these sums are
listed in appendix B, where the closed-form identities are (77chcr)–(77chcv)
and the semi closed-form ones are (77chcw)–(B). These semi closed-form
identities represent the relation between two single sums consisting of
rational functions and polygamma functions. By using these semi closed-form
identities, one is able to convert different single sums into the ones that we
refer to as unsimplifiable basis such as the sums (77cha)–(77chc). Note that
most of the summations in appendix A do not admit a closed-form evaluation
individually. To facilitate the cancellations among the summations, one has to
potentially first simplify these sums into semi closed-form ones via the
unsimplifiable basis. The closed-form results in this work become available by
observing the cancellations among the semi closed-form terms.

The simplification process of the summation forms (77chcm)–(77chcp) in
appendix A is outlined below. For an arbitrary $a$, the summation forms
(77chcm) and (77chco) only involve rational functions. Thus, these sums admit
semi closed-form expressions by using the identities listed in appendix B
along with the results (63) and

$\displaystyle\psi_{1}(mk)=\frac{1}{m^{2}}\sum_{i=0}^{m-1}\psi_{1}\left(\frac{i}{m}+k\right),\qquad
m\in\mathbb{Z^{+}}.$ (77bz)

The summations (77chcn) and (77chcp) may not be represented via the semi
closed-form terms. This fact naturally guides us to first consider the special
case of equal subsystem dimensions $a=n-m=0$, where the summations can be
presented in terms of semi-closed expressions. Specifically, the summation
representation (77chcn) is simplified by first using the identity

$\displaystyle\sum_{k=0}^{m}\frac{1}{\Gamma(k+1)\Gamma(c+k)\Gamma(m+1-k)\Gamma(m+b+1-k)}$
$\displaystyle=$
$\displaystyle\frac{\,{}_{2}F_{1}(-b-m,-m;c;1)}{\Gamma(c)\Gamma(m+1)\Gamma(b+m+1)}$
$\displaystyle=$
$\displaystyle\frac{\Gamma(b+c+2m)}{\Gamma(m+1)\Gamma(b+m+1)\Gamma(c+m)\Gamma(b+c+m)},\quad\Re(2m+c+b)>0$
(77cb)

and identities obtained by up to its second derivatives of $c$ and $b$,

$\displaystyle\sum_{k=0}^{m}\frac{\psi_{0}(m+b+1-k)}{\Gamma(k+1)\Gamma(c+k)\Gamma(m+1-k)\Gamma(m+b+1-k)}$
(77cc) $\displaystyle=$
$\displaystyle\frac{\Gamma(b+c+2m)(\psi_{0}(b+c+m)-\psi_{0}(b+c+2m)+\psi_{0}(b+m+1))}{\Gamma(m+1)\Gamma(b+m+1)\Gamma(c+m)\Gamma(b+c+m)}$
$\displaystyle\sum_{k=0}^{m}\frac{\psi_{0}(c+k)}{\Gamma(k+1)\Gamma(c+k)\Gamma(m+1-k)\Gamma(m+b+1-k)}$
(77cd) $\displaystyle=$
$\displaystyle\frac{\Gamma(b+c+2m)(\psi_{0}(b+c+m)-\psi_{0}(b+c+2m)+\psi_{0}(c+m))}{\Gamma(m+1)\Gamma(b+m+1)\Gamma(c+m)\Gamma(b+c+m)}$
$\displaystyle\sum_{k=0}^{m}\frac{\psi_{0}(c+k)\psi_{0}(m+b+1-k)}{\Gamma(k+1)\Gamma(c+k)\Gamma(m+1-k)\Gamma(m+b+1-k)}$
(77ce) $\displaystyle=$
$\displaystyle\frac{\Gamma(b+c+2m)}{\Gamma(m+1)\Gamma(b+m+1)\Gamma(c+m)\Gamma(b+c+m)}((\psi_{0}(b+c+m)-\psi_{0}(b+c+2m)+$
$\displaystyle\psi_{0}(b+m+1))(\psi_{0}(b+c+m)-\psi_{0}(b+c+2m)+\psi_{0}(c+m))-\psi_{1}(b+c+m)+$
$\displaystyle\psi_{1}(b+c+2m)).$

All the gamma functions cancel in pairs in the resulting summations, where the
remaining terms can be evaluated by the identities in appendix B that leads to
a semi closed-form representation. For the summation (77chcp), it directly
reduces to a sum of rational functions when $a=0$. This also leads to a semi
closed-form representation of (77chcp). With the above results, the integrals
$\mathrm{I_{A}}$ and $\mathrm{I_{B}}$ can now be obtained. For
$\mathrm{I_{A}}$, one has

$\displaystyle\mathrm{I_{A}}=$
$\displaystyle\sum_{k=1}^{n}\left(\frac{\psi_{0}(2k)}{k}-\frac{\psi_{0}(4k)}{2k}+\frac{\psi_{0}(4k)}{2k+1}\right)-\frac{24n^{2}-12n+1}{4(4n-1)}\psi_{1}(2n)-\frac{1}{4}\psi_{1}(n)+$
(77cf)
$\displaystyle(4n-1)\psi_{0}^{2}(4n)+2(1-4n)\psi_{0}(2n)\psi_{0}(4n)+\frac{1}{2}(8n-3)\psi_{0}^{2}(2n)+\psi_{0}(n)\psi_{0}(2n)-$
$\displaystyle\frac{1}{2}\psi_{0}^{2}(n)-\frac{16n^{3}+8n^{2}-1}{2n(2n+1)}\psi_{0}(4n)+\frac{1}{2}(8n+1)\psi_{0}(2n)-\left(\frac{1}{2n}+\ln
2\right)\psi_{0}(n)-$
$\displaystyle\frac{1}{2}\psi_{0}(n)\psi_{0}\\!\left(n+\frac{1}{2}\right)-\frac{2n+1}{2n}\psi_{0}\\!\left(n+\frac{1}{2}\right)+\frac{1}{2}\psi_{0}\\!\left(n+\frac{1}{4}\right)+\frac{n(5n-2)}{4n-1}\psi_{1}(1)+$
$\displaystyle\\!\left(\frac{1}{2}+\ln
2\right)\psi_{0}(1)+\frac{1}{2}\psi_{0}\\!\left(\frac{1}{2}\right)\psi_{0}(1)+\psi_{0}\\!\left(\frac{1}{2}\right)-\frac{1}{2}\psi_{0}\\!\left(\frac{1}{4}\right)-\frac{1}{n}\ln
2-n+2.$

Similarly, for $\mathrm{I_{B}}$, one arrives at

$\displaystyle\mathrm{I_{B}}=$
$\displaystyle\sum_{k=1}^{n}\left(\frac{\psi_{0}(2k)}{k}-\frac{\psi_{0}(4k)}{2k}+\frac{\psi_{0}(4k)}{2k+1}\right)+\frac{4n-1}{2}\psi_{1}(4n)-\frac{104n^{2}-60n+7}{8(4n-1)}\psi_{1}(2n)-$
(77cg)
$\displaystyle\frac{3}{8}\psi_{1}(n)+(4n-1)\psi_{0}^{2}(4n)-2(4n-1)\psi_{0}(2n)\psi_{0}(4n)+\frac{1}{2}(8n-3)\psi_{0}^{2}(2n)+\psi_{0}(n)\times$
$\displaystyle\psi_{0}(2n)-\frac{1}{2}\psi_{0}^{2}(n)+\frac{-16n^{3}-6n^{2}+n+1}{2n(2n+1)}\psi_{0}(4n)+4n\psi_{0}(2n)-\left(\frac{1}{2n}+\ln
2\right)\times$
$\displaystyle\psi_{0}(n)-\frac{1}{2}\psi_{0}(n)\psi_{0}\\!\left(n+\frac{1}{2}\right)-\frac{(2n+1)}{2n}\psi_{0}\\!\left(n+\frac{1}{2}\right)+\frac{1}{2}\psi_{0}\\!\left(n+\frac{1}{4}\right)+\frac{(5n-2)}{4n-1}\times$
$\displaystyle n\psi_{1}(1)+\left(\frac{1}{2}+\ln
2\right)\psi_{0}(1)+\frac{1}{2}\psi_{0}\\!\left(\frac{1}{2}\right)\psi_{0}(1)-\frac{1}{2}\psi_{0}\\!\left(\frac{1}{4}\right)+\psi_{0}\\!\left(\frac{1}{2}\right)-$
$\displaystyle\frac{1}{n}\ln 2-n+2.$

Now inserting the $\mathrm{I_{A}}$ expression (77cf) and the above
$\mathrm{I_{B}}$ expression (77cg) into (64), we observe substantial
cancellations among the terms of $\mathrm{I_{A}}-\mathrm{I_{B}}$. In
particular, the terms involving the three types of sums

$\displaystyle\sum_{k=1}^{n}\frac{\psi_{0}(2k)}{k}$ (77cha)
$\displaystyle\sum_{k=1}^{n}\frac{\psi_{0}(4k)}{k}$ (77chb)
$\displaystyle\sum_{k=1}^{n}\frac{\psi_{0}(4k)}{2k+1}$ (77chc)

cancel completely. The surviving terms give us the desired identity (19). This
completes the proof of proposition 19. By using the same approach as for the
equal subsystem dimension case $a=0$, we also compute the variance for cases
being $a=1,2,3$. Note that cancellation phenomenon also appears in other
statistic of entanglement measures, see for example [8, 17, 10, 11].

###  2.3 Computation of mean capacity

In this section, we compute the average capacity of entanglement in fermionic
Gaussian states. By definitions (18) and (33), computing the average capacity
requires only the one-point density (56) as

$\displaystyle\mathbb{E}\\!\left[C\right]=\mathrm{I_{C}}-\mathrm{I_{A}},$
(77chci)

where

$\displaystyle\mathrm{I_{C}}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\left(\frac{1+x}{2}\ln^{2}\left(\frac{1+x}{2}\right)+\frac{1-x}{2}\ln^{2}\left(\frac{1-x}{2}\right)\right)K(x,x)\,\mathrm{d}x,$
(77chcj)

and $\mathrm{I_{A}}$ has been defined in (67). Using the correlation kernel
(34) of one-point density, the integral $\mathrm{I_{C}}$ can be written as

$\displaystyle\mathrm{I_{C}}=\sum_{k=0}^{m-1}\frac{1}{h_{k}}\int_{-1}^{1}\frac{1+x}{2}\ln^{2}\left(\frac{1+x}{2}\right)\left(1-x^{2}\right)^{a}p_{k}^{2}(x)\,\mathrm{d}x.$
(77chck)

The above integral can be evaluated by taking twice derivatives of the
identity (61) with respect to $c$ before setting $c=a+1$. The resulting
summation of $\mathrm{I_{C}}$ is given by (77chcq) in appendix A. Similar to
$\mathrm{A_{1}}$, $\mathrm{I_{C}}$ also admits a semi-closed expression for
any $a$. To proceed further, we have to consider the case when $a=0$ due to
the intractability of $\mathrm{I_{A}}$ as discussed previously. Indeed, when
$a=0$, we immediately have

$\displaystyle\mathrm{I_{C}}=$
$\displaystyle\sum_{k=1}^{n}\left(\frac{\psi_{0}(2k)}{k}-\frac{\psi_{0}(4k)}{2k}+\frac{\psi_{0}(4k)}{2k+1}\right)-\frac{1}{4}(8n-3)\psi_{1}(2n)-\frac{3}{8}\psi_{1}(n)+(4n-1)\times$
(77chcl)
$\displaystyle\psi_{0}^{2}(4n)-2(4n-1)\psi_{0}(2n)\psi_{0}(4n)+\frac{1}{2}(8n-3)\psi_{0}^{2}(2n)-\frac{16n^{3}+8n^{2}-1}{2n(2n+1)}\psi_{0}(4n)+$
$\displaystyle\frac{\left(8n^{2}-3n-2\right)}{2n}\psi_{0}(2n)+\psi_{0}(n)+\frac{1}{2}\psi_{0}\\!\left(n+\frac{1}{4}\right)+\frac{1}{8}(16n-3)\psi_{1}(1)+\frac{1}{2}\psi_{0}^{2}(1)+$
$\displaystyle\frac{3}{2}\psi_{0}(1)-\frac{1}{2}\psi_{0}\\!\left(\frac{1}{4}\right)-2n+\frac{5}{2}.$

By inserting the $\mathrm{I_{A}}$ expression (77cf) and the above
$\mathrm{I_{C}}$ expression (77chcl) into (77chci), we arrive at the claimed
result (25) of the case $a=0$ after observing the cancellations. In the same
manner, the capacity formulas (26)–(28) that correspond to the cases $a=1,2,3$
are obtained.

##  3 Conclusion

In this work, we compute the two major second-order statistics of fermionic
Gaussian states – the variance of von Neumann entanglement entropy and the
average entanglement capacity. A key ingredient in obtaining these results is
based on representing the underlying integrals into appropriate summations as
well as the use of semi closed-form basis leading to a complete cancellation.
We also conjectured an explicit expression for the variance of von Neumann
entropy valid for any subsystem dimensions.

## Acknowledgments

The authors wish to thank Lucas Hackl, Santosh Kumar, and Nicholas Witte for
correspondence. The work of Lu Wei is supported in part by the U.S. National
Science Foundation ($\\#$2150486 and $\\#$2006612).

##  Appendix A Summation representations of integrals in $\mathrm{I_{A}}$,
$\mathrm{I_{B}}$, and $\mathrm{I_{C}}$

In this appendix, we list the summation representations of the integrals in
$\mathrm{I_{A}}$, $\mathrm{I_{B}}$, and $\mathrm{I_{C}}$.

$\displaystyle\mathrm{A_{1}}=$
$\displaystyle\sum_{k=0}^{m-1}2(2a+4k+1)\left(\rule{0.0pt}{21.33955pt}\sum_{j=2k-2}^{2k}\frac{(-1)^{j}(j+1)_{2}(a+j+1)_{2}}{\Gamma(2k-j+1)\Gamma(j-2k+3)(2a+j+2k+1)_{3}}\right.\times$
(77chcm)
$\displaystyle\bigg{(}(\psi_{0}(a+j+3)-\psi_{0}(2a+j+2k+4)-\psi_{0}(j-2k+3)+\psi_{0}(j+3))^{2}-$
$\displaystyle\psi_{1}(2a+j+2k+4)+\psi_{1}(a+j+3)-\psi_{1}(j-2k+3)+\psi_{1}(j+3)\bigg{)}+$
$\displaystyle\sum_{j=0}^{2k-3}\frac{2(j+1)_{2}(a+j+1)_{2}}{(2k-j-2)_{3}(2a+j+2k+1)_{3}}(\psi_{0}(2a+j+2k+4)-\psi_{0}(a+j+3)+$
$\displaystyle\\!\left.\psi_{0}(2k-j-2)-\psi_{0}(j+3))\rule{0.0pt}{21.33955pt}\right)$

$\displaystyle\mathrm{A_{2}}=$
$\displaystyle\sum_{k=0}^{m-1}\frac{(2a+4k+1)\Gamma(2k+1)\Gamma(2a+2k+1)}{\Gamma(2a+4k+4)}\left(\rule{0.0pt}{21.33955pt}\sum_{i=0}^{2k}\frac{2(i+1)(2k-i+1)}{\Gamma(i+1)\Gamma(a+i+1)}\times\right.$
(77chcn)
$\displaystyle\frac{\Gamma(a+2k+2)^{2}}{\Gamma(2k-i+1)\Gamma(a+2k-i+1)}((\psi_{0}(a+2k+2)-\psi_{0}(2a+4k+4)-\psi_{0}(2)+$
$\displaystyle\psi_{0}(2k-i+2))(\psi_{0}(a+2k+2)-\psi_{0}(2a+4k+4)+\psi_{0}(i+2)-\psi_{0}(2))-$
$\displaystyle\psi_{1}(2a+4k+4))-\sum_{j=0}^{2k}\frac{(j+1)\Gamma(a+2k+1)\Gamma(a+2k+3)}{\Gamma(j)\Gamma(a+j+1)\Gamma(2k-j+1)\Gamma(a-j+2k+1)}\times$
$\displaystyle((\psi_{0}(a+2k+1)-\psi_{0}(2a+4k+4)+\psi_{0}(2k-j+2)-\psi_{0}(1))(\psi_{0}(a+2k+3)-$
$\displaystyle\psi_{0}(2a+4k+4)+\psi_{0}(j+2)-\psi_{0}(3))-\psi_{1}(2a+4k+4))-\sum_{j=0}^{2k}\frac{(2k-j+1)}{\Gamma(j+1)}\times$
$\displaystyle\frac{\Gamma(a+2k+1)\Gamma(a+2k+3)}{\Gamma(a+j+1)\Gamma(2k-j)\Gamma(2k-j+a+1)}((\psi_{0}(a+2k+3)-\psi_{0}(2a+4k+4)$
$\displaystyle+\psi_{0}(2k-j+2)-\psi_{0}(3))(\psi_{0}(a+2k+1)-\psi_{0}(2a+4k+4)+\psi_{0}(j+2)-$
$\displaystyle\left.\psi_{0}(1))-\psi_{1}(2a+4k+4))\rule{0.0pt}{21.33955pt}\right)+\sum_{k=0}^{m-1}\frac{4(2a+4k+1)\Gamma(2k+1)\Gamma(2a+2k+1)}{\Gamma(2a+4k+4)}\times$
$\displaystyle\sum_{j=0}^{2k}\sum_{i=0}^{2k-j-2}\frac{(2k-i-j-1)(i+j+3)\Gamma(a-j+2k)\Gamma(a+j+2k+4)}{\Gamma(i+1)\Gamma(2k-i+1)\Gamma(a+i+j+3)\Gamma(a-i-j+2k-1)}\times$
$\displaystyle\frac{1}{(j+1)_{3}}(\psi_{0}(a+j+2k+4)-\psi_{0}(2a+4k+4)+\psi_{0}(i+j+4)-\psi_{0}(j+4))$

$\displaystyle\mathrm{B_{1}}=$
$\displaystyle\sum_{k=0}^{m-1}\left(\psi_{0}(a+2k)+\psi_{0}(2a+2k)-2\psi_{0}(2a+4k)-\frac{1}{2}\left(\frac{a}{a+2k+1}+\frac{a}{a+2k}+\right.\right.$
(77chco) $\displaystyle\\!\left.\left.\frac{2}{2a+4k+1}\right)+1\right)^{2}$

$\displaystyle\mathrm{B_{2}}=$
$\displaystyle\sum_{k=0}^{m-1}\sum_{j=1}^{m-k-1}\frac{\Gamma(2a+2k+1)\Gamma(2j+2k+1)(2a+4k+1)(2a+4j+4k+1)}{2\Gamma(2k+1)\Gamma(2a+2j+2k+1)j^{2}(2j-1)^{2}(2j+1)^{2}}\times$
(77chcp)
$\displaystyle\frac{\left(2a^{2}j+a^{2}+2aj^{2}+4ajk+3aj+4ak+a+2j^{2}+4jk+j+4k^{2}+2k\right)^{2}}{(a+j+2k)^{2}(a+j+2k+1)^{2}(2a+2j+4k+1)^{2}}$
$\displaystyle\mathrm{I_{C}}=$
$\displaystyle\\!\left(\psi_{0}(a+2)-\psi_{0}(2a+3)\right)^{2}+\psi_{1}(a+2)-\psi_{1}(2a+3)+\sum_{k=1}^{m-1}2(2a+4k+1)\times$
(77chcq)
$\displaystyle\left(\rule{0.0pt}{21.33955pt}\sum_{j=2k-1}^{2k}\frac{(-1)^{j}(j+1)(a+j+1)}{(2a+j+2k+1)_{2}}\bigg{(}(\psi_{0}(j+2)-\psi_{0}(2a+j+2k+3)+\right.$
$\displaystyle\psi_{0}(a+j+2)-\psi_{0}(j-2k+2))^{2}+\psi_{1}(a+j+2)-\psi_{1}(2a+j+2k+3)+$
$\displaystyle\psi_{1}(j+2)-\psi_{1}(j-2k+2)\bigg{)}+\sum_{j=0}^{2k-2}\frac{2(j+1)(a+j+1)}{(2k-j-1)_{2}(2a+j+2k+1)_{2}}\times$
$\displaystyle\\!\left.\left(\psi_{0}(a+j+2)-\psi_{0}(2a+j+2k+3)-\psi_{0}(2k-j-1)+\psi_{0}(j+2)\right)\rule{0.0pt}{21.33955pt}\right)$

##  Appendix B List of summation identities

In this appendix, we list the finite sum identities useful in simplifying the
summations in appendix A. Except for the identity (77chcv) that seems new, the
others can be found in [8, 10, 17, 11, 32]. Here, it is sufficient to set
$a,b\geq 0,a\neq b$ in identities (77chcr)–(77chct), (77chcw)–(77chcx), and
$a>m$ in (B). The derivation of the identity (77chcv) is provided in the end
of this appendix.

$\displaystyle\sum_{k=1}^{m}\psi_{0}(k+a)=(m+a)\psi_{0}(m+a+1)-a\psi_{0}(a+1)-m$
(77chcr)
$\sum_{k=1}^{m}\psi_{1}(k+a)=(m+a)\psi_{1}(m+a+1)-a\psi_{1}(a+1)+\psi_{0}(m+a+1)-\psi_{0}(a+1)$
(77chcs)
$\displaystyle\sum_{k=1}^{m}\frac{\psi_{0}(k+a)}{k+a}=\frac{1}{2}\left(\psi_{1}(m+a+1)-\psi_{1}(a+1)+\psi_{0}^{2}(m+a+1)-\psi_{0}^{2}(a+1)\right)$
(77chct)
$\displaystyle\sum_{k=1}^{m}\frac{\psi_{0}(m+1-k)}{k}=\psi_{0}^{2}(m+1)-\psi_{0}(1)\psi_{0}(m+1)+\psi_{1}(m+1)-\psi_{1}(1)$
(77chcu)
$\displaystyle\sum_{k=1}^{m}\frac{\psi_{0}(m+1+k)}{k}=\psi_{0}^{2}(m+1)-\psi_{0}(1)\psi_{0}(m+1)-\frac{1}{2}\psi_{1}(m+1)+\frac{\psi_{1}(1)}{2}$
(77chcv) $\displaystyle\sum_{k=1}^{m}\psi_{0}(k+a)\psi_{0}(k+b)$
$\displaystyle=$
$\displaystyle(b-a)\sum_{k=1}^{m-1}\frac{\psi_{0}(a+k)}{b+k}+(m+a)\psi_{0}(m+a)\psi_{0}(m+b)-a\times$
(77chcw)
$\displaystyle\psi_{0}(a+1)\psi_{0}(b+1)-(m+a-1)\psi_{0}(m+a)+a\psi_{0}(a+1)-$
$\displaystyle(m+b)\psi_{0}(m+b)+(b+1)\psi_{0}(b+1)+2m-2$
$\displaystyle\sum_{k=1}^{m}\frac{\psi_{0}(k+b)}{k+a}=$
$\displaystyle-\sum_{k=1}^{m}\frac{\psi_{0}(k+a)}{k+b}+\psi_{0}(m+a+1)\psi_{0}(m+b+1)-\psi_{0}(a+1)\times$
(77chcx)
$\displaystyle\psi_{0}(b+1)+\frac{1}{a-b}(\psi_{0}(m+a+1)-\psi_{0}(m+b+1)-\psi_{0}(a+1)+$
$\displaystyle\psi_{0}(b+1))$
$\displaystyle\sum_{k=1}^{m}\frac{\psi_{0}(a+1-k)}{k}=$
$\displaystyle-\sum_{k=1}^{m}\frac{\psi_{0}(k+a-m)}{k}+(\psi_{0}(a-m)+\psi_{0}(a+1))(\psi_{0}(m+1)-$
$\displaystyle\psi_{0}(1))+\frac{1}{2}\left((\psi_{0}(a-m)-\psi_{0}(a+1))^{2}+\psi_{1}(a+1)-\psi_{1}(a-m)\right)$

The detailed steps in deriving the identity (77chcv) are as follows.

$\displaystyle\sum_{k=1}^{m}\frac{\psi_{0}(m+1+k)}{k}$ (77chcz)
$\displaystyle=$
$\displaystyle\sum_{k=1}^{m}\sum_{l=1}^{m}\frac{1}{k(k+l)}+\sum_{k=1}^{m}\frac{\psi_{0}(k+1)}{k}$
(77chda) $\displaystyle=$ $\displaystyle
2\sum_{k=1}^{m}\frac{\psi_{0}(k+1)}{k}+\sum_{l=1}^{m}\frac{\psi_{0}(m+1)-\psi_{0}(1)}{l}-\sum_{l=1}^{m}\frac{\psi_{0}(m+1+l)}{l}.$
(77chdb)

The equality (77chda) is obtained by representing $\psi_{0}(m+1+k)$ in terms
of $\psi_{0}(k+1)$ and $m$ rational terms by using (16). The equality (77chdb)
is obtained by changing the summation order to first sum over $l$, where one
obtains a same summation as in (77chcz) except for a negative sign. The
remaining sums can be simplified that leads to the identity (77chcv).

## References

## References

  * [1] Lubkin E 1978 Entropy of an n-system from its correlation with a k-reservoir J. Math. Phys. 19 1028
  * [2] Page D N 1993 Average entropy of a subsystem Phys. Rev. Lett. 71 1291-4
  * [3] 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
  * [4] Sánchez-Ruiz J 1995 Simple proof of Page’s conjecture on the average entropy of a subsystem Phys. Rev.E 52 5653-5
  * [5] 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
  * [6] Giraud O 2007 Distribution of bipartite entanglement for random pure states 40 2793
  * [7] Vivo P, Pato M P and Oshanin G 2016 Random pure states: Quantifying bipartite entanglement beyond the linear statistics Phys. Rev.E 93 052106
  * [8] Wei L 2017 Proof of Vivo-Pato-Oshanin’s conjecture on the fluctuation of von Neumann entropy Phys. Rev.E 96 022106
  * [9] Wei L 2019 On the exact variance of Tsallis entanglement entropy in a random pure state Entropy 21 539
  * [10] Wei L 2020 Skewness of von Neumann entanglement entropy 53 075302
  * [11] Huang Y, Wei L and Collaku B 2021 Kurtosis of von Neumann entanglement entropy (DOI: 10.1088/1751-8121/ac367c) (also available at arXiv:2107.10978)
  * [12] Sommers H-J and Życzkowski K 2004 Statistical properties of random density matrices J. Phys. A: Math. Gen. 37 35
  * [13] Osipov V, Sommers H-J and Życzkowski K 2010 Random Bures mixed states and the distribution of their purity 43 055302
  * [14] Borot G and Nadal C 2012 Purity distribution for generalized random Bures mixed states 45 075209
  * [15] Sarkar A and Kumar S 2019 Bures-Hall ensemble: spectral densities and average entropies 52 295203
  * [16] Wei L 2020 Proof of Sarkar-Kumar conjectures on average entanglement entropies over the Bures-Hall ensemble 53 235203
  * [17] Wei L 2020 Exact variance of von Neumann entanglement entropy over the Bures-Hall measure Phys. Rev.E 102 062128
  * [18] Li S-H and Wei L 2021 Moments of quantum purity and biorthogonal polynomial recurrence 54 445204
  * [19] Bianchi E, Hackl L and Kieburg M 2021 The Page curve for fermionic Gaussian states Phys. Rev.B 103 L241118
  * [20] Bhattacharjee B, Nandy P and Pathak T 2021 Eigenstate capacity and Page curve in fermionic Gaussian states (arXiv:2109.00557)
  * [21] Hackl L and Jonsson R H 2019 Minimal energy cost of entanglement extraction Quantum 3 165
  * [22] Hackl L and Bianchi E 2021 Bosonic and fermionic Gaussian states from Kähler structures SciPost Phys. Core 4 025
  * [23] Bengtsson I and Życzkowski K 2006 Geometry of Quantum States: an Introduction to Quantum Entanglement (Cambridge: Cambridge University Press)
  * [24] Nandy P 2021 Capacity of entanglement in local operators 19
  * [25] de Boer J, Järvelä J and Keski-Vakkuri E 2019 Aspects of capacity of entanglement Phys. Rev.D 99 066012
  * [26] Okuyama K 2021 Capacity of entanglement in random pure state Phys. Lett.B 820 136600
  * [27] Brychkov Y A 2008 Handbook of Special Functions: Derivatives, Integrals, Series and Other Formulas (Boca Raton: CRC Press)
  * [28] Forrester P 2010 Log-gases and Random Matrices (Princeton: Princeton University Press)
  * [29] Erdélyi A, Magnus W, Oberhettinger F and Tricomi FG 1954 Tables of Integral Transforms. Vol. II (New York: McGraw-Hill)
  * [30] Szegő G 1975 Orthogonal Polynomials (Provindence: American Mathematical Society)
  * [31] Abramowitz M and Stegun I A 1972 Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (New York: Dover)
  * [32] Milgram M 2017 On some sums of digamma and polygamma functions (arXiv:0406338v3)


# Interactive Bayesian Generative Models for Abnormality Detection in
Vehicular Networks

Nobel J. William, Ali Krayani, Lucio Marcenaro, and Carlo Regazzoni DITEN,
University of Genoa, Italy
emails<EMAIL_ADDRESS><EMAIL_ADDRESS>{lucio.marcenaro<EMAIL_ADDRESS>
###### Abstract

The following paper proposes a novel Vehicle-to-Everything (V2X) network
abnormality detection scheme based on Bayesian generative models for enhanced
network self-awareness functionality at the Base station (BS). In the learning
phase, multi-modal data signals contrived by the vehicles’ integrated and
sensing module are imbued into data-driven Generalized Dynamic Bayesian
network (GDBN) models. Following that, during the testing phase, an
Interactive Modified Markov Jump Particle filter (IM-MJPF) is utilized to
forecast forthcoming network states and vehicle trajectories by leveraging the
assimilated semantics embedded in the coupled multi-GDBNs. This approach
involves learning statistically correlated association between evolving
trajectories and network communication links. Security and surveillance of
Internet of Vehicles (IOVs) links are performed online with high detection
probabilities by matching predicted with observed network connectivity maps
(graphs).

###### Index Terms:

Generative Bayesian models, Jammer detection, Vehicular networks, V2X,
Abnormality detection.

##  I Introduction

Vehicular networks such as Vehicle-to-everything (V2X), consisting the
vehicle-to-vehicle (V2V), vehicle-to-pedestrian (V2P), vehicle-to-
infrastructure (V2I), and vehicle-to-network (V2N) communication are envisaged
to be key components in Intelligent Transportation Systems (ITS) [1].
Integrated Sensing and Communication (ISAC) is expected to play a pivotal role
in the future V2X standard by enabling next-generation connected autonomous
vehicles (CAVs) [2]. Moreover, the inherent openness and dynamism of V2X
results in a restricted availability of precise Channel State Information
(CSI). Consequently, ensuring secure transmission at the Physical Layer (PHY)
poses a noteworthy challenge in ISAC-aided V2X systems. In order to address
these real-world conditions, the implementation of ISAC technology
necessitates the utilization of advanced machine learning techniques that
effectively integrate communication and sensing data. This integration enables
the development of ambient (situational) awareness and thus enabling cognitive
services at the global network level, to systematically mitigate potential
threats or detect abnormalities. In high-mobility V2X scenarios, it is crucial
for the transmitter or base station (BS) to possess the capability to
anticipate connectivity requirements beforehand. This predictive ability
allows for effective attack mitigation and resource allocation to ensure
compliance with the demanding QoS requirements for V2X communications.

Networked Sensing (NS) involves utilizing components of cellular networks,
such as Base Stations (BS) or Roadside Units (RSUs), along with vehicular
nodes to enhance the overall awareness of V2X networks [3]. This approach
provides diverse perspectives on the communication and sensing environment.
However, V2X networked sensing requires more attention to improve the
estimation of network topology and enhance resource management and network
configuration. Also, the integrity of V2I and V2V links is paramount for
effective resource management in dense V2X environments, as jamming attacks
can severely impact road vehicle safety and traffic efficiency. The V2I links
include collaborative information showing the vehicle’s location and other
details from active V2V links exchanging basic safety messages (BSMs) and
cooperative awareness messages (CAMs). This enables the receiving cognitive
BS/RSU to learn a semantic mapping between vehicular trajectory states and
their respective V2V network connectivity states. Unfortunately, there is
limited research in this area. However, the Bayesian method can aid in
inferring network structure even with unjustified measurement errors [4].
Furthermore, it is vital to anticipate how the network topology will change
over time to ensure the safety of V2X operations by using the rules learned
from past experiences. This will help in identifying potential issues like
jamming or spoofing attacks that could cause abnormalities (i.e., network
failure) and determine their root cause. By understanding the cause of the
problem, effective measures can be taken to prevent such issues from happening
again in the future.

Ensuring the safe operation of V2X applications requires detecting
abnormalities in the topology caused by sensor attacks or new mobility
motions, as well as network abnormalities caused by jamming or spoofing
attacks. Several studies have addressed these concerns. For instance, [5]
proposes a framework for detecting and classifying sensor attacks that use GPS
and LiDAR sensors. In [6], researchers analyze the impact of abnormal driver
behaviour on data transmission in VANET networks. In [7], a new approach to
driver profiling and abnormality detection is proposed, using LSTM for
modelling drivers. Additionally, [8] suggests a method for identifying jammers
that affect vehicular networks. However, previous research has shown that
multi-modal sensors, such as GPS, LiDAR, video cameras, and RF antennas,
typically work independently to achieve separate goals and do not collaborate
to improve V2X performance via abnormality detection, consequently enhancing
security and resource management. Moreover, a PHY layer security oriented
jammer detection scheme using Dynamic Bayesian network (DBN) was proposed in
[9] to detect attacks on V2I links, however, this work doesn’t provide a
paradigm for estimating V2X connectivity state at the global network level for
abnormality detection. Another related work in [10] discusses the use of
neural network-based Long Short-term Memory (LSTM) and Multilayer Perceptron
(MLP) system for anomaly detection in vehicular communications. However, these
deep learning techniques lack the ability to explain explicitly the causal
probabilistic relationships that result in network abnormalities.

In this work, inspired by [11], we propose a new method of combining vehicle
trajectory data collected by a video camera sensor with RF signals received
from each vehicle in order to detect any abnormalities in the V2X network
caused by jamming attacks. This cognitive service is designed to identify
network disruptions at the BS/RSU by combining multiple sensory data in a
joint probabilistic representation. The method uses a probabilistic
representation structured in a coupled Multi-Generative Dynamic Bayesian
Networks (M-GDBN) generated during training. It also utilizes real-time
vehicle trajectory measurements to predict the current and future V2X network
connectivity structure. The Interactive Modified Markov Jump Particle filter
(IM-MJPF) is used to make these predictions by inferring what the vehicles
could transmit based on their mobility pattern. The BS/RSU performs the
inference, considering the full or partial observability of both the vehicular
trajectories and the V2X network’s connectivity state/structure. The proposed
method enables the BS/RSU to jointly predict the evolution of both vehicles’
trajectories and RF signals and detect any abnormality causing a network
failure at the high abstraction level in the probabilistic representation. It
can also explain the cause behind the detected abnormalities and the sensor it
is coming from (e.g., video camera sensor or RF antenna) through the causal
probabilistic relationships encoded in the representation. In this study, we
conducted simulations to avoid the complexity of collecting data from multiple
vehicles during the initial training phase. This would have been difficult to
implement in real-world empirical analysis. We assumed that the Channel State
Information (CSI) reported to the BS was accurate because of the multi-agent
(vehicle) environment perception in our scenario. This integration enables
network-level cognitive services in the V2X context, specifically aimed at
mitigating anomalies.

##  II System Model and Problem Formulation

Figure 1: Illustration of the System model.

Assume $\mathrm{N}$ connected vehicles traversing a freeway within a single-
cell vehicular network, each transmitting it’s respective positional
information $p_{n,t}=[x_{n,t},y_{n,t}]$ where $\mathrm{n}\in\mathrm{N}$ at
each time instant $\mathrm{t}$ jointly with supplementary information
regarding it’s active V2V connections and CSI (i.e. channel quality, signal
strength, interference, and other relevant metrics) to the BS located at
$p_{BS}=[x_{BS},y_{BS}]$ via the $\mathrm{n^{th}}$ V2I link. Each vehicle
exchanges messages with other vincinal vehicles $n^{\prime}$ located at
$p_{n^{\prime},t}=[x_{n^{\prime},t},y_{n^{\prime},t}]$ on the freeway via the
$\mathrm{k^{th}}$ V2V links. The total bandwidth available is divided into
$\mathrm{N}$ physical resource blocks (PRBs) and each vehicle on the freeway
accesses a single PRB from the set of mutually orthogonal spaced frequency
bands denoted as $\mathcal{F}=\\{f_{1},f_{2},\dots,f_{N}\\}$ for its V2I link
(without loss of generality we assume $\mathrm{N}=\mathcal{F}$ here). The set
of V2V links $\mathcal{K}=\\{f_{1},f_{2},\dots,f_{K}\\}$ reuse the V2I up-link
in an interference free optimal manner as discussed in [12]. Moreover,
scheduling and interference management of V2V traffic is assisted by the BS
(in Mode 3) via the control signaling that assigns the resources used for V2V
communications in a dynamic manner. It is to be duely noted in this scenario,
the BS effectively assists in V2V links CSI management (not necessarily always
done in mode 3). The antagonist here is a constant (single) and periodic
(multi-attack) road side jammer/compromised RSU residing within the cell at
$p_{j}=[x_{j},y_{j}]$. The system model is illustrated as in Fig. 1. The
intelligent jammer injects interference at mulitple frequencies (PRBs) and
overloads the dynamic V2V links. The jamming detection problem (for each
vehicle) can be devised as the following binary hypothesis test:

$\left\\{\begin{array}[]{l}\mathcal{H}_{0}:\mathrm{z}_{t,k}=\mathrm{g}_{t,k}^{nn^{\prime}}\mathrm{x}_{t,k}+\mathrm{v}_{t,k},\\\
\mathcal{H}_{1}:\mathrm{z}_{t,k}=\mathrm{g}_{t,k}^{nn^{\prime}}\mathrm{x}_{t,k}+\mathrm{g}_{t,k}^{jn^{\prime}}\mathrm{x}_{t,k}^{j}+\mathrm{v}_{t,k},\\\
\end{array}\right.$ (1)

where $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ denote absence and presence of
jammer, respectively. $\mathrm{x}_{t,k}$ is the reference signal,
$\mathrm{v}_{t,k}$ is the additive white Gaussian noise (with variance
$\sigma^{2}_{n}$) and $\mathrm{z}_{t,k}$ is the received signal at the
individual vehicle transreceiver at time $\mathrm{t}$ over the
$\mathrm{k^{th}}$ V2V link. $\mathrm{g}_{t,k}^{nn^{\prime}}$ is the channel
gain between vehicles $n$ and $n^{\prime}$. $\mathrm{x}_{t,k}^{j}$ is the
jammer signal that overwhelms the vehicle radio receivers with channel power
gain $\mathrm{g}_{t,k}^{jn^{\prime}}$ at $n^{\prime}$ vehicle.
$d^{nn^{\prime}}_{t}=\sqrt{[(x_{n,t}-x_{n^{\prime},t})^{2}+(y_{n,t}-y_{n^{\prime},t})^{2}}$
is the distance between vehicle $n$ and vehicle $n^{\prime}$ at time $t$.

Mean while for the V2I links, the channel gain $g^{nBS}_{t,n}$ of the $n^{th}$
V2I link to the BS is $g^{nBS}_{t,n}=\alpha^{nBS}_{t,n}h^{nBS}_{t,n}$, where
$\alpha^{nBS}_{t,n}$ = $\mathrm{G\beta}d^{-\gamma}_{t,nBS}$ captures the
large-scale fading effects (path-loss and fading component). $\mathrm{G}$ is
the path-loss constant and $\beta$ is log normal shadow fading random
variable, $d^{nBS}_{t}=\sqrt{[(x_{n,t}-x_{BS})^{2}+(y_{n,t}-y_{BS})^{2}}$ is
the distance between $\mathrm{n^{th}}$ vehicle and the BS at time $t$.
$\gamma$ is power decay exponent and $h_{t,n}$ is the small-scale fading
component distributed according to $\mathcal{N}(0,1)$ as seen in [13]. To
accommodate the global network perspective, we reformulate on road vehicles as
nodes and the communication links between them at time $t$ as edges. We assume
unweighted and undirected edges here to represent existence of bidirectional
connection between vehicles $n$ and $n^{\prime}$. Thus, a discrete time
dynamic graph $\mathcal{G}_{t}=({V_{t},\mathcal{E}_{t}})$ with $V$ vertices
(nodes) and $\mathcal{E}$ edges as mapped in adjacency matrix
$\mathbf{A}\in\mathbb{R}^{V\times V}$ where $V=N$. The binary Hypothesis test
mentioned in (1) is modified to a graph perturbation scenario as follows:

$\left\\{\begin{array}[]{l}\mathcal{H}_{0}:\mathrm{Z}_{\mathcal{G}_{t}}=\mathrm{G}_{t}+\mathrm{v}_{t}\\\
\mathcal{H}_{1}:\mathrm{Z}_{\mathcal{G}_{t}}=\mathrm{f}(\mathrm{G}_{t})+\mathrm{v}_{t}\\\
\end{array}\right.$ (2)

where $\emph{f}(\mathrm{G}_{t})$ is a graph perturbation function informed by
the jammer, $\mathrm{v}_{t}$ is the process noise and
$\mathrm{Z}_{\mathcal{G}_{t}}$ is the observed graph.

##  III Proposed Scheme for Joint Network Structure Inference and Abnormality
Detection

###  III-A Environment Representation

We represent the physical and RF environment (sensory signals) onto a model
based data driven framework that internally embodies a Generative Bayesian
probabilistic model. Here, BS transceiver observes measurements $Z_{t}^{m}$
reported from individual vehicles (which includes position, estimated CSI and
activated V2V connections) on the freeway at each time instant $t$. The
individual physical states (positions) of $\mathrm{N}$ vehicles at time $t$
are mapped at a meta-cluster level termed as a Word $W^{p}_{t}$ within a
learned global odometry/positional dictionary $D^{p}$. Each Word comprises of
$N$ symbolised letters (aka clusters or superstates) $S^{p}_{n}$ with unique
mean $\mu$ and variance $\Sigma$. Similarly, the activated V2V connections at
each vehicle reported to the BS are contained in symbolized letter as
$S^{c}_{n}$. It may be pertinent to note that several symbols fuse to make a
letter here and is a function of the number of vehicles on the freeway.
Multiple letters ($S^{c}_{n}$) based on data received at the BS imbibe into
the Word $W^{c}_{t}$ representing the network connectivity state/map at time
$t$. A global dictionary $D^{c}$ containing multitude of Words $W^{c}$
expressing the universe of learned vehicular network configurations is
formulated during training. The dynamic $m$-th data signal evolution function
using Generalized State Space model at multiple hierarchical levels can be
expressed as:

$\mathrm{\tilde{W}^{m}}_{t}=\mathrm{{{f}^{1}}\hskip
1.0pt(\tilde{W}^{m}_{t-1})}+\tilde{w}^{m}_{t},$ (3)
$\mathrm{\tilde{S}^{m}}_{t,n}=\mathrm{{{f}^{2}}\hskip
1.0pt(\tilde{S}^{m}_{t-1,n})}+\tilde{w}^{m}_{t,n},$ (4)

where $m\in$ (positional, communication) and $\mathrm{f}(.)$ is a non-linear
function of the argument and $w^{m}_{t}$, $w^{m}_{t,n}$ encode the noise
affecting the processes at different hierarchical levels. The discrete
equations (3) and (4) illustrate the temporal evolution of Generalized Words
and Generalized letters respectively during the higher latent state
evolutionary process. (3) maps transitions from previous Generalized Word to
the current Word, embodying the dynamics of transitions between discrete-level
positional or network configurations over time. Moving down the hierarchical
network, the continuous state dynamic evolution equation (expressed in (5))
and the measurement model mapping the continuous state to the observation
space (in (6)) are defined as follows:

$\mathrm{\tilde{X}^{m}}_{t,n}=\mathrm{A\tilde{X}^{m}_{t-1,n}}+\mathrm{B\tilde{U}}_{\tilde{S}^{m}_{t,n}}+\tilde{w}^{m}_{t,n},$
(5)
$\mathrm{\tilde{Z}^{m}}_{t,n}=\mathrm{H\tilde{X}^{m}_{t,n}}+\tilde{w}^{m}_{t,n}.$
(6)

More specifically, the motion dynamics of each individual $n^{th}$ vehicle in
the cell or the V2V connections states are baked into (5). To account for
measurement uncertainty, the noise term $\tilde{w}^{m}_{t,n}$ is introduced.
$\mathrm{A}\in\mathbb{R}^{2d,2d}$ and $\mathrm{H}\in\mathbb{R}^{2d}$ are the
dynamic and control matrices, respectively.
$\mathrm{\tilde{U}}_{\tilde{S}^{m}_{t,n}}$ constitutes the control vector
which represents the dynamic rules of the signal’s temporal evolution. The
measurement model in (6) encapsulates the mapping from the observed sensory
signals $\mathrm{\tilde{Z}^{m}}_{t,n}$ to the underlying hidden state
$\mathrm{\tilde{X}^{m}_{t,n}}$ through the measurement matrix
$\mathrm{B}\in\mathbb{R}^{2d}$ where $d$ is data dimensionality and depends on
the features used for training such as XY-coordinates, velocities, relative
distances, RF IQ data, etc.

###  III-B Learning Dictionaries (containing multi-GDBNs)

The dynamic hierarchical models expressed as (3), (4), (5) and (6), provide us
with the scaffold to build the Generalized Dynamic Bayesian Network [14]
(refer to Fig. 2 ”blue-shaded region”) that embodies conditional dependencies
between latent variables and measured sensory states into a structure known as
‘probabilistic graphical model’ (PGM). The generative process model imbibed
into the PGM, explaining the sensory signals can be factorised as follows:

$\scriptsize\mathrm{P}(\mathrm{W}_{t}^{m},\mathrm{\tilde{S}}^{m}_{t,n},\mathrm{\tilde{X}}^{m}_{t,n},\mathrm{\tilde{Z}}^{m}_{t,n})=\mathrm{P}(\tilde{W}^{0}_{t})\mathrm{P}(\tilde{S}^{0}_{t,n})\mathrm{P}(\tilde{X}^{0}_{t,n})\\\
\prod_{t}P(\mathrm{W}_{t}^{m})\prod_{t,n}\mathrm{P}(\mathrm{\tilde{S}}^{m}_{t,n}|\mathrm{W}_{t}^{m})\mathrm{P}(\mathrm{\tilde{X}}^{m}_{t,n}|\mathrm{\tilde{S}}_{t,n}^{m})\mathrm{P}(\mathrm{\tilde{Z}}^{m}_{t,n}|\mathrm{\tilde{X}}^{m}_{t,n}),$
(7)

where $\mathrm{P}(\tilde{W}^{0}_{t})$, $\mathrm{P}(\tilde{S}^{0}_{t,n})$ and
$\mathrm{P}(\tilde{X}^{0}_{t,n})$ are initial prior distributions. The
likelihood is
$\mathrm{P}(\mathrm{\tilde{Z}}^{m}_{t,n}|\mathrm{\tilde{X}}^{m}_{t,n})$.

The transition probabilities depicting the hierarchical and temporal
relationships of the generalized state-space model are
$\mathrm{P}(\mathrm{\tilde{X}}^{m}_{t,n}|\mathrm{\tilde{S}}_{t,n}^{m})$ and
$\mathrm{P}(\mathrm{\tilde{S}}^{m}_{t,n}|\mathrm{W}_{t}^{m})$. The causality
relationship chain between the random hidden and observed variables
$\mathrm{{W}_{t}^{m}}$ $\rightarrow$ $\mathrm{\tilde{S_{t,n}^{m}}}$
$\rightarrow$ $\mathrm{\tilde{X_{t,n}^{m}}}$ $\rightarrow$
$\mathrm{\tilde{Z_{t,n}^{m}}}$ are imbued in (7).

The cognitive BS instantiates the perception of the vehicular cell network
with infantile knowledge (based on static assumptions regarding environmental
state evolution). Thus the BS predicts the vehicle trajectory or activated V2V
links at each vehicle by enacting a Null-force filter based on the model:
$\mathrm{\tilde{X}_{t,n}^{m}}$ = $\mathrm{A\tilde{X}^{m}_{t-1,n}}$\+
$\tilde{w}_{t,n}^{m}$. This model diverges from (5) due to the control vector
being null ($\mathrm{\tilde{U}_{{S}_{t,n}}^{m}}=0$). Hence, to explore the
dynamic rules informing the environmental states over time, new rules need to
be ascertained by exploiting generalized errors (GE) which is the difference
between the observation and prediction. The GEs projected at the different
hierarchical level are as follows: for the measurement space,
$\tilde{\varepsilon}_{Z_{t}}^{m}$ = $\mathrm{\tilde{Z}}_{t,n}^{m}$ \-
$\mathrm{H\tilde{X_{t,n}^{m}}}$ and for the generalized state space,
$\tilde{\varepsilon}_{{\mathrm{X}_{t,n}^{m}}}=\mathrm{H}^{-1}\mathrm{Z}_{t,n}^{m}-\mathrm{X}_{t,n}^{m}$.
Growing Neural Gas (GNG), an unsupervised clustering algorithm is employed on
the GEs and ensures dimensionality reduction of the environment
representational data as illustrated Fig. 4 into distinct identifiable symbols
which enable our approach to computationally tractable and intrepretable.
Assembling symbols to form letters which engender Words have dynamic and vivid
representational properties of the global state space. Furthermore, Fig. 5
illustrates the knowledge encoded about the dynamics at different scales of
this hierarchical representation. Clusters encoding the discrete variables are
symbolised as $\mathrm{\tilde{S}_{t}^{m}}$ = $({\mathrm{\tilde{S}_{1}^{m}}}$,
${\mathrm{\tilde{S}_{2}^{m}}}$, ${\mathrm{\tilde{S}_{3}^{m}}}$, …,
${\mathrm{\tilde{S}_{{\mathcal{L}_{n}}}^{m}}})$, where $\mathcal{L}_{n}^{m}$
denotes the number of clusters (aka letters) of the $m$ modality signal type
at $n^{th}$ vehicle and each cluster
$\mathrm{\tilde{S}_{{\mathcal{L}}_{n}^{m}}}\in$ $\mathrm{\tilde{S}^{m}}$
follows a Gaussian model characterised by
$\mathrm{\tilde{S}_{{\mathcal{L}_{n}}}^{m}}\sim\mathcal{N}(\tilde{\mu}_{\mathrm{\tilde{S}}_{\mathcal{L}_{n}^{m}}}=[\mu_{{\mathrm{\tilde{S}}_{\mathcal{L}_{n}^{m}}}},{\dot{\mu}}_{{\mathrm{\tilde{S}}_{\mathcal{L}_{n}}^{m}}}],\Sigma_{\mathrm{\tilde{S}}_{{\mathcal{L}}^{m}}})$.
Information on the dynamic state transitions (between letters) are encoded
into a time varying transition matrix $(\Pi_{\tau})$ by learning the time
mapped transition probabilities
$\pi_{i,j}=P(\mathrm{\tilde{S}}_{t}^{m}=i|\mathrm{\tilde{S}}_{t-1}^{m}=j,\tau)$
where $\tau$ refers to the time spent in cluster $j$ before shifting on to
cluster $i$. An additional sub-dictionary of letters $\mathcal{D}_{l}^{m}$ for
$m=c$ is formed to accommodate for the letter state transition mapping and
future network state tracking at each individual vehicle in the V2X network.

###  III-C Learning Coupled Multi-GDBN

The dictionary learning stage mentioned in the prior section is effectuated by
learning a letter subdictionary and the symbolized letters transition
probability model. Letters from subdictionaries of N vehicles are time
synchronised to form a global Word at every time instance $t$. Each global
Word here represents a configuration of a set of N vehicles on the road and
the corresponding active V2V network connectivity graph between the N
vehicles. Post fusion of Multi-GDBNs for each module $m$, we analyse time
synchronized joint multi-cluster interaction behaviour between vehicle
trajectory configurations and the respective Global V2X network which contains
the V2V and V2I links in this scenario. By tracking the interaction, we imply
the identification of multi-cluster activation across both modalities during a
viable experience. This learning and interaction is embodied into a Coupled
Multi-GDBN (CM-GDBN) composed of N GDBN representing each signal modality
(Trajectory and activated V2V connections). Thus, a stochastic coupling of the
activated hidden states across the Multi-GDBNs in the two modalities
$(\mathrm{\tilde{C}}_{t})\text{=}[\mathrm{\tilde{W}}_{t}^{p},\mathrm{\tilde{W}}_{t}^{c}]$
at time $t$ is performed. An Interactive Matrix
$\Phi\in\mathbb{R}^{W^{p},W^{c}}$ is generated encoding the activated Word
pairs which map into the firing clusters at the letter-level (refer to Fig.
5).

Figure 2: Cognitive Architecture: Interactive Coupled Multi-GDBN model.

$\Phi$ allows us to infer the activated V2V links at a given time $t$ from the
trajectory states (positions) of $N$ vehicles,

$\scriptstyle\boldsymbol{\Phi}=\begin{bmatrix}\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{1}^{c}|\mathrm{\tilde{W}}_{1}^{p})&\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{2}^{c}|\mathrm{\tilde{W}}_{1}^{p})&\dots&\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{{L{(W^{c})}}}|\mathrm{\tilde{W}}_{1}^{p})\\\
\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{1}^{c}|\mathrm{\tilde{W}}_{2}^{p})&\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{2}^{c}|\mathrm{\tilde{W}}_{2}^{p})&\dots&\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{{L{(W^{c})}}}|\mathrm{\tilde{W}}_{2}^{p})\\\
\vdots&\vdots&\ddots&\vdots\\\
\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{1}^{c}|\mathrm{\tilde{W}}_{1}^{p})&\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{2}^{c}|\mathrm{\tilde{W}}_{{L{(W^{p})}}})&\dots&\scriptstyle\mathrm{P}(\mathrm{\tilde{W}}_{{L{(W^{c})}}}|\mathrm{\tilde{W}}_{{L{(W^{p})}}})\end{bmatrix}$
(8)

where ${L}{(W^{{p}\text{ or }{c}})}$ represents the total number of words in
the Dictionary of each modality $m$.

###  III-D Joint Trajectory-Network Inference and Perception

To effectuate cognitive tasks and decisions, a Bayesian generative filter
named Interactive Modified Markov Jump filter (IM-MJPF) illustrated in Fig. 2
in detail is employed to access the knowledge embedded in CM-GDBN. The BS
starts inferring the network connectivity map of the V2X network based on
reported vehicular node positions in the cell using IM-MJPF. IM-MJPF is based
on the Modified Markov Jump filter (M-MJPF) in [14] howbeit with the
introduction of a bank of GDBNs (multi-GDBNs) at each informed modalities. It
is important to recollect that the M-MJPF hosts a Particle filter (PF) and
Kalman filter (KF) at the discrete and continuous level, respectively. Thus
aiding in multiple predictions that propagate in a feed-forward (top-down) and
feedback (bottom-up) manner which modify and enforce inherent beliefs about
the state of the vehicular configurations and corresponding on road active V2V
communication links. The PF contained in IM-MJPF proliferates a set of (P)
equally weighted $(1/L)$ particles, such that
$\mathrm{\tilde{W}}_{t}^{p},\mathrm{\tilde{S}}_{t,n}^{p}$
$\sim{{\pi({\tilde{W}_{t,n}^{p}})}}$, where ${\tilde{W}_{t,n}^{p}}$ positional
state of the global scene and $L$ being number of vehicles on the road. The BS
uses $\Phi$ in (8) to infer $({\mathrm{\tilde{W}}_{t}^{c}})$ the vehicular
network map at time $t$ and by further deduction into discrete superstate
variables $\mathrm{\tilde{S}}_{t,n}^{p}$. A PF is used to predict future
multiple superstate guided by top-down and inter-time predictions as mention
in (3).

The posterior probability is expressed as follows:
$\pi({\mathrm{\tilde{W}}_{t}^{p}},{\mathrm{\tilde{W}}_{t}^{c}}){=}P({\mathrm{\tilde{W}}_{t}^{c}}|{\mathrm{\tilde{W}}_{t-1}^{c}},{\mathrm{\tilde{W}}_{t-1}^{p}})$
and
$\pi({\mathrm{\tilde{W}}_{t}^{p}}){=}P({\mathrm{\tilde{W}}_{t}^{p}}|{\mathrm{\tilde{W}}_{t-1}^{p}},{\mathrm{\tilde{W}}_{t-1}^{c}})$
at the Metacluster (Word) level. The diagnostic (bottom-up update) message in
the network leading up to the meta-cluster level is as follows:

$\scriptsize\lambda\left({\mathrm{\tilde{W}}}_{t}^{c}\right)=\prod_{n=1}^{L}\lambda\left(\mathrm{\tilde{S}}_{t,n}^{p}\right)\prod_{n=1}^{L}P\left(\mathrm{\tilde{S}}_{t,n}^{p}\mid{\mathrm{\tilde{W}}}_{t,n}^{p}\right).$
(9)

At the superstate discrete level the feed-forward and feed-backward messages
are expressed as follows:

$\displaystyle\pi\left(\mathrm{\tilde{S}}_{t,n}^{p}\right)=P\left(\mathrm{\tilde{S}}_{t,n}^{p}\mid\mathrm{\tilde{S}}_{t-1,n}^{p},\ldots,\mathrm{\tilde{S}}_{t-1,L}^{p}\right),L\neq
n$ $\displaystyle\text{and
}\lambda\left(\mathrm{\tilde{S}}_{t,n}^{p}\right)=\lambda\left(\mathrm{\tilde{X}}_{t,n}^{p}\right)P\left(\mathrm{\tilde{X}}_{t,n}^{p}\mid\mathrm{\tilde{S}}_{t,n}^{p}\right).$

Posterior probabilities at the continuous state level $\tilde{X}$ mapped at
individual vehicles is defined by the expression:

$\pi\left(\mathrm{\tilde{X}}_{t,n}^{p}\right)=P\left(\mathrm{\tilde{X}}_{t,n}^{p},\mathrm{\tilde{S}}_{t,n}^{p}\mid\mathrm{\tilde{Z}}_{t-1,n}^{p}\right)\\\
=\int P\left(X_{t,n}^{p}\mid
X_{t-1,n}^{p},S_{t,n}^{p}\right)\lambda\left(X_{t-1,n}^{p}\right)dX_{t-1,n}^{p}$
(10)

and the diagnostic message is,

$\scriptsize\lambda\left(\mathrm{\tilde{X}}_{t,n}^{p}\right)=P\left(\mathrm{\tilde{Z}}_{t,n}^{p}\mid\mathrm{\tilde{X}}_{t,n}^{p}\right)\sim
N\left(\mu_{\mathrm{Z}_{n}^{p}},\Sigma_{{\mathrm{Z}}_{n}^{p}})\right).$ (11)

The posterior distribution is updated by the new sensory evidence
$\mathrm{\tilde{Z}_{t}^{m}}$ by purposing the diagnostic messages discussed
i.e. $\lambda\left(\mathrm{\tilde{X}}_{t,n}^{m}\right)$ ,
$\lambda\left(\mathrm{\tilde{S}}_{t,n}^{m}\right)$ and
$\lambda\left({\mathrm{\tilde{W}}}_{t}^{m}\right)$ at the State, Superstate
and Word level respectively.

###  III-E Vehicular Network Abnormality Detection

The BS surveys the current state of the vehicular network by verifying if the
V2V communications are under attack or in the optimal state based on predicted
road and network configurations. In this work we introduce a new abnormality
(Bayesian surprise) indicator at the Word level that calculates the similarity
between the predicted word and the observed word encoding the $m$ modules
(recall $m=p$ or $c$) based on the symmetric Kullback-Leibler-Divergence
Abnormality (KLDA) [15] defined as:

$\small{\Upsilon_{\tilde{W}_{t}^{m}}}=D_{KL}\left(\pi\left(\tilde{W}_{t}^{m}\right)\|\lambda\left(\tilde{W}_{t}^{m}\right)\right)+\\\
D_{KL}\left(\lambda\left(\tilde{W}_{t}^{m}\right)\|\pi\left(\tilde{W}_{t}^{m}\right)\right).$
(12)

KLDA is used to compare the predictive message $\pi(\tilde{W}_{t}^{m})$ and
diagnostic message $\lambda(\tilde{W}_{t}^{m})$ defined in the previous sub-
section which tracks abnormalities at the global V2X network. BS enacts to
verify if the current situation is normal or not i.e. the V2V links are under
attack by a jammer based on the hypothesis test:

$\left\\{\begin{array}[]{l}\mathcal{H}_{0}:\text{if
}{\Upsilon_{\tilde{W}_{t}^{m}}}<\xi\\\ \mathcal{H}_{1}:\text{if
}{\Upsilon_{\tilde{W}_{t}^{m}}}>\xi\\\ \end{array}\right.$ (13)

where $\xi$ =
$\mathbb{E}[{\Upsilon_{\tilde{W}_{t}^{m}}}]+\phi\sqrt{\mathbb{V}[{\Upsilon_{\tilde{W}_{t}^{m}}}]}$
(ideally $\phi=3$). Within $\xi$, ${\Upsilon_{\tilde{W}_{t}^{m}}}$ refers to
abnormality signals during training (i.e. under reference situation with no
jammer) (demonstrated as blue line in Fig. 8).

###  III-F Evaluation metrics

ROC (Receiver Operating Characteristic) curves are used in this study to
quantitatively assess the performance of our jamming detection scheme. ROC
curves here illustrates the trade-off between the true positive rate
(sensitivity) and the false positive rate (1 - specificity) for different
detection thresholds.

##  IV Simulation and Experimental Results

In this section we perform a performance test of our scheme to detect
vehicular network abnormalities using extensive simulations. Scenarios on the
freeway with multiple number of vehicles (where N = 4, 6, 10 and 14) are
studied using our detection scheme. We consider N vehicles maneuvering safely
in an interactive environment relaying their positions, CSI and active V2V
connections to BS. The vehicles trace a predefined trajectory with velocities
from the Interstate 80 Freeway Dataset [16] illustrated in Fig. 3. performing
different maneuvers on a freeway. The BS assigns one PRB to each vehicle in
the cell to realize the V2I link. In this scenario, the transmitted signal
carrying the vehicle’s status information and the jamming signal both utilize
Quadrature Phase Shift Keying (QPSK) modulation. The simulation is configured
with specific parameters: a carrier frequency of $2$ GHz, a bandwidth (BW) of
$1.4$ MHz, a cell radius of $500$ m, an RSU antenna height of $25$ m with a
gain of $8$ dBi, a receiver noise figure of $5$ dB, a vehicle antenna height
of $1.5$ m with a gain of $3$ dBi, a V2I transmit power of $23$ dBm and jammer
transmit power of $23$ dBm. Additionally, the simulation incorporates a
Signal-to-Noise Ratio (SNR) of $20$ dB, a path loss model based on
$(128.1+37.6\log{d})$, log-normal shadowing with a standard deviation of $8$
dB, and a fast fading channel following the Rayleigh distribution [11]. While
the standard minimum distance at which a vehicle-to-everything (V2X) nodes can
detect a signal depends on the specific technology and implementation being
used. In our system model, without loss of generality with respect to short-
range or cellular communications, for V2V transmit power of upto $23$ dBm is
converted to a vehicle to vehicle relative distance oriented relationship.
Therefore, we assume the adjacency matrix of graph $\mathcal{G}_{t}$ encodes
the V2V connectivity which again is explicitly mapped to the relative
distances between vehicles ${d_{ij}}$, where $i$ and $j$ indicate vehicle
index number $\in$ N. The maximum connectivity distance between $i^{th}$
vehicle to $j^{th}$ vehicle for all N vehicles on the freeway is defined by
scalar value $d_{k}$ (we assume $d_{k}$ $\leq 10m$ as connected) and thereby
transforming by abstraction the relative distance matrix to connectivity
adjacency matrix of the global vehicular network.

The multi-vehicle tracking capability of the BS is demonstrated via trajectory
prediction in Fig. 6 by inferring from variables that encode spatio-temporal
relationships between on road vehicles embedded inside CM-GDBN during the
learning phase. The spatio-temporal configurations of vehicles affect the
overlaying dynamic V2V connections persistent during an interaction/manoeuvre
experience. Consequently, the connectivity map between vehicles can be
estimated based on the spatio-temporal vehicular interactions as depicted in
Fig. 7.

The optimal decision making process at AI enabled cognitive BS is elucidated
in Fig. 8, where mismatches between expected and observed words trigger
regular abnormality significantly to cross a predefined threshold as defined
in (13). The abnormality tracking results for single and multi attack
demonstrate the effectiveness of the proposed framework and engender us to
interpret the observed phenomenon in more explainable fashion through causal
reasoning capability of PGMs. Moreover, in this study we performance test our
approach to detecting anomalies caused by the (two) jamming attack types,
single and multiple attack over different window intervals. ROC curves verify
the efficacy of our jamming detection scheme for the attack scenarios in Fig.
9 by learning the optimum environment and thereby embedding the knowledge into
an optimal dynamic Bayesian network. While the ROC performance seen here are
strong, it may be duly noted that as multiplicity of attacks increase ROC
performance remains adept at the KLDA level in comparison to random observed
abnormalities.

(a)

(b)

(c)

Figure 3: Examples of four different sets of vehicle configuration: a) N=6, b)
N=10, c) N=14 from NGSIM dataset.

(a)

(b)

Figure 4: Clustering trajectories using GNG: (a) Vehicle $4$, (b) Vehicle
$10$.

(a)

(b)

(c)

(d)

(e)

Figure 5: Learning stage: Letter-TM for (a) Normal and (b) Jammer. Word-TM for
(c) Normal and (d) Jammer and (e) Interaction Matrix for N = $4$ scenario.

(a)

(b)

Figure 6: Predicted Trajectories at the BS for different number of Vehicles:
(a) N = 10, (b) N = 14.

(a)

(b)

(c)

(d)

Figure 7: Predicted Network Graphs $\mathcal{(}G_{t})$ using Dictionaries: (a)
$4$ Vehicles, (b) $6$ Vehicles, (c) $10$ Vehicles, and (d) $14$ Vehicles.

(a)

(b)

(c)

Figure 8: BS tracking abnormality: (a) Trajectory Normality, (b) Abnormality
(single window attack), c) Abnormality (double window attack) at the Network
(Word) level for N = 10 scenario.

(a)

(b)

(c)

(d)

Figure 9: ROC (KLDA) Curves for Jammer detection: (a) 4 Vehicle Double-attack,
(b) 6 Vehicle Double-attack, (c) 10 Vehicle Double-attack, (d) 10 Vehicle
Single-attack.

##  V Conclusion

Our proposed approach uses Multi-CDBNs to imbibe new observations at scale
into its model over time, thus enabling incremental learning and inherent
explainability. This approach performs predictions based on past observed
data, which provides an interpretable representation of the model’s reasoning.
By predicting future V2X network states and structure through vehicle
trajectories on the freeway, discrepancies in vehicular connectivity can be
detected by matching the connectivity prediction with on-road radio
environment observations. Our framework uses a Bayesian generative mechanism
that effectuates probabilistic belief-based inferences of the current and
future network connectivity state, inspired by the Bayesian Brain hypothesis.
This approach has been validated successfully via multiple simulation tests.
In the future, we plan to perform attack localization and sensitivity analysis
tests for varying network sizes and channel conditions/impairments to evaluate
and further validate our proposed framework.

## References

  * [1] Study on LTE-Based V2X Services (Release 14), Tech. Specification Group Serv. Syst. Aspects (TSG SA).  Technical report, 3GPP TR 36.885, 2016.
  * [2] Yi Zhong et al.  Empowering the V2X Network by Integrated Sensing and Communications: Background, Design, Advances, and Opportunities.  IEEE Network, 36(4):54–60, 2022.
  * [3] Thorsten Wild, Volker Braun, and Harish Viswanathan.  Joint design of communication and sensing for beyond 5G and 6G systems.  IEEE Access, 9:30845–30857, 2021.
  * [4] Jean-Gabriel Young, George T Cantwell, and MEJ Newman.  Bayesian inference of network structure from unreliable data.  Journal of Complex Networks, 8(6):cnaa046, 2020.
  * [5] Mubeena Begum, Gunasekaran Raja, and Mohsen Guizani.  AI-based Sensor Attack Detection and Classification for Autonomous Vehicles in 6G-V2X Environment.  IEEE Transactions on Vehicular Technology, pages 1–10, 2023.
  * [6] Balasubramani S and D John Aravindhar.  Impact of Driver Behavior in VANET Data Transmission.  In 2022 6th International Conference on Intelligent Computing and Control Systems (ICICCS), pages 759–764, May 2022.
  * [7] Adrian Kłusek et al.  Driver Profiling by Using LSTM Networks with Kalman Filtering.  In 2018 IEEE Intelligent Vehicles Symposium (IV), pages 1983–1988, June 2018.
  * [8] Nalam Venkata Abhishek and Mohan Gurusamy.  JaDe: Low Power Jamming Detection Using Machine Learning in Vehicular Networks.  IEEE Wireless Communications Letters, 10(10):2210–2214, Oct 2021.
  * [9] Ali Krayani, Nobel J. William, Lucio Marcenaro, and Carlo Regazzoni.  Jammer Detection in Vehicular V2X Networks.  In 2022 Microwave Mediterranean Symposium (MMS), pages 1–5, 2022.
  * [10] Michele Russo, Maxime Labonne, Alexis Olivereau, and Mohammad Rmayti.  Anomaly detection in vehicle-to-infrastructure communications.  In 2018 IEEE 87th Vehicular Technology Conference (VTC Spring), pages 1–6. IEEE, 2018.
  * [11] Ali Krayani, Gabriele Barabino, Lucio Marcenaro, and Carlo Regazzoni.  Integrated Sensing and Communication for Joint GPS Spoofing and Jamming Detection in Vehicular V2X Networks.  In 2023 IEEE Wireless Communications and Networking Conference (WCNC), pages 1–7, 2023.
  * [12] Le Liang, Shijie Xie, Geoffrey Ye Li, Zhi Ding, and Xingxing Yu.  Graph-based resource sharing in vehicular communication.  IEEE Transactions on Wireless Communications, 17(7):4579–4592, 2018.
  * [13] Chongtao Guo, Le Liang, and Geoffrey Ye Li.  Resource allocation for high-reliability low-latency vehicular communications with packet retransmission.  IEEE Transactions on Vehicular Technology, 68(7):6219–6230, 2019.
  * [14] Ali Krayani, Atm Shafiul Alam, Lucio Marcenaro, Arumugam Nallanathan, and Carlo Regazzoni.  Automatic Jamming Signal Classification in Cognitive UAV Radios.  IEEE Transactions on Vehicular Technology, 71(12):12972–12988, 2022.
  * [15] Ali Krayani, Mohamad Baydoun, Lucio Marcenaro, Atm S. Alam, and Carlo Regazzoni.  Self-Learning Bayesian Generative Models for Jammer Detection in Cognitive-UAV-Radios.  In GLOBECOM 2020 - 2020 IEEE Global Communications Conference, pages 1–7, Dec 2020.
  * [16] Federal Highway Administration (FHWA).  I-80 Freeway Dataset.  www.fhwa.dot.gov/publications/research/operations/06137/, 2006.