Entanglement criteria for the bosonic and fermionic induced ensembles

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.


Introduction
In Quantum Mechanics, quantum states are modeled by vectors in a complex Hilbert space. In the presence of an environment, one uses the formalism of density matrices B Stephane Dartois stephane.dartois@unimelb.edu.au Ion Nechita ion.nechita@univ-tlse3.fr Adrian Tanasa ntanasa@u-bordeaux.fr 1 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) [28], the Bures measure (motivated by statistical considerations) [15], random matrix product states (motivated by condensed matter theory) [13], and random graph states (which encode a given entanglement structure) [9]. 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 antisymmetric 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 [28], 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 [17,24] of these random density matrices in the large N limit. For the bipartite fermionic ensemble, we find any 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 (for two bosons) is computed in Theorem 5.1, which we state informally here: Theorem 1.1 Let ρ N be a N 2 × N 2 element from the bosonic ensemble. The spectrum of ρ 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 [10,25] 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 [1]. 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 Sect. 2, we recall some basic definitions from quantum information theory and from the theory of symmetric operators. In Sects. 3 and 4, we introduce, respectively, the fermionic and the bosonic induced ensembles. In Sect. 5, we state the main result on the asymptotic spectrum of a random bosonic quantum state. In Sect. 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 Sect. 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 Sect. 12 for the proof of the main result. We conclude in Sect. 13, where we also lay out some directions for future research.

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 modeled by density matrices. Each quantum system comes with a complex Hilbert space H, which we shall consider here to be finite dimensional: H ≈ C N ; we say that the quantum system has N degrees of freedom. A density matrix is a positive semidefinite operator acting on H, with unit trace: ρ ∈ M N (C), ρ ≥ 0, Tr ρ = 1.
The Hilbert space corresponding to two quantum systems is obtained by taking the tensor product of the individual Hilbert spaces: A bipartite quantum state ρ AB is said to be separable if it can be written as a convex mixture of product states: for density matrices α i , β 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 [14]. For this reason, several computationally efficient necessary or sufficient conditions for entanglement have been developed [18,Section VI.B].
In this paper, we shall be concerned with the most famous entanglement criterion, the positive partial transpose criterion [17,24]. (We shall also briefly discuss the realignment criterion in Sect. 4.) The starting observation is that the partial transposition of a separable state is positive semidefinite: Hence, if the partial transposition of a quantum state ρ 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 [18, 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 th tensor power H ⊗r of H. A natural basis of H ⊗r is given by the set where the e i 's form a basis of H, and [N ] denotes the set {1, 2, . . . , N }. Let σ ∈ S r , where we denote by S r the symmetric group on r elements. We define the action of σ on H ⊗r by extending linearly its action on basis elements, so that We now define two orthogonal projectors where ε(σ ) is the signature of the permutation σ . P s,r are projectors on the symmetric subspace Sym r H and P a,r are projectors on the antisymmetric subspace (or r th exterior power) ∧ r H. We define the r th symmetric power (resp. the r th exterior power) of H as the subspace P s,r (H ⊗r ) (resp. P a,r (H ⊗r )). We recall The orthogonal projection on the symmetric subspace can be written as [16] For example, in the bipartite case, we have P s,2 = 1 2 (I + F) and P a,2 = 1 2 where F is the flip (or swap) operator F x ⊗ y = y ⊗ x while I is the identity permutation and is represented by a N 2 × 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 )/2 and P a,2 A = (A − A )/2. Note that this corresponds to symmetrizing and antisymmetrizing matrices; in the complex case, we do not obtain, however, selfadjoint 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 [19] 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 ≤ r ≤ k): Proof We use the integral representation from Eq. (4):

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 [13,23,[27][28][29]. 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 [28]. The general setting is as follows. Consider a random pure state |ψ ∈ H ⊗ H E , uniformly distributed on the unit sphere the tensor product between the Hilbert space of interest H and an auxiliary Hilbert space H E , called the environment. Define a random density matrix ρ := Tr E |ψ ψ| by tracing out the environment. In this way, we obtain a random matrix acting on H; the dimension of the environment, dim H E , is a parameter of the model. It was recognized in [23,26] that the induced measure can be understood as a normalized Wishart measure: where Y = GG * , for G : H → H E a random Gaussian matrix (an element of the Ginibre ensemble). In this paper, we shall employ several times the Wick-Isserlis formula, which allows one to compute Gaussian integrals in terms of covariances [20]; see also [8] for the graphical interpretation used in this paper. 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 ∈ ∧ r H ⊗ H E . We identify the space of antisymmetric tensors where P a,r is defined in Sect. 2 to be the projector on the antisymmetric subspace of H ⊗r . We use this identification to construct the antisymmetric random Gaussian tensor from a complex random standard Gaussian tensor G ∈ H ⊗r ⊗ 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 := Tr E (G a G * a ). Note that, using the identification from (5), we can also write A = P a,r 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 H.
The fermionic induced ensemble of parameters (N , M, r ) is the set of random bipartite density matrices of size N r defined by where A is an antisymmetrized Wishart random matrix from Eq. (6). Matrices from the fermionic induced ensemble are supported on the antisymmetric subspace of (C N ) ⊗r .
Let us now comment on our choice of model. As already detailed in 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 Rb (which are fermions), or states of photons or 4 He atoms (which are bosons, whose ensemble we introduce later in sect. 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 H ⊗r , while no symmetry is imposed on the tensor product with the environment system 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 We shall consider an asymptotic regime where both N , M → ∞ in such way that M = cN 2 for some fixed constant c ∈ (0, ∞). This choice of scaling is the standard one needed for the convergence of the standard Wishart random matrices to the Marcenko-Pastur distribution: with a = (1 − √ c) 2 and b = (1 + √ c) 2 . The limiting eigenvalue distribution of the properly normalized random matrix A is the same as that of a non-symmetrized Wishart matrix. Below, E denotes the expectation functional with respect to the distribution of the random matrix W .
Proof After restricting it to its support (∧ 2 (C N )), the random matrix A is indeed a Wishart matrix of parameters (N a [2], M). The result follows from the standard theorem of Marcenko and Pastur, see [22], [ [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 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 ∼ cN 2 ). Before we commence, recall that the maximally entangled state is the unit norm vector where {|i } i∈ [N ] is the canonical orthonormal basis of C N . We also write ω = | | for the corresponding density matrix.

Proposition 3.3 The average of the partial transpose A is given by
The eigenvalues of EA are as follows: Proof The proof is a direct computation. Using Wick-Isserlis theorem, we compute |i j ji|.
Applying the partial transpose to the flip operator, we deduce that Writing P a,2 = (I − F)/2 yields the formula for EA . The exact form of the spectrum follows from the fact that ω is a rank one projection.
We take away from the above proof that the maximally entangled vector | is an eigenvector of EA with a negative eigenvalue for all N ≥ 2, c > 0, asymptotically of the form − c 2 N 3 + O(N 2 ). The next theorem shows that the same holds asymptotically for the random matrix A .
then lim N →∞ Eδ N = 0. In particular, the random variable δ N converges in probability toward zero as N → ∞.

Proof
We start with a more explicit formula of δ N : Computing the average, we obtain, thanks to Proposition 3.3, We then use Wick-Isserlis theorem to compute E |(A ) 2 | , and we find an expression involving the partial trace Tr 1 P a,2 , We compute Tr 1 (P a,2 ) = 1 2 (N − 1)I N from which we deduce showing that Eδ N → 0. We then use Markov's inequality to bound the probability that δ N is positive. Indeed, for all a > 0, we have concluding the proof.
These two results are here sufficient to conclude that the random matrix A is not positive semidefinite, as indeed it has a negative eigenvalue at scale N 3 . In particular, this implies that the random matrix A is entangled; this is in agreement with the following deterministic result, stating that any positive semidefinite operator supported on the antisymmetric subspace is entangled. Proposition 3.5 Let X be a bipartite operator of size N 2 such that X ≥ 0. Then, where P a is the projection on the antisymmetric subspace of C N ⊗ C N . In particular, all operators supported on P a are not PPT.

Proof
The result is a simple linear algebra computation: where we have used We have seen that any (random) matrix supported on the antisymmetric subspace is not PPT; hence, it is entangled. This is in contrast to the result of Aubrun [3] 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 ∧ 2 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.

The bosonic induced ensemble
In this section, we shall consider a setting analog to sect. 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 ∈ Sym r (H) ⊗ H E . In practice, we shall use the identification where we recall that P s is the orthogonal projection on the symmetric subspace of H ⊗ H. The tensor factor 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 H in contact with an environment 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 = Tr E (G s G * s ). Again, using the identification from (15), we can also write where G is a random (complex) Gaussian vector in the Hilbert space H ⊗r ⊗ H E ; see Fig. 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.
where W is the symmetrized Wishart random matrix from (16). Matrices from the bosonic induced ensemble are supported on the symmetric subspace of (C N ) ⊗r .
Mimicking the construction of Sect. 3, we set the dimension of the Hilbert spaces as follows 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 antisymmetrized matrix.

Proposition 4.2
The random matrix W converges in moments, as N → ∞, M ∼ cN 2 , to a Marcenko-Pastur distribution of parameter 2c: Proof After restricting it to its support (Sym 2 (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 Marcenko and Pastur, see [22], [6, Theorem 3.6], or [11,Proposition 2.4]. The parameter of the limiting distribution is given by

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 M = (α, γ ) where α is a permutation of S p and γ = (1 2 3 · · · p) is the full cycle: where F(M) is the number of faces of M and V • (M) its number of white vertices. We refer the reader to [11,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 α, or, equivalently, over planar bipartite maps, or dessins d'enfants, M with one black vertex. The limiting moments are values of the Narayana polynomial:

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 × N 2 random matrix, supported on the symmetric sub- Its partial transpose is defined as the action of the transposition operator on the second tensor factor: The next theorem is the main result of our paper, providing a description of the spectrum of W in the large N limit. Interestingly, W has a large eigenvalue of order N 3 , and N 2 − 1 eigenvalues of order N 2 .

converges in moments to a semicircular distribution of mean c/2 and variance c/4.
We recall here that the semicircular distribution with average m and variance σ 2 is given by Let us record here an important corollary for the theory of quantum information. The partial transposition criterion [17,24] states that a separable (i.e., non-entangled) density matrix ρ has a positive semidefinite partial transposition: ρ ≥ 0. This fact is mostly used in the reverse direction, as an entanglement criterion: Given a density matrix σ , having a partial transposition with at least one negative eigenvalue (σ 0) is entangled. We have thus the following result, regarding the bosonic induced ensemble defined in Definition 4.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 ρ N does not exclude the absence of negative outliers. In the standard (non-symmetrized) Wishart case, Aubrun [3] 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 behavior of outliers of ρ N ; at the same time, in this work, we cannot study them directly because of the large eigenvalue of W , 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 Sect. 12. We start, in the spirit of Sect. 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 Sect. 3, we compute the average of the partially transposed matrix W 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 is the canonical orthonormal basis of C N .

Proposition 5.4 The average of the partial transpose W is given by
The eigenvalues of EW are as follows Proof The proof is a direct computation. Using Wick-Isserlis theorem, we compute EW = M 2 (I N 2 + F). Applying the partial transposition on the flip operator |i j ji|, Writing P s = (I + F)/2 yields the formula for EW . The spectrum follows from the fact that | | is a rank one projection.
We note that | is an eigenvector of P s with eigenvalue (N + 1)/2. Hence, the average EW has an eigenvalue M(N + 1)/2 ∼ cN 3 /2 with eigenvector | . The fact that the average of the random matrix W has a large eigenvalue does not imply, in general, that W itself also does. Due to the concentration of measure phenomenon, we show next that this is indeed the case, asymptotically.

Theorem 5.5 The maximally entangled vector | is an approximate eigenvector of the random matrix N −3 W with corresponding approximate eigenvalue c/2: If
then lim N →∞ Eε N = 0. In particular, the random variable ε N converges in probability toward zero as N → ∞.
Proof Manipulating the expression of ε N , we have Computing the average and using Proposition 5.4, we have We compute E |(W ) 2 | using the Wick-Isserlis theorem and find, using Lemma 2.1, 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 [10,25]. The realignment criterion states that any separable bipartite density matrix ρ satisfies the following inequality: where · 1 denotes the Schatten 1-norm or the nuclear norm (i.e., the sum of the singular values of a matrix), while ρ R denotes the realignment (or the reshuffling, see [7, Chapter 10.2]) of a matrix ρ, defined algebraically by or graphically in Fig. 2.
The key observation here is that the realignment and the partial transposition of a matrix are related as follows: In the case of a (non-normalized) random matrix W from the induced bosonic ensemble, entanglement follows from the inequality W 1 > Tr W . Note that, asymptotically as N → ∞ and M ∼ cN 2 , the trace behaves as Tr W ∼ (c/2)N 4 , while Theorem 5.1 states that 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: But the above happens precisely when the probability measure SC c/2, √ 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 0 = 4 [3]) is strictly larger than the threshold for the realignment criterion (c R 0 = (8/3π) 2 ≈ 0.72 [2]), proving that, in some sense, the former criterion is stronger than the latter. 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.
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 → j} for i, j ∈ V to better picture the ordering. Possibly non-distinct elements in the pairs mean the digraphs can have loops.
Note that incoming edges of a vertex i are the edges of the form {a → i}, while outgoing edges are the edges of the form {i → 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 with j k = |{covers of G in k 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.
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 [12, Definition 2.1] and below; the case we are interested in is the Eulerian case). These skein relations read graphically 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 [21] that we show here. We give a proof below, using the graphical skein relations.
Proof We have using the skein relations described above at the cut vertex v 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 choices 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, Therefore, 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 (20) leads to 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 [1], defined here according to [5]. For a (non-oriented) graph 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.
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 = . . . a . . . b . . . a . . . b . . ., see [1] 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) In [1], 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
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 , . . . , 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. (26)); hence, deg J (G i ) ≤ 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 We then use the above bound to show the next lemma. Proof Assume that a connected digraph G has only k < p out of p vertices which are cut vertices. Then using skein relations (19), 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 and G . 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, we have that deg J (G ), deg J (G ) ≤ p. We also know from equation (24) 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 σ 1 ,σ 2 , for σ 1 , σ 2 ∈ S p , the 2-in/2-out digraph on p vertices whose (oriented) edge (multi-)set is where we used the same notation for the permutations σ 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 σ 1 , σ 2 , . . . , σ k ∈ S p such that 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.
Assume γ = (1, 2, 3, . . . , n) ∈ S n is the full cycle permutation, then In particular, Proof The proof follows from the fact that G γ,γ −1 is of the form depicted in Fig. 3. Then it is sufficient to notice that at each vertex of G γ,γ −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 δ 2,k in the above formula). 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 n k choice of states leading to a cycle cover with k cycles of G γ,γ −1 . It follows that the degree of J (G γ,γ −1 ) is n. This last result can also easily be obtained by constructing an interlace graph H γ,γ −1 from an Eulerian cycle C γ,γ −1 of G γ,γ −1 . It is easy to exhibit such an Eulerian cycle: C γ,γ −1 = 1, 2, 3, . . . , n, 1, n, . . . , 3, 2. A pair (1, q) is clearly interlaced for any q = 2, . . . , n while any other pair (a, b) for a and b = 1 is not interlaced. Thus, H γ,γ −1 is the star graph with the vertex 1 at its center, and thus, max X minor of H dim 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).
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 α and β, once the moves f (i) have been applied at each vertex i ∈ [p], see (21).

Moments of W 0
In this section, we shall express the moments of the random matrix W (see Fig. 4) as a sum of circuit counting polynomials, allowing us to determine their asymptotic behavior. 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.
We first establish the exact form of the moments of W in terms of J polynomials.

Proposition 7.1
The moments of the random matrix W ∈ M N 2 (C) from Eq. (17) are given by: Proof The proof is a rather simple application of the Wick formula, in its tensor network incarnation [8]. We need to evaluate the expectation of the trace of the product of p copies of the matrix from Fig. 4. By the Wick theorem, the result is a sum over permutations α ∈ S p of diagrams obtained by connecting the ith G box to the α(i)th G * box Consider now the 2-in/2-out digraph G γ α,γ −1 α having edges (see (28) for the general case) {i → γ (α(i))} i∈ p {i → γ −1 (α(i))} i∈ 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 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 α = id and α = (12): where, for the case α = id, we have used Proposition 6.9 with n = 2.
To analyze larger moments, we need the following technical result.
We describe next the asymptotic behavior of the larger moments of W .
In conclusion, all the terms with α = id are subdominant with respect to the term α = id, finishing the proof.
To finish this section, note that the asymptotic behavior from (35), corresponding to the moments p ≥ 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).

t-channel random matrix and graphical representation
In this section, we introduce the t-channel 1 random matrix that we denote W t . We do so as it is better suited for our aims in Sects. 9, 10 and 11. We define W t component-wise by and we provide a diagrammatic representation of the above definition below In the above diagrammatic representation, the red arc represents the symmetrizer P s . The black square represents the complex conjugateḠ of the Ginibre matrix G of Sect. 4 while the white square represents G itself.Ḡ 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: 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 C| By the projector property, we have P 2 = P . We want to study the sequence of moments E(Tr(Q p )) of the random matrix In order to use graphical methods for tensor networks evaluation (see [4] for an overview of tensor networks techniques and ideas), we will use the following relation This is easily shown by using the fact that F 2 = I , F W t F = W * t and [P , F] = 0. These moments have a graphical representation as ladder diagrams introduced in Sect. 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 , 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 Sect. 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 Sects. 9 and 10.
We denote P | {1, . . . , p} an interval partition of {1, . . . , p}. That is, P is a set of subsets P i ⊆ {1, . . . , p} of P such that i P i = P and P i are subintervals of {1, . . . , p} seen cyclically. For instance,  Expanding the projectors P appearing in the definition of Q, we have 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 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 (43) is Tr((W ) p ). We use this remark later to deduce the diagrammatics of E Tr((W ) p ) .

Diagrammatics for the moments of Q and W t
We start by considering the diagram representation of expression involving W t in the expansion (42). They are easily obtained by stacking building blocks of equations (37),(38). We have If p ∈ 2N + 1 is odd, we have the slightly twisted representation 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 α ∈ S p . Indeed, since G is normally distributed, with vanishing pseudo-variance, Wick pairings only pair instances of G with instances ofḠ. We represent such a pairing by adding oriented dotted edges from black squares to white ones in diagrams of the type appearing in equations (44) and (45). Each Wick pairing can be indexed by a permutation in α ∈ S p . Indeed the black box number i is paired with the white box j if and only if α(i) = j. This translates diagrammatically as a dotted edge between black square i and white square j, oriented from black to white. See examples in Fig. 6. 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 α(i) = j, then the symmetrizer associated with 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 two 2 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 detail this evaluation combinatorially. Since each vertex represents twice a symmetrizer P s = 1 2 (I + F), we can evaluate the contraction by deciding for a state 3 that is the assignation of either the identity I : 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 1 2 p , with p the number of vertices. This last factor just takes into account the normalization factor 1 2 of the symmetrizer.
The digraph resulting from the Wick pairing α 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 σ 1 = γ α, σ 2 = γ −1 α with γ = (1, 2, . . . , p). Indeed, a vertex labeled i corresponds, by definition, to a pair (i, α(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 α(i). Note that the edges exiting from vertex i of the digraph are adjacent to the black boxes of the gray pairs γ (α(i)) and γ −1 (α(i)) that is in the digraph they are adjacent to the vertices γ (α(i)), γ −1 (α(i) oriented as i → γ (α(i)) and i → γ −1 (α(i)). Hence, σ 1 = γ α, σ 2 = γ −1 α. See the local construction of the vertices in Fig. 7.
We follow a similar train of thoughts for the terms in the sum from equation (42). Indeed, the product (W t ) |P f ,i | mod 2 (W * t W t ) |P f ,i |/2 is also represented diagrammatically by stacking building blocks of equations (37) and (38). The Fig. 7 On this figure, we show how the vertices of the digraph are constructed from the pairing α by identifying the black box i with the white box α(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 (Color figure online) only difference being that the vector | has the following representation Hence, we have the diagrammatic below 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, and to each P i ∈ P | {1, . . . , p} we associate the bijection C P i : P i → {1, . . . , |P i |} defined by where P i (q) is the q th element of P i . Then we define for each P i the permutation γ a,P i = C −1 P i γ a,|P i | C P i and γ b,P i = C −1 P i γ b,|P i | C P i . Finally, we construct yet another two permutations, a,P , b,P := i γ a,P i , i γ b,P i . We have several figures to illustrate the relations with the diagrams. See Fig. 8 for γ a,|P i | and γ b,|P i | . See also Fig. 9 for an example of permutations a,P , b, p . Then we have as a consequence of Wick theorem that: where we recall that i |P f ,i | = 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 α(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 α(i) while these edges are adjacent to the black boxes a,P (α(i)), b,P (α(i)), which correspond to the vertices a,P (α(i)) and b,P (α(i)) in the digraph. This is sufficient to conclude that we have the relation (53).
We summarize the content of the above discussion by the proposition below: Reformulating the above proposition and using Remark 8.1, we can recover the moments of the random matrix W from Proposition 7.1 We use the shorthand notation, for α ∈ S p w(α, γ ) := N 2#α c #α 2 p J (G γ α,γ −1 α ).

Circuit polynomial techniques for E 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 α contributing to the asymptotics of the random matrices W and Q.
We have the following bound Proposition 10.1 Let α ∈ S p , then there exists R α, f > 0, that can depend on α and on f but not on N with K (G a,P f α, b,P f α ) the number of connected components of the corresponding digraph. Moreover, we have: Proof We now come to the second bound in the case of 2-in/2-out digraphs G a,P f α, b,P f α . The main difference lies in their number of connected component. Indeed, a,P f b,P f ∈ a,P f α, b,P f α acts transitively on the sets P f ,i separately, but is not transitive on {1, . . . , p}. Thus, the number of connected components is bounded from above by |P f |. Hence, we have Again, using the relation between the circuit polynomial and the interlace polynomial on each connected component of G a,P f α, b,P f α , we obtain We also show the two technical Lemmas 10.2 and 10.3. They prove useful in the last section, to determine the limiting empirical distribution of eigenvalues of Q, where # i α is the number of cycles of length i of α. If K (G γ α,γ −1 α ) = 2, Similarly, let e be a positive integer such that 2#α + p + K (G a,P α, b,P α ) = 2 p + e . Then, In particular, if α is not allowed to have fixed points, then e ≤ 0. If we have equality e = 0, then # i α = 0 for all i ≥ 3.
Proof We first focus on e. We denote # i α the number of cycles of length i of α. We have the following relations From those, we deduce We have two cases, first assume that K (G γ α,γ −1 α ) = 1. Then, a consequence of proposition 7.2 is that α is allowed to have fixed points; thus, Now, assume that K (G γ α,γ −1 α ) = 2. A consequence of proposition 7.2 is that α is not allowed to have fixed point (as those are not parity changing). More generally, α is not allowed to have cycles of odd length. Therefore, The above inequality is saturated when α 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 . We have where the second and last term of the right-hand side is bounded from above by 0.
As we understand from the above proposition, fixed points and cycles of length 2 of α play an important role. Indeed, the scaling of a given Wick pairing α can exceed N 2 p+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 α ∈ S p has only cycles of length 2 and is parity changing (i.e., α(i) mod 2 = i mod 2). Assume also that there exists two cycles of α of the form (a, i + 1) and We use the parity changing assumption together with the constraint that α 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 proved from similar considerations. However, the above form of the lemma will be sufficient for our purposes.
Proof The parity changing assumption ensures that G γ α,γ −1 α 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 α, the isomorphism ϕ is set to ϕ(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 γ α,γ −1 α 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 α had more crossings), this means vertices a and i cannot become cut vertices by changing α while keeping the cycles (a, i + 1) and (i, b) fixed.
Denote G (a) γ α,γ −1 α the connected component containing vertex a. We focus on this connected component for now. By Lemma 6.7, we have deg J (G (a) γ α,γ −1 α ; N ) ≤ 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. (19) on a we are left with two graphs G and G such that Both G and G have a maximum of p/2 − 1 vertices. Therefore, from Lemma 6.6 By the isomorphism argument, we have the same bound on the degree of J (G 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.

Tensor network evaluation of E 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 . In particular, we do not expand the projectors P appearing in the expression of Q. In that case, for each Wick pairing α we obtain a contraction of two types of tensors, P s : C N ⊗ C N → C N ⊗ C N as before and P : C N ⊗ C N → C N ⊗ C N which is the projector on the complement of span{| } introduced earlier in Eq. (39).
We work here with a graphical representation that is very similar to the ones introduced in previous paragraphs. Here, we represent Tr(Q p ) directly, without expanding. We do so by stacking the building blocks introduced in (37), (38) together with an These building blocks are stacked in between (37) and (38) blocks to represent the operator P . The diagram thus obtained is a tensor network representing Tr(Q n ) To compute E(Tr(Q p )) we proceed as before, with each Wick pairing identifying black box i with white box α(i) and leading to a vertex labeled i. This vertex represents a P s operator. However, the resulting oriented graph also contains -boxes that represents P . This leads to a slightly more complicated connection pattern. We show in Fig. 11 the local structure around a vertex i. The only change is that now vertices are only connected to -boxes. Locally it is simple to check that vertex i is connected to -boxes i, γ (i) via ingoing edges, while it is connected to -boxes α(i), γ (α(i)) via outgoing edges. This is due to the fact that the white box α(i) (of which adjacent edges are outgoing) is connected toboxes α(i), γ (α(i)) while the black box α(i) (of which adjacent edges are ingoing) is connected to -boxes i and γ (i). In the rest of this paper, we denote by T α the tensor network (and its evaluation) corresponding to a pairing α constructed in the way described above.
Using this tensor network construction, we have Note that we have the relation (75) In the next paragraph, we prove reduction relations for the evaluation of such tensor networks at fixed points, and simple transpositions of α. 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 σ = c 1 c 2 . . . c #σ ∈ S p . We denote by σ ÷ c q := c 1 c 2 . . .ĉ q . . . c #σ the permutation acting on the set {1, . . . , p} \ {i : i ∈ c q }. We have the following reduction property for tensor network T α such that α have the corresponding feature Proposition 11.1 Let α ∈ S p . Then, • if α as a fixed point, that we denote i, then whereα is the permutation satisfying α =α(i) in cycle notation. In particular, if α has several fixed points, consider α R the reduced permutation acting on {1, . . . , p}\ Stab(α), then one has • if α contains a cycle of the form (i, i + 1), then where (N ) = 1 4 (N 2 + 2(N − 1 N )).
Proof We proceed to the proof by diagrammatic manipulation of the tensor network. Fixed points i of α can be treated locally by realizing that the corresponding tensor network always has the same structure around vertex i (see the left-hand side of equation (79)) where we used the fact thatΩ Together with the fact that P 2 = P , this shows that The second statement can be obtained by looking at what happens locally on the vertices i and i + 1 when α has a cycle of the form (i, i + 1). One has the following diagrammatic relation obtained by using the definition of P . The leftmost term of the right-hand side of equation (82) leads after further expanding the two black vertices while the rightmost term of the right-hand side of equation (82) can also be extended further. Together with using again the relation (80), we have Putting everything together, we obtain Thanks to this proposition, we can consider permutations without fixed points by just remembering that for every α ∈ S p , there is a reduced permutation α R without fixed points such that equation (76) 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 (75), this is due to cancelations between terms in the sum over words f . The tensor networks formulation allows to formulate these cancelations compactly. One particular case is if α ∈ S p is the identity permutation. Then we have Another particular case is when α ∈ 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 α which is in NC 2 ( p − 2). Using this property, we can recursively reduce all the cycles of such a permutation and we obtain

Proof of the main result
We now have all the tools we need to come to the proof of the main result (Theorem 5.1). This proof follows the steps: (1) Determining the limiting spectrum of Q. This is Theorem 12. Limiting spectrum of Q. We tackle the first step of our three steps plan of proof. Assume that the permutation α do not contain fixed points. By virtue of Lemma 10.2, we know that e ≤ 2 and e ≤ 0 for such α. Moreover, these inequalities can be saturated if α has cycles of length 2. This implies that there exists a positive constant R, independent of N , which can be attained if α 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 Indeed, we have, graphically, Use now the fact that P s is a projection, and thus, Tr α,α (P s ) = (Tr P s ) #α = N [2] #α .
Interlacing of eigenvalues and asymptotic of W . 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 We use this theorem to recover the behavior of the limiting eigenvalues of W from the limiting spectrum of Q and the large eigenvalue of W . 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 behaves, asymptotically as N → ∞, as c 2 N 3 . To this end, we shall compare it with two quantities: the overlap |W | and wih the moments of the random matrix W .
On the one hand, we have (we denote λ max = λ 1 for clarity), for all p ≥ 2, establishing the first claim. In order to prove the second part of the statement, concerning the lower N 2 − 1 eigenvalues of W , we shall make use of Theorems 12.1 and 12.3. Indeed, if we denote by μ 1 ≥ · · · ≥ μ N 2 −1 the (nonzero) eigenvalues of the matrix N −2 Q and by λ 1 ≥ λ 2 ≥ · · · ≥ λ N 2 −1 ≥ λ N 2 those of N −2 W , by Cauchy's interlacing theorem, we have On the other hand, we know by Theorem 12.1 that the empirical distribution of the μ eigenvalues converges to the desired shifted semicircle distribution. Hence, in order to conclude, we need to control the smallest eigenvalue λ N 2 and show that it does not contribute asymptotically to moments. We have

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. antisymmetric 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 nonsymmetric case discussed in [3]. 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 toward 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 [3]; 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 ≥ 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 antisymmetric) 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.