# Entanglement criteria for the
bosonic and fermionic induced ensembles

S. Dartois<EMAIL_ADDRESS>Univ. Bordeaux, LaBRI, UMR 5800,
F-33400 Talence, France  ,  Ion Nechita<EMAIL_ADDRESS>Laboratoire
de Physique Théorique, Université de Toulouse, CNRS, UPS, France.  and  A.
Tanasa<EMAIL_ADDRESS>Univ. Bordeaux, LaBRI, UMR 5800, F-33400
Talence, France and DFT, H. Hulubei Nat. Inst. Nucl. Engineering, P.O.Box
MG-6, 077125 Magurele, Romania

###### Abstract.

We introduce the bosonic and fermionic ensembles of density matrices and study
their entanglement. In the fermionic case, we show that random bipartite
fermionic density matrices have non-positive partial transposition, hence they
are typically entangled. The similar analysis in the bosonic case is more
delicate, due to a large positive outlier eigenvalue. We compute the
asymptotic ratio between the size of the environment and the size of the
system Hilbert space for which random bipartite bosonic density matrices fail
the PPT criterion, being thus entangled. We also relate moment computations
for tensor-symmetric random matrices to evaluations of the circuit-counting
and interlace graph polynomials for directed graphs.

Keywords:  random fermionic and bosonic states, random tensors and matrices,
entanglement, circuit polynomial

###### Contents

  1. 1 Introduction
  2. 2 Quantum states, entanglement, symmetry
  3. 3 The fermionic induced ensemble
  4. 4 The bosonic induced ensemble
  5. 5 The partial transposition — statement of main results
  6. 6 The circuit counting graph polynomial
  7. 7 Moments of $W^{\Gamma}$
  8. 8 $t$-channel random matrix and graphical representation
  9. 9 Diagrammatics for the moments of $Q$ and $W_{t}$
  10. 10 Circuit polynomial techniques for $\mathbb{E}\operatorname{Tr}(Q^{p})$
  11. 11 Tensor network evaluation of $\mathbb{E}\operatorname{Tr}(Q^{p})$
  12. 12 Proof of the main result
  13. 13 Conclusion

##  1\. Introduction

In Quantum Mechanics, quantum states are modelled by vectors in a complex
Hilbert space. In the presence of an environment, one uses the formalism of
density matrices to describe quantum states. The dimension of the
corresponding Hilbert space grows exponentially with the number of quantum
particles one wants to describe, hence it is very difficult to precisely
characterize the interesting physical properties of arbitrary density
matrices. One resorts then to understanding the properties of _typical quantum
states_ , in the sense that one would like to find properties which hold with
large probability, for some suitable (natural from a physical perspective)
probability measure on the set of density matrices.

Several such natural probability measures have been considered in the
literature: the induced measure (closely related to the Wishart ensemble from
random matrix theory) [ŻS01], the Bures measure (motivated by statistical
considerations) [Hal98], random matrix product states (motivated by condensed
matter theory) [GdOZ10], and random graph states (which encode a given
entanglement structure) [CNŻ10]. Our first contribution in this work is to
introduce density matrix ensembles for bipartite quantum systems which have a
given _symmetry_. We consider respectively, the symmetric and the anti-
symmetric subspace of a tensor product of Hilbert spaces, and introduce random
density matrices supported on these subspaces. Since the construction is very
similar to that of the induced ensemble [ŻS01], we call these the _bosonic_ ,
respectively the _fermionic induced ensembles_ of mixed quantum states.

We then study the entanglement properties of states from these ensembles,
focusing on the bipartite case (two fermions and two bosons). We analyze the
spectrum of the _partial transposition_ [HHH96, Per96] of these random density
matrices in the large $N$ limit. For the bipartite fermionic ensemble, we find
(see Corollary 3.5) that a typical fermionic mixed state is entangled. This is
due to the presence of a large negative eigenvalue in the spectrum of the
partial transposition of these states. In the bosonic case, the situation is
more complicated, since the symmetry of the state is responsible for a large
positive eigenvalue of the partial transposition. This outlier eigenvalue
makes the study of the spectrum of the partial transposition more delicate.
The asymptotic spectrum of the partial transposition of a random bosonic
density matrix is computed in Theorem 5.1, which we state informally here:

######  Theorem 1.1.

Let $\rho_{N}$ be a $N^{2}\times N^{2}$ element from the bosonic ensemble. The
spectrum of $\rho_{N}$ contains, in the large $N$ limit:

  * •

a large, _outlier_ , positive eigenvalue, of order $1/N$

  * •

a _bulk_ , containing $N^{2}-1$ eigenvalues, of order $1/N^{2}$, whose
empirical distribution follows, when properly normalized, a _shifted
semicircle distribution_.

The bulk of the spectrum is characterized by the shape parameter of the
bosonic induced ensemble, and we observe a _threshold phenomenon_ : the
limiting probability distribution of the bulk eigenvalues has negative support
if and only if the shape parameter $c$ is smaller than 4.

We discuss the implications of this result for the entanglement of typical
bosonic states in Corollary 5.2: when the shape parameter of the bosonic
induced ensemble is smaller than $4$, a typical symmetric mixed state will be
entangled. We also show that the _realignment criterion_ [CW02, Rud03] gives
exactly the same entanglement threshold as the positive partial transposition
(PPT) criterion.

We analyze the asymptotic spectrum of the random matrices we consider (and
their partial transposition) using the moment method. Due to the imposed
symmetry of the models, direct computations are cumbersome; we put forward a
novel connection between such moment computations and a well-known graph
polynomial: the _circuit counting polynomial_ of directed graphs. In
combinatorics and graph theory, this polynomial has received a lot of
attention recently, in relation to the _interlace polynomial_ [ABS04]. We use
this connection in our derivations to bound the asymptotic moments of the
random matrices we are considering, showcasing the importance of the newly
discovered parallel between moment computations and graph polynomials.

The paper is organized as follows. In Section 2 we recall some basic
definitions from quantum information theory and from the theory of symmetric
operators. In Sections 3 and 4 we introduce, respectively, the fermionic and
the bosonic induced ensembles. In Section 5 we state the main result on the
asymptotic spectrum of a random bosonic quantum state. In Section 6 we recall
the main facts on the circuit counting graph polynomial, relating it to
generalized traces of the symmetric projection. This connection is used in
Section 7 to express the moments of the partial transposition of a random
bosonic density matrix to evaluations of the circuit counting graph
polynomial. Sections 8, 9, 10, 11 contain technical results used in Section 12
for the proof of the main result. We conclude in Section 13, where we also lay
out some directions for future research.

##  2\. Quantum states, entanglement, symmetry

In this section we recall the basic notions from quantum information theory
which motivate our work, and then gather some basic mathematical facts
regarding symmetric operators needed later. We start with a brief overview of
quantum states, the notion of quantum entanglement, and the partial
transposition entanglement criterion.

In quantum mechanics, quantum states are modelled by _density matrices_. Each
quantum system comes with a complex Hilbert space $\mathcal{H}$, which we
shall consider here to be finite dimensional:
$\mathcal{H}\approx\mathbb{C}^{N}$; we say that the quantum system has $N$
degrees of freedom. A density matrix is a positive semidefinite operator
acting on $\mathcal{H}$, with unit trace:

$\rho\in\mathcal{M}_{N}(\mathbb{C}),\,\rho\geq 0,\,\operatorname{Tr}\rho=1.$

The Hilbert space corresponding to two quantum systems is obtained by taking
the _tensor product_ of the individual Hilbert spaces:

$\mathcal{H}_{AB}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}.$

A bipartite quantum state $\rho_{AB}$ is said to be _separable_ if it can be
written as a convex mixture of product states:

$\rho_{AB}=\sum_{i=1}^{k}p_{i}\alpha_{i}\otimes\beta_{i},$

for density matrices $\alpha_{i}$, $\beta_{i}$, and probabilities $p_{i}$.
Non-separable quantum states are called _entangled_. Deciding whether a given
bipartite density matrix is separable or entangled is a computationally hard
problem [Gur04]. For this reason, several computationally efficient necessary
or sufficient conditions for entanglement have been developed [HHHH09, Section
VI.B].

In this paper, we shall be concerned with the most famous entanglement
criterion, the _positive partial transpose_ criterion [Per96, HHH96] (we shall
also briefly discuss the realignment criterion in Section 4). The starting
observation is that the partial transposition of a separable state is positive
semidefinite:

$\rho_{AB}^{\Gamma}:=[\operatorname{id}\otimes\operatorname{Tr}](\rho_{AB})=\sum_{i}\alpha_{i}\otimes\beta_{i}^{\top}\geq
0.$

Hence, if the partial transposition of a quantum state $\rho_{AB}$ is _not_
positive semidefinite, the state must be entangled. This criterion is known to
be exact only if the total Hilbert space dimension is at most 6; we refer the
reader to [HHHH09, Section VI.B.1] for further properties.

We shall now discuss several mathematical considerations regarding symmetry in
quantum information theory; more precisely, we shall investigate the symmetric
subspace of a tensor product of Hilbert spaces, and the related operators.
Consider the $r^{\textrm{th}}$ tensor power $\mathcal{H}^{\otimes r}$ of
$\mathcal{H}$. A natural basis of $\mathcal{H}^{\otimes r}$ is given by the
set $\\{e_{i_{1}}\otimes\ldots\otimes
e_{i_{r}}\,:\,i_{1},\ldots,i_{r}\in[N]\\}$ where the $e_{i}$’s form a basis of
$\mathcal{H}$, and $[N]$ denotes the set $\\{1,2,\ldots,N\\}$. Let $\sigma\in
S_{r}$. We define the action of $\sigma$ on $\mathcal{H}^{\otimes r}$ by
extending linearly its action on basis elements, so that,

(1) $\sigma(e_{i_{1}}\otimes\ldots\otimes
e_{i_{r}})=e_{i_{\sigma(1)}}\otimes\ldots\otimes e_{i_{\sigma(r)}}.$

We now define two orthogonal projectors

(2) $P_{s,r}=\frac{1}{r!}\sum_{\sigma\in S_{r}}\sigma\quad\textrm{ and }\quad
P_{a,r}=\frac{1}{r!}\sum_{\sigma\in S_{r}}\varepsilon(\sigma)\sigma$

where $\varepsilon(\sigma)$ is the signature of the permutation $\sigma$.
$P_{s,r}$ are projectors on the symmetric subspace
$\mathrm{Sym}_{r}\mathcal{H}$ and $P_{a,r}$ are projectors on the
antisymmetric subspace (or $r^{\textrm{th}}$ exterior power)
$\wedge^{r}\mathcal{H}$. We define the $r^{\textrm{th}}$ symmetric power
(resp. the $r^{\textrm{th}}$ exterior power) of $\mathcal{H}$ as the subspace
$P_{s,r}(\mathcal{H}^{\otimes r})$ (resp. $P_{a,r}(\mathcal{H}^{\otimes r})$).
We recall

(3) $N_{s}[r]:=\mathrm{dim}\ \mathrm{Sym}_{r}\mathcal{H}=\binom{N+r-1}{r}\
\mathrm{and}\ N_{a}[r]:=\mathrm{dim}\ \wedge^{r}\mathcal{H}=\binom{N}{r}.$

The orthogonal projection on the symmetric subspace can be written as [Har13]

(4) $P_{s,r}=\binom{N+r-1}{r}\int_{\|x\|=1}(xx^{*})^{\otimes r}\mathrm{d}x,$

For example, in the bipartite case, we have

$P_{s,2}=\frac{1}{2}(I+F)\ \textrm{and }P_{a,2}=\frac{1}{2}(I-F)$

where $F$ is the _flip_ (or _swap_) operator $Fx\otimes y=y\otimes x$ while
$I$ is the identity permutation and is represented by a $N^{2}\times N^{2}$
identity matrix. Bipartite tensors can be identified with matrices in a
straightforward way, and, in that case, we have $P_{s,2}A=(A+A^{\top})/2$ and
$P_{a,2}A=(A-A^{\top})/2$. Note that this corresponds to symmetrizing and
anti-symmetrizing matrices; in the complex case, we do not obtain however
self-adjoint or skew-adjoint matrix, since we are not taking complex
conjugates of the entries. Since most of our work will use these projectors
for $r=2$, we use a simplified notation for this case and denote
$P_{s}=P_{s,2}$ and $P_{a}=P_{a,2}$.

We state below a computation which will be useful later (see [Hol18] for a
related problem, stated in representation theoretic language).

######  Lemma 2.1.

The partial trace of the projection on the symmetric subspace can be computed
as follows: (below, $0\leq r\leq k$):

$[\operatorname{id}_{r}\otimes\operatorname{Tr}_{k-r}]P_{s,k}=\frac{N[k]}{N[r]}P_{s,r}.$

###### Proof.

We use the integral representation from eq. (4):

$\displaystyle[\operatorname{id}_{r}\otimes\operatorname{Tr}_{k-r}]P_{s,k}$
$\displaystyle=N[k]\int_{\|x\|=1}[\operatorname{id}_{r}\otimes\operatorname{Tr}_{k-r}]\left((xx^{*})^{\otimes
k}\right)\mathrm{d}x,$ $\displaystyle=N[k]\int_{\|x\|=1}(xx^{*})^{\otimes
r}\mathrm{d}x$ $\displaystyle=\frac{N[k]}{N[r]}P_{s,r}.$

∎

##  3\. The fermionic induced ensemble

In this section we introduce the first random matrix model for quantum states
that we shall study in this paper, that is the matrix model for random
fermionic states. Random quantum states have played an important role in
quantum physics and quantum information theory, being used to study properties
of typical quantum systems and also as a source of examples of states with
interesting properties.

Ensembles of random density matrices have been considered by several authors,
having different physical and mathematical motivations [ŻS01, ŻS03, Nec07,
GdOZ10, ŻPNC11]. The obvious candidate, the _Lebesgue measure_ on the convex
body of density matrices, is part of a one-parameter family of probability
measures, called the _induced ensembles_ [ŻS01]. The general setting is as
follows. Consider a random pure state
$\ket{\psi}\in\mathcal{H}\otimes\mathcal{H}_{E}$, uniformly distributed on the
unit sphere the tensor product between the Hilbert space of interest
$\mathcal{H}$ and an auxiliary Hilbert space $\mathcal{H}_{E}$, called the
_environment_. Define a random density matrix

$\rho:=\operatorname{Tr}_{E}|\psi\rangle\langle\psi|$

by tracing out the environment. In this way, we obtain a random matrix acting
on $\mathcal{H}$; the dimension of the environment, $\dim\mathcal{H}_{E}$, is
a parameter of the model. It was recognized in [SŻ04, Nec07] that the induced
measure can be understood as a _normalized Wishart measure_ :

$\rho=\frac{Y}{\operatorname{Tr}Y},$

where $Y=GG^{*}$, for $G:\mathcal{H}\to\mathcal{H}_{E}$ a random Gaussian
matrix (an element of the _Ginibre ensemble_).

We start by introducing the random matrix model for fermionic quantum states
that we call the _fermionic induced ensemble_. We then show the entanglement
results for this model of random fermionic states.

In practice, we study an antisymmetric random Gaussian tensor
$G_{a}\in\wedge^{r}\mathcal{H}\otimes\mathcal{H}_{E}$. We identify the space
of antisymmetric tensors

(5)
$\wedge^{r}\mathcal{H}\cong\textrm{Ran}(P_{a,r})\subseteq\mathcal{H}^{\otimes
r},$

where $P_{a,r}$ is defined in Section 2 to be the projector on the
antisymmetric subspace of $\mathcal{H}^{\otimes r}$. We use this
identification to construct the antisymmetric random Gaussian tensor from a
complex random standard Gaussian tensor $G\in\mathcal{H}^{\otimes
r}\otimes\mathcal{H}_{E}$, that is a tensor whose entries are standard normal
random variables. We will be interested in the marginals obtained by tracing
out the environmental degrees of freedom. We study the random matrix
$A:=\operatorname{Tr}_{E}(G_{a}G_{a}^{*})$. Note that, using the
identification from (5), we can also write

(6) $A=P_{a,r}\operatorname{Tr}_{E}(GG^{*})P_{a,r},$

where $G$ is the complex standard Gaussian tensor. In the language of quantum
information theory, normalizing matrices $A$ yields random density matrices
having an _antisymmetry_ constraint, reflecting their fermionic character on
the $r$ copies of $\mathcal{H}$.

######  Definition 3.1.

The _fermionic induced ensemble_ of parameters $(N,M,r)$ is the set of random
bipartite density matrices of size $N^{r}$ defined by

(7) $\rho^{(A)}=\frac{A}{\operatorname{Tr}A},$

where $A$ is an antisymmetrized Wishart random matrix from Eq. (6). Matrices
from the fermionic induced ensemble are supported on the antisymmetric
subspace of $(\mathbb{C}^{N})^{\otimes r}$.

Let us now comment on our choice of model. As already detailed in the
introduction, our main motivation is that these types of states are more
realistic than any (possibly non-symmetric) random states, since quite often
in quantum experiments we consider states of electrons or say
${}^{87}\mathrm{Rb}$ (which are fermions), or states of photons or
${}^{4}\mathrm{He}$ atoms (which are bosons, whose ensemble we introduce later
in section 4). The corresponding degrees of freedom must have
antisymmetric/symmetric wave functions, since they are fermions/bosons, while
the environment has no reason to have any kind of symmetries. This explains
why we impose the symmetry condition on the tensor product
$\mathcal{H}^{\otimes r}$, while no symmetry is imposed on the tensor product
with the environment system $\mathcal{H}_{E}$.

In the rest of this section, we set $r=2$. This is for technical reasons. If
in the fermionic case, we could rather easily extend our results to general
$r$, it does not appear feasible with our current techniques to study the
bosonic case for general $r$. In order to stay coherent all along the paper we
decided to present our results in the $r=2$ case for both fermionic and
bosonic induced ensembles.

We set the dimension of the Hilbert spaces as follows

$\displaystyle\dim\mathcal{H}$ $\displaystyle=:N$
$\displaystyle\dim\mathcal{H}_{E}$ $\displaystyle=:M.$

We shall consider an asymptotic regime where both $N,M\to\infty$ in such way
that $M=cN^{2}$ for some fixed constant $c\in(0,\infty)$. This choice of
scaling is the standard one needed for the convergence of the standard Wishart
random matrices to the Marc̆enko-Pastur distribution:

(8)
$\mathrm{d}\mathrm{MP}_{c}=\max(1-c,0)\delta_{0}+\frac{\sqrt{(b-x)(x-a)}}{2\pi
x}\;\mathbf{1}_{[a,b]}(x)\,\mathrm{d}x,$

with $a=(1-\sqrt{c})^{2}$ and $b=(1+\sqrt{c})^{2}$.

The limiting eigenvalue distribution of the properly normalized random matrix
$A$ is the same as that of a non-symmetrized Wishart matrix.

######  Proposition 3.2.

The random matrix $A$ converges in moments, as $N\to\infty$, $M\sim cN^{2}$,
to a Marc̆enko-Pastur distribution of parameter $2c$:

$\forall p\geq
1,\qquad\lim_{N\to\infty}\mathbb{E}\frac{1}{N_{a}[2]}\operatorname{Tr}\left[\left(\frac{A}{N_{a}[2]}\right)^{p}\right]=\int
x^{p}\mathrm{d}\mathrm{MP}_{2c}(x).$

###### Proof.

After restricting it to its support ($\mathrm{Sym}_{2}(\mathbb{C}^{N})$), the
random matrix $W$ is indeed a Wishart matrix of parameters $(N_{a}[2],M)$. The
result follows from the standard theorem of Marc̆enko and Pastur, see [MP67],
[BS10, Theorem 3.6], or [DLN20, Proposition 2.4]. The parameter of the
limiting distribution is given by

$\lim_{N\to\infty}\frac{M}{N_{a}[2]}=2c.$

∎

In order to understand the entanglement properties of the fermionic random
states we just introduced, we proceed by first studying the average of their
partial transposition. We show that $A^{\Gamma}$ has, on average, a negative
eigenvalue, on a larger scale than the rest of the spectrum, for all $c>0$ (in
the asymptotic regime where $M\sim cN^{2}$). Before we commence, recall that
the _maximally entangled state_ is the unit norm vector

$\ket{\Omega}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{ii}\in\mathbb{C}^{N}\otimes\mathbb{C}^{N},$

where $\\{\ket{i}\\}_{i\in[N]}$ is the canonical orthonormal basis of
$\mathbb{C}^{N}$. We also write $\omega=|\Omega\rangle\langle\Omega|$ for the
corresponding density matrix.

######  Proposition 3.3.

The average of the partial transpose $A^{\Gamma}$ is given by

(9) $\mathbb{E}A^{\Gamma}=MP_{a,2}=\frac{M}{2}I_{N^{2}}-\frac{MN}{2}\omega.$

The eigenvalues of $\mathbb{E}A^{\Gamma}$ are as follows:

(10)
$\mathrm{spec}(\mathbb{E}A^{\Gamma})=\begin{cases}-\frac{M(N-1)}{2}&\textrm{
with multiplicity }1\\\ \frac{M}{2}&\textrm{with multiplicity
}N^{2}-1.\end{cases}$

###### Proof.

The proof is a direct computation. Using Wick-Isserlis theorem we compute
$\mathbb{E}A=\frac{M}{2}(I_{N^{2}}-F)$, where $F$ is the flip (or swap)
operator

$F=\sum_{i,j=1}^{N}|ij\rangle\langle ji|.$

Applying the partial transpose to the flip operator, we deduce that

$F^{\Gamma}=\sum_{i,j=1}^{N}(|ij\rangle\langle
ji|)^{\Gamma}=\sum_{i,j=1}^{N}|ii\rangle\langle
jj|=N|\Omega\rangle\langle\Omega|=N\omega.$

Writing $P_{a,2}=(I-F)/2$ yields the formula for $\mathbb{E}A^{\Gamma}$. The
exact form of the spectrum follows from the fact that $\omega$ is a rank-one
projection. ∎

We take away from the above proof that the maximally entangled vector
$\lvert\Omega\rangle$ is an eigenvector of $\mathbb{E}A^{\Gamma}$ with a
negative eigenvalue for all $N\geq 2,c>0$, asymptotically of the form
$-\frac{c}{2}N^{3}+O(N^{2})$. The next theorem shows that, the same holds
asymptotically for the random matrix $A^{\Gamma}$.

######  Theorem 3.4.

The maximally entangled vector $\lvert\Omega\rangle$ is an approximate
eigenvector of the random matrix $N^{-3}A^{\Gamma}$ with corresponding
approximate eigenvalue $-c/2$: if

(11)
$\delta_{N}:=\|N^{-3}A^{\Gamma}\lvert\Omega\rangle+c/2\lvert\Omega\rangle\|^{2},$

then $\lim_{N\rightarrow\infty}\mathbb{E}\delta_{N}=0$. In particular, the
random variable $\delta_{N}$ converges in probability towards zero as
$N\rightarrow\infty$.

###### Proof.

We start with a more explicit formula of $\delta_{N}$:

(12)
$\delta_{N}=N^{-6}\langle\Omega\rvert(A^{\Gamma})^{2}\lvert\Omega\rangle+cN^{-3}\langle\Omega\rvert
A^{\Gamma}\lvert\Omega\rangle+\frac{c^{2}}{4}.$

Computing the average, we obtain, thanks to Proposition 3.3,

(13)
$\mathbb{E}\delta_{N}=N^{-6}\mathbb{E}\langle\Omega\rvert(A^{\Gamma})^{2}\lvert\Omega\rangle-\frac{c^{2}}{4}+o(1).$

We then use Wick-Isserlis theorem to compute
$\mathbb{E}\langle\Omega\rvert(A^{\Gamma})^{2}\lvert\Omega\rangle$, and we
find an expression involving the partial trace $\operatorname{Tr}_{1}P_{a,2}$,

$\displaystyle\mathbb{E}\langle\Omega\rvert(A^{\Gamma})^{2}\lvert\Omega\rangle$
$\displaystyle=M^{2}\langle\Omega|(P_{a,2}^{\Gamma})^{2}|\Omega\rangle+\frac{M}{N}\operatorname{Tr}\left[(\operatorname{Tr}_{1}P_{a,2})^{2}\right]$
$\displaystyle=\frac{M^{2}}{N}\operatorname{Tr}\left[(\operatorname{Tr}_{1}P_{a,2})^{2}\right]+\frac{M}{N}\operatorname{Tr}\left[(\operatorname{Tr}_{1}P_{a,2})^{2}\right]$

We compute $\operatorname{Tr}_{1}(P_{a,2})=\frac{1}{2}(N-1)I_{N}$ from which
we deduce

$\displaystyle\mathbb{E}\langle\Omega\rvert(A^{\Gamma})^{2}\lvert\Omega\rangle$
$\displaystyle=\left[\frac{c^{2}}{4}N^{6}+\frac{c}{4}N^{4}\right]\left(1-\frac{2}{N}+\frac{1}{N^{2}}\right)=\frac{c^{2}}{4}N^{6}+O(N^{5}),$

showing that $\mathbb{E}\delta_{N}\to 0$. We then use Markov’s inequality to
bound the probability that $\delta_{N}$ is positive. Indeed, for all $a>0$, we
have

(14) $\mathbb{P}(\delta_{N}\geq a)\leq\frac{1}{a}\mathbb{E}\delta_{N}\to 0,$

concluding the proof. ∎

These two results are here sufficient to conclude that the random matrix
$A^{\Gamma}$ is not positive semi-definite, as indeed it has a negative
eigenvalue at scale $N^{3}$. As a consequence, the corresponding normalized
random state $\rho^{(A)}$ is entangled. More precisely, we have the corollary

######  Corollary 3.5.

Consider a sequence $(\rho^{(A)}_{N})_{N\geq 2}$ of density matrices from the
fermionic induced ensemble of parameters $(N,M,2)$, with $M=cN^{2}$. Then, for
all $c>0$,

(15) $\lim_{N\rightarrow\infty}\mathbb{P}[(\rho^{(A)}_{N})^{\Gamma}\geq 0]=0,$

hence,

(16) $\lim_{N\rightarrow\infty}\mathbb{P}[\rho^{(A)}_{N}\textrm{ is
entangled}\,]=1.$

This is in contrast with the result of Aubrun [Aub12] who considered random
distinguishable states. In his framework, it was shown that the analog of the
above statement is valid if and only if the environmental Hilbert space is
large enough (i.e. for $c>4$). What we show here is that for states of
fermions, typical states are always entangled, this does not depend on the
size of the environment as long as it is commensurable with the two particle
Hilbert space $\wedge^{2}\mathcal{H}$ size. In the next part we study the same
problem for random states of bosons. Bosonic states prove to be considerably
more difficult to study and the rest of this paper is devoted to them.

##  4\. The bosonic induced ensemble

In this section, we shall consider a setting analog to the one of the previous
section 3, only now the system Hilbert space has a _symmetric tensor product
structure_. In other words, we shall start from a random Gaussian tensor
$G_{s}\in\mathrm{Sym}_{r}(\mathcal{H})\otimes\mathcal{H}_{E}$. In practice, we
shall use the identification

(17)
$\mathrm{Sym}_{r}(\mathcal{H})\cong\operatorname{Ran}(P_{s,r})\subseteq\mathcal{H}^{\otimes
r},$

where we recall that $P_{s}$ is the orthogonal projection on the symmetric
subspace of $\mathcal{H}\otimes\mathcal{H}$. The tensor factor
$\mathcal{H}_{E}$ is again the environment Hilbert space. These random vectors
$G_{s}$ can be seen as (un-normalized) bosonic states on two copies of
$\mathcal{H}$ in contact with an environment $\mathcal{H}_{E}$. Similarly, we
will be interested in the marginals obtained by tracing out the environmental
degrees of freedom. We study the corresponding random matrix
$W=\operatorname{Tr}_{E}(G_{s}G_{s}^{*})$. Again, using the identification
from (17), we can also write

(18) $W=P_{s,r}\operatorname{Tr}_{E}(GG^{*})P_{s,r},$

where $G$ is a random (complex) Gaussian vector in the Hilbert space
$\mathcal{H}^{\otimes r}\otimes\mathcal{H}_{E}$; see Figure 1 for a pictorial
representation of the random matrix $W$ at $r=2$. Normalizing matrices $W$
yields random density matrices which have a _symmetry_ constraint.

######  Definition 4.1.

The _bosonic induced ensemble_ of parameters $(N,M,r)$ is the set of random
bipartite density matrices of size $N^{r}$ defined by

$\rho=\frac{W}{\operatorname{Tr}W},$

where $W$ is the symmetrized Wishart random matrix from (18). Matrices from
the bosonic induced ensemble are supported on the symmetric subspace of
$(\mathbb{C}^{N})^{\otimes r}$.

Figure 1. A graphical representation of the symmetrized Wishart matrix $W$
that we shall study in this work. The grey boxes surrounding the matrix
$GG^{*}$ are symmetrizers.

Mimicking the construction of Section 3, we set the dimension of the Hilbert
spaces as follows

$\displaystyle\dim\mathcal{H}$ $\displaystyle=:N$
$\displaystyle\dim\mathcal{H}_{E}$ $\displaystyle:=M.$

Doing so, the limiting eigenvalue distribution of the properly normalized
random matrix $W$ is the same as that of a non-symmetrized Wishart and also an
anti-symmetrized matrix.

######  Proposition 4.2.

The random matrix $W$ converges in moments, as $N\to\infty$, $M\sim cN^{2}$,
to a Marc̆enko-Pastur distribution of parameter $2c$:

$\forall p\geq
1,\qquad\lim_{N\to\infty}\mathbb{E}\frac{1}{N_{s}[2]}\operatorname{Tr}\left[\left(\frac{W}{N_{s}[2]}\right)^{p}\right]=\int
x^{p}\mathrm{d}\mathrm{MP}_{2c}(x).$

###### Proof.

After restricting it to its support ($\mathrm{Sym}_{2}(\mathbb{C}^{N})$), the
random matrix $W$ is indeed a Wishart matrix of parameters $(N_{s}[2],M)$. The
result follows from the standard theorem of Marc̆enko and Pastur, see [MP67],
[BS10, Theorem 3.6], or [DLN20, Proposition 2.4]. The parameter of the
limiting distribution is given by

$\lim_{N\to\infty}\frac{M}{N_{s}[2]}=2c.$

∎

######  Remark 4.3.

The exact moments of order $p$ of the random matrix $W$ can be computed in
different manners, either as a sum over permutations, or as a sum over
bipartite combinatorial maps with one black vertex and $p$ edges
$\mathcal{M}=(\alpha,\gamma)$ where $\alpha$ is a permutation of $S_{p}$ and
$\gamma=(1\,2\,3\,\cdots\,p)$ is the full cycle:

$\displaystyle\mathbb{E}\operatorname{Tr}\left[\left(\frac{W}{N_{s}[2]}\right)^{p}\right]$
$\displaystyle=\sum_{\alpha\in
S_{p}}N_{s}[2]^{\\#(\gamma^{-1}\alpha)}M^{\\#\alpha}$
$\displaystyle=\sum_{\mathcal{M}}N_{s}[2]^{F(\mathcal{M})}M^{V_{\circ}(M)},$

where $F(\mathcal{M})$ is the number of faces of $\mathcal{M}$ and
$V_{\circ}(\mathcal{M})$ its number of white vertices. We refer the reader to
[DLN20, Section 2] for the dictionary between these two approaches.

######  Remark 4.4.

The asymptotic moments from Proposition 4.2 correspond to sums over _geodesic
permutations_ $\alpha$, or, equivalently, over _planar bipartite maps, or
dessins d’enfants,_ $\mathcal{M}$ with one black vertex. The limiting moments
are values of the Narayana polynomial:

$\mathrm{Nar}(p,2c)=\int
x^{p}\mathrm{d}\mathrm{MP}_{2c}(x)=\sum_{k=0}^{p-1}\frac{1}{(k+1)(2c)^{k}}\binom{p-1}{k}\binom{p}{k}.$

##  5\. The partial transposition — statement of main results

In this section we state the main result of this work, the characterization of
the limiting spectrum of the _partial transposition_ of the random matrix $W$,
in the large $N$ limit. Recall that $W$ is a $N^{2}\times N^{2}$ random
matrix, supported on the symmetric subspace
$\mathrm{Sym}_{2}(\mathbb{C}^{N})\subseteq\mathbb{C}^{N}\otimes\mathbb{C}^{N}$.
Its partial transpose is defined as the action of the transposition operator
on the second tensor factor:

(19) $W^{\Gamma}=[\mathrm{id}\otimes\mathrm{transp}](W).$

The next theorem is the main result of our paper, providing a description of
the spectrum of $W^{\Gamma}$ in the large $N$ limit. Interestingly,
$W^{\Gamma}$ has a large eigenvalue of order $N^{3}$, and $N^{2}-1$
eigenvalues of order $N^{2}$.

######  Theorem 5.1.

Let $\lambda_{1}\geq\lambda_{2}\geq\ldots\lambda_{N^{2}}$ be the eigenvalues
of the random matrix $W^{\Gamma}$ from Eq. (19). Then, in the asymptotic
regime $N\to\infty$ and $M\sim cN^{2}$,

  * •

$N^{-3}\lambda_{1}\rightarrow\frac{c}{2}$ in probability

  * •

The empirical distribution $\sum_{i=2}^{N^{2}}\delta_{N^{-2}\lambda_{i}}$
converges in moments to a semi-circular distribution of mean $c/2$ and
variance $c/4$.

We recall here that the _semicircular distribution_ with average $m$ and
variance $\sigma^{2}$ is given by

$\mathrm{d}\mathrm{SC}_{m,\sigma}=\frac{1}{2\pi\sigma}\sqrt{4\sigma^{2}-(x-m)^{2}}\mathbf{1}_{[-2\sigma+m,2\sigma+m]}(x)\,\mathrm{d}x.$

Let us record here an important corollary for the theory of quantum
information. The _partial transposition criterion_ [Per96, HHH96] states that
a _separable_ (i.e. non-entangled) density matrix $\rho$ has a positive
semidefinite partial transposition: $\rho^{\Gamma}\geq 0$. This fact is mostly
used in the reverse direction, as an entanglement criterion: given a density
matrix $\sigma$ having a partial transposition with at least one negative
eigenvalue ($\sigma^{\Gamma}\ngeq 0$) is entangled. We have thus the following
result, regarding the bosonic induced ensemble defined in Definition 4.1.

######  Corollary 5.2.

Consider a sequence $(\rho_{N})_{N\geq 1}$ of density matrices from the
bosonic induced ensemble of parameters $(N,M,2)$, with $M\sim cN^{2}$ as
$N\to\infty$. If $c<4$, then

$\lim_{N\to\infty}\mathbb{P}[\rho_{N}^{\Gamma}\geq 0]=0.$

In particular,

$\lim_{N\to\infty}\mathbb{P}[\rho_{N}\text{ is entangled }]=1.$

######  Remark 5.3.

Note that the other regime, where $c>4$, is much more involved, since the
absence of negative support of the limiting measure of $\rho_{N}^{\Gamma}$
does not exclude the absence of negative outliers. In the standard (non-
symmetrized) Wishart case, Aubrun [Aub12] computed large moments of the
partial transposition in order to achieve this conclusion. We leave these
considerations open. Indeed, in our case, considerations on large moments
would not allow us to conclude. This is due to our use of the Cauchy
interlacing theorem. The study of large moments of the compressed matrix $Q$,
introduced later, does not allow us to retain enough information on the
behaviour of outliers of $\rho_{N}^{\Gamma}$; at the same time, in this work,
we cannot study them directly because of the large eigenvalue of $W^{\Gamma}$,
as will be seen later.

The proof of the main theorem is rather involved, and we shall proceed in
several steps in the following sections, culminating with the final proof
given in Section 12. We start, in the spirit of Section 3, with an analysis of
the average state and of the outlier eigenvalue. However, in the bosonic case,
this will not be sufficient to make a statement about entanglement. Hence the
need of a more detailed study of the spectrum, provided by our main result,
Theorem 5.1.

Following the approach of Section 3, we compute the average of the partially
transposed matrix $W^{\Gamma}$, and show that already on average, there is an
outlier, _positive_ , eigenvalue on a scale larger than the rest of the
spectrum. Recall that the _maximally entangled state_ is the unit norm vector

$\ket{\Omega}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{ii}\in\mathbb{C}^{N}\otimes\mathbb{C}^{N},$

where $\\{\ket{i}\\}_{i\in[N]}$ is the canonical orthonormal basis of
$\mathbb{C}^{N}$.

######  Proposition 5.4.

The average of the partial transpose $W^{\Gamma}$ is given by

$\mathbb{E}W^{\Gamma}=MP_{s}^{\Gamma}=\frac{M}{2}I_{N^{2}}+\frac{MN}{2}|\Omega\rangle\langle\Omega|.$

The eigenvalues of $\mathbb{E}W^{\Gamma}$ are as follows

$\operatorname{spec}(\mathbb{E}W^{\Gamma})=\begin{cases}\frac{M(N+1)}{2}\qquad&\text{
with multiplicity $1$}\\\ \frac{M}{2}\qquad&\text{ with multiplicity
$N^{2}-1$.}\end{cases}$

###### Proof.

The proof is a direct computation. Using Wick-Isserlis theorem we compute
$\mathbb{E}W=\frac{M}{2}(I_{N^{2}}+F)$. Applying the partial transposition on
the flip operator

$F=\sum_{i,j=1}^{N}|ij\rangle\langle ji|,$

we deduce that

$F^{\Gamma}=\sum_{i,j=1}^{N}(|ij\rangle\langle
ji|)^{\Gamma}=\sum_{i,j=1}^{N}|ii\rangle\langle
jj|=N|\Omega\rangle\langle\Omega|.$

Writing $P_{s}=(I+F)/2$ yields the formula for $\mathbb{E}W^{\Gamma}$. The
spectrum follows from the fact that $|\Omega\rangle\langle\Omega|$ is a rank-
one projection. ∎

We note that $\lvert\Omega\rangle$ is an eigenvector of $P_{s}^{\Gamma}$ with
eigenvalue $(N+1)/2$. Hence, the average $\mathbb{E}W^{\Gamma}$ has an
eigenvalue $M(N+1)/2\sim cN^{3}/2$ with eigenvector $\ket{\Omega}$. We show
next that the same holds, asymptotically, for $W^{\Gamma}$.

######  Theorem 5.5.

The maximally entangled vector $\ket{\Omega}$ is an approximate eigenvector of
the random matrix $N^{-3}W^{\Gamma}$ with corresponding approximate eigenvalue
$c/2$: if

$\varepsilon_{N}:=\|N^{-3}W^{\Gamma}\ket{\Omega}-c/2\ket{\Omega}\|^{2},$

then $\lim_{N\to\infty}\mathbb{E}\varepsilon_{N}=0$. In particular, the random
variable $\varepsilon_{N}$ converges in probability towards zero as
$N\to\infty$.

###### Proof.

Manipulating the expression of $\varepsilon_{N}$ we have

$\varepsilon_{N}=N^{-6}\langle\Omega|(W^{\Gamma})^{2}|\Omega\rangle-
cN^{-3}\langle\Omega|W^{\Gamma}|\Omega\rangle+c^{2}/4.$

Computing the average and using Proposition 5.4 we have

$\mathbb{E}\varepsilon_{N}=N^{-6}\mathbb{E}\langle\Omega|(W^{\Gamma})^{2}|\Omega\rangle-c^{2}/4+o(1).$

We compute $\mathbb{E}\langle\Omega|(W^{\Gamma})^{2}|\Omega\rangle$ using the
Wick-Isserlis theorem and find, using Lemma 2.1,

$\displaystyle\mathbb{E}\langle\Omega|(W^{\Gamma})^{2}|\Omega\rangle$
$\displaystyle=M^{2}\langle\Omega|(P_{s}^{\Gamma})^{2}|\Omega\rangle+\frac{M}{N}\operatorname{Tr}\left[(\operatorname{Tr}_{1}P_{s})^{2}\right]$
$\displaystyle=\frac{M^{2}}{N}\operatorname{Tr}\left[(\operatorname{Tr}_{1}P_{s})^{2}\right]+\frac{M}{N}\operatorname{Tr}\left[(\operatorname{Tr}_{1}P_{s})^{2}\right]$
$\displaystyle=\left[N^{6}\frac{c^{2}}{4}+N^{4}\frac{c}{4}\right]\left(1+\frac{2}{N}+\frac{1}{N^{2}}\right)=\frac{c^{2}}{4}N^{6}+O(N^{5}),$

establishing the first claim. The convergence in probability follows by
Markov’s inequality. ∎

Let us finish this section by addressing another important entanglement
criterion: _realignment_ [CW02, Rud03]. The realignment criterion states that
any separable bipartite density matrix $\rho$ satisfies the following
inequality:

$\|\rho^{R}\|_{1}\leq 1,$

where $\|\cdot\|_{1}$ denotes the Schatten 1-norm or the nuclear norm (i.e.
the sum of the singular values of a matrix), while $\rho^{R}$ denotes the
realignment (or the reshuffling, see [BZ06, Chapter 10.2]) of a matrix $\rho$,
defined algebraically by

$\langle ij|\rho^{R}|kl\rangle=\langle ik|\rho|jl\rangle,$

or graphically in Figure 2.

Figure 2. The realignment of a bipartite matrix $\rho$.

The key observation here is that the realignment and the partial transposition
of a matrix are related as follows:

$\rho^{R}=\left(\rho^{\Gamma}\right)\cdot
F\quad\implies\quad\|\rho^{R}\|_{1}=\|\rho^{\Gamma}\|_{1}.$

In the case of a (non-normalized) random matrix $W$ from the induced bosonic
ensemble, entanglement follows from the inequality
$\|W^{\Gamma}\|_{1}>\operatorname{Tr}W$. Note that, asymptotically as
$N\to\infty$ and $M\sim cN^{2}$, the trace behaves as
$\operatorname{Tr}W\sim(c/2)N^{4}$, while Theorem 5.1 states that

$\|W^{\Gamma}\|_{1}\sim(c/2)N^{3}+N^{4}\int|x|\,\mathrm{d}\mathrm{SC}_{c/2,\sqrt{c}/2}(x).$

We note that the large outlier eigenvalue does not contribute asymptotically
to the value of the Schatten 1-norm. Hence, in order for the realignment
criterion to detect entanglement in the random bosonic matrix $W$, the
following inequality needs to hold:

$\int|x|\,\mathrm{d}\mathrm{SC}_{c/2,\sqrt{c}/2}(x)>\frac{c}{2}=\int
x\,\mathrm{d}\mathrm{SC}_{c/2,\sqrt{c}/2}(x).$

But the above happens precisely when the probability measure
$\mathrm{SC}_{c/2,\sqrt{c}/2}$ has negative support, which is also the
condition for the partial transpose criterion to detect the entanglement of
$W$. Hence, we conclude that, in the case of the induced bosonic ensemble, the
partial transpose and the realignment criteria have _identical thresholds_
($c_{0}=4$). This situation is to be contrasted with the case of the usual
induced ensemble, where the threshold for the partial transpose criterion
($c^{\Gamma}_{0}=4$ [Aub12]) is strictly larger than the threshold for the
realignment criterion ($c_{0}^{R}=(8/3\pi)^{2}\approx 0.72$ [AN12]), proving
that, in some sense, the former criterion is stronger than the latter.

##  6\. The circuit counting graph polynomial

In this subsection we discuss some basic facts about the _circuit counting
polynomial_ of directed graphs. Of special importance to us is Theorem 6.10,
relating the generalized bipartite trace of the symmetric projector
$P_{s}=P_{s,2}$ to the circuit counting polynomial $J$ of an associated graph.

We start by recalling the definition of a digraph.

######  Definition 6.1.

A digraph is a pair $G=(V,A)$ where $V$ is a set of vertices and $A$ is a
multi-set of ordered pairs of (possibly non distinct) elements of $V$. The
fact that $A$ is a multi-set allows for multiple edges.

In this paper, we will denote when needed the elements of $A$ as
$\\{i\rightarrow j\\}$ for $i,j\in V$ to better picture the ordering. Possibly
non-distinct elements in the pairs means the digraphs can have loops.
Note that incoming edges of a vertex $i$ are the edges of the form
$\\{a\rightarrow i\\}$, while outgoing edges are the edges of the form
$\\{i\rightarrow a\\}$. A loop edge is both incoming and outgoing for the
vertex it is adjacent to. A digraph is said to be $k$-in/$k$-out if each
vertex has $k$ incoming edges and $k$ outgoing edges.

We introduce the circuit counting polynomial $J$ of a digraph $G$ as

(20) $J(G;x)=\sum_{k\geq 1}j_{k}x^{k},$

with $j_{k}=|\\{\textrm{covers of }G\textrm{ in }k\textrm{ cycles}\\}|$. Note
that the cycles must follow the edge orientations. The circuit counting
polynomial possesses several properties useful to us. As we will be interested
only in the case of digraphs that are $2$-in/$2$-out, we state these
properties only in the case of $2$-in/$2$-out digraphs. The circuit counting
polynomial has the following trivial multiplicativity property.

######  Lemma 6.2.

Let $G$ be a 2-in/2-out digraph having connected components
$G_{1},G_{2},\ldots,G_{k}$ (which are connected 2-in/2-out digraphs). Then

$J(G;x)=\prod_{i=1}^{k}J(G_{i};x).$

Next, we consider skein relations. The reason is that a cycle cover of $G$
corresponds to a choice of state for each vertex of $G$ (as defined in [EM04,
Definition 2.1] and below; the case we are interested in is the Eulerian
case). These skein relations read graphically

(21) $J\left(\raisebox{-45.5244pt}{\includegraphics[scale={0.4}]{J-before-
move.pdf}};x\right)\ =\
J\left(\raisebox{-45.5244pt}{\includegraphics[scale={0.4}]{J-after-
move-1.pdf}};x\right)\ +\
J\left(\raisebox{-45.5244pt}{\includegraphics[scale={0.4}]{J-after-
move-2.pdf}};x\right).$

Additionally, when evaluated on connected graphs, it is of maximal degree for
graphs that have only cut vertices, where we recall that a cut vertex is a
vertex whose removal disconnects the graph. This fact is a simple extension of
[LV83], that we show here. We give a proof below, using the graphical skein
relations.

######  Theorem 6.3.

Let $G$ be a $2$-in/$2$-out connected digraph and $v$ be a cut vertex. Denote
by $G\backslash v$ the _connected_ digraph obtained from $G$ by erasing the
vertex $v$ and reconnecting the edges in an orientation preserving way. Then

(22) $J(G;x)=(x+1)J(G\backslash v;x).$

###### Proof.

We have using the skein relations described above at the cut vertex $v$

(23) $J\left(\raisebox{-28.45274pt}{\includegraphics[scale={0.4}]{J-cut-
vertex.pdf}};x\right)\ =\
J\left(\raisebox{-28.45274pt}{\includegraphics[scale={0.4}]{J-cut-vertex-
move1.pdf}};x\right)\ +\
J\left(\raisebox{-28.45274pt}{\includegraphics[scale={0.4}]{J-cut-vertex-
move2.pdf}};x\right)$

where the gray areas represent the rest of the graph. Since $v$ is a cut
vertex these gray areas are connected only through $v$. The blue dashed lines
inside the gray parts show how a cycle that would pass through the edges
entering and exiting the gray parts would behave. Notice in particular that
due to the $2$-in/$2$-out constraint, a cycle going through one edge entering
a gray part has to go out of this gray part through the edge exiting it. Now
notice that one of the choice of state for $v$ (the leftmost term of the right
hand side) leads to a disconnected graph that has one more cycle than the
other choice of state that leads to a connected graph (rightmost term of the
right hand side). Thus we have, graphically,

(24) $J\left(\raisebox{-28.45274pt}{\includegraphics[scale={0.4}]{J-cut-
vertex-
move1.pdf}};x\right)=xJ\left(\raisebox{-28.45274pt}{\includegraphics[scale={0.4}]{J-cut-
vertex-move2.pdf}};x\right).$

Therefore,

(25) $J\left(\raisebox{-28.45274pt}{\includegraphics[scale={0.4}]{J-cut-
vertex.pdf}};x\right)\ =\
(x+1)J\left(\raisebox{-28.45274pt}{\includegraphics[scale={0.4}]{J-cut-vertex-
move2.pdf}};x\right).$

∎

Let $G$ be a $2$-in/$2$-out digraph with $p$ vertices such that if every
vertex is a cut vertex. Then repeated use of equation (22) leads to

(26) $J(G;x)=(x+1)^{p}J(L;x)=x(x+1)^{p},$

where $L$ is the trivial digraph made of an edge looping on itself with no
vertices. We later show that digraphs with only cut vertices are the ones
maximizing the degree of the circuit counting polynomial.

In the case of $2$-in/$2$-digraphs, the circuit counting polynomial is related
to the _interlace polynomial_ [ABS04], defined here according to [BH11]. For a
(non oriented) graph $H=(V,E)$, the interlace polynomial $q(H;x)$ is

(27) $q(H;x)=\sum_{\begin{subarray}{c}S\subseteq
V\end{subarray}}(x-1)^{\textrm{dim }\textrm{ker }A_{H[S]}},$

where $H[S]$ is the sub-graph induced by the subset $S$ of vertices of $H$ and
$A_{H[S]}$ denotes the adjacency matrix of $H[S]$. Note that $\textrm{dim
}\textrm{ker }A_{H[S]}$ is often called the nullity of $H[S]$. Note also that
the degree of the interlace polynomial of a graph is bounded by its number of
vertices:

(28) $\deg q(H;x)\leq|H|.$

We use repeatedly the following property. For connected 2-in/2-out digraphs,
the circuit polynomial is related to the interlace polynomial $q(x)$ of a
corresponding interlace graph via the following theorem.

######  Theorem 6.4 ([ABS04, Theorem 24]).

Consider $G$ an Eulerian $2$-in/$2$-out digraph, $C$ any Euler circuit of $G$
and $H=H(C)$ its interlace graph. Then

(29) $J(G;x)=xq(H;x+1)=xm(G;x+1)$

with $q$ the interlace polynomial, and $m$ the Martin polynomial.

Let us recall a bit of terminology. An Euler circuit $C$ of a connected
digraph $G$ is a circuit passing through every edge of $G$ exactly once. For a
$2$-in/$2$-out digraph $G$ the interlace graph $H(C)$ of one of its Euler
circuits $C$ is the (non oriented!) graph defined as having the same vertex
set than $G$ and having an edge between vertices $a$ and $b$ if those are
interlaced in $C$. That is, they appear in the following way $C=\ldots a\ldots
b\ldots a\ldots b\ldots$, see [ABS04] for these definitions. It is easy to
deduce that connected $2$-in/$2$-out digraphs $G$ having an Euler circuit with
no interlace pairs have an interlace graph $H(C)$ maximizing the degree of
$q(H(C);x)$ since $H(C)$ is empty (thus $\textrm{dim }\textrm{ker
}A_{H(C)[S]}=|S|$ for all $S\subseteq V$).

In [ABS04] the following result is established.

######  Lemma 6.5.

If a $2$-in/$2$-out digraph $G$ has an Euler circuit $C$ with no interlace
pairs, then $G$ has only one Euler circuit.

A simple consequence is that $2$-in/$2$-out digraphs $G$ maximizing the degree
of the interlace polynomial of their interlace graphs have a unique Euler
circuit, see also Lemma 6.7. We can now give a very useful bound for the
degree of the circuit counting polynomial. We denote by $K(G)$ the number of
connected components of the graph $G$; here, we do not care about the
orientation of the edges in $G$ when defining connected components.

######  Lemma 6.6.

Given a 2-in/2-out digraph $G$ on $n$ vertices, we have

$\deg J(G;x)\leq n+K(G).$

###### Proof.

We shall use the relation between the circuit counting polynomial and the
interlace polynomial from Theorem 6.4 on the connected components
$G_{1},G_{2},\ldots,G_{K(G)}$ of $G$. The degree of the interlace polynomial
of a connected $2$-in/$2$-out digraph with $p$ vertices is bounded by $p$ (see
Eq. (28)), hence $\deg J(G_{i})\leq n_{i}+1$, where $n_{i}$ is the number of
vertices of $G_{i}$. Using the multiplicativity property of $J$ from Lemma
6.2, we have

$\deg J(G;x)=\sum_{i=1}^{K(G)}\deg
J(G_{i};x)\leq\sum_{i=1}^{K(G)}n_{i}+1=n+K(G).$

∎

We then use the above bound to show the next lemma.

######  Lemma 6.7.

Among all connected digraphs with a given number of vertices, the ones
maximizing the degree of the circuit-counting polynomial $J$ are precisely the
ones having only cut vertices.

###### Proof.

Assume that a connected digraph $G$ has only $k<p$ out of $p$ vertices which
are cut vertices. Then using skein relations (21), we can pick a vertex that
is not a cut vertex of $G$ and reduce it. We obtain two connected
$2$-in/$2$-out digraphs with $p-1$ vertices $G^{\prime}$ and
$G^{\prime\prime}$. Indeed, they are connected because the vertex we chose to
reduce is not a cut-vertex, hence none of the moves disconnect the digraph.
Then using the bound of Lemma 6.6 above we have that $\textrm{deg}\
J(G^{\prime}),\textrm{deg}\ J(G^{\prime\prime})\leq p$. We also know from
equation (26) that connected $2$-in/$2$-out digraphs with $p$ vertices have
circuit counting polynomial of degree $p+1$, proving the claim. ∎

We introduce a practical notation for the rest of this paper. We denote
$G_{\sigma_{1},\sigma_{2}}$, for $\sigma_{1},\sigma_{2}\in S_{p}$, the
$2$-in/$2$-out digraph on $p$ vertices whose (oriented) edge (multi-)set is

(30)
$E=\\{i\rightarrow\sigma_{1}(i)\,:\,i\in[p]\\}\sqcup\\{i\rightarrow\sigma_{2}(i)\,:\,i\in[p]\\}.$

Note in particular that there is an edge of the form ${i\rightarrow k}$ if
either the matrix of the permutation $\sigma_{1}$ has element $(k,i)$ equal to
$1$ or the matrix of the permutation $\sigma_{2}$ has element $(k,i)$ equal to
$1$. If there are two edges of the form ${i\rightarrow k}$ then both matrices
of the permutations $\sigma_{1},\sigma_{2}$ have elements $(k,i)$ equal to
$1$. Finally, there is no edge from $i$ to $k$ if and only if both
$\sigma_{1},\sigma_{2}$ have elements $(k,i)$ equal to $0$. These facts
together with a moment of reflection reveal that the adjacency matrix
$A_{G_{\sigma_{1},\sigma_{2}}}$ of the $2$-in/$2$-out digraph is the sum of
the two matrices of the permutations $\sigma_{1},\sigma_{2}$

(31) $A_{G_{\sigma_{1},\sigma_{2}}}=\sigma_{1}+\sigma_{2},$

where we used the same notation for the permutations $\sigma_{1,2}$ and their
matrix representations. In general, we also have the reciprocal, that holds in
a more general setting,

######  Proposition 6.8.

Given a $k$-in/$k$-out digraph $G$ with $p$ vertices, there exists $k$
permutations $\sigma_{1},\sigma_{2},\ldots,\sigma_{k}\in S_{p}$ such that

(32) $A_{G}=\sum_{i=1}^{k}\sigma_{i}$

where $A_{G}$ is the adjacency matrix of the digraph $G$.

###### Proof.

First notice that due to the regularity of the digraph $G$, $A_{G}$ is, up-to-
normalization, a bistochatic matrix. Indeed, the $k$-in/$k$-out property
implies that every row and every column of $A_{G}$ must sum to $k$. From this
remark, the proposition is a trivial consequence of the Birkhoff-von Neumann
algorithm applied to $A_{G}$. ∎

This proposition implies that by choosing to index $2$-in/$2$-out digraphs
with two permutations we did not restrict the set of digraphs we are looking
at. However note that the choice of permutations, given a digraph, is not
unique and there may be different collections of permutations representing the
same digraph. It is possible to describe classes of permutations that lead to
equivalent digraphs, but this will be of no use to us.

Figure 3. Structure of a digraph $G_{\gamma,\gamma^{-1}}$.

Assume $\gamma=(1,2,3,\ldots,n)\in S_{n}$ is the full cycle permutation, then

######  Proposition 6.9.

Consider the $2$-in/$2$-out digraph $G_{\gamma,\gamma^{-1}}$, where $\gamma\in
S_{n}$ is the full cycle. The number of covers of $G_{\gamma,\gamma^{-1}}$
having $k$ cycles is given by

(33) $j_{k}=\binom{n}{k}+\delta_{2,k}.$

In particular,

(34) $\mathrm{deg}\,J(G_{\gamma,\gamma^{-1}})=n.$

###### Proof.

The proof follows from the fact that $G_{\gamma,\gamma^{-1}}$ is of the form
depicted in Fig 3. Then it is sufficient to notice that at each vertex of
$G_{\gamma,\gamma^{-1}}$ we have the choice to either come back to where we
came from (vertex state $1$) or to keep going in the same direction (vertex
state $2$). The case for which we decide to keep going at all vertices (i.e.
we pick state $2$ at all vertices) leads to two cycles (responsible for the
$\delta_{2,k}$ in the above formula). While if we decide to have one vertex
with state $1$ and the others with state $2$ then we get one cycle. From that,
a moment of reflection reveals that each time we assign an additional vertex
the state $1$ we create an additional cycle. Thus there are $\binom{n}{k}$
choice of states leading to a cycle cover with $k$ cycles of
$G_{\gamma,\gamma^{-1}}$.

It follows that the degree of $J(G_{\gamma,\gamma^{-1}})$ is $n$. This last
result can also easily be obtained by constructing an interlace graph
$H_{\gamma,\gamma^{-1}}$ from an Eulerian cycle $C_{\gamma,\gamma^{-1}}$ of
$G_{\gamma,\gamma^{-1}}$. It is easy to exhibit such an Eulerian cycle:
$C_{\gamma,\gamma^{-1}}=1,2,3,\ldots,n,1,n,\ldots,3,2$. A pair $(1,q)$ is
clearly interlaced for any $q=2,\ldots,n$ while any other pair $(a,b)$ for
$a\textrm{ and }b\neq 1$ is not interlaced. Thus $H_{\gamma,\gamma^{-1}}$ is
the star graph with the vertex $1$ at its center and thus
$\underset{{X\textrm{ minor of }H}}{\textrm{max }}\textrm{dim }\textrm{ker
}X=n$ the maximizing minor being the graph obtained from $H$ by removing
vertex $1$. ∎

Finally, we now relate the circuit counting graph polynomial $J$ to a
generalized trace of the projector on the symmetric subspace $P_{s}$ (see also
Proposition 9.1 for a more general result obtained using the formalism of
tensor networks).

######  Theorem 6.10.

For all permutations $\alpha,\beta\in S_{p}$, we have

$\operatorname{Tr}_{\alpha,\beta}(P_{s},\ldots,P_{s})=2^{-p}J(G_{\alpha,\beta},N).$

###### Proof.

We start from $P_{s}=(I+F)/2$ and expand the left hand side of the equation in
the statement as a sum over all possible choices of terms:

$\operatorname{Tr}_{\alpha,\beta}(P_{s},\ldots,P_{s})=2^{-p}\sum_{f:[p]\to\\{I,F\\}}\operatorname{Tr}_{\alpha,\beta}\left(\bigotimes_{i=1}^{p}f(i)\right).$

Above, the generalized trace of the right hand side, in diagrammatic notation,
is a collection of loops evaluating to $N$ raised to the power of the number
of loops. Clearly, each such loop corresponds to a circuit in the digraph $G$
given by the permutations $\alpha$ and $\beta$, once the moves $f(i)$ have
been applied at each vertex $i\in[p]$, see (23). ∎

##  7\. Moments of $W^{\Gamma}$

In this section, we shall express the moments of the random matrix
$W^{\Gamma}$ (see Figure 4) as a sum of circuit counting polynomials, allowing
us to determine their asymptotic behaviour. On the way, we shall establish
several useful properties of the aforementioned polynomials, which shall also
be useful in the later sections. The main result, Theorem 7.3 showcases the
power of the relation between the generalized bipartite traces of $P_{s}$ and
the circuit counting polynomial from Theorem 6.10.

Figure 4. The partial transposition $W^{\Gamma}$ of a symmetrized Wishart
matrix.

We first establish the exact form of the moments of $W^{\Gamma}$ in terms of
$J$ polynomials.

######  Proposition 7.1.

The moments of the random matrix
$W^{\Gamma}\in\mathcal{M}_{N^{2}}(\mathbb{C})$ from Eq. (19) are given by:

(35) $\forall p\geq
1,\qquad\mathbb{E}\operatorname{Tr}\left((W^{\Gamma})^{p}\right)=2^{-p}\sum_{\alpha\in
S_{p}}M^{\\#\alpha}J(G_{\gamma\alpha,\gamma^{-1}\alpha};N).$

###### Proof.

The proof is a rather simple application of the Wick formula, in its tensor
network incarnation [CN11]. We need to evaluate the expectation of the trace
of the product of $p$ copies of the matrix from Figure 4. By the Wick theorem,
the result is a sum over permutations $\alpha\in S_{p}$ of diagrams obtained
by connecting the $i$-th $G$ box to the $\alpha(i)$-th $G^{*}$ box

Consider now the 2-in/2-out digraph $G_{\gamma\alpha,\gamma^{-1}\alpha}$
having edges (see (30) for the general case)

(36) $\\{i\mapsto\gamma(\alpha(i))\\}_{i\in
p}\sqcup\\{i\mapsto\gamma^{-1}(\alpha(i))\\}_{i\in p}.$

The formula in the statement follows now from Theorem 6.10. ∎

Let us compute next the exact form of the first two moments. For the first
moment, we have

$\mathbb{E}\operatorname{Tr}W^{\Gamma}=\frac{M}{2}J(\,\includegraphics[valign={c}]{G-1-1.pdf}\,;N)=\frac{MN(N+1)}{2}\sim\frac{c}{2}N^{4}.$

Note that this computation is consistent with the result from Proposition 5.4.

For the second moment ($p=2$), we have to sum over the cases
$\alpha=\operatorname{id}$ and $\alpha=(12)$:

$\displaystyle\mathbb{E}\operatorname{Tr}\left((W^{\Gamma})^{2}\right)$
$\displaystyle=\frac{M^{2}}{4}J(\,\includegraphics[valign={c}]{G-t-t.pdf}\,;N)+\frac{M}{4}J(\,\includegraphics[valign={c}]{G-2-2.pdf}\,;N)$
$\displaystyle=\frac{M^{2}}{4}2N(N+1)+\frac{M}{4}N^{2}(N+1)^{2}\sim\left(\frac{c^{2}}{2}+\frac{c}{4}\right)N^{6},$

where, for the case $\alpha=\operatorname{id}$, we have used Proposition 6.9
with $n=2$.

To analyze larger moments, we need the following technical result.

######  Proposition 7.2.

For a given permutation $\alpha\in S_{p}$, consider the 2-in/2-out digraph
$G_{\gamma\alpha,\gamma^{-1}\alpha}$ defined in Eq. (36). Then,
$G_{\gamma\alpha,\gamma^{-1}\alpha}$ has either one or two connected
components. It has two connected components if and only if the permutation
$\alpha$ is _parity changing_ : for all $i\in[p]$, $\alpha(i)\neq
i(\mathrm{mod}\,2)$.

###### Proof.

In the non-oriented graph $G_{\gamma\alpha,\gamma^{-1}\alpha}$, we have the
following property: for any permutation
$\sigma\in\langle\gamma\alpha,\gamma^{-1}\alpha\rangle$, and every $i\in[p]$,
there is a path between the vertices $i$ and $\sigma(i)$. Note also that
$\gamma^{2}\in\langle\gamma\alpha,\gamma^{-1}\alpha\rangle$. If $p$ is odd,
$\gamma^{2}$ is a full cycle, hence $G_{\gamma\alpha,\gamma^{-1}\alpha}$ is
connected. If $p$ is even, then $\gamma^{2}$ acts transitively on
$\\{2,4,\ldots,p\\}$ and $\\{1,3,\ldots,p-1\\}$ separately. Thus
$K(G_{\gamma\alpha,\gamma^{-1}\alpha})=1$ or $2$, depending on $\alpha$. Since
$\gamma$ is a permutation that changes the parity,
$K(G_{\gamma\alpha,\gamma^{-1}\alpha})=2$ if and only if $\alpha$ also changes
the parity. This concludes the proof. ∎

We describe next the asymptotic behaviour of the larger moments of
$W^{\Gamma}$.

######  Theorem 7.3.

For $p\geq 3$, in the limit $N\to\infty$, the moments of $W^{\Gamma}$
asymptotically behave as

(37)
$\mathbb{E}\operatorname{Tr}\left((W^{\Gamma})^{p}\right)=\frac{c^{p}}{2^{p}}N^{3p}+O(N^{3p-1}).$

###### Proof.

We start from the sum over permutations from (35), and we distinguish three
cases, as follows:

  * •

$\\#\alpha=p$, i.e. $\alpha=\operatorname{id}$. We know from Proposition 6.9
that $J(G_{\gamma,\gamma^{-1}};N)\sim N^{p}$, and thus the term
$\alpha=\operatorname{id}$ behaves as $(c/2)^{p}N^{3p}$.

  * •

$\\#\alpha=p-1$, i.e. $\alpha$ is a transposition. According to Lemma 6.6,
such a term is bounded by a polynomial whose leading term in $N$ has exponent

$2\\#\alpha+p+K(G_{\gamma\alpha,\gamma^{-1}\alpha})=3p-2+K(G_{\gamma\alpha,\gamma^{-1}\alpha}).$

Since $\alpha$ is a transposition and $p\geq 3$, $\alpha$ must have a fixed
point, so it cannot be parity changing. Hence, using Proposition 7.2,
$K(G_{\gamma\alpha,\gamma^{-1}\alpha})=1$, proving that such terms are
subdominant with respect to the previous case $\alpha=\operatorname{id}$.

  * •

$\\#\alpha\leq p-2$. Reasoning as in the previous case, these terms are
clearly subdominant, since $K(G_{\gamma\alpha,\gamma^{-1}\alpha})\leq 2$.

In conclusion, all the terms with $\alpha\neq\operatorname{id}$ are
subdominant with respect to the term $\alpha=\operatorname{id}$, finishing the
proof. ∎

To finish this section, note that the asymptotic behaviour from (37),
corresponding to the moments $p\geq 3$, does not match the ones for $p=1,2$.
This is a signature of the presence of an outlier, an eigenvalue on a larger
scale, which is the phenomenon described in our main result, Theorem 5.1.

##  8\. $t$-channel random matrix and graphical representation

In this section we introduce the _$t$ -channel_111in reference to the $s$,
$t$, $u$ channels of particle physics, as the diagrammatic representation is
similar. In this spirit $W$ corresponds to the $s$-channel. random matrix that
we denote $W_{t}$. We do so as it is better suited for our aims in Sections 9,
10 and 11. We define $W_{t}$ component-wise by

(38) $(W_{t})_{ij,kl}:=(W)_{jl,ki}$

and we provide a diagrammatic representation of the above definition below

(39)
$(W_{t})_{ij,kl}=\raisebox{-31.29802pt}{\includegraphics[scale={0.6}]{W_St.pdf}}.$

In the above diagrammatic representation the red arc represents the
symmetrizer $P_{s}$. The black square represents the complex conjugate
$\bar{G}$ of the Ginibre matrix $G$ of section 4 while the white square
represents $G$ itself. $\bar{G}$ is seen as a linear form on the Hilbert
space, hence input vectors, which is why we display the corresponding tensor
legs with ingoing edges. For similar reason we display tensor legs of $G$ with
outgoing edges. The letter $E$ denotes the environmental Hilbert space that is
traced out. Note also that we have the following equalities:

(40)
$(W_{t}^{*})_{ij,kl}=\raisebox{-31.29802pt}{\includegraphics[scale={0.6}]{W_St-
star.pdf}}$

We prefer to use black and white squares here to simplify the representations
since we will use it intensively in the coming sections.

We introduce the projector on the complement of $\mathbb{C}\ket{\Omega}$

(41) $P_{\overline{\Omega}}=I-\lvert\Omega\rangle\langle\Omega\rvert.$

By the projector property we have
$P_{\overline{\Omega}}^{2}=P_{\overline{\Omega}}$. We want to study the
sequence of moments $\mathbb{E}(\operatorname{Tr}(Q^{p}))$ of the random
matrix

(42) $Q:=P_{\overline{\Omega}}W^{\Gamma}P_{\overline{\Omega}}.$

In order to use graphical methods for tensor networks evaluation (see [BB17]
for an overview of tensor networks techniques and ideas), we will use the
following relation

(43)
$\operatorname{Tr}(Q^{p})=\operatorname{Tr}((FP_{\overline{\Omega}}W_{t}P_{\overline{\Omega}})^{p\textrm{
mod
}2}(P_{\overline{\Omega}}W_{t}^{*}P_{\overline{\Omega}}W_{t}P_{\overline{\Omega}})^{\lfloor
p/2\rfloor}).$

This is easily shown by using the fact that $F^{2}=I$, $FW_{t}F=W_{t}^{*}$ and
$[P_{\overline{\Omega}},F]=0$. These moments have a graphical representation
as ladder diagrams introduced in the next Section 9, and their Wick expansion
produces terms that are indexed by tensor networks that are quotients of the
ladder diagram by the action of Wick pairings. When there is no projectors
$P_{\overline{\Omega}}$, the tensor network evaluates to the circuit counting
polynomial of the tensor network graph. However, due to the presence of the
projector, there are two types of tensors appearing in the tensor network and
this will give a different answer. This is the route we follow in Section 11.
Note however that if we expand all projectors then we obtain, for each term of
the expansion, a $2$-in/$2$-out digraph which represents a tensor network
evaluating to the circuit counting polynomial of the graph. This is the method
we use in this section and in the next Sections 9 and 10.

We denote $P\models\\{1,\ldots,p\\}$ an _interval partition_ of
$\\{1,\ldots,p\\}$. That is, $P$ is a set of subsets
$P_{i}\subseteq\\{1,\ldots,p\\}$ of $P$ such that $\bigsqcup_{i}P_{i}=P$ and
$P_{i}$ are sub-intervals of $\\{1,\ldots,p\\}$ seen cyclically. For instance,

  * •

$p=4$, $P=\\{\\{1,2\\},\\{3,4\\}\\},\ P_{1}=\\{1,2\\},\ P_{2}=\\{3,4\\}$

  * •

$p=4$, $P=\\{\\{2,3\\},\\{4,1\\}\\},\ P_{1}=\\{2,3\\},\ P_{2}=\\{4,1\\}$

  * •

$p=4$, $P=\\{\\{2\\},\\{3\\},\\{4,1\\}\\},\ P_{1}=\\{2\\},\ P_{2}=\\{3\\},\
P_{3}=\\{4,1\\}$.

  * •

$p=3$, $P=\\{\\{1,2,3\\}\\},\ P_{1}=\\{1,2,3\\}$,

a non-example is given by the partition

  * •

$p=5$, $P=\\{\\{1,2,4\\},\\{3,5\\}\\}$.

Assume we have a word $f\in\\{0,1\\}^{p}$. Assuming there exists $i$ such that
$f(i)=1$, we associate to $f$ a subdivision of $P_{f}\models\\{1,\ldots,p\\}$
by putting bars between elements $i-1,i\in\\{1,\ldots,p\\}$ if and only if
$f(i)=1$. $P_{j}$ is the interval defined as the set of elements between the
$(j-1)^{\mathrm{th}}$ and $j^{\mathrm{th}}$ bar.

Figure 5. Examples of a interval subdivisions of $\\{1,\ldots,4\\}$ associated
to words $f\in\\{0,1\\}^{4}$ using the bars construction.

See Fig. 5 for a pictorial description.

Expanding the projectors $P_{\overline{\Omega}}$ appearing in the definition
of $Q$, we have

(44) $\operatorname{Tr}(Q^{p})=\operatorname{Tr}((FW_{t})^{p\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\lfloor p/2\rfloor})\\\
+\sum_{\begin{subarray}{c}f\in\\{0,1\\}^{p}\\\ \exists
i,f(i)=1\end{subarray}}\left(\frac{-1}{N}\right)^{|P_{f}|}\prod_{P_{f,i}\in
P}N\langle\Omega|(W_{t})^{|P_{f,i}|\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\left\lfloor|P_{f,i}|/2\right\rfloor}\lvert\Omega\rangle$

This expansion forms the basis for the expansion of the moments of $Q$ as sums
of circuit counting polynomials of $2$-in/$2$-out digraphs. The term of the
form

(45) $N\langle\Omega|(W_{t})^{|P_{f,i}|\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\left\lfloor|P_{f,i}|/2\right\rfloor}\lvert\Omega\rangle$

can be represented diagrammatically, and we will use this diagrammatic
representation to understand the Wick expansion of the moments of $Q$.

######  Remark 8.1.

Note also that the first term of the right hand side of the equation (45)
above is $\operatorname{Tr}((W^{\Gamma})^{p})$. We use this remark later to
deduce the diagrammatics of
$\mathbb{E}\left(\operatorname{Tr}((W^{\Gamma})^{p})\right)$.

##  9\. Diagrammatics for the moments of $Q$ and $W_{t}$

We start by considering the diagram representation of expression involving
$W_{t}$ in the expansion (44). They are easily obtained by stacking building
blocks of equations (39),(40). We have

(46) $\forall p\in 2\mathbb{N},\ \operatorname{Tr}((FW_{t})^{p\textrm{ mod
}2})(W_{t}^{*}W_{t})^{\lfloor
p/2\rfloor}=\operatorname{Tr}((W_{t}^{*}W_{t})^{p/2})\\\
=\raisebox{-31.29802pt}{\includegraphics[scale={0.5}]{trace-to-ladder-
even.pdf}}\ $

If $p\in 2\mathbb{N}+1$ is odd, we have the slightly twisted representation

(47) $\operatorname{Tr}((FW_{t})(W_{t}^{*}W_{t})^{\lfloor p/2\rfloor})\ =\
\raisebox{-34.1433pt}{\includegraphics[scale={0.5}]{trace-to-ladder-odd.pdf}}$

The gray boxes highlight the copies of $W_{t}$ (resp. $W_{t}^{*}$) which
appear in the product in the trace. We number these copies as shown on the
graphical representation. The black and white squares inherit the numbering of
the gray box they lie in. This diagrammatic allows us to describe the Wick
pairings as permutations $\alpha\in S_{p}$. Indeed, since $G$ is normally
distributed, with vanishing pseudo-variance, Wick pairings only pair instances
of $G$ with instances of $\bar{G}$. We represent such a pairing by adding
oriented dotted edges from black squares to white ones in diagrams of the type
appearing in equations (46) and (47). Each Wick pairing can be indexed by a
permutation in $\alpha\in S_{p}$. Indeed the black box number $i$ is paired
with the white box $j$ if and only if $\alpha(i)=j$. This translates
diagrammatically as a dotted edge between black square $i$ and white square
$j$, oriented from black to white. See examples on Fig. 6.

Figure 6. Two examples of Wick pairings indexed by permutations $\alpha$. The
dotted lines allows to track how the random states are paired together.

Now we realize that each such pairing tells us about the composition of
symmetrizers (denoted as red circle arcs), whose inputs and outputs are kept
track of using the orientations of the black edges in the diagrams. In fact,
these Wick pairings tell us how to match the indices of these symmetrizers
together. In particular, note that if $\alpha(i)=j$, then the symmetrizer
associated to the black box $i$ is composed (in the usual operator sense) with
the symmetrizer associated with the white box $j$. Since the symmetrizer is a
projector, we are left with a unique symmetrizer for each dotted edge whose
inputs are the black edges ending on the black box $i$ and whose outputs are
the black edges exiting the white box $j$. Since a symmetrizer symmetrizes
over both the inputs and the outputs it is not necessary to keep track of the
difference between a first input (resp. output) and a second input (resp.
output). Thus, the corresponding contraction of symmetrizers is indexed by a
directed graph (digraph) whose vertices are $2$-in/$2$-out and represent
two222the factor of $2$ is conventional. Its purpose is to recover the
definition of the $J$ polynomial later. times $P_{s}$. Such a digraph can be
seen as a tensor network for the tensor $P_{s}$ The value of the tensor
network is just the evaluation of the corresponding contraction. We explain in
details this evaluation combinatorially. Since each vertex represents twice a
symmetrizer $P_{s}=\frac{1}{2}(I+F)$, we can evaluate the contraction by
deciding for a _state_ 333we use here the terminology that is familiar to the
graph polynomial literature. Note though that vertex states should not be
confused with quantum states. We hope that the context will be clear enough to
avoid confusion., that is the assignation of either the identity
$I:\mathbb{C}^{N}\otimes\mathbb{C}^{N}\rightarrow\mathbb{C}^{N}\otimes\mathbb{C}^{N}$
or the flip
$F:\mathbb{C}^{N}\otimes\mathbb{C}^{N}\rightarrow\mathbb{C}^{N}\otimes\mathbb{C}^{N}$,
at each vertex. The possible states $\\{S_{1},S_{2}\\}$ are

(48)
$\raisebox{-17.07164pt}{\includegraphics[scale={0.7}]{state_no_state.pdf}}\longrightarrow\raisebox{-34.1433pt}{\includegraphics[scale={0.7}]{state_S_1.pdf}},\quad\raisebox{-34.1433pt}{\includegraphics[scale={0.7}]{state_S_2.pdf}}.$

For each state assignation to all the vertices of the graph, we obtain a
collection of cycles, called a cycle cover of the graph. Each of these cycles
are weighted by $N$. Hence the weight of a particular state assignation is $N$
raised to the number of resulting cycles. Summing over all possible state
assignation will result in the circuit counting polynomial of the graph. Hence
the weight of the tensor network we constructed is just the circuit counting
polynomial of the underlying graph times $\left(\frac{1}{2}\right)^{p}$, with
$p$ the number of vertices. This last factor just takes into account the
normalization factor $\frac{1}{2}$ of the symmetrizer.

The digraph resulting from the Wick pairing $\alpha$ can be indexed by the
data of two permutations that also allows us to construct directly the
adjacency matrix of the digraph.

We describe this permutation representation of $2$-in/$2$-out digraphs
obtained from a Wick pairing. These digraphs are indexed by
$\sigma_{1}=\gamma\alpha,\ \sigma_{2}=\gamma^{-1}\alpha$ with
$\gamma=(1,2,\ldots,p)$. Indeed, a vertex labeled $i$ corresponds, by
definition, to a pair $(i,\alpha(i))$. The edges entering this vertex are the
edges entering the black square of the pair $i$ in the ladder graph while the
edges exiting this vertex are the edges exiting the white square of the gray
pair $\alpha(i)$. Note that the edges exiting from vertex $i$ of the digraph,
are adjacent to the black boxes of the gray pairs $\gamma(\alpha(i))$ and
$\gamma^{-1}(\alpha(i))$, that is in the digraph they are adjacent to the
vertices $\gamma(\alpha(i))$, $\gamma^{-1}(\alpha(i)$ oriented as
$i\rightarrow\gamma(\alpha(i))$ and $i\rightarrow\gamma^{-1}(\alpha(i))$.
Hence, $\sigma_{1}=\gamma\alpha,\ \sigma_{2}=\gamma^{-1}\alpha$. See the local
construction of the vertices on Fig. 7.

Figure 7. On this figure we show how the vertices of the digraph are
constructed from the pairing $\alpha$ by identifying the black box $i$ with
the white box $\alpha(i)$. We forget the environment (green) lines in the
digraph. We highlighted the edges adjacent to the black and white boxes to be
identified in orange and blue with their orientations so that it is easier to
spot them in the digraph.

We follow a similar train of thoughts for the terms in the sum

(49) $\sum_{\begin{subarray}{c}f\in\\{0,1\\}^{p}\\\ \exists
i,f(i)=1\end{subarray}}\left(\frac{-1}{N}\right)^{|P_{f}|}\prod_{P_{f,i}\in
P}N\langle\Omega|(W_{t})^{|P_{f,i}|\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\left\lfloor|P_{f,i}|/2\right\rfloor}\lvert\Omega\rangle$

from equation (44). Indeed, the product $(W_{t})^{|P_{f,i}|\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\left\lfloor|P_{f,i}|/2\right\rfloor}$ is also
represented diagrammatically by stacking building blocks of equations (39) and
(40). The only difference being that the vector $\lvert\Omega\rangle$ has the
following representation

(50) $\lvert\Omega\rangle=\frac{1}{\sqrt{N}}\
\raisebox{-8.53581pt}{\includegraphics[scale={0.6}]{omega-vector.pdf}}$

Hence, we have the diagrammatic below

(51) $N\langle\Omega\rvert(W_{t})^{|P_{f,i}|\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\left\lfloor|P_{f,i}|/2\right\rfloor}\lvert\Omega\rangle=\
\raisebox{-34.1433pt}{\includegraphics[scale={0.48}]{overlaps-general.pdf}}$

We introduce permutations that we use to give a description of the digraphs
appearing in the computation of the above sum. We start with the following
permutations on $l$ elements, whose definition differs in the odd and even
cases,

(52) $\displaystyle\gamma_{a,l}$
$\displaystyle=\begin{cases}(1)(23)(45)\ldots(2k,2k+1),\ \textrm{if
}l=2k+1\textrm{ is odd}\\\ (1)(23)(45)\ldots(2k-2,2k-1)(2k),\ \textrm{if
}l=2k\textrm{ is even}.\end{cases}$ (53) $\displaystyle\gamma_{b,l}$
$\displaystyle=\begin{cases}(12)(34)\ldots(2k-1,2k)(2k+1),\ \textrm{if
}l=2k+1\textrm{ is odd}\\\ (12)(34)\ldots(2k-1,2k),\ \textrm{if }l=2k\textrm{
is even}\end{cases}$

and to each $P_{i}\in P\models\\{1,\ldots,p\\}$ we associate the bijection
$C_{P_{i}}:P_{i}\rightarrow\\{1,\ldots,\lvert P_{i}\rvert\\}$ defined by

(54) $C_{P_{i}}(P_{i}(q))=q$

where $P_{i}(q)$ is the $q^{\textrm{th}}$ element of $P_{i}$. Then we define
for each $P_{i}$ the permutation
$\gamma_{a,P_{i}}=C_{P_{i}}^{-1}\gamma_{a,\lvert P_{i}\rvert}C_{P_{i}}$ and
$\gamma_{b,P_{i}}=C_{P_{i}}^{-1}\gamma_{b,\lvert P_{i}\rvert}C_{P_{i}}$.
Finally, we construct yet another two permutations,
$\Gamma_{a,P},\,\Gamma_{b,P}:=\prod_{i}\gamma_{a,P_{i}},\,\prod_{i}\gamma_{b,P_{i}}$.
We have several figures to illustrate the relations with the diagrams. See
Fig. 8 for $\gamma_{a,\lvert P_{i}\rvert}$ and $\gamma_{b,\lvert
P_{i}\lvert}$. See also Fig. 9 for an example of permutations
$\Gamma_{a,P},\Gamma_{b,p}$.

Figure 8. Illustration of the relation between the diagrams and the
permutations $\gamma_{a,|P_{i}|}$ and $\gamma_{b,|P_{i}|}$. On the topmost
picture: the even $|P_{i}|$ case. On the bottom-most picture: the odd
$|P_{i}|$ case.

Then we have as a consequence of Wick theorem that:

(55) $\mathbb{E}\left(\prod_{P_{f,i}\in
P}N\langle\Omega|(W_{t})^{|P_{f,i}|\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\left\lfloor|P_{f,i}|/2\right\rfloor}\lvert\Omega\rangle\right)=\sum_{\alpha\in
S_{p}}N^{2\\#\alpha}\frac{c^{\\#\alpha}}{2^{p}}J(G_{\Gamma_{a,P}\alpha,\Gamma_{b,P}\alpha},N)$

where we recall that $\sum_{i}\lvert P_{f,i}\rvert=p$.

In the same way than in the previous case, the Wick theorem identifies black
box $i$ with white box $j$ for any pairing such that $\alpha(i)=j$. For the
same reason than before, it leads to a $2$-in/$2$-out digraph. However, now
the edges exiting a vertex $i$ of the digraph are the edges exiting the white
box $\alpha(i)$ while these edges are adjacent to the black boxes
$\Gamma_{a,P}(\alpha(i)),\Gamma_{b,P}(\alpha(i))$, which correspond to the
vertices $\Gamma_{a,P}(\alpha(i))$ and $\Gamma_{b,P}(\alpha(i))$ in the
digraph. This is sufficient to conclude that we have the relation (55).

We summarize the content of the above discussion by the proposition below:

######  Proposition 9.1.

As a consequence of Wick theorem we have,

(56) $\mathbb{E}\left(\operatorname{Tr}((FW_{t}^{p\textrm{ mod
}2})(W_{t}^{*}W_{t})^{\lfloor p/2\rfloor})\right)=\sum_{\alpha\in
S_{p}}N^{2\\#\alpha}\frac{c^{\\#\alpha}}{2^{p}}J(G_{\gamma\alpha,\gamma^{-1}\alpha};N)$

and

(57) $\sum_{\begin{subarray}{c}f\in\\{0,1\\}^{p}\\\ \exists
i,f(i)=1\end{subarray}}\left(\frac{-1}{N}\right)^{|P_{f}|}\prod_{P_{f,i}\in
P}N\langle\Omega|(W_{t})^{|P_{f,i}|\textrm{ mod
}2}(W_{t}^{*}W_{t})^{\left\lfloor|P_{f,i}|/2\right\rfloor}\lvert\Omega\rangle=\\\
\sum_{\begin{subarray}{c}f\in\\{0,1\\}^{p}\\\ \exists
i,f(i)=1\end{subarray}}\sum_{\alpha\in
S_{p}}\left(-1\right)^{|P_{f}|}N^{2\\#\alpha-\lvert
P_{f}\rvert}\frac{c^{\\#\alpha}}{2^{p}}J(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha};N).$

Reformulating the above proposition and using the Remark 8.1, we can recover
the moments of the random matrix $W^{\Gamma}$ from Proposition 7.1

(58)
$\mathbb{E}\left(\operatorname{Tr}\left((W^{\Gamma})^{p}\right)\right)=\sum_{\alpha\in
S_{p}}N^{2\\#\alpha}\frac{c^{\\#\alpha}}{2^{p}}J(G_{\gamma\alpha,\gamma^{-1}\alpha};N).$

We use the short-hand notation, for $\alpha\in S_{p}$

(59)
$w(\alpha,\gamma):=N^{2\\#\alpha}\frac{c^{\\#\alpha}}{2^{p}}J(G_{\gamma\alpha,\gamma^{-1}\alpha}).$

Figure 9. Example of permutations $\Gamma_{a,P}$, $\Gamma_{b,P}$ for
$P=\\{\\{2,3\\},\\{1,4\\}\\}$.

##  10\. Circuit polynomial techniques for
$\mathbb{E}\operatorname{Tr}(Q^{p})$

In this section we prove a sequence of technical bounds and combinatorial
properties that we use to determine the set of permutations $\alpha$
contributing to the asymptotics of the random matrices $W^{\Gamma}$ and $Q$.

We have the following bound

######  Proposition 10.1.

Let $\alpha\in S_{p}$, then there exists $R_{\alpha,f}>0$, that can depend on
$\alpha$ and on $f$ but not on $N$

(60) $N^{\\#2\alpha-\lvert
P_{f}\rvert}J(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha};N)\leq
R_{\alpha,f}N^{2\\#\alpha-\lvert
P_{f}\rvert+p+K(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha})},$

with $K(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha})$ the number of
connected components of the corresponding digraph. Moreover, we have:

$K(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha})-\lvert P_{f}\rvert\leq
0.$

###### Proof.

We now come to the second bound in the case of $2$-in/$2$-out digraphs
$G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha}$. The main difference lies
in their number of connected component. Indeed,
$\Gamma_{a,P_{f}}\Gamma_{b,P_{f}}\in\langle\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha\rangle$
acts transitively on the sets $P_{f,i}$ separately, but is not transitive on
$\\{1,\ldots,p\\}$. Thus the number of connected components is bounded from
above by $\lvert P_{f}\rvert$. Hence we have

(61) $K(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha})-\lvert
P_{f}\rvert\leq 0.$

Again, using the relation between the circuit polynomial and the interlace
polynomial on each connected component of
$G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha}$, we obtain

(62) $\textrm{deg }J(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha};N)\leq
p+K(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha}).$

∎

We also show the two next technical Lemma 10.2 and 10.3. They prove useful in
the last section, to determine the limiting empirical distribution of
eigenvalues of $Q$,

######  Lemma 10.2.

Let $\alpha\in S_{p}$ such that
$2\\#\alpha+p+K(G_{\gamma\alpha,\gamma^{-1}\alpha})=2p+e$ with $e$ a positive
integer. Then, if $K(G_{\gamma\alpha,\gamma^{-1}\alpha})=1$

(63) $e=(\\#_{1}\alpha+1)+\sum_{i\geq
3}(2-i)\\#_{i}\alpha\leq\\#_{1}\alpha+1,$

where $\\#_{i}\alpha$ is the number of cycles of length $i$ of $\alpha$. If
$K(G_{\gamma\alpha,\gamma^{-1}\alpha})=2$,

(64) $e=2+\sum_{k\geq 2}2(1-k)C_{2k}\leq 2.$

Similarly, let $e^{\prime}$ be a positive integer such that
$2\\#\alpha+p+K(G_{\Gamma_{a,P}\alpha,\Gamma_{b,P}\alpha})=2p+e^{\prime}$.
Then,

(65) $e^{\prime}=\\#_{1}\alpha+\sum_{i\geq
3}(2-i)\\#_{i}\alpha+\left(K(G_{\Gamma_{a,P}\alpha,\Gamma_{b,P}\alpha})-\lvert
P\rvert\right).$

In particular, if $\alpha$ is not allowed to have fixed points, then
$e^{\prime}\leq 0$. If we have equality $e^{\prime}=0$, then $\\#_{i}\alpha=0$
for all $i\geq 3$.

###### Proof.

We first focus on $e$. We denote $\\#_{i}\alpha$ the number of cycles of
length $i$ of $\alpha$. We have the following relations

(66) $p=\sum_{i\geq
1}i\\#_{i}\alpha\quad\textrm{and}\quad\\#\alpha=\sum_{i\geq 1}\\#_{i}\alpha.$

From those we deduce

(67) $e=\sum_{i\geq
1}(2-i)\\#_{i}\alpha+K(G_{\gamma\alpha,\gamma^{-1}\alpha}).$

We have two cases, first assume that
$K(G_{\gamma\alpha,\gamma^{-1}\alpha})=1$. Then, a consequence of proposition
7.2 is that $\alpha$ is allowed to have fixed points, thus

(68) $e=(\\#_{1}\alpha+1)+\sum_{i\geq
3}(2-i)\\#_{i}\alpha\leq\\#_{1}\alpha+1.$

Now, assume that $K(G_{\gamma\alpha,\gamma^{-1}\alpha})=2$. A consequence of
proposition 7.2 is that $\alpha$ is not allowed to have fixed point (as those
are not parity changing). More generally, $\alpha$ is not allowed to have
cycles of odd length. Therefore,

(69) $e=2+\sum_{k\geq 2}2(1-k)\\#_{2k}\alpha\leq 2.$

The above inequality is saturated when $\alpha$ has only cycles of length $2$
as any longer cycles will contribute negatively to the sum of the right hand
side. This proves the two first statements of the lemma.
We now come to $e^{\prime}$. We have

(70) $e^{\prime}=\\#_{1}\alpha+\sum_{i\geq
3}(2-i)\\#_{i}\alpha+\left(K(G_{\Gamma_{a,P}\alpha,\Gamma_{b,P}\alpha})-\lvert
P\rvert\right),$

where the second and last term of the right hand side is bounded from above by
$0$. Hence, $e^{\prime}\leq\\#_{1}\alpha$. In particular, if $\alpha$ does not
have fixed points, $e^{\prime}\leq 0$ with equality if $\\#_{i}\alpha=0\
\forall i\geq 3$ and $\alpha(P_{i})=P_{i},\ \forall P_{i}\in P$. ∎

As we understand from the above proposition, fixed points and cycles of length
$2$ of $\alpha$ play an important role. Indeed, the scaling of a given Wick
pairing $\alpha$ can exceed $N^{2p+2}$ if and only if there are enough fixed
points, and not too many cycles of length greater or equal to $3$. The lemma
below is an important step to the proof that we discuss fully in the next
sections.

######  Lemma 10.3.

Assume $\alpha\in S_{p}$ has only cycles of length $2$ and is parity changing
(i.e. $\alpha(i)\textrm{ mod }2\neq i\textrm{ mod }2$). Assume also that there
exists two cycles of $\alpha$ of the form $(a,i+1)$ and $(i,b)$ for
$a<i<i+1<b$. Then, $\textrm{deg }J(G_{\gamma\alpha,\gamma^{-1}\alpha};N)\leq
p$.

We use the parity changing assumption together with the constraint that
$\alpha$ has only length $2$ cycles only to simplify the proof and avoid
having to consider several cases. Indeed, a stronger form of this result can
easily be proven from similar considerations. However, the above form of the
lemma will be sufficient for our purposes.

###### Proof.

The parity changing assumption ensures that
$G_{\gamma\alpha,\gamma^{-1}\alpha}$ has $2$ connected components. We note
that these two connected components must be isomorphic graphs. Indeed, we can
exhibit an isomorphism. For each cycle of the form $(k,l)$ of $\alpha$, the
isomorphism $\varphi$ is set to $\varphi(k)=l$. Hence, we can focus on one
connected component only and bound the degree of the circuit counting
polynomial of this connected component. We use Fig. 10 to show the structure
of $G_{\gamma\alpha,\gamma^{-1}\alpha}$ around vertices $a,b,i,i+1$. From this
figure we note that neither vertex $a$ nor vertex $i$ are cut vertices.
Moreover, since the figure we show is minimally connected (gray parts inside
each connected components could be further connected if $\alpha$ had more
crossings), this means vertices $a$ and $i$ cannot become cut-vertices by
changing $\alpha$ while keeping the cycles $(a,i+1)$ and $(i,b)$ fixed.

Denote $G^{(a)}_{\gamma\alpha,\gamma^{-1}\alpha}$ the connected component
containing vertex $a$. We focus on this connected component for now. By Lemma
6.7, we have $\textrm{deg }J(G_{\gamma\alpha,\gamma^{-1}\alpha}^{(a)};N)\leq
p/2$. Indeed, since for instance $a$ is not a cut-vertex we know that once we
reduced vertex $a$ by summing over its two possible states using the skein
relations of Eq. (21) on $a$ we are left with two graphs $G^{\prime}$ and
$G^{\prime\prime}$ such that

(71)
$J(G^{(a)}_{\gamma\alpha,\gamma\alpha^{-1}};N)=J(G^{\prime};N)+J(G^{\prime\prime};N).$

Both $G^{\prime}$ and $G^{\prime\prime}$ have a maximum of $p/2-1$ vertices.
Therefore, from Lemma 6.6

(72) $\textrm{deg
}J(G^{(a)}_{\gamma\alpha,\gamma\alpha^{-1}};N)=\textrm{max}(\textrm{deg
}J(G^{\prime};N),\textrm{deg }J(G^{\prime\prime};N))\leq p/2.$

By the isomorphism argument, we have the same bound on the degree of
$J(G^{(b)}_{\gamma\alpha,\gamma^{-1},\alpha};N)$, hence

(73) $\textrm{deg }J(G_{\gamma\alpha,\gamma^{-1}\alpha})\leq p.$

∎

In the next section, we show that the contribution of fixed points is smaller
than the previous bounds would let us believe. In order to study fixed points
and cycles of length $2$ we introduce more general tensor networks.

Figure 10. Local structure of the graph $G_{\gamma\alpha,\gamma^{-1}\alpha}$
for $\alpha$ satisfying the constraints of Lemma 10.3. The axe of symmetry
corresponds to the ladder axe of symmetry. The figure shows the smallest
connexity case. Indeed, the gray part of each connected component could be
further connected if $\alpha$ have more crossings. For readability we color
the edges differently depending on which vertices they are adjacent to. The
color change in the middle if the edges are shared between vertices $a$ and
$i$ or between vertices $b$ and $i+1$.

##  11\. Tensor network evaluation of $\mathbb{E}\operatorname{Tr}(Q^{p})$

In this section, we express the moments of the random matrix $Q$ as a sum over
evaluations of tensor networks made of the tensors $P_{s}$ and
$P_{\overline{\Omega}}$. In particular, we do not expand the projectors
$P_{\overline{\Omega}}$ appearing in the expression of $Q$. In that case, for
each Wick pairing $\alpha$ we obtain a contraction of two types of tensors,
$P_{s}:\mathbb{C}^{N}\otimes\mathbb{C}^{N}\rightarrow\mathbb{C}^{N}\otimes\mathbb{C}^{N}$
as before and
$P_{\overline{\Omega}}:\mathbb{C}^{N}\otimes\mathbb{C}^{N}\rightarrow\mathbb{C}^{N}\otimes\mathbb{C}^{N}$
which is the projector on the complement of
$\textrm{span}\\{\lvert\Omega\rangle\\}$ introduced earlier in Eq. (41).

We work here with a graphical representation that is very similar to the ones
introduced in previous paragraphs. Here, we represent
$\operatorname{Tr}(Q^{p})$ directly, without expanding. We do so by stacking
the building blocks introduced in (39), (40) together with an additional
building block representing $P_{\overline{\Omega}}$

(74)
$P_{\overline{\Omega}}=I-\lvert\Omega\rangle\langle\Omega\rvert=\raisebox{-17.07164pt}{\includegraphics[scale={0.8}]{proj_comp_omega.pdf}}$

These building blocks are stacked in between (39), (40) blocks to represent
the operator $P_{\overline{\Omega}}$. The diagram thus obtained is a tensor
network representing $\operatorname{Tr}(Q^{n})$

(75)
$\operatorname{Tr}(Q^{n})=\raisebox{-34.1433pt}{\includegraphics[scale={0.50}]{trace-
Q-p.pdf}}$

To compute $\mathbb{E}(\operatorname{Tr}(Q^{p}))$ we proceed as before, with
each Wick pairing identifying black box $i$ with white box $\alpha(i)$ and
leading to a vertex labeled $i$. This vertex represents a $P_{s}$ operator.
However, the resulting oriented graph also contains $\overline{\Omega}$-boxes
that represents $P_{\overline{\Omega}}$. This leads to a slightly more
complicated connection pattern. We show on Fig. 11 the local structure around
a vertex $i$.

Figure 11. Local structure of a tensor network given a Wick pairing $\alpha$.
Left: the Wick pairing on the diagram representing
$\mathbb{E}(\operatorname{Tr}(Q^{p}))$. Right: the corresponding local
structure of the tensor network of $P_{s}$ and $P_{\overline{\Omega}}$
operators.

The only change is that now vertices are only connected to
$\overline{\Omega}$-boxes. Locally it is simple to check that vertex $i$ is
connected to $\overline{\Omega}$-boxes $i$, $\gamma(i)$ via ingoing edges,
while it is connected to $\overline{\Omega}$-boxes $\alpha(i)$,
$\gamma(\alpha(i))$ via outgoing edges. This is due to the fact that the white
box $\alpha(i)$ (of which adjacent edges are outgoing) is connected to
$\overline{\Omega}$-boxes $\alpha(i)$, $\gamma(\alpha(i))$ while the black box
$\alpha(i)$ (of which adjacent edges are ingoing) is connected to
$\overline{\Omega}$-boxes $i$ and $\gamma(i)$. In the rest of this paper we
denote by $T_{\alpha}$ the tensor network (and its evaluation) corresponding
to a pairing $\alpha$ constructed in the way described above.

Using this tensor network construction, we have

(76) $\mathbb{E}\operatorname{Tr}(Q^{p})=\sum_{\alpha\in
S_{p}}M^{\\#\alpha}T_{\alpha}.$

Note that we have the relation

(77)
$T_{\alpha}=\frac{1}{2^{p}}\left(J(G_{\gamma\alpha,\gamma^{-1}\alpha};N)+\sum_{\begin{subarray}{c}f\in\\{0,1\\}^{p}\\\
\exists i,f(i)=1\end{subarray}}\left(-1\right)^{|P_{f}|}N^{-\lvert
P_{f}\rvert}J(G_{\Gamma_{a,P_{f}}\alpha,\Gamma_{b,P_{f}}\alpha};N)\right).$

In the next paragraph, we prove reduction relations for the evaluation of such
tensor networks at fixed points, and simple transpositions of $\alpha$. This
will allow us to consider permutations that do not contain fixed points and
only have cycles of length $2$.

Reduction property for tensor network evaluation

Le $c_{q}$ be a cycle of a permutation $\sigma=c_{1}c_{2}\ldots
c_{\\#\sigma}\in S_{p}$. We denote by $\sigma\div
c_{q}:=c_{1}c_{2}\ldots\hat{c_{q}}\ldots c_{\\#\sigma}$ the permutation acting
on the set $\\{1,\ldots,p\\}\setminus\\{i\,:\,i\in c_{q}\\}$. We have the
following reduction property for tensor network $T_{\alpha}$ such that
$\alpha$ have the corresponding feature

######  Proposition 11.1.

Let $\alpha\in S_{p}$. Then,

  * •

if $\alpha$ as a fixed point, that we denote $i$, then

(78) $T_{\alpha}=\frac{1}{2}T_{\alpha\div(i)},$

where $\tilde{\alpha}$ is the permutation satisfying
$\alpha=\tilde{\alpha}(i)$ in cycle notation. In particular, if $\alpha$ has
several fixed points, consider $\alpha_{R}$ the reduced permutation acting on
$\\{1,\ldots,p\\}\backslash\operatorname{Stab}(\alpha)$, then one has

(79)
$T_{\alpha}=\left(\frac{1}{2}\right)^{|\operatorname{Stab}(\alpha)|}T_{\alpha_{R}}$

  * •

if $\alpha$ contains a cycle of the form $(i,i+1)$, then

(80) $T_{\alpha}=\Lambda(N)T_{\alpha\div(i,i+1)},$

where $\Lambda(N)=\frac{1}{4}(N^{2}+2(N-\frac{1}{N}))$.

Note that the special cases where $\alpha=(1)$ and $\alpha=(12)$ are dealt
with by setting
$T_{\alpha\div(1)}=\operatorname{Tr}(P_{\overline{\Omega}})=N^{2}-1$ for the
first case and
$T_{\alpha\div(12)}=\operatorname{Tr}(P_{\overline{\Omega}})=N^{2}-1$ for the
second case. For convenience, we denote
$T_{\emptyset}:=\operatorname{Tr}(P_{\overline{\Omega}})$.

###### Proof.

We proceed to the proof by diagrammatic manipulation of the tensor network.
Fixed points $i$ of $\alpha$ can be treated locally by realizing that the
corresponding tensor network always have the same structure around vertex $i$
(see the left hand side of equation (81))

(81) $\displaystyle T_{\alpha}\ =\
\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{Fixed-points-red-
lhs.pdf}}\ $ $\displaystyle=\ \frac{1}{2}\
\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{Fixed-points-red-
rhs-1.pdf}}+\frac{1}{2}\
\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{Fixed-points-red-
rhs-2.pdf}}\ =\frac{1}{2}\
\raisebox{-56.9055pt}{\includegraphics[scale={0.4}]{Fixed-points-red-
rhs-1.pdf}}$

where we used the fact that

(82) $\raisebox{-22.76219pt}{\includegraphics[scale={0.7}]{PcompOmega.pdf}}\
=\ \sqrt{N}P_{\overline{\Omega}}\lvert\Omega\rangle=0.$

Together with the fact that $P_{\overline{\Omega}}^{2}=P_{\overline{\Omega}}$
his shows that

(83) $T_{\alpha}=\frac{1}{2}T_{\alpha\div(i)}.$

The second statement can be obtained by looking at what happens locally on the
vertices $i$ and $i+1$ when $\alpha$ has a cycle of the form $(i,i+1)$. One
has the following diagrammatic relation

(84) $T_{\alpha}\ =\
\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{i-i+1-red-lhs.pdf}}\ =\
\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{i-i+1-red-
rhs-1a.pdf}}-\frac{1}{N}\
\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{i-i+1-red-rhs-1b.pdf}}$

obtained by using the definition of $P_{\overline{\Omega}}$. The leftmost term
of the right hand side of equation (84) leads after further expanding the two
black vertices

(85) $\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{i-i+1-red-
rhs-1a.pdf}}=\frac{1}{4}(N+1)^{2}\
\raisebox{-56.9055pt}{\includegraphics[scale={0.4}]{Fixed-points-red-
rhs-1.pdf}}$

while the rightmost term of the right hand side of equation (84) can also be
extended further. Together with using again the relation (82) we have

(86) $\frac{1}{N}\
\raisebox{-56.9055pt}{\includegraphics[scale={0.42}]{i-i+1-red-rhs-1b.pdf}}\
=\ \frac{1}{4}\left(\frac{N+1}{N}+\frac{1}{N}\right)\
\raisebox{-56.9055pt}{\includegraphics[scale={0.4}]{Fixed-points-red-
rhs-1.pdf}}.$

Putting everything together we obtain

(87) $T_{\alpha}=\frac{1}{4}\left((N+1)^{2}-\frac{N+1}{N}-\frac{1}{N}\right)\
\raisebox{-56.9055pt}{\includegraphics[scale={0.4}]{Fixed-points-red-
rhs-1.pdf}}$

hence $\Lambda(N)=\frac{1}{4}(N^{2}+2(N-\frac{1}{N}))$. ∎

Thanks to this proposition we can consider permutations without fixed points
by just remembering that for every $\alpha\in S_{p}$, there is a reduced
permutation $\alpha_{R}$ without fixed points such that equation (78) is
satisfied. In particular it is important to note that by getting rid of the
fixed points, we did not produce any factor of $N$. In terms of the right hand
side of (77), this is due to cancellations between terms in the sum over words
$f$. The tensor networks formulation allows to formulate these cancellations
compactly. One particular case is if $\alpha\in S_{p}$ is the identity
permutation. Then we have

(88) $T_{\alpha=id}=\frac{1}{2^{p}}(N^{2}-1).$

Another particular case is when $\alpha\in\textrm{NC}_{2}(p)$, in that case it
is always possible to find a cycle of the form $(i,i+1)$ and reducing it leads
to another permutation $\alpha^{\prime}$ which is in $\textrm{NC}_{2}(p-2)$.
Using this property we can recursively reduce all the cycles of such a
permutation and we obtain

(89)
$T_{\alpha\in\textrm{NC}_{2}(p)}=\Lambda(N)^{p/2}(N^{2}-1)=\frac{1}{2^{p}}N^{p+2}+O(N^{p}).$

##  12\. Proof of the main result

We now have all the tools we need to come to the proof of the main theorem
5.1. This proof follows the steps:

  1. (1)

Determining the limiting spectrum of $Q$. This is Theorem 12.1. It is proved
by nested applications of Proposition 11.1, Lemma 10.2 and 10.3.

  2. (2)

Showing a localization result for the overlap
$\left\langle\Omega\left\lvert\frac{W^{\Gamma}}{N^{3}}\right\rvert\Omega\right\rangle$
in Proposition 12.2.

  3. (3)

Proving the main theorem by using Proposition 12.2 and Theorem 7.3 to control
the largest and smallest eigenvalues of $W^{\Gamma}$ and using Cauchy’s
interlacing theorem on $\frac{W^{\Gamma}}{N^{3}}$ and $Q$ to obtain the
convergence of the empirical distribution of the $N^{2}-1$ smallest
eigenvalues of $W^{\Gamma}$.

Limiting spectrum of $Q$. We tackle the first step of our three steps plan of
proof. Assume that the permutation $\alpha$ do not contain fixed points. By
virtue of Lemma 10.2, we know that $e\leq 2$ and $e^{\prime}\leq 0$ for such
$\alpha$. Moreover, these inequalities can be saturated if $\alpha$ has cycles
of length $2$. This implies that there exists a positive constant $R$,
independent of $N$,

(90) $\lvert T_{\alpha}\rvert\leq RN^{2p+2}$

which can be attained if $\alpha$ has only length $2$ cycles. This remark
allows us to prove the following theorem

######  Theorem 12.1.

The limiting moments of the matrix $Q$ are given by

(91)
$\frac{1}{N^{2p+2}}\mathbb{E}\left(\operatorname{Tr}(Q^{p})\right)=\frac{1}{2^{p}}\sum_{\alpha\in\textrm{NC}_{1,2}(p)}c^{\\#\alpha}.$

This completely determines the limiting law of eigenvalues of $Q$ as being a
shifted semi-circular law of mean $c/2$ and variance $c/4$. In particular, the
limiting support is $[-\sqrt{c}+c/2,\sqrt{c}+c/2]$ which is contained in the
positive real line if and only if $c\geq 4$.

###### Proof.

We start with the relation

(92) $\mathbb{E}\left(\operatorname{Tr}(Q^{p})\right)=\sum_{\alpha\in
S_{p}}N^{2\\#\alpha}c^{\\#\alpha}T_{\alpha}.$

Thanks to Proposition 11.1 we can reduce to permutations $\alpha^{\prime}$
that do not contain fixed points at the cost of a factor $\frac{1}{2}$, that
is

(93)
$\mathbb{E}\left(\operatorname{Tr}(Q^{p})\right)=\sum_{q=0}^{p}\frac{1}{2^{q}}\binom{p}{q}\sum_{\alpha^{\prime}\in
S_{p-q}}N^{2q+2\\#\alpha^{\prime}}c^{q+\\#\alpha^{\prime}}T_{\alpha^{\prime}}.$

According to Lemma 10.2, we have

(94) $\sum_{q=0}^{p}\frac{1}{2^{q}}\binom{p}{q}\sum_{\alpha^{\prime}\in
S_{p-q}}N^{2q+2\\#\alpha^{\prime}}c^{q+\\#\alpha^{\prime}}T_{\alpha^{\prime}}=\sum_{q=0}^{p}\frac{1}{2^{q}}\binom{p}{q}\sum_{\begin{subarray}{c}\alpha^{\prime}\in
S_{p-q}\\\
\\#_{2}\alpha^{\prime}=\frac{p-q}{2}\end{subarray}}N^{2q+2\\#\alpha^{\prime}}c^{q+\\#\alpha^{\prime}}T_{\alpha^{\prime}}+O(N^{2p+1})$

as indeed all the cycles of $\alpha^{\prime}$ must be of length $2$ since
$\alpha^{\prime}$ is not allowed to have fixed points. Making use of Lemma
10.3 and the relation of equation (77), we can further reduce the sum

(95) $\sum_{\begin{subarray}{c}\alpha^{\prime}\in S_{p-q}\\\
\\#_{2}\alpha^{\prime}=\frac{p-q}{2}\end{subarray}}N^{2q+2\\#\alpha^{\prime}}c^{q+\\#\alpha^{\prime}}T_{\alpha^{\prime}}=\sum_{\alpha^{\prime}\in\textrm{NC}_{2}(p-q)}N^{2q+2\\#\alpha^{\prime}}c^{q+\frac{p-q}{2}}T_{\alpha^{\prime}}+O(N^{2p+1})$

because if $\alpha$ only have cycles of length $2$ and
$\alpha\notin\textrm{NC}_{2}(p-q)$ then $\alpha$ can be reduced to a non
trivial permutation $\tilde{\alpha}$ that do not contain cycles of the form
$(i,i+1)$ by using Proposition 11.1 and recursively reducing cycles of this
type. Then, $\tilde{\alpha}$ must have two cycles of the form $(a,i+1),(i,b)$
with $a<i<i+1<b$. This allows us to use Lemma 10.3 to bound the contribution
of $T_{\tilde{\alpha}}$. Finally, using the consequence of Proposition 11.1
for permutations in $\textrm{NC}_{2}(k)$, we have that

(96)
$\displaystyle\sum_{\alpha^{\prime}\in\textrm{NC}_{2}(p-q)}N^{2q+2\\#\alpha^{\prime}}c^{q+\frac{p-q}{2}}T_{\alpha^{\prime}}$
$\displaystyle=\textrm{Cat}_{\frac{p-q}{2}}N^{2q+2\\#\alpha^{\prime}}c^{q+\frac{p-q}{2}}\Lambda(N)^{\frac{p-q}{2}}T_{\emptyset}$
(97)
$\displaystyle=\frac{1}{2^{p-q}}c^{q+\frac{p-q}{2}}\textrm{Cat}_{\frac{p-q}{2}}N^{2p}(N^{2}-1),$

provided $p-q$ is even. Otherwise the sum is zero because there are no terms
in the sum. Putting these nested sums together, we have that

(98)
$\displaystyle\mathbb{E}(\operatorname{Tr}(Q^{p}))=\frac{1}{2^{p}}\sum_{q=0}^{p}\binom{p}{q}\mathbf{1}_{p-q=0\textrm{
mod }2}c^{q+\frac{p-q}{2}}\textrm{Cat}_{\frac{p-q}{2}}N^{2p+2}+O(N^{2p+1}).$

This sum rewrites as

(99)
$\mathbb{E}(\operatorname{Tr}(Q^{p}))=\frac{N^{2p+2}}{2^{p}}\sum_{\alpha\in\textrm{NC}_{1,2}(p)}c^{\\#\alpha}+O(N^{2p+1}).$

∎

Localization of the overlap. We now come to the second step of our three steps
plan. We state and prove the localization result for the overlap below.

######  Proposition 12.2.

For all $p\geq 1$,

$\mathbb{E}\left\langle\Omega\left|\frac{W^{\Gamma}}{N^{3}}\right|\Omega\right\rangle^{p}=N^{-4p}\sum_{\alpha\in
S_{p}}M^{\\#\alpha}N[2]^{\\#\alpha}.$

In particular, we have

$\mathbb{E}\left\langle\Omega\left|\frac{W^{\Gamma}}{N^{3}}\right|\Omega\right\rangle=\frac{c}{2},\qquad\qquad\operatorname{Var}\left\langle\Omega\left|\frac{W^{\Gamma}}{N^{3}}\right|\Omega\right\rangle=\frac{MN[2]}{N^{8}}\sim\frac{c}{2}N^{-4},$

and

$\lim_{N\to\infty}\mathbb{E}\left\langle\Omega\left|\frac{W^{\Gamma}}{N^{3}}\right|\Omega\right\rangle^{p}=\left(\frac{c}{2}\right)^{p}.$

###### Proof.

To show the first statement, use the Wick formula to write

$\mathbb{E}\left\langle\Omega\left|\frac{W^{\Gamma}}{N^{3}}\right|\Omega\right\rangle^{p}=N^{-4p}\sum_{\alpha\in
S_{p}}M^{\\#\alpha}\operatorname{Tr}_{\alpha,\alpha}(P_{s}).$

Indeed, we have, graphically,

Use now the fact that that $P_{s}$ is a projection, and thus

$\operatorname{Tr}_{\alpha,\alpha}(P_{s})=(\operatorname{Tr}P_{s})^{\\#\alpha}=N[2]^{\\#\alpha}.$

∎

Interlacing of eigenvalues and asymptotic of $W^{\Gamma}$. We are now ready to
prove the main result of the paper. As a reminder, let us start by stating a
classical result, namely Cauchy’s interlacing theorem

######  Theorem 12.3 (Cauchy’s interlacing theorem).

Let $A$ be a Hermitian operator on a Hilbert space $\mathcal{H}$, and let $B$
be its compression to a codimension $k$ subspace $\mathcal{N}$. Then for
$j=1,2,\ldots,n-k$

(100) $\lambda_{j}(A)\geq\lambda_{j}(B)\geq\lambda_{j+k}(A).$

We use this theorem to recover the behaviour of the limiting eigenvalues of
$W^{\Gamma}$ from the limiting spectrum of $Q$ and the large eigenvalue of
$W^{\Gamma}$. We now come to the proof of theorem 5.1.

###### Proof of the Theorem 5.1.

Let us first show that the largest eigenvalue of the matrix $W^{\Gamma}$
behaves, asymptotically as $N\to\infty$, as $\frac{c}{2}N^{3}$. To this end,
we shall compare it with two quantities: the overlap
$\langle\Omega|W^{\Gamma}|\Omega\rangle$ and wih the moments of the random
matrix $W^{\Gamma}$.

On the one hand, we have (we denote $\lambda_{\max}=\lambda_{1}$ for clarity),
for all $p\geq 2$,

$\operatorname{Tr}\left(\frac{W^{\Gamma}}{N^{3}}\right)^{2p}\geq\lambda_{\max}^{2p}\left(\frac{W^{\Gamma}}{N^{3}}\right)\geq\left\langle\Omega\left\lvert\frac{W^{\Gamma}}{N^{3}}\right\rvert\Omega\right\rangle^{2p}.$

Using Theorem 7.3 and Proposition 12.2, we conclude that

$\lim_{N\to\infty}\lambda_{\max}^{2p}\left(\frac{W^{\Gamma}}{N^{3}}\right)=\left(\frac{c}{2}\right)^{2p}.$

By Markov’s inequality, we have thus

$\lim_{N\to\infty}\lambda_{\max}^{2}\left(\frac{W^{\Gamma}}{N^{3}}\right)=\left(\frac{c}{2}\right)^{2}$

in probability. But since $\operatorname{Tr}W^{\Gamma}=\operatorname{Tr}W\geq
0$, $\lambda_{\max}(W^{\Gamma})\geq 0$ almost surely, so we also have, in
probability,

$\lim_{N\to\infty}\lambda_{\max}\left(\frac{W^{\Gamma}}{N^{3}}\right)=\frac{c}{2},$

establishing the first claim.

In order to prove the second part of the statement, concerning the lower
$N^{2}-1$ eigenvalues of $W^{\Gamma}$, we shall make use of Theorems 12.1 and
12.3. Indeed, if we denote by $\mu_{1}\geq\cdots\geq\mu_{N^{2}-1}$ the (non-
zero) eigenvalues of the matrix $N^{-2}Q$ and by
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{N^{2}-1}\geq\lambda_{N^{2}}$
those of $N^{-2}W^{\Gamma}$, by Cauchy’s interlacing theorem, we have

$\lambda_{1}\geq\mu_{1}\geq\lambda_{2}\geq\mu_{2}\geq\cdots\geq\lambda_{N^{2}-1}\geq\mu_{N^{2}-1}\geq\lambda_{N^{2}}.$

On the other hand, we know by Theorem 12.1 that the empirical distribution of
the $\mu$ eigenvalues converges to the desired shifted semicircle
distribution. Hence, in order to conclude, we need to control the smallest
eigenvalue $\lambda_{N^{2}}$ and show that it does not contribute
asymptotically to moments. We have

$\displaystyle\mathbb{E}\lambda_{\min}^{4}\left(\frac{W^{\Gamma}}{N^{3}}\right)$
$\displaystyle\leq\mathbb{E}\operatorname{Tr}\left(\frac{W^{\Gamma}}{N^{3}}\right)^{4}-\mathbb{E}\lambda_{\max}^{4}\left(\frac{W^{\Gamma}}{N^{3}}\right)$
$\displaystyle\leq\mathbb{E}\operatorname{Tr}\left(\frac{W^{\Gamma}}{N^{3}}\right)^{4}-\mathbb{E}\left\langle\Omega\left|\frac{W^{\Gamma}}{N^{3}}\right|\Omega\right\rangle^{4}\leq
o(1)+\left(\frac{c}{2}\right)^{4}-\left(\frac{c}{2}\right)^{4}=o(1).$

Hence, $\mathbb{E}|\lambda_{\min}(W^{\Gamma})|=o(N^{3})$, which implies that
$\mathbb{E}|\lambda_{N^{2}}|=o(N)$, completing the proof. ∎

##  13\. Conclusion

In this paper we have introduced the ensembles of bosonic and fermionic
density matrices, which are random quantum states supported on the symmetric,
resp. anti-symmetric subspaces of a tensor product space. We have then studied
the entanglement of typical states from these ensembles, focusing on the
partial transposition criterion.

In the fermionic case, we have found that the partial transposition of such a
random density matrix has a large negative eigenvalue, hence a typical
fermionic bipartite density matrix is entangled. The bosonic case is more
delicate, since the symmetry condition imposes a large positive eigenvalue for
the partial transposition, so a finer analysis of the spectrum is required.

We use the moment method to completely describe the spectrum of the partial
transposition of a symmetric random density matrix. We find that the bulk of
the spectrum follows a shifted semicircular distribution at scale $N^{2}$,
similarly to the non-symmetric case discussed in [Aub12]. We also find a large
outlier positive eigenvalue, on the scale $N^{3}$, which is a signature of the
symmetry. Our method relies on establishing a connection between the moments
of symmetric random matrices and the circuit counting graph polynomial. This
allows us to characterize the asymptotic ratio between the system size and the
environment size for which the PPT criterion is strong enough to certify
entanglement.

The current paper is a first step in the study of the entanglement properties
of the bosonic and fermionic ensembles of quantum states introduced here,
several questions being left for future work. In the bipartite case (which is
the main focus of this paper), the question of the convergence of the smallest
eigenvalue of the partial transposition towards the left edge of the spectrum
is left open. Probably, this would require the detailed analysis of large
moments of the corresponding random matrix, similarly to the method used in
[Aub12]; such an analysis is hindered in the current situation by the presence
of the outlier eigenvalue on a larger scale.

We have not addressed the multipartite case, $r\geq 3$, mainly because the
tools used in this work (connection to graph polynomials) are not efficient
enough in the general case. However, first investigations of this general case
suggest that an interesting phenomenon is at play. Indeed, we expect the
spectrum of the partial transpose to split at several scales, such that the
empirical eigenvalues distribution has several ($r-1$) connected dense parts
with each part appearing at a specific scale. We expect that the rigorous
study of this general case will make heavy use of representation theory.

Another interesting direction for future study is the de Finetti theorem for
the bosonic ensemble. This would require analyzing the separability property
of the partial trace of a symmetric random density matrix in the regime where
the parameter $r$ is large, and the Hilbert space dimension $N$ is fixed. This
asymptotic regime usually requires completely different techniques in random
matrix theory.

In conclusion, we have set up and studied basic entanglement properties of
symmetric (and anti-symmetric) random density matrices. Several questions
regarding these ensembles are left open, and we hope that our work will
generate interest in these ensembles of random matrices / tensors, and their
connection to combinatorics and graph polynomials.

### Acknowledgments.

S. Dartois has been partially supported by the Australian Research Council
grant DP170102028, the French ANR project ANR-18-CE47-0010 (QUDATA). I.
Nechita has been supported by the ANR project ESQuisses, grant number
ANR-20-CE47-0014-01. Both S. Dartois and I. Nechita were supported on this
project by the ANR project “Investissements d’avenir” ANR-11-LABX-0040. A.
Tanasa has been partially supported by the ANR-20-CE48-0018 “3DMaps” grant.

## References

  * [ABS04] Richard Arratia, Béla Bollobás, and Gregory B Sorkin.  The interlace polynomial of a graph.  Journal of Combinatorial Theory, Series B, 92(2):199–233, 2004\.
  * [AN12] Guillaume Aubrun and Ion Nechita.  Realigning random states.  Journal of Mathematical Physics, 53(10):102210, 2012.
  * [Aub12] Guillaume Aubrun.  Partial transposition of random states and non-centered semicircular distributions.  Random Matrices: Theory and Applications, 1(02):1250001, 2012.
  * [BB17] Jacob Biamonte and Ville Bergholm.  Tensor networks in a nutshell.  arXiv preprint arXiv:1708.00006, 2017.
  * [BH11] Robert Brijder and Hendrik Jan Hoogeboom.  Nullity invariance for pivot and the interlace polynomial.  Linear algebra and its applications, 435(2):277–288, 2011.
  * [BS10] Zhidong Bai and Jack W Silverstein.  Spectral analysis of large dimensional random matrices, volume 20.  Springer, 2010.
  * [BZ06] Ingemar Bengtsson and Karol Zyczkowski.  Geometry of quantum states: an introduction to quantum entanglement.  Cambridge University Press, 2006.
  * [CN11] Benoît Collins and Ion Nechita.  Gaussianization and eigenvalue statistics for random quantum channels (iii).  The Annals of Applied Probability, pages 1136–1179, 2011.
  * [CNŻ10] Benoît Collins, Ion Nechita, and Karol Życzkowski.  Random graph states, maximal flow and fuss–catalan distributions.  Journal of Physics A: Mathematical and Theoretical, 43(27):275303, 2010.
  * [CW02] Kai Chen and Ling-An Wu.  A matrix realignment method for recognizing entanglement.  arXiv preprint quant-ph/0205017, 2002.
  * [DLN20] Stephane Dartois, Luca Lionni, and Ion Nechita.  The joint distribution of the marginals of multipartite random quantum states.  Random Matrices: Theory and Applications, 9(03):2050010, 2020.
  * [EM04] Joanna A Ellis-Monaghan.  Identities for circuit partition polynomials, with applications to the tutte polynomial.  Advances in Applied Mathematics, 32(1-2):188–197, 2004.
  * [GdOZ10] Silvano Garnerone, Thiago R de Oliveira, and Paolo Zanardi.  Typicality in random matrix product states.  Physical Review A, 81(3):032336, 2010.
  * [Gur04] Leonid Gurvits.  Classical complexity and quantum entanglement.  Journal of Computer and System Sciences, 69(3):448–484, 2004.
  * [Hal98] Michael JW Hall.  Random quantum correlations and density operator distributions.  Physics Letters A, 242(3):123–129, 1998.
  * [Har13] Aram W Harrow.  The church of the symmetric subspace.  arXiv preprint arXiv:1308.6595, 2013.
  * [HHH96] Michał Horodecki, Paweł Horodecki, and Ryszard Horodecki.  Separability of mixed states: necessary and sufficient conditions.  Physics Letters A, 223(1):1–8, 1996.
  * [HHHH09] Ryszard Horodecki, Paweł Horodecki, Michał Horodecki, and Karol Horodecki.  Quantum entanglement.  Reviews of Modern Physics, 81(2):865, 2009.
  * [Hol18] Randall R Holmes.  Partial trace of a full symmetrizer.  arXiv preprint arXiv:1805.09721, 2018.
  * [LV83] Michel Las Vergnas.  Le polynôme de martin d’un graphe eulérien.  In North-Holland Mathematics Studies, volume 75, pages 397–411. Elsevier, 1983.
  * [MP67] Vladimir A Marčenko and Leonid Andreevich Pastur.  Distribution of eigenvalues for some sets of random matrices.  Sbornik: Mathematics, 1(4):457–483, 1967.
  * [Nec07] Ion Nechita.  Asymptotics of random density matrices.  Annales Henri Poincaré, 8(8):1521–1538, 2007.
  * [Per96] Asher Peres.  Separability criterion for density matrices.  Physical Review Letters, 77(8):1413, 1996.
  * [Rud03] Oliver Rudolph.  On the cross norm criterion for separability.  Journal of Physics A: Mathematical and General, 36(21):5825, 2003\.
  * [SŻ04] Hans-Jürgen Sommers and Karol Życzkowski.  Statistical properties of random density matrices.  Journal of Physics A: Mathematical and General, 37(35):8457, 2004\.
  * [ŻPNC11] Karol Życzkowski, Karol A Penson, Ion Nechita, and Benoit Collins.  Generating random density matrices.  Journal of Mathematical Physics, 52(6):062201, 2011.
  * [ŻS01] Karol Życzkowski and Hans-Jürgen Sommers.  Induced measures in the space of mixed quantum states.  Journal of Physics A: Mathematical and General, 34(35):7111, 2001\.
  * [ŻS03] Karol Życzkowski and Hans-Jürgen Sommers.  Hilbert-schmidt volume of the set of mixed quantum states.  Journal of Physics A: Mathematical and General, 36(39):10115, 2003\.