Stephane Dartois
email: stephane.dartois@unimelb.edu.au

Luca Lionni
email: luca.lionni@yukawa.kyoto-u.ac.jp

ON THE JOINT DISTRIBUTION OF THE MARGINALS OF MULTIPARTITE RANDOM QUANTUM STATES

Keywords: Wishart ensemble, random quantum state, random tensor, marginal of a quantum state, combinatorial map, planar map

We study the joint distribution of the set of all marginals of a random Wishart matrix acting on a tensor product Hilbert space. We compute the limiting free mixed cumulants of the marginals, and we show that in the balanced asymptotical regime, the marginals are asymptotically free. We connect the matrix integrals relevant to the study of operators on tensor product spaces with the corresponding classes of combinatorial maps, for which we develop the combinatorial machinery necessary for the asymptotic study. Finally, we present some applications to the theory of random quantum states in quantum information theory.

1. Introduction 2. The limiting eigenvalue distribution of random density matrices 2.1. Random density matrices 2.2. The limiting eigenvalue distribution 2.3. A combinatorial map version of the proof 2.4. Free cumulants and maps 3. The four-partite case: joint distribution of two marginals 3.1. Exact expression for the moments 3.2. The balanced asymptotical regime 3.3. Comparing with the Fuss-Catalan matrix model 3.4. The unbalanced asymptotical regime 4. The general multipartite case 4.1. Exact expression for the moments 4.2. The balanced asymptotical regime 4.3. The unbalanced asymptotical regime References

Introduction

The Wishart ensemble was historically the first probability distribution on matrices which was studied [Wis28]. The main motivation for Wishart was statistics; later, Wigner [Wig55] modelled complex, analytically intractable Hamiltonians in nuclear physics with Hermitian random matrices, and the field of Random Matrix Theory [START_REF] Madan | Random matrices[END_REF][START_REF] Greg W Anderson | An introduction to random matrices[END_REF] was born. Nowadays, random matrices play a significant role in many sub-fields of mathematics, such as operator algebras, combinatorics and algebraic geometry, integrable systems and partial differential equations, as well as in other disciplines such as theoretical physics or telecommunication. In this paper, we are motivated by a recent application of random matrices to Quantum Information Theory [START_REF] Michael | Quantum computation and quantum information[END_REF].

The mathematical formalism of Quantum Information Theory is constructed upon the central notion of quantum states, also known in the physical literature as density matrices. These are positive semidefinite d × d complex matrices, normalized to have unit trace; here, d is the number of degrees of freedom of the quantum system under consideration. Density matrices model "open quantum systems", that is quantum systems which interact with an environment (which, most of the times, is too complicated to be taken under consideration). Isolated systems, (which are called "closed") are modeled traditionally by unit vectors in C d , which we choose to identify with the rank-one projections on the corresponding vector space; these projections are the extremal points of the convex set of density matrices.

One might want to study random quantum states for several different reasons. Foremost, we would like to understand what are the typical mathematical (or physical, or even informationtheoretical) properties of a typical state, where typical should be understood as randomly distributed with respect to some natural (or physically relevant) probability distribution. Another reason one would like to understand random density matrices comes from the empirical observation that, in situations where explicit examples satisfying some desired properties are hard to come by, one should simply pick the sought-for object at random; in many cases, with large probability, the random sample will have the desired properties. Random quantum states (and random quantum channels) have been a valuable source of (counter-)examples in Quantum Information Theory (see, e.g. the recent review paper [START_REF] Collins | Random matrix techniques in quantum information theory[END_REF]).

There is a large literature on random density matrices and their applications to quantum information theory, most of it focusing on spectral properties of one random matrix. In particular, the focus was on random states of single quantum systems and bipartite quantum systems (mostly related to the study of entanglement). In this work, we tackle a fundamentally different question:

Given a multipartite random density operator, what is the joint probability distribution of its different marginals?

Recall that, for a quantum state of a multi-partite system, the marginal of a subset of systems is the partial trace with respect to the complementary set of systems. We shall study the question above in different settings, of increasing generality, first with just 4-partite states, and then for quantum states with arbitrarily many subsystems. Depending on the relative sizes of the subsystems and on their rates of asymptotic growth, we shall exhibit two types of behavior. In one situation, where the growth rates of the system dimensions are the same, we prove that the whole set of (balanced) marginals are asymptotically free, meaning, in broad terms, that they behave like independent random matrices (although they might share one or more subsystems). In a different regime, where the dimensions of some of the subsystems are being kept fixed, we do not have asymptotic freeness, but we provide exact formulas for the limiting joint free cumulants, in terms of the types of marginals involved. We state next, informally, two of the main results of this paper, corresponding to the asymptotic regimes described above; for the more general and precise results, we refer the reader to Theorem 4.13 and, respectively, Theorem 3.16.

Theorem. Let ρ (N ) ∈ M N 2r (C) be a random pure state of 2r-partite quantum system, where each subsystem is N -dimensional. Then, the 2r r marginals ρ (N )

S , with S ⊆ {1, 2, . . . , 2r} a set of cardinality r, are asymptotically free: their (rescaled) joint distribution converges in moments to that of 2r r free Marčenko-Pastur elements. Equivalently, their limiting joint distribution is the same as that of 2r r independent copies of, say, ρ

(N )
{1,2,...,r} . In the unbalanced case, we have, informally, for 4-partite systems, the following result.

Theorem. Let ρ N ∈ M N m 2 M (C) be a random pure state of 4-partite quantum system H ABCD , where dim H A = N , dim H B = dim H C = m and dim H D = M ∼ cN . Assuming m, c are fixed constants, the (rescaled) joint distribution of the marginals (ρ AB , ρ BC ) converges, in moments, as N → ∞, to a pair of non-commutative random variables (x AB , x AC ) having the following free cumulants:

κ(x f (1) , x f (2) , . . . , x f (p) ) = cm -alt(f ) , where f ∈ {AB, AC} p is an arbitrary word in the letters AB, AC, and alt(f ) is the number of different consecutive values of f , counted cyclically: alt(f ) := |{a : f (a) = f (a + 1)}|, where f (p + 1) := f (1).

The contribution of our paper is threefold. First, we introduce, in full detail, the notions of combinatorial maps relevant for the random matrix computations we perform, and we extend them to the case of matrices having a tensor product structure. Our presentation starts at a basic level, gradually adding layers of complexity, and can be used by readers with a quantum information background as an introduction to the subject. We develop the necessary combinatorial techniques to deal with the types of maps appearing in our study (combinatorial maps with vertices of two colors and edges of 2r colors). Secondly, we contribute to the theory of random matrices and free probability by computing the limiting distribution of a family of random matrices obtained as marginals of a unique random object. Although, globally, the random matrix model is standard (Wishart matrices), taking (intersecting) marginals (i.e. partial traces) and considering their joint distribution is new; for this reason, in order to emphasize the importance of the tensor product structure of the Hilbert space, we shall call the random matrices we study Wishart tensors. We prove asymptotic freeness in the balanced case, in a very general setting, and obtain the limiting free cumulants in the unbalanced setting; the explicit form of the free cumulants (see the second informal theorem above) is very interesting, involving a parameter counting the number of different consecutive letters appearing in the respective word. Thirdly, from the point of view of quantum information theory, our study shows that the marginals of a random pure quantum state behave independently in the balanced case and in the large N limit: the moment statistics of the whole set of marginals are the same as if the marginals were independent. The situation is different in the unbalanced case: there is a strong correlation between, say, the marginals ρ AB and ρ AC of a pure random 4-partite quantum state ψ ABCD when dim H A dim H B = dim H C . Since the main focus of our paper is on random quantum states over Hilbert spaces with a tensor product structure, let us give now a brief survey of the literature on the subject, emphasizing the point of contact with our work. Physicists started working on ensembles of quantum states in the early '90s, when Page computed the average entropy of entanglement of a bipartite random pure state [Pag93]. The study of probability measures induced by metrics in the one-party case was initiated by Hall in [START_REF] Michael | Random quantum correlations and density operator distributions[END_REF] and developed by Sommers and Życzkowski in [ ŻS01, ŻS03, SZ03, S Ż04]; see also [OS Ż10] for the limiting eigenvalue distribution of the Bures ensemble. In the multi-partite case, the study of random tensors in quantum information theory was initiated [AHH12], where superpositions of random product states were investigated. Later, specific models of randomness were studied in [PGVWC07] in the case of random matrix product states and in [CN Ż10, CN Ż13] for random graph states. In [CDKW14], Christandl, Doran, Kousidis, and Walter studied the joint distribution of all the single particle marginals of a multi-partite quantum state in a very general setting, allowing for different distributions of the global state, and making use of the Duistermaat-Heckman measures from Lie theory. The current work is, to our knowledge, the first instance where the question of the joint distribution of the possibly overlapping marginals of a random quantum state is considered; we do so in the simplest framework, that of the Wishart ensemble.

In the framework of quantum information theory, this corresponds to considering a random pure quantum state on a multipartite Hilbert space, and tracing out some of the subsystems to obtain the marginals.

The paper is organized as follows. In Section 2 we recall some well-known results about Wishart matrices and random density matrices, which can be seen as marginals of bipartite Wishart tensors. This simple situation is also the occasion to introduce the machinery of combinatorial maps. We provide two proofs of the classical Marčenko-Pastur theorem, one using the language of permutations and their metric properties, and another one using combinatorial maps; the reader can see from this example how the two approaches mirror each other. In Section 3, we study in full detail the case of 4-partite Wishart tensors. We consider two different asymptotical regimes: a balanced regime, where the dimension of all the spaces are equal, and an unbalanced regime, where two of the four spaces have fixed dimension. We compute the limiting joint distribution of the two 2-marginals in both regimes: in the balanced case, we show that the marginals are asymptotically free, while in the unbalanced case we compute the (non-trivial) limiting mixed free cumulants. Finally, in Section 4, we study the general multipartite case. In the balanced case, we show again that the marginals are asymptotically free; in the other asymptotic regimes, we only have partial results: we list the different regimes that take place, but leave their detailed description for future work.

Acknowledgments. L.L. is a JSPS International Research Fellow. The work of S.D. was partially supported by the Australian Research Council grant DP170102028. I.N.'s research has been supported by the ANR projects StoQ (grant number ANR-14-CE25-0003-01) and NEXT (grant number ANR-10-LABX-0037-NEXT), and by the PHC Sakura program (grant number 38615VA). I.N. also acknowledges the hospitality of the Technische Universität München. The authors would like to thank the Institut Henri Poincaré in Paris for its hospitality and for hosting the trimester on "Analysis in Quantum Information Theory", during which part of this work was undertaken. S.D. and I.N. would also like to thank the organizers of the "QUATR-17" conference in Skoltech/Moscow, and especially Leonid Chekhov, for bringing together researchers in random tensor theory and quantum information theory.

The limiting eigenvalue distribution of random density matrices

In this section, we discuss the different ensembles of random density matrices from the literature, focusing on the induced ensemble, which will be the one we shall study in the later sections. We also compute the limiting eigenvalue distribution of (rescaled) random density matrices, using two different formalisms: an algebraic one, emphasizing the role of permutations, and a combinatorial one, featuring the theory of combinatorial maps. Although the two proofs given will be equivalent, we present both in full detail in order to introduce the main objects and to prepare the reader for the more complicated situations discussed in the later sections.

2.1. Random density matrices. To start, let us fix some notation. Density matrices with N degrees of freedom are represented by unit trace, positive semidefinite N × N matrices:

M 1,+ N := {ρ ∈ M N (C) : ρ ≥ 0 and Tr ρ = 1}
. This is a convex body, whose extreme points are rank one projections which we identify (up to a phase) with vectors x ∈ C N , x = 1, called pure states.

We consider first the canonical distribution on pure states, that is the Lebesgue measure on the unit sphere of C N . Integrating polynomials in the state's coordinates with respect to this measure is quite straightforward, see e.g. [START_REF] Gerald B Folland | Real analysis: modern techniques and their applications[END_REF]. The idea is to relate the spherical integral to a Gaussian one with the help of a change of variable to the polar coordinates, and then use Wick's (or Isserlis' [Iss18]) formula to evaluate the Gaussian integral; we recall this result next.

Proposition 2.1. Let X 1 , . . . , X k be a k-tuple of random variables having a joint (complex) Gauss-

ian distribution. If k is odd, then E[X 1 • • • X k ] = 0. If k = 2l is even, then E[X 1 • • • X k ] = p={{i 1 ,j 1 },...,{i l ,j l }} pairing of {1,...,k} l s=1 E[X is X js ].
Let us move now to ensembles on the whole set of N × N density matrices, M 1,+ N . Again, there is a natural candidate here, the normalization of the Lebesgue measure on the ambient space. It turns out however, that this measure is just a specialization of a 1-parameter family of probability distributions, called the induced measures. Introduced by Życzkowski and Sommers in [ ŻS01], these measures have the advantage of being interesting and natural both from the physical and the mathematical perspectives. Let us start with some motivating consideration from quantum physics. Assume the physical system we are interested in (which has N degrees of freedom) is not isolated, but coupled to an environment, having M degrees of freedom. In most physical applications, the environment is big and inaccessible, so we choose not to model it; in other words, if Ψ ∈ C N ⊗ C M is the (pure) quantum state describing jointly the system and the environment, we only have access to the state of the system ρ = [id N ⊗ Tr M ](ΨΨ * ) ∈ M 1,+ N . In the equation above, we assume that the vector Ψ is normalized, Ψ = 1. The main idea of [ ŻS01] is to consider Ψ uniformly distributed on the unit sphere of C N ⊗ C M ; this leads to the following definition.

Definition 2.2. The induced measure of parameters (N, M ) is the image measure of the uniform probability distribution on the unit sphere of C N ⊗ C M through the map

C N ⊗ C M → M 1,+ N Ψ → [id N ⊗ Tr M ](ΨΨ * ),
where ΨΨ * denotes the self-adjoint rank one projection on CΨ.

Importantly, the uniform (Lebesgue) probability measure on the convex body M 1,+ N is exactly the induced measure with parameters (N, N ) [ ŻS01, Section 3].

A detailed mathematical analysis of the induced measures defined above was performed in [Nec07], where it was emphasized that quantum states distributed along the induced measures are just normalized Wishart random matrices. To make this observation more precise, let us briefly remind the reader the definition of the Wishart ensemble (we refer the reader to [HP00, Section 4], [BS10, Chapter 3], [HT03] for detailed treatments of the Wishart ensemble from a random matrix theory perspective, and to the excellent book [AS17, Section 6.2.3] for a quantum information theory point of view). Let X ∈ M N ×M (C) be a Ginibre random matrix, that is a matrix having i.i.d. complex standard Gaussian entries (no symmetry is assumed here). For completeness, we rewrite formally what we mean by standard i.i.d. complex Gaussian entries. The entries of X form a set of N M complex numbers x i,j each with density

1 2iπ e -|x i,j | 2 dx i,j dx i,j , (1)
which rewrites in terms of the real part r i,j and imaginary part s i,j of x i,j as 1 π e -(r 2 i,j +s 2 i,j ) dr i,j ds i,j .

(2)

A Wishart matrix of parameters (N, M ) is defined as W = XX * . The relation to random density matrices from the induced ensemble has been made mathematically rigorous in [Nec07, Lemma 1].

Proposition 2.3. Let W be a random Wishart matrix of parameters (N, M ). Then,

ρ = W Tr W (3)
is a random density matrix distributed along the induced measure from Definition 2.2 of parameter (N, M ).

Although the normalization by the trace is a highly non-trivial (and non-linear) operation, in practice, for large random matrices, it does not pose technical difficulties. There are two reasons for this: first, in equation (3), the trace Tr W and the normalized density matrix ρ are independent random variables; this fact is similar to the result in classical probability which states that the norm and the direction of a (standard) Gaussian vector are independent random variables. The second reason which allows us to deal in a simple manner with the trace normalization is that the trace of a Wishart random matrix is a chi-squared random variable and thus concentrates very well around its average

E Tr W = N M (see [AS17, Exercise 6.40]) ∀t > 0, P [| Tr W -N M | > tN M ] ≤ 2 exp - t 2 N M 2 + 4t/3 . (4)
2.2. The limiting eigenvalue distribution. We compute in this section the eigenvalue distribution of Wishart matrices, and thus of random density matrices, in the large N limit. We shall work with the simpler model of Wishart matrices, and then translate the results to quantum states in Corollary 2.7. We shall present two proofs of the well-known convergence to the Marčenko-Pastur distribution, one using permutations and the other one using combinatorial maps; this will be the occasion to introduce these two proof techniques and to familiarize the reader with the main objects appearing in the respective theories.

In the limit of large matrix dimension (N → ∞), the behavior of Wishart matrices depends on the asymptotic ratio M/N . The most common situation is when M/N → c ∈ (0, ∞), in which case the matrix converges to the well-known Marčenko-Pastur distribution [MP67]. Although this result is straightforward and very well-known, we provide a self-contained proof in order to compare the approach of this section with the one in the next section.

Proposition 2.4. Let W N be a sequence of random Wishart matrices of parameters (N, M N ), where M N is an integer sequence with the property that

M N ∼ cN as N → ∞, where c ∈ (0, ∞) is a constant. The sequence N -1 W N converges, in moments, towards the Marčenko-Pastur distribution ∀p ≥ 1, lim N →∞ E 1 N Tr W N N p = x p dMP c (x),
where

dMP c = max(1 -c, 0)δ 0 + (b -x)(x -a) 2πx 1 [a,b] (x) dx, ( 5
)
with a = (1 - √ c) 2 and b = (1 + √ c) 2 .
We plot the density of the Marčenko-Pastur distribution, along with Monte-Carlo simulations in Figure 1; for other regimes, see [AS17, Theorem 6.27]. The mass term for c < 1 is easily explained by the fact that in this case M < N and the rank of the Proof. We shall use the method of moments, i.e. compute Gaussian expectations of the form

N × N matrix W is M < N thus W shall have N -M ∼ N →∞ N (1 -c) null eigenvalues.
E Tr(W p ) = (N,M ) Tr (XX * ) p dµ(X), (6)
with dµ(X) = 1 (2iπ) N M e -Tr(XX * ) dX * dX, and dX * dX = 1≤i≤N 1≤j≤M dx i,j dx i,j .

The proof consists of three steps. First, we show the exact formula

EN -1 Tr (N -1 W N ) p = α∈Sp N #(γα) M #α , (7)
where γ ∈ S p is the full cycle permutation γ = (1, 2, 3, . . . p) and #α denotes the number of cycles of the permutation α. Note that we dropped the dependence on N of the parameter M , in order to keep the notation light; the reader should keep in mind that M = M N is a function of N which grows as M ∼ cN . The second step, Lemma 2.5, will consist in analyzing the dominating terms in (7). We will show that the surviving permutations are in bijection with non-crossing partitions, recovering the moments of the Marčenko-Pastur distribution lim

N →∞ EN -1 Tr (N -1 W N ) p = α∈NC(p) c #α .
Using Voiculescu's R-transform, we compute in a third step the Cauchy transform of the probability measure having the moments above, and then, by Stieltjes inversion, we recover the exact expression of the density (5).

First step. To show (7), we need to perform the integration on the left-hand-side with the help of the Wick formula from Proposition 2.1. We are going to use a graphical reading of the Wick formula introduced in [CN11]. In this framework, matrices (and more generally, tensors) are represented by boxes having decorations corresponding to the vector spaces the matrix is acting on. The decorations have two attributes: shape, distinguishing vector spaces of various dimensions, and shading, distinguishing primal (filled symbols) from dual (empty symbols) spaces. We depict in Figure 2 (from left to right) the diagram of a Wishart matrix W = XX * , then, in the center, the same diagram, with X * replaced by the transpose of X, and then the diagram for the second moment Tr(W 2 ). In order to establish the formula (7), we need to apply the Wick formula to the quantity Tr(W p ), for an arbitrary p ≥ 1. In [CN11, Section 3.3] (see also [START_REF] Collins | Random matrix techniques in quantum information theory[END_REF]Section III

.C]), it has been shown that X X * X X X X X X 1 2
Figure 2. Diagrams for Wishart matrices. On the left, the diagram for W = XX * , with X ∈ M N ×M . The Hilbert space C N is depicted by round decorations, while C M is depicted by square decorations. The tensor contraction between the two square decorations corresponds to the matrix product X • X * . In the center panel, we have the same diagram, after replacing X * by ( X) ; notice that taking the transposition amounts to inverting the shading of the decorations of the X box. In the last panel, we depict the scalar Tr(W 2 ) = Tr(XX * XX * ). computing a Gaussian expectation can be done in a graphical way, as follows. Given a diagram containing boxes X and X corresponding to random matrices (or tensors) with i.i.d. standard complex Gaussian entries, if the number of X-boxes is different than the number of X-boxes, the expectation (over the randomness in X) is zero. If one has, say, p X-boxes and p X-boxes, the expectation of the diagram is a sum indexed by permutations α ∈ S p , where the terms are obtained by deleting the X and the X boxes, and connecting the corresponding attached decorations with the permutation α: the decorations of the i-th X-box are to be connected to the corresponding decorations of the α(i)-th X box. We would like to warn the reader at this point that the above convention is opposite from the one used in [CN11, Section 3.3], where α was connecting X-boxes to X-boxes. It turns out that the current convention makes the connection between the Wick graphical calculus and the theory of combinatorial maps more transparent, justifying our choice.

As an example, in our moment problem, the two diagrams appearing when computing E Tr(XX * XX * ) are depicted in Figure 3. Notice that the new diagrams are entirely made out of loops, so their (scalar) values are given by N a M b , where a (resp. b) is the number of loops corresponding to C N , that is to round decorations (resp. C M , i.e. square decorations). Figure 3. The two diagrams corresponding to the graphical application of the Wick formula for E Tr(W 2 ), see Figure 2, rightmost panels. On the left, the diagram corresponding to the identity permutation (the blue wires connect the decorations of each X-box to the corresponding ones of the X-box directly following it). On the right is the diagram corresponding to the permutation α = (12). The two diagrams contain only loops, so their values are N M 2 ,respectively N 2 M .

Moving to the general case of an arbitrary p, we write E Tr(W p ) = α∈Sp D α , where D α is the diagram obtained by deleting the p X-and X-boxes and by connecting the corresponding decorations according to the permutation α. It is clear that D α consists only of loops corresponding to the Hilbert spaces C N and C M ; it follows that in order to evaluate such a diagram, one has to count the number of loops of each type. Let us start by counting the loops attached to square decorations, each giving a contribution of M . Note that in the original diagram (before taking the expectation), the square decoration of the i-th X-box is connected to the square decoration of the X-box belonging to the same, i-th, group. It is then easy to see that each distinct cycle of the permutation α gives rise to a loop, then the number of M -loops is #α, thus giving a total contribution of M #α . The same reasoning can be applied when counting the contribution of loops attached to round decorations (which correspond to the Hilbert space C N ), with one difference: in the initial wiring of the diagram (before taking the expectation), the i-th X-box is connected to the (i -1)-th X-box (where the subtraction operation is understood cyclically, modulo p). In other words, the initial wiring is given by the full-cycle permutation γ, with γ(i) = i + 1; note that we are numbering the boxes 1, 2, . . . , p from right to left. A similar combinatorial argument shows that the number of loops is, in this case, #(γα), for a final contribution of N #(γα) . We conclude that D α = N #(γα) M #α , proving (7).

Second step. We now move to the second step of the proof, which is computing the limit N → ∞ of the moment formula (7).

Lemma 2.5. The asymptotic moments of the (normalized) random matrices W N are given by a sum over non-crossing partitions

lim N →∞ EN -1 Tr (N -1 W N ) p = α∈NC(p) c # α. ( 8
)
Proof of Lemma 2.5. Since we need to find the dominating terms in the sum, we have to maximize the function S p α → #α + #(γα). The following lemma contains the key combinatorial insight which allows us to perform this task. The result below is contained in [START_REF] Biane | Some properties of crossings and partitions[END_REF] (see also [NS06, Lecture 23] for a textbook presentation).

Lemma 2.6. For a permutation α ∈ S p , let |α| denote the minimum number of transpositions that multiply to α; |α| is called the length of the permutation α and satisfies the relations

|α| + #α = p |α| = |α -1 | |αβ| = |βα|
for all permutations α, β ∈ S p . The mapping d(α, β) := |α -1 β| defines a distance on S p . For two fixed permutations α, β ∈ S p , the permutations χ ∈ S p saturating the triangle inequality

d(α, χ) + d(χ, β) ≥ d(α, β)
are called geodesic; we write α -χ -β. The set of geodesic permutations between the identity permutation id and the full cycle permutation γ ∈ S p is in bijection with the set of non-crossing partitions NC(p): a partition α ∈ NC(p) encodes the cycle structure of χ, and the elements inside a given cycle of χ have the same cyclic ordering as in γ.

Using the lemma above and the asymptotic relation M N ∼ cN , the exponent of N in the general term of (7) can be bounded as follows

#α + #(γα) = p -|α| + p -|γα| = 2p -(|α| + |α -1 γ -1 |) (9) ≤ 2p -|γ -1 | = p + 1,
where we have used the triangle inequality. The permutations α saturating this inequality are precisely the geodesic ones, i.e. the ones satisfying id -α -γ. These are in bijection with noncrossing partitions α ∈ NC(p), and one has #α = #α, where the # notation denotes at the same time the number of cycles of a permutation α and the number of blocks of the corresponding non-crossing partition α. We have shown

E Tr(W p ) ∼ N p+1 α∈NC(p) c # α.
Third step. The general term in the sum in the right hand side of (8) is a multiplicative function over the blocks of non-crossing partition α: the contribution of each cycle is c, independently of the length of the block. We have thus identified the free cumulants of the limiting distribution of the random matrices W N : κ n = c, for all n ≥ 1 (we refer the reader to [NS06, Lecture 11] for the definition and the basic properties of free cumulants). In order to obtain the density of the probability distribution having the moments above, we use Voiculescu's R-transform machinery. We have

R(z) = ∞ n=0 κ n+1 z n = ∞ n=0 cz n = c 1 -z .
The Cauchy transform and the R-transform are related by the implicit equation G(R(z)+1/z) = z, which we solve for G, obtaining

G(w) = 1 -c + w -(w -a)(w -b) 2w .
Above, we have chosen the solution of the second degree equation in G such that G(w) ∼ w -1 as w → ∞, and we have set a = (1 -√ c) 2 and b = (1 + √ c) 2 . Note that the function G can have poles only at w = 0, with residue max(1 -c, 0), explaining the atom at 0, when 0 < c < 1. The expression for the density is obtained using the Stieltjes inversion formula

dMP c dx = - 1 π lim ε→0 G(x + iε).
For random density matrices, one has to simply take into account the trace normalization: if ρ N is a sequence of random density matrices from the induced ensemble with parameters (N, M N ), we can write ρ N = W N / Tr W N , for W N a sequence of random Wishart matrices of parameters (N, M N ). We have then the following corollary.

Corollary 2.7. Let ρ N ∈ M 1,+ N be a sequence of random density matrices from the induced ensemble of parameters (N, M N ), where M N is an integer sequence with the property that M N ∼ cN as N → ∞, with c ∈ (0, ∞) a constant. The sequence cN ρ N converges, in moments, towards the Marčenko-Pastur distribution MP c .

Proof. Write

cN ρ N = cN W N Tr W N = W N N • cN 2 Tr W N .
From equation (4), it follows that cN 2 /(Tr W N ) converges, almost surely, to 1. Together with Proposition 2.4, this proves the claim.

Note that one can prove much stronger statements of convergence than the ones we cited; importantly, one can show that the largest eigenvalue of (properly normalized) Wishart and random density matrices converges, almost surely, towards the right edge of the support of the limiting Marčenko-Pastur distribution

[BSY88]. Remark 2.8. Note that if W ∈ M N 1 N 2 (C) is a Wishart tensor of parameters (N 1 N 2 , M ), then W = Tr 2 (W ) is a Wishart matrix of parameters (N 1 , N 2 M ). Indeed, if W = XX * , with X ∈ M N 1 N 2 ×M (C) a Gaussian matrix, then W = Y Y * , where Y ∈ M N 1 ×N 2 M (C)
is the matrix obtained by "reshaping" X into a matrix of appropriate dimensions. This equivalence comes from the fact that both the partial trace and the matrix multiplication correspond to tensor contractions. In the random density matrix picture, we have, for a random vector

Ψ ∈ C N 1 ⊗ C N 2 ⊗ C M , ρ 1 = [id N 1 ⊗ Tr N 2 ](ρ 23 ) = [id N 1 ⊗ Tr N 2 ⊗ Tr M ](ΨΨ * ).

A combinatorial map version of the proof.

In this subsection we re-prove the formulas (8) and (9) using combinatorial map methods instead of the results on distances between permutations (Lemma 2.6). Let us first provide a few definitions. Definition 2.9. A connected labeled bicolored combinatorial map, or simply a bicolored map, is a triplet M = (E, σ • , σ • ) where

• E is a set of edges labeled from 1 to p • σ • and σ • are permutations on E,

• the group σ • , σ • generated by σ • and σ • acts transitively on E.

We call white vertices (resp. black vertices, resp. f aces) the elements of the unique decomposition of σ

• (resp. σ • , resp. σ • σ • ) in disjoint cycles. A bicolored map M = (E, σ • , σ • ) satisfies the following classical result: #σ • + #σ • + #(σ • σ • ) -p = 2 -2g(M), ( 10
)
where g(M) is a non-negative integer. A map can also be seen as a graph drawn on a twodimensional surface, up to continuous deformations of the edges: each disjoint cycle of σ • (resp. σ • ) is a white (resp. black) point (vertex) of that surface, each element of E is an arc between a black and a white vertex, and the ordering of the cycles of σ • , and σ • define an ordering of the half-edges incident to the vertices. The graph is said to be cellularly embedded (or simply embedded, in the context of this paper) on that surface if the connected components of the complement of the graph on that surface are homeomorphic to discs, in which case these components correspond to the faces, the disjoint cycles of σ • σ • . The genus (number of holes) of a surface on which the underlying graph of a map can be embedded is precisely the integer g(M) in (10). A genus 0 map can therefore be drawn on the plane without crossing. Finally, the transitivity condition forces the combinatorial map to be connected. With the notations of the proof of Prop. 2.4, we define a bicolored map on E = {1, . . . , p} by setting σ • := γ, and σ • := α. As illustrated in Fig. 4, in this representation, a single black vertex of valency p is associated to Tr(W p ), each matrix W corresponding to an edge incident to it. The trace induces a cyclic counter-clockwise ordering of these half-edges around the vertex 1 (the matrices are labeled growingly from right to left, and the corresponding half-edges are labeled growingly when going counter-clockwise around the black vertex). The free index corresponding to the matrix X (resp X * ) is associated to the left (resp. right) side of the half-edge XX * . Each Wick pairing induces a permutation σ • which defines the white vertices. The maps corresponding to the examples of Figure 3 are shown in Figure 5.

Tr(W 4 ) = N a,b,c,d=1 W a,b W b,c W c,d W d,a 1 2 3 4 ↔ W c,d W a,b W b,c W d,a
Notation 2.10. Throughout the text, we will denote M p the set of labeled connected bicolored maps with a single black vertex and p edges, and M g p its subset of maps having genus g.

1 As will appear in the following, the moments are expressed in terms of labeled maps. This is because the trace corresponding to the only black vertex does not come with a factor 1/p!. The two elements of M 2 = M 0 2 are shown in Fig. 5, and an example of map in M 0 7 is shown on the left of Fig. 6. Denoting F (M) the number of faces of a map M, and V (M) its total number of vertices, we can therefore restate the key relation (7) as follows:

E Tr W p = M∈Mp N F (M) M V (M)-1 . ( 11
)
The genus (10) of a connected map in M p writes

2 -2g(M) = V (M) + F (M) -p, (12)
so that, if we assume M ∼ cN ,

E Tr W p = (1 + o(1)) M∈Mp N p+1-2g(M) c V (M)-1 . ( 13
)
Since g ≥ 0, the bound (9) on the number of cycles #α+#(γα) is now obtained using the positivity of the genus instead of the distances between permutations, and the leading terms in N correspond to planar maps. Therefore, with the Notation 2.10, lim

N →+∞ 1 N p+1 E Tr W p = M∈M 0 p c V (M)-1 . (14)
It is easily seen by studying the permutations σ • = γ and σ • = α, that elements of M 0 p with V (M) vertices correspond to non-crossing partitions with V (M) -1 disjoint cycles. Another simple way of seeing this is to notice that maps in M 0 p with V (M) vertices and F (M) faces are in bijection with bicolored labeled plane trees with p edges, V (M) -1 white vertices, and F (M) black vertices (plane trees are all bicolored), as explained below. Such trees are themselves in bijection with non-crossing partitions with V (M) -1 disjoint cycles throughout a dual mapping.

Starting from a map in M 0 p , a vertex is added in each face and an edge is added between each corner of each white vertex and the newly added vertices (the ordering of edges around vertices is given by the clockwise ordering of appearance of corners around faces)2 . This procedure is illustrated in Fig. 6. It is a particular case of Tutte's bijection for bicolored maps [Tut75] and is known to be bijective. The map obtained from vertices, and its number of edges is the number of corners of M around white vertices, which is p, and is therefore a tree.

M ∈ M 0 p is connected, has V (M )-1+F (M ) = p+1
In the case where c = 1, lim N →+∞ 1 N p+1 E Tr W p is just the cardinal of M 0 p , namely the number of plane trees with p edges, which is well-known [HPT64] to be the Catalan number

C p = 1 p+1 2p p , lim N →+∞ 1 N p+1 E Tr W p c=1 = C p . (15)
At finite N and for c = 1, the contribution of order N p+1-2g is given by the number of maps in M p whose genus is g. The very same transformation which mapped maps in M 0 p to plane trees can be applied. It maps elements of M g p to bicolored maps of genus g with a single face. Maps with a single face are called unicellular maps [Cha10]. Conversely, the inverse transformation maps bicolored unicellular maps of genus g to genus g bicolored maps with a single black vertex. Unicellular maps of genus g are the object of many studies, and they can be counted exactly using various recursive formulas (see for instance the Harer-Zagier formula [START_REF] Harer | The euler characteristic of the moduli space of curves[END_REF], Lehman-Walsh formula [WL72], Goupil-Schaeffer formula [GS98], and Chapuy's formula [START_REF] Chapuy | A new combinatorial identity for unicellular maps, via a direct bijective approach[END_REF]CFF13,[START_REF] Carrell | Simple recurrence formulas to count maps on orientable surfaces[END_REF]). Using such results, (13) can be expressed exactly for any p.

2.4.

Free cumulants and maps. We elaborate here on the relation between moments and free cumulants in the case of a Wishart distribution. As was pointed out in the proof of Proposition 2.4 (third step), the free cumulants of the limiting distribution are κ n = c for all n ≥ 1. We consider again the expectation value of the traces of powers of the matrix W

E Tr(W p ) = N i 1 ,...,ip=1 E(W i 1 ip . . . W i 3 i 2 W i 2 i 1 ), (16)
= N i 1 ,...,ip=1 J {1,...,p} Jr∈J K(W Jr ), (17)
where the second sum is taken over all partitions J of {1, . . . , p}, and we define

K(W Jr ) = K({W i k+1 i k } k∈Jr k∈Zp ) (18)
as the classical cumulants of the family of random variables

{W i k+1 i k } k∈Jr k∈Zp .
The claim (17) can be obtained from the definition of the generating function of classical cumulants. Depending on conventions, this cumulant generating function can be defined either as the logarithm of the characteristic function or of the moment generating function of the corresponding probability density. For reference on classical cumulants, see for instance [NS06, Lecture 11, Appendix], [Nov14] or [Luk70]. As moments, these classical cumulants can be represented graphically. In the mathematical-physics literature, the graphical objects representing them are often described as connected combinatorial maps (or ribbon graphs) with half-edges or connected opened combinatorial maps, see Remark 2.13 for more details. In our case, as is shown in the next paragraphs, they will be represented as usual combinatorial maps with one white vertex 3 and one black vertex. The fact that combinatorial maps with one white vertex and one black vertex represent classical cumulants is the translation of the decomposition of classical cumulants in terms of classical moments for some random variables X 1 , . . . , X s , which is the cumulant-moment formula in classical 3 The fact that there is only one white vertex is the direct translation of the "connected" conditions introduced in the mathematical-physics literature.

probability [NS06, Definition 11.28], making use of the Möbius function on the lattice of partitions:

K(X |Jr| , . . . , X 1 ) = I {1,...,|Jr|} (-1) |I|-1 (|I| -1)! Im∈I E   k∈Im X k   , (19)
for the random variables X k = W i k+1 i k . The moments E k∈Im X k are then computed using Wick's theorem (Prop. 2.1). As detailed previously, these Wick pairings can be represented graphically, either using X-and X * -boxes, or using the combinatorial map representation. In this last representation, identifying W i k+1 i k with the integer k ∈ I m , the Wick pairings induce a permutation α m on the set I m (see the proof of Prop. 2.4). Each cycle of α m contributes with N c . As before, we can define a bicolored map in M |Im| by setting σ • = α m , and σ • = γ m , where γ m is the cycle (x 1 , . . . , x |Im| ), in which x 1 < x 2 < . . . < x |Im| are all the elements of I m .

The difference with what has been done before, where we were considering the expectation of a trace, is that here, the indices {i k , i k+1 } are not summed. This means that the faces of the map (I m , α m , γ m ) (the disjoint cycles of γ m α m ) do not factorize in the expression of E k∈Im X k , and instead we have a Kronecker delta for each corner of each white vertex:

E   k∈Im W i k+1 i k   = αm∈S |Im| N c (#αm) k∈Im δ i k ,i αm(k)+1 . (20)
Furthermore, because of the factors (-1) |I|-1 (|I| -1)! in the expression of the cumulants (19), and because the products of Kronecker deltas factorize, the terms corresponding to maps which have more than one white vertex will be canceled by products of terms corresponding to maps with fewer white vertices. The only non-vanishing terms are those for which I is reduced to {1, . . . , |J r |}, and the cumulants can therefore be expressed in terms of cycles of J r of length |J r |, or equivalently using a graphical expansion over bicolored maps with two vertices,

K(W Jr ) = N c M∈M |Jr | V (M)=2 k∈Jr δ i k ,i α M (k)+1 , (21)
where we have denoted α M the cycle σ • for the map M.

We do not prove these statements here, however we give what we hope is an illuminating example in the present case. Consider the cumulant of three matrix elements

W i 6 i 5 , W i 4 i 3 , W i 2 i 1 , K(W i 6 i 5 , W i 4 i 3 , W i 2 i 1 ) = E(W i 6 i 5 W i 4 i 3 W i 2 i 1 ) -E(W i 6 i 5 W i 4 i 3 )E(W i 2 i 1 ) -E(W i 6 i 5 )E(W i 4 i 3 W i 2 i 1 ) -E(W i 6 i 5 W i 2 i 1 )E(W i 4 i 3 ) + 2E(W i 6 i 5 )E(W i 4 i 3 )E(W i 2 i 1 ). (22)
The graphical representation of the Wick pairings of equation ( 22) is shown in Fig. 7. The dashed edges indicate which X, X * are paired together. In order to simplify (and compactify) the drawings, the vertices of empty types represent the X-boxes, while the filled vertices represent the X * -boxes. In this representation, it is the disconnected parts of the graphical expansion of E(W i 6 i 5 W i 4 i 3 W i 2 i 1 ) which are suppressed by the remaining terms of the right hand side of ( 22), that are all disconnected in nature. The weights associated to the only remaining terms are N c δ i 1 i 6 δ i 2 i 3 δ i 4 i 5 and N c δ i 1 i 4 δ i 2 i 5 δ i 3 i 6 . In Fig. 8, we represent the same expansion using the representation in terms of combinatorial maps.

+ + + + + + - - - +2 K(W i6i5 , W i4i3 , W i2i1 ) = + i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 Figure 7. Graphical expansion of the cumulant K(W i 6 i 5 , W i 4 i 3 , W i 2 i 1
) from its moments decomposition and the Wick theorem. Identical terms are shown in the same color.

i 4 i 3 i 5 i 6 i 1 i 2 + + + + + + + - - - +2 i 4 i 3 i 2 i 1 i 6 i 5 K(W i6i5 , W i4i3 , W i2i1 ) = i 1 i 2 i 6 i 5 i 4 i 3 i 3 i 4 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 + i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 2 i 1 i 6 i 5 i 4 i 3 i 5 i 6 i 1 i 2 i 3 i 4 i 2 i 1 i 6 i 5 i 1 i 2 i 6 i 5 i 4 i 3 Figure 8
. This is the translation of the previous Figure 7 in the language of combinatorial maps. Every corner on a white vertex carries a δ i k ,i k , where i k and i k are the incident indices on the drawing. The only graphs which do not cancel are the ones corresponding to N c δ

i 1 i 6 δ i 2 i 3 δ i 4 i 5 and N c δ i 1 i 4 δ i 2 i 5 δ i 3 i 6 (shown in black).
Using the expression (21) in the classical moment-cumulant formula ( 16) above, we have

E Tr(W p ) = N i 1 ,...,ip=1 J {1,...,p} Jr∈J N c M∈M |Jr | V (M)=2 k∈Jr δ i k ,i α M (k)+1 (23) = J {1,...,p} N c |J | N i 1 ,...,ip=1 Jr∈J M∈M |Jr | V (M)=2 k∈Jr δ i k ,i α M (k)+1 . ( 24
)
We can rewrite the terms on the second line in terms of bicolored maps in M p . The partition J divides the edges incident to its black vertex in groups J r which are all linked to the same white vertex, one for each group. For each one of these groups, a sum is taken over all possible maps with two vertices, of the product of Kronecker deltas over the corners around the white vertex.

Any choice of |J | two-vertices maps provides a different map in M p . Denoting M p [J ] the subset of M p of maps for which for every J r ∈ J , (α M ) |Jr is a unique cycle of length |J r |, we can therefore rewrite this sum as

E Tr(W p ) = J {1,...,p} N c |J | N i 1 ,...,ip=1 M∈Mp[J ] k∈{1,...,p} δ i k ,i α M (k)+1 . (25)
Now, we have precisely one index per corner around the black vertex, and we still have one Kronecker delta per corner around a white vertex, so the sum over i 1 , . . . , i p yields a factor N F (M) , so that

E Tr(W p ) = J {1,...,p} N c |J | M∈Mp[J ] N F (M) = M∈Mp N c V (M)-1 N F (M) , (26)
and we indeed recover (11), and therefore the results of Lemma 2.5 which state that the leading terms correspond to the planar maps. However, we can now pull back the planarity result to the level of the cumulants. From the calculation above, we see that the cumulants are the two-vertices maps which are then re-arranged (their black vertices are crushed into a single one, and the edges are naturally ordered around this single vertex as they are labeled from 1 to p). In order for the resulting map to contribute to the large N limit, it has to be planar, thus the collection of two-vertices maps we start from all have to be planar. Indeed it is a classical result that the genus of a map is greater or equal to that of a submap (see e.g. [MT01, Prop. 4.1.5, Sec. 4.1, p. 102]). We reprove it in our case in Prop. 2.11. If all the two-vertices maps are planar, in addition, the partition J has to be non-crossing for the map to be planar.

Proposition 2.11. Consider a map M ∈ M p , and the submap M obtained by keeping the black vertex, a single white vertex, and all the edges between them. If M is non-planar, then M is also non-planar.

Proof. Suppose that g(M ) > 0, and consider the set of edges of the combinatorial map that do not belong to M . Among these edges, we first remove all the ones adjacent to univalent white vertices, and we remove the isolated white vertices created in this process. The obtained map has the same genus as the initial map. For each of the remaining edge we proceed as follows. If the edge is adjacent to only one face and to a white vertex which is not univalent, then the map has a non zero genus and the proposition is true. If the edge is adjacent to two different faces, withdrawing it does not change the genus. We do so, and then repeat this process. At each step, we may create univalent white vertices, which we remove, together with the isolated white vertices thus created. Either we find during the process an edge which is not a bridge and is adjacent to only one face (and the proposition is proved), or we have erased all the edges that did not belong to M , without changing the genus of the map. We are thus left with a map with non vanishing genus g(M) = g(M ) > 0.

We can now claim that the free cumulant of order p is represented by the unique planar map with one white vertex, one black vertex and p edges. The interest of this remark lies in the fact that it gives a heuristic for reading freeness directly out of the combinatorial maps expression of moments (intuition which will be used in Sections 3.2 and 4.2). Indeed looking at moments of alternating products of different random matrices, one expects to express these moments using colored combinatorial maps 4 (where the colors label different types of edges indexing Wick pairing between different matrices). If one believes 5 that white vertices with adjacent edges of different colors correspond to mixed cumulants, it is easy to guess if in the large N limit the different matrices converge to different free elements, just by confronting the weights of the white vertices with only one adjacent color to those of white vertices with strictly more than one adjacent color.

Remark 2.12. We notice here that the following expectation values of products of matrix elements such that i k = i q , ∀k = q (indices are fixed and not summed)

E(W i 1 ip . . . W i 3 i 2 W i 2 i 1 ) = M = N c , (27)
thus are equal up to a rescaling in N to the free cumulants κ n = c. A similar fact has been used in the study of Hermitian (formal) matrix integrals and simple combinatorial maps in [START_REF] Borot | Simple maps, Hurwitz numbers, and Topological Recursion[END_REF]. In this reference a combinatorial maps interpretation of the relation between moments and classical cumulants is given by bijective means valid at all orders in the large N asymptotic expansion. In particular the bijective method is valid at the first order of the expansion and thus can be interpreted in the context of free probability and the corresponding free cumulants.

Remark 2.13. The cumulants can also be expressed in terms of the white vertices with incidentpending half-edges, whose both sides are labeled by pairs (i k , i k+1 ), while the half-edges themselves are labeled by the first index of the pair k ∈ I r . In this representation, the moments are obtained by gluing the cumulants around the black vertex while respecting the natural order of natural integers, and the free cumulants are the connected terms, for which the indices of the half-edges satisfy a non-crossing condition.

The four-partite case: joint distribution of two marginals

In this section, we describe the joint distribution of two marginals, ρ AB and ρ AC , of a random pure quantum state ρ ABCD = |ψ ψ| ABCD . We prove all results both using metric properties of the permutation groups and the combinatorial maps approach with the underlying discrete geometry, in parallel. As before, we state our results for Wishart tensors; the corresponding results for random density matrices can be obtained by renormalizing the Wishart matrices (see Remark 3.9, as well as the very last paragraph of this section).

Let X = X ABCD ∈ H A ⊗ H B ⊗ H C ⊗ H D be a random Gaussian tensor, i.e. a random tensor with i.i.d. standard complex Gaussian entries. We are going to study here the joint distribution of the two marginals

W AB := [id ⊗ id ⊗ Tr ⊗ Tr](XX * )
and

W AC := [id ⊗ Tr ⊗ id ⊗ Tr](XX * ), (28)
in some asymptotic regime, where the dimension of some of the Hilbert spaces grow to infinity. In order to study mixed moments of the two random matrices W AB and W AC , we are going to assume from now on that dim

H B = dim H C .
We consider the following two asymptotic regimes:

• the balanced regime: dim H B = dim H C = N → ∞, dim H D / dim H A → c, for some constant c ∈ (0, ∞); • the unbalanced regime: dim H B = dim H C = m, dim H A ∼ N , dim H D ∼ cN , N → ∞, for
some constants m ≥ 1 and c ∈ (0, ∞). As special cases of the balanced asymptotic regime, we shall emphasize the following two sub-cases (see Remark 3.7):

• dim

H B = dim H C = N , dim H A ∼ c 1 N , dim H D ∼ c 4 N , N → ∞, for some constants c 1,4 ∈ (0, ∞); • dim H B = dim H C = N , dim H A = k, dim H D = l, N → ∞, for some integer constants k, l ≥ 1.
The balanced regime will be studied in Section 3.2, while the unbalanced regime will be studied in Section 3.4. We first compute the exact combinatorial expression for the moments.

Exact expression for the moments.

In this section, we prove a non-asymptotical result, describing the mixed moments of the marginals W AB , W AC as a combinatorial sum. This result is the starting point for the asymptotic computations in Sections 3.2 and 3.4. In order to describe the mixed moments of the two marginals, we write

W f := ←- 1≤a≤p W f (a) ,
for a word f ∈ {AB, AC} p of length p in the letters (AB), (AC). Here again, the matrices are labeled from right to left. We also denote N X := dim H X , for X = A, B, C, D; note that we have

N B = N C =: N B,C .
Theorem 3.1. The mixed moments of the two marginals W AB , W AC defined in (28) are given by the following exact formula:

E Tr W f = α∈Sp N #(γα) A N #α D N L(f,α) B,C , (29)
where γ = (123 . . . p) is the full cycle permutation and L(f, α) is a combinatorial function described in equations (30)-(33).

Proof. The proof is a standard application of the graphical Wick formalism from [CN11], similar to the proof of the moment formula in Proposition 2.4. First, consider the diagram associated to the trace of the operator W f (see Figure 9 for a simple example): each of the p X (or X) boxes has 4 decorations, corresponding to the 4 Hilbert spaces H A,B,C,D . The sum over permutations α ∈ S p comes from the application of the Wick formula, and the different factors in the general term count the different types of loops, as follows:

• There are #(γα) loops associated to the Hilbert space C N A , since the initial wiring of the boxes is given by the cyclic permutation γ, corresponding to matrix multiplication. The total number of loops is the number of cycles of the permutation obtained as the product of the permutation describing the initial wiring (here, γ, connecting X boxes to X boxes) and the permutation coming from the graphical Wick formula (α, connecting X boxes to X boxes). • There are #α loops associated to the Hilbert space C N D , since the initial wiring of the boxes is given by the identity, corresponding to the partial trace over the fourth tensor factor (which appears in all the W f (i) ). • We denote by L(f, α) the number of loops corresponding to the second (B) and the third (C) tensor factors. This quantity is less obvious to evaluate, since it depends in a non-trivial way on both the specific word f and on the permutation α. These loops contribute each a factor of N B,C .

We discuss next the exponent L(f, α), counting the number of loops associated to systems B and C. As in the other cases, we have

L(f, α) = #(γ f α), ( 30
)
where γf , α ∈ S 2p encode, respectively, the initial and the Wick wirings of the B, C decorations. In order to define properly these permutations, let us relabel the 2p decorations of the X boxes by (1, B), (2, B), . . . , (p, B), (1, C), (2, C), . . . , (p, C) to keep track of the system they refer to (as usual, the numbering is increasing from right to left). Then, using cyclical notation (i.e. p + 1 ≡ 1), we have

γf (a, B) =      (a, B) if f (a) = AC (a + 1, B) if f (a) = AB and f (a + 1) = AB (a + 1, C) if f (a) = AB and f (a + 1) = AC (31) γf (a, C) =      (a, C) if f (a) = AB (a + 1, B) if f (a) = AC and f (a + 1) = AB (a + 1, C) if f (a) = AC and f (a + 1) = AC (32) α(i, Z) = (α(i), Z), for Z = B, C. (33)
As a direct application of (29), we have, for the example in Fig. 9,

E Tr(W AB W AC ) = N A N 3 B,C N 2 D + N 2 A N B,C N D ,
where the two terms in the right hand side correspond respectively to the identity permutation and to the transposition (12).

We shall now present a different point of view on Theorem 3.1, using combinatorial maps. First let us introduces some important (and somehow non-standard) notation for integer intervals: for integers a, b, we denote

a, b := {a, a + 1, . . . , b -1, b}. (34)
Orbits. Before going to the combinatorial maps version of this result, we provide a more formal definition of the objects we count, in terms of the permutation. If I ⊂ {A, B, C, D}, we denote

[p] I := i∈I (a, i) a∈ 1,p , (35)
so that for instance, [p] B,C = (a, B) a∈ 1,p (a, C) a∈ 1,p , where a runs through the set of edges6 . In some sense, [p] I is the set of all potential colored edges with color set I that can be obtained from the p edges of a map.

Definition 3.2 (Orbits). We call orbits for f and α, the cycles respectively induced by γα, γf α, and α on the sets [p] A , [p] B,C , and [p] D . 7

Of course, these are in bijection with the loops in the box representation. Note that the orbits involving the colors A and D are already understood from Section 2 (corresponding to the bipartite quantum system case), and only the orbits involving the colors B and C present a new behavior.

Combinatorial maps. As just mentioned, as far as the colors A and D are concerned, the behavior is unchanged from Section 2. So, naturally, if we define the map M = (γ, α) ∈ M p , the orbits corresponding to the color A are simply the faces of M, while there is one orbit of color D per white vertex of the map. Colors A and D therefore contribute with a factor N

F (M) A N V (M)-1 D .
We now focus on the orbits involving colors B and C. We can provide each edge a ∈ 1, p of the map with a color j(a) ∈ {B, C}: the non traced color in {B, C} for the matrix number a of the word f (i.e. the color different from A in f (a)). For instance, the moment E(Tr W AB W AC ) of the trace of Fig. 9 involves the two maps in M 2 (Fig. 5), but in addition, one edge carries the color B, the other the color C. We consider the set of the colored edges,

[p] B,C f = a, j(a) a∈ 1,p . (36)
It is the subset of [p] B,C corresponding to the colors which are not traced. Some orbits contain at least one element of [p] B,C f , or equivalently, follow at least one edge of the map M. In the map language, when such an orbit reaches a corner on the black vertex after having followed a colored edge a, j(a) , a ∈ 1, p and j(a) ∈ {B, C}, it will go to the colored edge a + 1, j(a + 1) . Apart from the information on the color, the orbit locally behaves, around the black vertex, as the face of the map M, which goes from a to a + 1.

However, the orbits and the faces of the map do not behave alike around white vertices. Indeed, arriving from a colored edge (a, j(a)) at a white vertex v • , the orbit goes to the first edge around v • counterclockwise that also carries the color j(a), namely α q (a), where q = min s ∈ N * | j(a) = j(α s (a)) .

(37)

If the integer q is greater than one, the behavior differs from that of the face, which goes to α(a).

A simple way of solving this issue is to locally split the white vertex into two vertices, one for each color B and C, as illustrated in Fig. 10. Note that this operation has no effect if all the edges Performing this operation on all the white vertices, we obtain a map M f ∈ M p , with the same black vertex γ as M, but with two permutations α B and α C , which are simply the restrictions of α to the two sets of edges that respectively carry the colors B and C. This defines a tricolored map M f = (σ • , α B , α C ), where σ • = γ, and such that edges only link white vertices and edges of 7 Where we canonically identified [p] A and [p] D with 1, p . color B or C (i.e. it is also a bicolored map (σ • , α B α C ) in the sense of Def. 2.9). Because if j(a) is e.g. B, α B (a) = α q (a) defined in (37), the faces of M f are precisely the orbits, and therefore the colors B and C contribute with a factor N

F (M f ) B,C
. And lastly, some orbits do not follow any edge of the map. They correspond to the restriction of γf α to the set

[p] B,C \ [p] B,C f , (38)
and in the map language, to white vertices of M whose incident edges all share the same color. Such white vertices are those left invariant by the operation of Fig. 10, and therefore their number is just

2V (M) -V (M f ) -1.
(39) Indeed, 2(V (M) -1) -(V (M f ) -1) counts one for each white vertex with only adjacent edges of a single color, while it counts zero for each white vertex with both adjacent colors.

We have therefore the following formulation of Thm. 3.1.

Theorem 3.1 in the maps formulation. For f a word of length p in the alphabet {AB, AC}, the mixed moments of the two marginals W AB , W AC are expressed exactly using a sum over combinatorial maps

E Tr(W f ) = M∈Mp N F (M) A N V (M)-1 D N L(f, M) B,C , (40)
where L(f, M) coincides with L(f, α) and can be expressed as

L(f, M) = F (M f ) + 2V (M) -V (M f ) -1. (41)
A bound on the exponent L. We establish now some key results that will be used in the following sections, regarding the range of the functional L(f, •). Note that the other exponents in (29) or (40) have been dealt with in Lemma 2.6 in terms of permutations, where it has been shown that |α| + |α -1 γ -1 | ≥ |γ -1 | = p -1, with equality if the permutation α is geodesic with respect to γ -1 , i.e. it is associated to a noncrossing partition (the cycle structure of α is non-crossing and inside each cycle the elements have the same cyclic ordering as γ -1 , or in terms of maps, using the genus (12). The next results are a generalization of the above fact, which corresponds to the situation where f is constant. Proposition 3.3. For f a word of length p in the alphabet {AB, AC} and α ∈ S p , we have the bound L(f, M) ≤ p + 1, (42) with equality iff the map M is planar, and all the edges of M incident to a common white vertex have the same color.

In terms of the permutation α, L(f, α) = p + 1 with equality iff α is geodesic and α ≤ ker f (for the usual partial order on partitions), where ker f is the partition having (at most) two blocks, f -1 ({AB}) and f -1 ({AC}).

Proof in the map language. The relation of the Euler characteristic for the map M f writes

2 -2g(M f ) = V (M f ) -p + F (M f ). ( 43
)
we can therefore rewrite (41) as

L(f, M) = p + 1 -2g(M f ) -2∆ f (M), (44)
where we have denoted ∆ f (M) the number of white vertices reached by both colors B and C,

∆ f (M) = V (M f ) -V (M). ( 45
)
In particular, this vanishes iff in M, all the edges incident to a common white vertex have the same color, and is positive otherwise. Therefore, L(f,

M) = p + 1 iff ∆ f (M) = g(M f ) = 0. Moreover, if ∆ f (M) = 0, the maps M f and M coincide, so that L(f, M) = p + 1 iff ∆ f (M) = g(M) = 0.
Note that, using the again the genus formula for M, we obtain the following exact expression for the moments

E Tr(W f ) = M∈Mp N p+2-2g(M)-V (M) A N V (M)-1 D N p+1-2g(M f )-2∆ f (M) B,C . ( 46
)
Proof in terms of permutations. Let us first restate the proposition in the form of a lemma. To recover Prop. 3.3, it suffices to use the fact that

#(γ f α) + |γ f α| = 2p. ( 47
)
Lemma 3.4. For any f ∈ {AB, AC} p and α ∈ S p , we have

|γ f α| ≥ p -1,
with equality iff α is geodesic and α ≤ ker f (for the usual partial order on partitions), where ker f is the partition having (at most) two blocks, f -1 ({AB}) and f -1 ({AC}).

Proof. Let us start by re-writing the permutation γf ∈ S 2p (we are using the labeling of (35), Indeed, one can check formula (48) by comparing it with equations (31) and (32). We now apply [CMSS07, Lemma 5.5] to the permutations γf , α ∈ S 2p :

[p] B,C = {(1, B), . . . , (p, B), (1, C), . . . , (p, C)}) γf = δ C γ B ⊕ id C δ C , (48)
|γ f α| + |γ f | + |α| ≥ 2|γ f ∨ α|, (49)
where the join operation ∨ on the right hand side should be understood as acting on the partitions induced by the cycles of the two permutations (the join of two partitions is their least upper bound). The notation | • | is extended to partitions as |π| := 2p -#π, for a partition π of [2p], where #π denotes the number of blocks of π. Since γf has one cycle of length p and p fixed points, and α = α B ⊕ α C (using the same notation as above), the inequality becomes

|γ f α| ≥ 4p -2#(γ f ∨ α) -(p -1) -(2p -2#α) = p + 1 + 2(#α -#(γ f ∨ α)).
We claim that #(γ f ∨ α) = 1 + pure(α, f ), where pure(α, f ) denotes the number of cycles of α on which the restriction of f is constant, i.e. in the map language, the white vertices whose incident edges all carry the same colors, so that pure(α, f ) is given by (39). Indeed, the blocks of the partition γf ∨ α correspond to the blocks of α B and α C merged by the block of length p of γf , together with each block of α B,C matching only fixed points of γf . It is now clear that the latter are exactly blocks of α B on which f ≡ AC or blocks of α C on which f ≡ AB; we conclude that there are exactly pure(α, f ) of those, proving the claim. Hence,

|γ f α| ≥ p -1 + 2(#α -pure(α, f )) ≥ p -1,
proving the inequality.

Let us now characterize the permutations α which saturate the inequality (assuming f is fixed). First, the last inequality above should be an equality, so α = pure(α, f ); in other words, f should be constant on the blocks of α, which is the condition α ≤ ker f appearing in the statement of the lemma. Moreover, we can easily see that this condition is equivalent to the fact that the permutations α and δ C commute. Indeed, by direct computation, we can see that, for any a ∈ 1, p and any Z = B, C, we have

[δ C αδ C ](a, Z) = (α(a), Z) if f (a) = f (α(a)) (α(a), Z) if f (a) = f (α(a)), (50)
where we denote by Z the complement of Z in {B, C}. Hence, δ C αδ C = α iff f is constant on the cycles of α, which is the claimed statement. Using (48), we have

|γ f α| = |δ C γ B ⊕ id C δ C α| = | γ B ⊕ id C α| = |(γα) B ⊕ α C | = |α| + |αγ| ≥ p -1,
with equality iff α is geodesic, proving the other equality condition in the statement and finishing the proof.

Reformulation of the exponent L. In this paragraph, we give a more intuitive expression for the exponent L. We define alt(f, α) as the total number of changes of colors around the cycles of α

alt(f, α) = {a ∈ 1, p : f (a) = f (α(a))} . (51)
In the map language, alt(f, M) is the total number of corners around white vertices of M whose incident edges a and α(a) have different colors j(a) = j(α(a)).

Proposition 3.5. For any f ∈ {AB, AC} p and α ∈ S p , where α is geodesic (or equivalently, where M = (γ, α) is planar), we have

L(f, α) = p + 1 -alt(f, α). ( 52
)
We also provide two proofs, first in terms of permutations, and then in terms of maps. In the following, for M = (γ, α), we identify L(f, α) and L(f, M), as well as alt(f, α) and alt(f, M). We stress that (52) is true only for a geodesic permutation α, or equivalently a planar map M. In general, for non-planar maps, 2(g

(M f ) + ∆ f (M)) = alt(f, M). For instance, if ∆ f (M) = 0, M f = M, so that 2(g(M f ) + ∆ f (M)) = 2g(M), but alt(f, M) = 0.
Proof in terms of permutations. Let us assume that α is geodesic (w.r.t. γ -1 , i.e. it satisfies |α|

+ |αγ| = |γ -1 | = p -1) and prove |γ f α| = p -1 + alt(f, α), (53)
which proves the proposition using (47). We write, as before,

|γ f α| = |δ C γ B ⊕ id C δ C α| = |δ C αδ C (γ B ⊕ id C )| ≤ |δ C αδ C α-1 | + |α(γ B ⊕ id C )| = |α| + |αγ| + |δ C αδ C α-1 | = p -1 + alt(f, α),
where we have used the triangle inequality and the fact that

δ C αδ C α-1 = a : f (a) =f (α -1 (a)) (a, B)(a, C) ,
which follows easily from (50). To conclude, we need to show that the triangle inequality used above is saturated, which is equivalent to the three permutations

δ C αδ C -α -(γ B ) -1 ⊕ id C
lying on a geodesic. Since δ C αδ C α-1 is a product of disjoint transpositions, this is in turn equivalent to the fact that, for all a such that f (a) = f (α -1 (a)), the elements (a, B) and (a, C) are contained in the same cycle of the permutation ε := δ C αδ C (γ B ⊕ id C ). This permutation acts in the following way:

(a, B) → (α(a + 1), B) if f (α(a + 1)) = f (a + 1) (α(a + 1), C) if f (α(a + 1)) = f (a + 1) (a, C) → (α(a), C) if f (α(a)) = f (a) (α(a), B) if f (α(a)) = f (a).
In other words, on the B-level, Proof in the map language. We assume that the map M ∈ M p is planar, and prove by induction that

alt(f, M) = 2g(M f ) + 2∆ f (M). ( 54
)
If v • is a white vertex and v • is the only black vertex, we denote M v• the submap obtained from M by keeping only v • and v • and the edges linking them. We can apply the operation of Fig. 10 to the white vertex of M v• , thus obtaining a map M v• f . Because the map M is planar, the map M f can be constructed by recursively inserting the maps M v• f in the corners in the appropriate way (inverse operation of that shown in Fig. 17), where v • spans the white vertices of M. Therefore, any operation on the edges of a given M v• f does not affect the faces incident to other

M v •
f , and the genus of M f is the sum of the genera of all the

M v• f g(M f ) = v•∈M white vertex g(M v• f ). ( 55
)
Note that this is not true for a non-planar map M. We can therefore study what happens locally, for a single M v• . The second step is to notice that if two consecutive edges around a v • have the same color, then removing one or the other will not affect the genus g(M v• f ). We can therefore consider that there are no two consecutive edges of the same color. The case with 3 edges of each color is shown in Fig. 12 a) below.

B C B C B C a) B C B C B C b) B C B C C c) B C B C d) Figure 12. a) A six edges M v• . b-d) Recursive operations on M v• f .
Firstly, in the case where M v• is made of only two edges of two different colors, we indeed have 2 = alt(f,

M v• ) = 2g(M v• f ) + 2∆ f (M v• ) = 2 • 0 + 2 • (3 -2) = 2.
If now there are more than two edges, it is easy to see that a single face visits all the edges twice in M v• f (see the dotted face in Fig. 12 b)). Furthermore, these edges are not bridges, and therefore, deleting any edge in M v• f , the genus decreases by one (there is one more face and one less edge), and the number of corners incident to edges of different colors around v • in M v• decreases by two. In the resulting map (Fig. 12 c)), two edges have two incident faces. Deleting one of them, the genus does not vary, nor the number of corners incident to edges of different colors, and we recover a map with less edges and with the same property (Fig. 12 d)). By induction we deduce that the number of corners incident to edges of different colors around v • in M v• is twice the genus of M v• f plus two. Summing over white vertices (55), we obtain the sought relation alt(f, M) = 2g(M f ) + 2∆ f (M).

The balanced asymptotical regime.

We focus in this section on the asymptotic regime where N := dim H B = dim H C → ∞, which we call balanced; for the case when the size of the B, C subsystems stays bounded, see Section 3.4. In this regime, where H B,C grow, it turns out that the asymptotic behavior of H A,D is not so important, as long as the ratio dim H D / dim H A converges to a positive constant as N → ∞. We present next the main result of this section, and discuss several particular asymptotic scenarios later.

Theorem 3.6. Let X N ∈ C N A ⊗ C N ⊗ C N ⊗ C N D
be a sequence of random Gaussian tensors, where N A,D are arbitrary functions of N satisfying N D ∼ cN A as N → ∞, for some constant c ∈ (0, ∞). Then, the normalized marginals

(N -1 A N -1 W (N ) AB , N -1 A N -1 W (N )
AC ) defined in (28) converge in distribution, as N → ∞, to a pair of identically distributed and free elements (x AB , x AC ), where x AB and x AC have a MP c distribution. Equivalently, for any word in the two marginals f ∈ {AB, AC} p , we have

lim N →∞ E(N A N ) -p-1 Tr ←- 1≤i≤p W (N ) f (i) = α∈NC(p), α≤ker f c #α , ( 56
)
where ker f is the partition having two blocks corresponding to the occurrences of AB (resp. AC) in the word f . This rewrites in terms of planar bicolored maps with one white vertex as

lim N →∞ E(N A N ) -p-1 Tr ←- 1≤i≤p W (N ) f (i) = M∈M 0 p , ∆ f (M)=0 c V (M)-1 , (57)
where the sum is restricted to maps whose white vertices only have incident edges of the same color, and we recall that M 0 p is the subset of the elements of M p of vanishing genus. Proof. We first prove formula (56). Starting from the exact moment formula of Theorem 3.1, let us analyze the contribution of N , through its exponent L(f, α). From Proposition 3.3, L(f, α) ≤ p + 1, with equality iff α is geodesic and α ≤ ker f . Hence, we have

E Tr W (N ) f = (1 + o(1))N p+1 α∈NC(p), α≤ker f N #(γα) A N #α D = (1 + o(1))(N A N ) p+1 α∈NC(p), α≤ker f N D N A #α = (1 + o(1))(N A N ) p+1 α∈NC(p), α≤ker f c #α ,
proving the claimed formula. Above, we have used the key fact that, for the surviving α terms (i.e. the permutations which are geodesic w.r.t. γ -1 ), #(γα) = p + 1 -#α.

We now show how the moment formula (56) implies the main claim. Using the moment-cumulant formula [NS06, Proposition 11.4], one can read the asymptotic free cumulants directly off the moment formula:

lim N →∞ E(N A N ) -p-1 Tr ←- 1≤i≤p W (N ) f (i) = α∈NC(p) b block of α c1 f is constant on b .
Hence, mixed cumulants vanish (implying freeness, see [NS06, Theorem 11.16]), and the distribution of the limiting variables x AB , x AC is Marčenko-Pastur of parameter c, ending the proof.

Remark 3.7. As a special case of the result above, one can consider the case where all the Hilbert spaces have, up to constants, the same dimension dim

H A = c 1 N , dim H D = c 4 N and dim H B = dim H C = N .
Then, the normalized marginals (c 1 N 2 ) -1 (W AB , W AC ) converge in moments, as N → ∞, towards two free elements having MP c 4 /c 1 distribution. The multi-partite equivalent of this result will be considered in Section 4.2. Similarly, when the subsystems A and D have fixed dimension dim

H A = k, dim H D = l and dim H B = dim H C = N the normalized marginals (kN ) -1 (W AB , W AC ) converge in moments, as N → ∞, towards two free elements having MP l/k distribution.
Remark 3.8. One can interpret the asymptotical freeness of the two marginals W AB and W AC in the following way: the two marginals behave as if they come from independent random tensors X and Y :

W AB = [id ⊗ id ⊗ Tr ⊗ Tr](XX * ) and WAC = [id ⊗ Tr ⊗ id ⊗ Tr](Y Y * ).
Indeed, for the random matrices above, the conclusion of the theorem above follows from the very general asymptotic freeness results of Voiculescu: the random matrices W AB and WAC are independent and unitarily invariant, and they converge to Marčenko-Pastur elements. One can understand this parallel using the fact that, in the asymptotical regime under consideration here, the amount of fresh randomness (the Hilbert spaces H B,C ) is growing. This situation is to be contrasted with the behavior of the marginals in the fixed B, C regime discussed in Section 3.4.

Remark 3.9. The remark above has interesting applications to quantum information theory. As the quantum marginals ρ AB and ρ AC are rescaled versions of the Wishart matrices W AB,AC , the same asymptotic freeness result holds, with a different scaling (precisely, it is the random matrices N D N ρ AB , N D N ρ AC which are asymptotically free). Thus, the previous remark implies that, in the large N limit, the quantum marginals ρ AB and ρ AC "forget" that they are marginals of the same quantum state |ψ ABCD and behave like independent random density matrices from the induced ensemble with parameter c = lim N D /N A . In particular, these marginals become uncorrelated asymptotically, the intuition for this fact being that the amount of "fresh randomness" from the systems H B and H C , which have dimension growing to infinity, is enough to erase the correlations from system H A , and this independently on the ratio dim H A vs. dim H B,C .

As an application of Theorem 3.6, let us consider the product of the two marginals W AB and W AC , or, to be exact, its self-adjoint version

P := W 1/2 AB W AC W 1/2
AB . Applying Theorem 3.6, the random matrix P converges in moments to the element x

1/2 AB x AC x 1/2
AB , where x AB , x AC are two free elements having distribution MP c . In free probability theory, the (self-adjoint) multiplication operation of free elements is known as the free multiplicative convolution (denoted by ), see [NS06, Lecture 14]. In our case, we are interested in the probability measure MP c MP c = MP 2 c . Exact formulas for the densities of the above distributions have been computed in [P Ż11] in the case c = 1 and in [DC14] in the general case. We compare Monte Carlo simulations to the exact densities in Figure 13. We discuss next another formulation of the asymptotic moment formula (56), useful in practice when one has to evaluate specific mixed moments in the AB/AC marginals. Proposition 3.10. In the same setting as Theorem 3.6 with c = 1, the asymptotical mixed moments of the marginals W AB , W AC are indexed by arbitrary words f ∈ {AB, AC} p , or, equivalently, by two integer vectors r, s:

←- 1≤i≤p W f (i) = W rm AB W sm AC • • • W r 2 AB W s 2 AC W r 1 AB W s 1 AC .
Asymptotically, we have Proof. The result follows from Theorem 3.6 and from the formula for the moments of a product of free random variables [NS06, Theorem 14.4]: denoting by x AB and x AC the limits in distribution of W AB , resp. W AC and by tr the expectation in the non-commutative probability space where x AB and x AC live, we have tr(

lim N →∞ (N A N ) -p-1 E Tr W rm AB W sm AC • • • W r 1 AB W s 1 AC = σ,π∈NC(m), σ≤π
x rq AB x sq AC • • • x r 1 AB x s 1 AC ) = σ∈NC(q) tr σ (x rq AB , . . . , x r 1 AB )κ σ Kr (x sq AC , . . . , x s 1 AC ) = σ∈NC(q)
Cat σ (r)

π ≤σ Kr tr π (x sm AC , . . . , x s 1 AC ) Mob(π , σ Kr ) = σ≤π Cat σ (r) Cat π Kr (s) Mob(σ, π),
where we have used π = π Kr .

As an application of the proposition above, we give below the explicit moments in the two simplest cases. When q = 1, writing p = r 1 + s 1 , we have just one term in the sum (corresponding to σ = π = ) lim

N →∞ (c 1 N 2 ) -p-1 E Tr W r 1 AB W s 1 AC = Cat r 1 Cat s 1 . When q = 2
, write similarly p = r 1 + s 1 + r 2 + s 2 ; this time, the sum contains three terms, corresponding respectively to σ = π = , {σ = , π = }, and σ = π = :

lim N →∞ (c 1 N 2 ) -p-1 E Tr W r 2 AB W s 2 AC W r 1 AB W s 1 AC = Cat r 1 Cat r 2 Cat s 1 +s 2 -Cat r 1 Cat r 2 Cat s 1 Cat s 2 + Cat r 1 +r 2 Cat s 1 Cat s 2 .
Remark 3.11. The number of terms in equation (58) is given by the number of pairs (σ, π) ∈ NC(m) 2 such that σ ≤ π. The set of all such pairs is known as the set of 2-chains (or intervals) in the lattice of non-crossing partitions and has been enumerated by Kreweras in [Kre72]: their number is given by the Fuss-Catalan numbers of order 2 (sequence A001764 in [START_REF] Sloane | The on-line encyclopedia of integer sequences[END_REF])

FC 2 (q) := 1 2q + 1 3q q .
Remarkably, these numbers are also the q-th moments of the free multiplicative square of the Marčenko-Pastur distribution of parameter 1:

∀q ≥ 1, x q dMP 2 1 (x) = 1 2q + 1 3q q .
Bijection with trees. At the end of Section 2 (Fig. 6), we described a bijective mapping between planar bicolored maps in M 0 p . Naturally, this bijection still applies in the present context, as the sum in (57) still involves elements of M 0 p . However we now have additional coloring information, which will translate into a coloring of the trees. The edges carry a color in {B, C}, but because the maps involved in (57) satisfy ∆ f (M) = 0, all the edges incident to the same white vertex share the same color. We color each white vertex with the color of the incident edges. We add a vertex in every face, and new edges linking it to the corners of the face, so that the result remains planar, and delete the initial edges and the black vertex. A tree is always bicolored, but here we see that the vertices of the corresponding tree are partitioned into two sets: white vertices, and vertices carrying a color B or C, such that edges may only link white vertices to colored vertices. Furthermore, the labeling of the edges now translates into a labeling of the colored corners of the tree: starting on corner 1 and following the face around the tree counter-clockwisely, we alternatively encounter corners incident to white vertices, and corners incident to colored vertices, labeled with p, then p -1, . . . until returning to corner 1 (the ordering is reversed with respect to the initial bicolored map, because of the dual bijective mapping).

All the labeled colored trees are not images of the colored bicolored maps: the corner corresponding to the edge number a should be incident to a vertex of color j(a). Therefore, starting from corner 1 and going around the tree in the clockwise direction, we should encounter corners incident to vertices of color j(1), then j(2), then j(3), and so on.

We denote T f the set of tricolored trees with white vertices, and vertices of color B or C, and with a such that the colors of the colored vertices encountered in the counterclockwise face are (j(1), j(p), j(p -1), . . . j(1)), which is also the reverse word f (adding color A on every vertex). In the case where c = 1, the quantity lim

N →∞ E(N A N ) -p-1 Tr W rq AB W sq AC • • • W r 1 AB W s 1 AC of Prop.
3.10 therefore counts the number of tricolored trees in T p . With the notations of Prop. 3.10, we therefore have:

|T f | = σ,π∈NC(q), σ≤π Cat σ (r)Cat π Kr (s) Mob(σ, π). ( 59
)
In particular, this is always a non-negative quantity.

Remark that because of the tree structure, we can easily find recursive relations for the quantity |T f |. 8 In order to write the recursive relations, we rather denote 9

W d 1 ,...,dq u 1 ,...,uq = |T f | = lim N →∞ (N A N ) -p-1 E Tr W u 1 AB W d 1 AC • • • W uq AB W dq AC . ( 60
)
8 Note that these relations can also be found from the Schwinger-Dyson equations applied to the matrix formulation of the moments. 9 For practical reasons, the labeling is reversed with respect to our usual convention.

We find:

W d 1 ,...,dq+1 u 1 ,...,uq = q-1 k=1 d k -1 s=0 W d 1 ,...,s u 1 ,...,u k W d k+1 ,...,dq+d k -s u k+1 ,...,uq + dq-1 s=0 W d 1 ,...,s u 1 ,...,uq C dq-s (61) W d 1 ,...,dq u 1 ,...,uq+1 = q k=1 u k -1 s=0 W d 1 ,...,d k-1 ,dq u 1 ,...,u k-1 ,s W d k ,d k+1 ,...,d q-1 uq+u k -s,u k+1 ,...,u q-1 + uq-1 s=0 W d 1 ,...,dq u 1 ,...,s C uq-s . (62)
However, solving these relations directly is a difficult task, and the solution is found considerably more easily using the techniques of the proof of Prop. 3.10.

3.3.

Comparing with the Fuss-Catalan matrix model.

In this subsection, we would like to compare the matrix model discussed above in the balanced regime where c 1 = c 4 = 1 (i.e.

H A = H B = H C = H D = C N )
with another one having the same asymptotic moments (the Fuss-Catalan numbers), the so-called free Bessel laws of parameter 2 from [START_REF] Banica | Free bessel laws[END_REF]. More precisely, the latter matrix model is given by Q = X 1 X 2 X * 2 X * 1 , where X 1,2 are i.i.d. N × N complex Gaussian random matrices. The exact moments of the random matrix Q are given by

E Tr Q p = α 1 ,α 2 ∈Sp N #(γα 1 α 2 )+#α 1 +#α 2 . ( 63
)
Indeed, in our usual representation in maps, the sum is taken over bicolored maps with one black vertex, such that white vertices only have incident edges of the same color. In that sense, white vertices inherit the color 1 or 2 of their incident edges. If α 1 (resp. α 2 ) is the permutation whose cycles encode the vertices of color 1 (resp. 2), then we have one orbit for each white vertex (so #α 1 + #α 2 ), and one for each face of the map (so #(γα 1 α 2 )). We wish to compare the expression (63) with the exact moments (46) of the model studied in Section 3.1, which we report here for N A = N B = N C = N D , and for

P = W 1/2 AB W AC W 1/2
AB : denoting f 0 the length 2p word AB, AC, AB, AC, . . . (there are no two consecutive edges of the same color),

E Tr P p = N 2(2p+1) M∈M 2p N -2g(M)-2g(M f 0 )-2∆ f 0 (M) . ( 64
)
Using that the Euler characteristics of the map writes F (M) + V (M) = 2p + 2 -2g(M), we notice that (63) is very similar to (64), provided that we impose the condition

∆ f 0 (M) = 0. ( 65
)
In that case, M = M f 0 , so that g(M f 0 ) = g(M). We rewrite

E Tr Q p = N 2p+1 M∈M 2p ∆ f 0 (M)=0 N -2g(M) . ( 66
)
As a consequence, the asymptotic moments of the two matrix models are identical, but the lower orders are different:

EN -2 Tr(N -4 P ) = 1 + N -2 EN -2 Tr(N -4 P ) 2 = 3 + 8N -2 + 8N -4 + 5N -6 EN -2 Tr(N -4 P ) 3 = 12 + 54N -2 + 135N -4 + 278N -6 + 170N -8 + 71N -10 , EN -1 Tr(N -2 Q) = 1 EN -1 Tr(N -2 Q) 2 = 3 + N -2 EN -1 Tr(N -2 Q) 3 = 12 + 21N -2 + 3N -4 . g = 0 ∆ f 0 = 0 g f 0 = 0 N 0 g = 0 ∆ f 0 = 2 g f 0 = 0 N -4 2× g = 0 ∆ f 0 = 1 g f 0 = 1 N -4 g = 1 ∆ f 0 = 1 g f 0 = 0 N -4 4× g = 1 ∆ f 0 = 0 g f 0 = 1 N -4 g = 1 ∆ f 0 = 1 g f 0 = 1 N -6 4× 4× g = 0 ∆ f 0 = 1 g f 0 = 0 N -2
Figure 15. Maps contributing to EN -2 Tr(N -4 P ) 2 . We have denoted g f 0 = g(M f 0 ). The edges labeled 1 is always the upper left one (blue). The only maps contributing to EN -1 Tr(N -2 Q) 2 are those with ∆ f 0 = 0, however in that case, the contribution in N is only given by the genus g of the map.

The coefficients of the lower orders in N of P are bigger than the ones of Q. This originates from the fact that in the Q case, the lower orders come from maps in M 2p with ∆ f 0 (M) = 0 with nontrivial genus, while in the case of P , lower orders are obtained either from higher genus combinatorial maps 10 or combinatorial maps containing both white vertices adjacent to both type AB and type AC edges (i.e. with non-vanishing ∆ f 0 ). In Figure 15, we show all the maps contributing to EN -2 Tr(N -4 P ) 2 . The values of g, ∆ f 0 and g f 0 = g(M f 0 ) are shown. For instance the maps in the upper right box of Fig. 15 are planar but contain one white vertex adjacent to one edge of type AB and one edge of type AC, and therefore contribute to EN -2 Tr(N -4 P ) 2 at order N -2 . The maps contributing to EN -1 Tr(N -2 Q) 2 are those for which ∆ f 0 = 0. Their contribution to

EN -1 Tr(N -2 Q) 2 is N -2g , but their contribution to EN -2 Tr(N -4 P ) 2 is N -4g , because g(M f 0 ) = g(M).
Keeping in mind the aim of comparing the two matrix models mixed moments with the moments of P we start by defining the following operation on combinatorial maps Definition 3.12. Let M 1 , M 2 ∈ M(p) × M(p ) be two combinatorial maps with respectively p and p edges. Both maps have a labeling of the edges. We define the gluing convolution as

: M(p) × M(p ) → M(p + p ) (67) (M 1 , M 2 ) -→ M ( 68
)
where M is the empty map if p = p and is otherwise obtained from M 1 and M 2 by stacking their black vertices one onto the other, in such way that the edges of M 2 slip into the corners of the black vertex of M 1 and edges of M 1 and M 2 alternate around the black vertex of the newly created map M. In order to select a unique way to perform this operation, we ask that the edge with label i in M 1 is followed by the edge labeled i in M 2 when following the ordering of the edges around the black vertex of M. Graphically one obtains the local construction shown on Fig. 16.

10 It is also important to notice that in the P case the genus appears in two ways in the exponent of N -once for the combinatorial map M and once for the combinatorial map M f 0 . This has to be taken into account when comparing the Q case against the P case.

G

1 2 p p-1 M 1 G' 1 2 p-1 p M 2 = G G' M 1 M 2 Figure 16
. Gluing convolution of two maps. The gluing convolution is made locally around the black vertex. The rest of the maps stays untouched.

We also define the splitting of a vertex.

Definition 3.13. A vertex-splitting is a local move on a vertex with at least two corners of a combinatorial map. It is performed by choosing two corners of the considered vertex and splitting the vertex along a straight line between these two corners.

→ → →

Figure 17. A vertex-splitting move applied to a white vertex.

Notice that the vertex splitting move can be understood in the language of permutations. In this language, a p-valent white vertex is a cycle of length p, (α c ) of α = c∈α (α c ). Such a cycle writes (α c ) = (a 1 a 2 . . . a p ). Choosing two corners of a white vertex amounts to picking two elements a k , a k in (α c ) such that a k = a k . This corresponds to the choice of the two corners of the white vertex located between edges a k , a k+1 and a k , a k +1 . These elements appear in the cycle (α c ) = (a 1 a 2 . . . a k a k+1 . . . a k a k +1 a p ). The splitting of a white vertex is just the composition of (α c ) with the transposition (a k a k ), as (α c )(a k a k ) = (α c )(α c ) where (α c ) = (a 1 a 2 . . . a k a k +1 . . . a p ) and (α c ) = (a k+1 . . . a k ) forms the two new white vertices.

We consider the partial order on M 0 p defined as follows, for M, M ∈ M 0 p two planar combinatorial maps, we say that M ≤ M if and only if there exists a finite sequence of maps {M i } F i=0 such that M 0 = M , M F = M, and M i+1 can be obtained from M i by applying a vertex-splitting move on one white vertex of M i .

We now define the Tutte dual of a map. It is a particular case of one of the Tutte bijections for bicolored maps [Tut75], as was the bijection presented in Section 2, Fig. 6. Though it can be defined for more general sets of maps we assume M ∈ M p , Definition 3.14. The Tutte dual T (M) of M ∈ M p is obtained as follows. In every face of M we apply the following rules

• We draw a white vertex inside the face under consideration.

• We draw a new edge between each corner adjacent to the black vertex inside this face and the newly created white vertex inside this face. Once we have followed this procedure for every face of M, we erase all initial edges and white vertices. We notice that the Tutte dual preserves the genus of the map, that is to say that g(M) = g(T (M)). Remark that the only difference with the bijection presented in Section 2, Fig. 6, is that the new edges are added between the new vertices and the black vertex, instead of the white vertices. A consequence of this choice is that after the bijection of Fig. 6, one no longer has elements of M p , while we do in the present case, thus the name "dual". We report the reader to [Tut75] for the general bijection encoding both cases. We then have the following proposition Proposition 3.15. Let M 1 , M 2 be two planar maps in M 0 p (i.e. g(M 1 ) = g(M 2 ) = 0). Then g(M 1 M 2 ) = 0 if and only if M 2 ≤ T (M 1 ).

Proof. First notice that by construction of T

(M 1 ), g(M 1 M 2 ) = 0 if M 2 = T (M 1 ).
Then when performing a vertex-splitting move on a white vertex of M 2 in M 1 M 2 the number of vertices is raised by one V → V = V + 1, the number of edges stays the same, and the number of faces is decreased by one F → F = F -1 as a consequence the genus stays constant under such a move. Thus we have M 2 ≤ T (M 1 ) ⇒ g(M 1 M 2 ) = 0. Assume now that M 2 ≤ T (M 1 ). Therefore there exists a vertex in M 2 adjacent to two edges of M 2 in M 1 M 2 whose starting points are at corners that belong to two different faces of M 1 . Thus edges of M 1 and M 2 cross in M 1 M 2 , which is equivalent to saying that g(M 1 M 2 ) = 0. Thus we have

M 2 ≤ T (M 1 ) ⇐ g(M 1 M 2 ).
If we come back to the definition of the Tutte dual we notice that it is is equivalent to the Kreweras complementation. This similarity is pictured on the Fig. 18, where the graphical representation of non-crossing partition has been borrowed from [NS06]. Note however that the Tutte dual is defined for every maps not just the planar ones (the latter corresponding to non-crossing partitions). We can write

E Tr P p = N 4p+2 M∈M 2p N -2g(M)-2(g(M f 0 )+∆ f 0 (M)) (69) E Tr Q p = N 2p+1 M 1 ,M 2 ∈M(p) N -2g(M 1 M 2 ) (70) = N 2p+1 M∈M 2p : ∆ f 0 (M)=0 N -2g(M) , (71)
The proof of the equations (70), (71) is straightforward. From (70) we see the connection with the formula (63), while from (71) we see why the moments of P are larger than those of Q: the terms in the sum for the moments of Q (71) are a subset of the terms in (69). Indeed, we only have

M 2p | ∆ f 0 =0 = M p M p ,
N →∞ 1 N 2p+1 E Tr Q p = M 1 ∈Mp M 2 ≤T (M 1 ) 1 , (72)
which mimics, in the language of maps, the expression obtained in [NS06, Lecture 14, Theorem 14.4] with cumulants equal to one in the language of permutations. This translates the fact that in the large N limit the moments of Q are the moments of the multiplicative convolution of two Marčenko-Pastur laws.

3.4. The unbalanced asymptotical regime.

In this section, we study the asymptotical regime where the Hilbert spaces H B and H C have fixed dimension, while the dimensions of H A and H D grow to infinity. To be more precise, we assume in this section that

• dim H B = dim H C = m, for some positive integer constant m; • dim H A = N , dim H D = cN , for some constant c ∈ (0, ∞), where N → ∞.
Contrary to the results proven in Section 3.2, in this setting, the random matrices W AB and W AC are no longer asymptotically free. One can understand this fact, stated precisely in the theorem below, by noticing that the "shared randomness" between the two random matrices (dim H A = N → ∞) is much larger than the "fresh randomness" (dim

H B = dim H C = m, fixed).
Theorem 3.16. In the asymptotical regime described above, the pairs of random matrices (mN ) -1 W AB , (mN ) -1 W AC converge in distribution, as N → ∞, to a pair of non-commutative random variables (x AB , x AC ) having the following free cumulants:

κ(x f (1) , x f (2) , . . . , x f (p) ) = cm -alt(f ) , (73)
where f ∈ {AB, AC} p is an arbitrary word in the letters AB, AC, and alt(f ) is the number of different consecutive values of f , counted cyclically: 11 The extension of the operation to sets of maps is understood in the straightforward way as the set of all maps with 2p edges that can be obtained by making the gluing convolution of two maps with p edges.

alt(f ) := |{a : f (a) = f (a + 1)}|,
where f (p + 1) := f (1). Equivalently, for any word in the two marginals f ∈ {AB, AC} p , we have

E Tr W f = (1 + o(1))(mN ) p+1 α∈NC(p) c #α m -alt(f,α) , (74)
where alt(f, α) has been defined in (51) as

alt(f, α) = {a ∈ 1, p : f (a) = f (α(a))} .
Proof. We begin by analyzing the exact moment formula of Theorem 3.1:

E Tr W f = α∈Sp N #(γα) cN #α m L(f,α) .
Note that since m is fixed, the dominating terms correspond to permutations α maximizing the exponent #α + #(αγ); these permutations have been shown before to be exactly the non-crossing ones, so we have

E Tr W f = (1 + o(1))N p+1 α∈NC(p) c #α m L(f,α) .
For α non-crossing, we know from Proposition 3.5 that L(f, α) = p+1-alt(f, α), and the conclusion concerning the cumulants follows from Speicher's moment-cumulant formula.

Remark 3.17. In the degenerate case m = 1 (i.e. there are no B and C systems), one has that W AB = W AC = W A is a Wishart matrix of parameters (N, cN ), and one recovers the result from the classical Marčenko-Pastur setting, Proposition 2.4: the free cumulants of x AB = x AC are all equal to c (see equation (8) and the comments following it).

Remark 3.18. If, after taking the limit N → ∞ in the theorem above, one takes the limit m → ∞, we recover the result of Theorem 3.6. Indeed, the only free cumulants surviving are the ones with alt(f ) = 0 (i.e. mixed cumulants vanish), proving that x AB and x AC are asymptotically free.

Remark 3.19. Elaborating on the preceding remarks, we notice that the eigenvalues distribution of the matrix

P = W 1/2 AB W AC W 1/2
AB for m = 1 is the distribution of the square of the eigenvalues of a Wishart matrix, the corresponding density, that we denote dρ m=1,c (x) writes dρ m=1,c (x) = dMP c ( √ x), while if m → ∞ we end up with dρ m=∞,c (x) = dMP 2 c (x). This fact does not depend on the rate at which m is sent to infinity with respect to the N limit. Thus the parameter m allows one to canonically interpolate between the two different distributions dMP c ( √ x) and dMP 2 c (x), the former being the distribution of the square of eigenvalues of a Wishart random matrix, while the latter is the free multiplicative convolution of two Marčenko-Pastur distribution. Notice also, that one can also consider the (non-selfadjoint) matrix P = W AB W AC which has the same moments and eigenvalues to reach the same conclusion. Finally, we plot Monte-Carlo simulations of the eigenvalues of P versus the square (resp. the free multiplicative square) of a Marčenko-Pastur distribution in Figure 19.

We record below some low mixed moments in the variables x AB , x AC (we denote by tr the expectation in the non-commutative probability space where these variables live): Let us now explore the consequences of the formulas for the mixed free cumulants from this section to quantum information theory. As in Section 3.2, the two quantum marginals ρ AB and ρ AC , when properly rescaled, converge in moments, jointly, to a pair of non-commutative random variables having free cumulants as in (73): N D N (ρ AB , ρ AC ) → (x AB , x AC ). This allows one to compute the asymptotic value of any correlation function involving ρ AB and ρ AC . For example, the rescaled overlap between the matrices converges to lim

tr(x AB x AC ) = c 2 + c m 2 tr(x AB x AB x AC ) = c 3 + c 2 + 2c 2 m 2 + c m 2 tr(x AB x AB x AC x AC ) = c 4 + 2c 3 + c 2 + 4c 3 m 2 + 4c 2 m 2 + c m 2 + c 2 m 4 tr(x AB x AC x AB x AC ) = c 4 + 2c 3 + 4c 3 m 2 + 4c 2 m 2 + 2c 2 m 4 + c m 4 .
N →∞ N A N ρ AB , ρ AC = 1 c 2 κ(x AB , x AC ) = 1 + 1 cm 2 .
The computation above should be compared with the similar limit from the balanced regime of Section 3.2, where the two marginals were uncorrelated: lim N →∞ N A N ρ AB , ρ AC = 0.

The general multipartite case

In this section we consider the more general situation of a random Wishart tensor defined on a Hilbert space which is factorized in an arbitrary number of factors. The section consists of three parts: we first derive the general, non-asymptotic mixed moment formula, and then consider two asymptotic regimes: the balanced regime, where all tensor factors have the same dimension, and the unbalanced regime, where some of the tensor factors (the ones corresponding to the "moving legs") are being kept fixed. In the balanced case, we prove that the marginals are asymptotically free (Proposition 4.6 in the Wishart setting and Theorem 4.13 in quantum information language), while in the unbalanced case, we show in an example that it is not possible to factorize the expression of the mixed cumulant functions over the cycles of the non-crossing partitions. Indeed they depend more finely on the structure of the non-crossing partitions. This implies that we cannot give an expression for the mixed free cumulants. However, we expect that this situation can be dealt with in the framework of free probability with amalgamation. Such results will be presented in a following paper.

We consider complex tensor X of size

N 1 × N 2 × • • • × N n , an un-normalized quantum state in the said Hilbert space C N 1 ⊗ • • • ⊗ C N 2 .
The density matrix of the corresponding pure state is X ⊗ X * , the (un-normalized) unit rank projection on the space CX. For a given set I ⊂ {1, . . . , n}, we denote I = {1, . . . , n} \ I, and define the reduced density matrix as the tensor X. I X obtained by summing, for each i ∈ I, the index of position i of X with the index of position i of X, X. I X = [id I ⊗ Tr I ](XX * ). This partial contraction of two tensors can also be understood as the matrix [X. I X], whose first (resp. second) sub-index of position j ∈ I is the free-index of position j of X (resp. X). For instance, for n = 4, choosing I = {3, 4},

[X • 3,4 X] i 1 ,i 2 ; i 1 ,i 2 = N 3 i 3 =1 N 4 i 4 =1 X i 1 ,i 2 ,i 3 ,i 4 Xi 1 ,i 2 ,i 3 ,i 4 . ( 75
)
There is a canonical one-to-one correspondence which maps I to {1, . . . , |I|} while preserving the ordering of natural integers. We denote S |I| the set of permutations of |I| elements. In the following, we implicitly make use of these canonical bijections when saying that a permutation σ ∈ S |I| acts on I and has I of same cardinality as an image. For instance if I = {A, B, C} and I = {A, C, E}, the identity id :

I → I is understood as the map A → A, B → C, C → E.
Given I ∈ N, a permutation σ ∈ S I , and two matrices P and Q, whose two indices have I sub-indices, we define the product P • σ Q as the twisted contraction

P • σ Q i 1 ,••• ,i I ; i 1 ,••• ,i I = j 1 ,••• ,j I j 1 ,••• ,j I I b=1 δ j σ(b) j b P i 1 ,••• ,i I ; j 1 ,•••j I Q j 1 ,••• ,j I ; i 1 ,••• ,i I . (76)
We define the associated trace Tr • σ accordingly. We are interested in computing expectations of the form

E Tr • σp [X. Ip X] • σ p-1 . . . • σ 3 [X. I 3 X] • σ 2 [X. I 2 X] • σ 1 [X. I 1 X], (77)
for some integer p, some non-necessarily distinct sets I a which all have the same number of elements |I a | = I, and some permutations σ a ∈ S I , (with our convention, σ a : I a → I a+1 , and I p+1 = I 1 ). We denote W I = [X. I X], so that the objects under focus are rewritten as

E Tr σ W f ({N i }) = E Tr • σp W Ip • σ p-1 . . . • σ 3 W I 3 • σ 2 W I 2 • σ 1 W I 1 , (78)
where we respectively denoted f and σ the ordered lists f = [I 1 , . . . , I p ] of p subsets of 1, n , and σ = [σ 1 , . . . , σ p ] of p permutations in S I . Note that depending on the {I a } and the {σ a }, all choices are not possible for {N 1 , . . . , N n }. In general the above defined object is always meaningful when

N 1 = • • • = N n ,
but interesting cases can be considered for specific {I a } and {σ a }.

To highlight this, we separate for a given f , the colors which are traced for every W Ia = [X. Ia X], which, without loss of generality, we can suppose to be the colors from r to n, {r, . . . , n} = p a=1

I a . (79)
We also suppose that the colors from 1 to l < r are the common fixed points of all the {σ a },

∀i ∈ 1, l , ∀a ∈ 1, p , σ a (i) = i, (80)
where we use •, • to denote integer intervals, see (34). For each a ∈ 1, p , among the colors l + 1, r -1 , the colors I a \ r, n = I a ∩ l + 1, r -1 are traced, and σ a only acts non-trivially on the colors J a = I a \ 1, l = I a ∩ l + 1, r -1 , (81) whose cardinal we denote k = I -l.

(82) Note that the sets J a generalize the color j(a) of Section 3.1. In the 4-partite case ABCD we considered in Section 3.1, we had l = 1 (which corresponded to color A), r = 4 (which corresponded to color D), and k = 1. The only non-trivial action of the permutations σ was either J a = {2} ≡ {B} or J a = {3} ≡ {C}, color which was denoted by j(a). The notations in the general case are illustrated in Fig. 20.

n 1 l l + 1 r -1 r X X * n 1 l l + 1 r -1 r X X * n 1 l l + 1 r -1 r X X * Figure 20.
A typical sample of an expectation (78). The colors linked by red edges belong to the corresponding J a s, while the blue edges link colors in the I a \ 1, l .

Before moving on to the moment computation, let us point out that the data defining the moment (σ, f ) could be replaced by a single list of ordered subsets of 1, n . 4.1. Exact expression for the moments. Proposition 4.1. We suppose that

N l+1 = • • • = N r-1 = N J 12
. Then, with the previous notations,

E Tr σ W f ({N i }) = α∈Sp l i=1 N #(γα) i n j=r N #α j N L(f ,σ,α) J , ( 83
)
in which γ = (12 • • • p), and L(f , σ, α) is a combinatorial function defined in (91).

Proof. As in Section 3.1, the moments will be expressed as a sum over Wick wirings α ∈ S p -or equivalently over maps in M p -of some weight. Each loop in the box representation contributes to this weight with a factor N i . We want to describe these loops as orbits as was done in Def. 3.2, i.e. in terms of the cycles of some permutations acting on the product of r -l -1 copies of 1, p , of the form [p] i = {(1, i), . . . , (p, i)} for each color i in l + 1, r -1 .

As detailed several times along this paper, to every permutation α ∈ S p corresponds a combinatorial map M with one black vertex and #α white vertices, and whose p labeled edges are disposed from 1 to p counterclockwise around the black vertex and correspond to matrices

W Ia = [X. Ia X].
The edge labeled a therefore carries the set I a of I colors, l of which belong to every edge.

Let us take a closer look at the orbits in the map formulation. As was previously the case in Section 3.1, an orbit which has color i, when it arrives on a white vertex, leaves this vertex on the next edge also carrying the color i, counterclockwise. When an orbit of color i arrives on a black vertex from the edge labeled a, it goes to the following edge around that vertex counterclockwise, labeled a + 1, but changes color to σ a (i).

For the colors 1 to l, these are actually the usual faces of the map, F (M), as was the case for color A in Section 3.1. The colors 1, l therefore contribute with a factor l i=1 N

F (M) i = l i=1 N #(γα) i
. We may therefore as well forget these colors, and label the edges with the sets J a instead of I a (as 12 We put independent Ni's for the colors i ∈ 1, n which are fixed points of all the initial wirings or which are always traced, and put a common NJ for the others. We stress however that more general cases might possibly be considered, if the colors can be separated into two sets if l + 1, r -1 = K1 K2, such that the support of any orbit is either included in K1, or in K2. We may then choose a different NJ 1 and NJ 2 for colors in K1 and K2.

was done for j(a) in Sec. 3.1).

The colors which belong to the sets I a will be taken care of further, and we now focus on the remaining colors, which belong to at least one set J a . In this general case, their behavior might differ from the particular case previously treated in Section 3.1. Indeed, for colors in r + 1, l -1 , a given edge might appear several times on the same orbit (at most k times). Therefore the orbits cannot in general be defined as the cycles of a permutation of the p edges, or equivalently the faces of a combinatorial map without colors, which was the key point in Section 3.1. In order to bypass this difficulty, we define a labeling for pk copies of the edges, one per each couple (a, i), where a ∈ 1, p labels the edge, and i ∈ J a is a color which is neither traced, nor a fixed point of all the σ b . On these pk elements, we define the following permutation

Γ f ,σ : (a, i) → (a + 1, σ a (i)). ( 84
)
We label each cycle of the permutation starting from the smallest color of the edge of smallest label in the cycle. For instance, if J 1 = {B, C}, J 2 = {B, D}, J 3 = {C, D}, σ 1 = id, σ 2 = id, and σ 3 is the transposition in S 2 , then

Γ f ,σ = (1, B), (2, B), (3, C), (1, C), (2, D), (3, D) , (85)
and if J 1 = {B, C}, J 2 = {B, D}, J 3 = {B, C}, J 4 = {C, D}, σ 1 = id, σ 2 = id, and σ 3 and σ 4 are the transposition in S 2 , then

Γ f ,σ = (1, B), (2, B), (3, B), (4, D) (1, C), (2, D), (3, C), (4, C) . (86)
For colors in l +1, r -1 , the behavior of the faces around white vertices is similar to Section 3.1, and the duplication operation of Figure 10 generalizes locally: performing the following duplication operation (illustrated for I = 3) on every white vertex does not change locally the incident external orbit, and removes the color conditions on the white vertices. This would precisely create one copy of an edge for each one of the kp couples (a, i). More precisely, if i ∈ J a , denoting α i (a) the first edge following a around the white endpoint of a and containing color i, α i (a) = α q (a), where q = min{s ∈ N

* | i ∈ J α s (a) }, (87)
(q depends on α, i, and a) then the permutation defining the resulting white vertices is

α f : (a, i) → (α i (a), i). (88)
We can now define a (non-necessarily connected) combinatorial map M f ,σ = (Γ f ,σ , α f ) from the two permutations Γ f ,σ and α f , acting on 1, p × 1, k . See the examples in Fig. 22 and23.

Among the orbits which act on colors in l + 1, n , one has to take into account the orbits which are entirely included in the white vertices of M. The number of such orbits contained in a given vertex is l, plus the number of colors in l + 1, r -1 which do not appear on any edge. The total number of the latter can be expressed as

(r -l -1) V (M) -1 -V (M f ,σ ) -V • M f ,σ )), (89)
where we have denoted the number of black vertices of M f ,σ as

V • M f ,σ ) = V • (f , σ) = #Γ f ,σ .
This is an integer between 1 and k which does not depend on α. These orbits therefore contribute with a factor

n i=r N V (M)-1 i N (r-l-1)(V (M)-1)-(V (M f ,σ )-V•(f ,σ)) J
, where we recall that 13 V (M) -1 = #α.

The remaining orbits run along at least one edge. Because the permutations Γ f ,σ and α f encode precisely the way they behave locally around the black and the white vertices respectively, their total number is precisely given by the number of faces of the map M f ,σ , 13 Notice also that if one wants to make explicit the analogy between M and M f ,σ in this expression, one can rewrite it as n i=r N

F (M f ,σ ) = #(Γ f ,σ α f ),
V (M)-1 i N (r-l-1)(V (M)-V•(M))-(V (M f ,σ )-V•(f ,σ)) J or n i=r N V (M)-1 i N (r-l-1)(V (M)-#γ)-(V (M f ,σ )-#Γ f ,σ ) J
hence contributing with a factor N

F (M f ,σ ) J
. Putting all of this together, we have shown that

E Tr σ W f ({N i }) = M∈Mp l i=1 N F (M) i n j=r N V (M)-1 j N L(f ,σ,M) J , (90)
where we have defined

L(f , σ, M) := F (M f ,σ ) + (r -l -1) V (M) -1 -V (M f ,σ ) -V • M f ,σ )). ( 91
)
Identifying L(f , σ, α) = L(f , σ, (γ, α)), this concludes the proof.

In the above proof we defined the companion map M f ,σ to M. This companion map is a generalization of the previous companion map M f introduced in section 3.1. Since this is an important object, we recall its construction in the following definition Definition 4.2. For all integer p, we associate to the triplet (M, f , σ), with M a combinatorial map in M p , f the list of colors of the edges of M and σ the list of permutations labeling the corners of the black vertex of M, a combinatorial map M f ,σ called unfolded map from the following data:

• The set E f := (a, i) a∈ 1,p , i∈Ja is the set of edges.

• The permutations α f and Γ f ,σ , where α f is the permutation that defines the white vertices after local duplication, and Γ f ,σ defines the black vertices of M f ,σ . Denoting

α i (a) = α q (a), where q = min{s ∈ N * | i ∈ J α s (a) },
the permutation α f is defined as

α f : E f → E f (a, i) → (α i (a), i),
and Γ f ,σ is defined as

Γ f ,σ : E f → E f (a, i) → (a + 1, σ a (i)).
If we compare with Definition 2.9 then we notice that we do not require that the group α f , Γ f ,σ acts transitively on E f , this is because M f ,σ can be disconnected. Remark also that each white vertex of M is duplicated into k = |J a | white vertices of M f ,σ if all the incident edges are labeled by the same set J a ; and into strictly more than k white vertices if at most k -1 colors are common to all the edges. We therefore define the quantity

∆ f ,σ (M) := V (M f ,σ ) -V • (f , σ) -k(V (M) -1), (92)
which vanishes if and only if all the edges incident to a given white vertex of M share the same color set J a , and is positive otherwise.

Moreover, if the number of black vertices is V • (f , σ), the number of connected components K of a given map M f ,σ is an integer in {1, . . . , V • (f , σ)}, so that we define

Σ(M f ,σ ) = V • (f , σ) -K(M f ,σ ), (93)
which is an integer between 0 and V • (f , σ) -1.

Theorem 4.3. With the previous notations, the mixed moments of the marginals {W Ia } are expressed exactly using a sum over combinatorial maps

E Tr σ W f = M∈Mp l i=1 N 2+p-V (M)-2g(M) i n j=r N V (M)-1 j N kp+V•(f ,σ)+(r-l-2k-1)(V (M)-1)-L(f ,σ,M) J , (94)
in which we have denoted

L(f , σ, M) = 2 g(M f ,σ ) + ∆ f ,σ (M) + Σ(M f ,σ ) , (95)
where g(M f ,σ ) ≥ 0 is the genus of the map M f ,σ , and ∆ f ,σ (M) ≥ 0 and Σ(M f ,σ ) ≥ 0 have been defined above.

Before starting with the proof, notice that the quantity ∆ f ,σ (M) seems to penalize the mixed cumulants in the large N J regime, if we trust the heuristic relating monochromatic white vertices (i.e. white vertices incident to edges of the same color) to non-mixed cumulants. This is an indication for freeness at large N J .

Proof. The only thing we need to prove is that

L(f , σ, α) = kp + V • (f , σ) + (r -l -2k -1)(V (M) - 1) -L(f , σ, M), which follows from the Euler characteristics of the map M f ,σ , 2K(M f ,σ ) -2g(M f ,σ ) = F (M f ,σ ) -kp + V (M f ,σ ). ( 96
)
We conclude using Proposition 4.1 and the definitions of ∆ f ,σ (M) and Σ(M f ,σ ).

Remark 4.4. In the case where k = 0, the result above degenerates, and we recover the classical Marčenko-Pastur result from Proposition 2.4 (see also equation (14) for the combinatorial map approach). Indeed, in this case r = l + 1, and the factor with N J is trivially equal to 1; moreover, the quantities σ and f are trivial and do not play any role. The result reads

E Tr W n = M∈Mp l i=1 N 2+p-V (M)-2g(M) i n j=r N V (M)-1 j .
If we denote N tot := l i=1 N i and consider the scaling where N tot := n j=r N j ∼ cN tot , we obtain

E Tr W n = (1 + o(1))N 1+p tot M∈Mp N -2g(M) tot c V (M)-1 .
Assuming moreover N tot → ∞, the limit selects planar maps, and we obtain (14):

lim

Ntot→∞ E 1 N tot Tr W n N tot p = M∈M 0 p c V (M)-1 .
Remark 4.5. The first non-trivial case corresponds to k = 1, where just one of the n legs of the tensors can "move around". In this case, it is clear that the permutations σ do not play any role, and f is just a function from 1, p to the set of tensor legs, selecting for each tensor the leg which "moves around". We have f : 1, p → r ], and, since σ is trivial, V • (f ) = 1. Assuming, to keep things simple, that all the vector spaces have dimension N , the result of Theorem 4.3 reads

E Tr W f = M∈Mp N (r-1)[2+p-V (M)-2g(M)] N V (M)-1 N p+V•(f )+(r-2)(V (M)-1)-L(f ,M) = N r(p+1) M∈Mp N -2(r-1)g(M)-L(f ,M) .
In order to find the dominating contributions, we have to identify the maps M ∈ M p which cancel both g(M) and L(f , M). The first condition corresponds to M being planar, while the second one is equivalent to the cancellation of the three quantities g(M f ) ,∆ f (M), Σ(M f ) from (95). We focus on the ∆ quantity: it follows from (92) that this quantity is zero iff all the edges incident to a white vertex in M have the same color given by f ; in other words, the partition introduced by the white vertices of M on 1, p has to be smaller than the partition f introduces on the same set.

We claim now that the conditions g(M) = 0 and ∆ f (M) = 0 imply that the other two quantities appearing with a negative sign in the exponent of N cancel. Indeed, Σ(M f ) = 0, since V • (f ) = 1.

Since ∆ f (M) = 0, it also follows that M f = M, and thus g(M f ) = 0. To summarize, we reach the same conclusion as in Theorem 3.6, see equation (57):

lim N →∞ N -r(p+1) E Tr W f = |M 0,f p |,
where M 0,f p := {M ∈ M p : M is planar and f assigns the same color to the edges incident to the white vertices of M}.

The balanced asymptotical regime.

As a first step, we focus on the limit where N 1 , . . . , N n all grow to infinity, while the ratio between N 1 . . . N l and N r . . . N n converges to a fixed constant at infinity. More precisely, we can consider for instance ∀i ∈ 1, l ,

N i = N J = N while ∀j ∈ [[r, n]], lim N →∞ N j /N = c 1 n-r+1 or pick j ∈ [[r, n]] such that N j ∼ cN at infinity, while N i = N for i = j.
In that case, Theorem 4.3 writes

E Tr σ W f = N (k+l)(p-2)+l+n+V•(f ,σ) M∈Mp h(N ) V (M)-1 N -2lg(M)-L(f ,σ,M)-(2k+2l-n)(V (M)-2) , (97)
where h(N ) is an arbitrary function such that h(N ) → c when N → +∞. Notice that the factor N V•(f ,σ) is out of the sum since, as already emphasized earlier, it only depends on f , σ. The terms g(M), L(f , σ, M) and V (M) -2 are all non-negative, so we have three cases, depending on the sign of 2k + 2l -n.

If 2k + 2l -n < 0, the maps M which survive in the large N limit are those which satisfy g(M) = 0 (if l = 0), L(f , σ, M) = 0, and V (M) = p + 1.

(

) 98
Since there is a unique solution to this system, given by the only tree in M p , in that case, in the large N limit, lim

N →∞ 1 N (n-k-l)p+l+V•(f ,σ) E Tr σ W f = c p . ( 99
)
This case is therefore trivial, and we do not consider it. In the following, we first consider the case 2k + 2l -n = 0 in Subsection 4.2.1, and then the case 2k + 2l -n > 0 in Subsection 4.2.2.

4.2.1.

Case where half the colors are traced.

In this section we study the case where 2k + 2l -n vanishes, or equivalently

I = l + k = n 2 , (100)
the number of colors n is even, and for each edge, precisely half the colors are traced. The equation (97) rewrites as

E Tr σ W f = N (k+l)p+l+V•(f ,σ) M∈Mp h(N ) V (M)-1 N -2lg(M)-L(f ,σ,M) , (101)
where here h(N ) is an arbitrary function whose limit at infinity is c.

Large N limit, case l = 0. In the large N limit, the sum restricts to planar maps, which are such that edges incident to a given white vertex have the same set of colors I a , and for which the map M f ,σ is planar and has as many connected components as black vertices. We define the set of "free" maps

M (0)free f ,σ := M ∈ M p | g(M) = 0, g(M f ,σ ) = 0, ∆ f ,σ (M) = 0, and Σ(M f ,σ ) = 0 . (102)
The following result is just a consequence of Theorem 4.3.

Proposition 4.6. In the large N 1 , . . . , N n regime, if l + k = n 2 and with the previous notations,

lim N →∞ 1 N n 2 p+l+V•(f ,σ) E Tr σ W f = M∈M (0)free f ,σ c V (M)-1 . (103)
Note that since V • (f , σ) might be smaller than k, the moments corresponding to σ and f might have a weaker scaling in N than the usual N n 2 (p+1) .

General case including l = 0. In the case where l = 0, the map M is a priori non necessarily planar in the large N limit. In this paragraph, we show that maps with vanishing L are in fact necessarily planar. As a consequence, the system defining the "free" maps in (102) can in fact be reduced to

M (0)free f ,σ = M ∈ M p | g(M f ,σ ) = 0, ∆ f ,σ (M) = 0, and Σ(M f ,σ ) = 0 , (104)
thus extending the preceding result to the case l = 0. Indeed, we have the following result

Proposition 4.7. A map M ∈ M p in a triplet (M, f , σ) which satisfies g(M f ,σ ) = 0 and ∆ f ,σ (M) = 0, (105)
is necessarily planar: g(M) = 0.

Proof. Relying on the following Lemma 4.8, the proposition is proven easily by noticing that, if g(M f ,σ ) = 0, the genus g(M) is a non-negative integer which is smaller than one.

Lemma 4.8. If ∆ f ,σ (M) = 0, we have the following bound for the number of faces of the combinatorial map M f ,σ ,

F (M f ,σ ) ≤ kF (M), and g(M) ≤ g(M f ,σ ) k + k -1 k . (106)
Proof. Because ∆ f ,σ (M) = 0, the orbits, though different from the faces, cover precisely the faces of M. The difference is that if a face starts on an edge a ∈ 1, p with a color i ∈ J a , after going around M, it might come back for the first time on the same side of the edge a for a color j = i.

In that case, it would follow the face of M again a certain number of times, until it reaches again the same side of the edge a, for the same color i. An orbit, which is projected to a face of M f ,σ , is therefore a multiple of a face of M. Each side of each edge a is used exactly once per color in J a . If all the orbits go exactly once around the map, then each face of M corresponds to k orbits (and thus faces of M f ,σ ), as they never pass twice on the same side of an edge, and an edge has k colors which are not traced. If not, there are less orbits, so that the first inequality is proven.

To prove the second inequality, we consider the Euler characteristics of M, 2g

(M) = 2-V (M)+ p -F (M). Since ∆ f ,σ (M) = 0, we replace V (M f ,σ ) -V • (f , σ) = k(V (M) -1), 2kg(M) = 2k -kV (M) + kp -kF (M) (107) = k -V (M f ,σ ) -V • (f , σ) + kp -kF (M). ( 108
) Now using the Euler characteristics of M f ,σ , kp -V (M f ,σ ) = 2g(M f ,σ ) + F (M f ,σ ) -2K(M f ,σ ), we find that 2kg(M) = k + 2g(M f ,σ ) + F (M f ,σ ) -2K(M f ,σ ) + V • (f , σ) -kF (M). ( 109
)
Using the first inequality of (106), we are left with

2kg(M) ≤ k + 2g(M f ,σ ) -2K(M f ,σ ) + V • (f , σ). (110)
We conclude using that K(M f ,σ ) ≥ 1 and V • (f , σ) ≤ k.

Remark 4.9. A consequence of Proposition 4.7 is that keeping the N i finite for i ∈ 1, l does not lead to more interesting behavior. Indeed, a map with L(f , σ, M) = 0 also has vanishing genus, so that the large N limit also selects all the maps in M (0)free f ,σ

, and we recover the results of the present case.

Application to the asymptotic freeness of permuted marginals of random Gaussian tensors. We now consider an application of Theorem 4.3 to a new situation, where we consider the set of all possible marginals of the tensor X ⊗ X * , with all possible permutations acting on half of the tensor legs. More precisely, let n be even, and put, for an (n/2)-subset I of 1, n , W I = [id I ⊗ Tr 1,n \I ](XX * ), where X ∈ (C N ) ⊗n is a Gaussian tensor; note that we assume here, for the sake of simplicity, that all the Hilbert space dimensions are N . The general case of different Hilbert space dimensions, say N i ∼ c i N for some constants c i > 0 and N → ∞, can be easily obtained.

Let us introduce, for every permutation π ∈ S n/2 ,

W I,π = P π W I P -1 π ∈ M N n/2 (C) (111)
the matrix obtained by permuting the tensor legs of W I according to the permutation π. For example, in the case n = 4, there are 12 possible matrices W I,π . We present in Figure 24 This is therefore a particular case of the large N 1 , . . . , N n regime, for specific permutations σ a ∈ σ, which factorize as

∀a ∈ 1, p , σ a = π -1 a+1 π a . (112)
Note that the permutations σ a might share some fixed points (recall that we denote by l the number of these common fixed points, and we also put k = n/2 -l). Proposition 4.6 applies, and the maps that contribute to the large N limit are a subset of M (0)free f ,σ

: in particular, they are planar, and the edges that are attached to the same white vertices all have the same set of colors J a . Moreover, the hypothesis (111) has stronger consequences. Indeed, we have the following lemmas.

Lemma 4.10. If the permutations in σ factorize as in (112), for any M ∈ M p , the map M f ,σ has precisely k black vertices. As a consequence, if Σ(M f ,σ ) = 0, the maps M f ,σ have k connected components.

Proof. In the combinatorial maps setting, the picture is the following: when arriving to an edge, the color of an orbit changes by multiplying with π -1 a , and when leaving the edge, one multiplies with the inverse permutation π a . Formally, since the black vertices of M f ,σ correspond to the cycles of the permutation Γ f ,σ from (84), we show that Γ f ,σ has precisely k cycles. Indeed, we have (a, i)

Γ f ,σ ---→ (a + 1, π -1 a+1 π a (i)) Γ f ,σ ---→ (a + 2, π -1 a+2 π a+1 π -1 a+1 π a (i) π -1 a+2 πa(i) ) Γ f ,σ ---→ • • • Γ f ,σ ---→ (a + p a , π -1 a+p π a (i) i ),
for all i ∈ J a , proving the claim (the addition on the first coordinate is done modulo p).

Lemma 4.11. If M is such that ∆ f ,σ (M) = Σ(M f ,σ ) = 0, and the permutations in σ factorize as in (112), the edges incident to a given white vertex are all labeled by the same permutation π.

Proof. Consider two edges a and b incident to the same white vertex in M. As ∆ f ,σ (M) = 0, they have the same set of colors J a = J b . Take a color i ∈ J a . In the unfolded map M f ,σ , the white vertex corresponding to the color i ∈ J a has attached to it, among others, the edges (a, i) and (b, i). The edge (a, i) has its other end attached to the black vertex corresponding to the cycle containing (a, i) in the permutation Γ f ,σ . In order for the map M f ,σ to have exactly k connected components (see Lemma 4.10), the edge (b, i) must be connected to the same black vertex (otherwise, the aforementioned white vertex would be connected to two different black vertices, decreasing the number of connected components). However, from the proof of Lemma 4.10, we see that the black vertex has attached to it the following edges: {(c, π -1 c π a (i))} c∈ 1,p . Hence, we must have π -1 b π a (i) = i for all i, which is the claim. Lemma 4.12. If M is such that ∆ f ,σ (M) = Σ(M f ,σ ) = 0, and the permutations in σ factorize as in (112), the unfolded map M f ,σ consists of k copies of the original map M (when discarding the edge coloring).

Proof. From the proof of Lemma 4.10, M f ,σ has k connected components, and each one of them has one black vertex with p incident edges, ordered from 1 to p. Discarding the edge coloring, the permutation defining the black vertex of each one of the connected components is therefore γ. These p edges are attached to white vertices in M f ,σ . Now since all the edges a incident to some white vertex v of M have the same sets J a , each white vertex v of M corresponds to k = |J a | white vertices {v i } i∈Ja in the unfolded map M f ,σ , with the same edges (a, α(a), α 2 (a), • • • ) attached to it, in the same order. Therefore, discarding the edge coloring, the permutation defining the white vertices of a given connected component of M f ,σ is just α. This proves that each connected component of M f ,σ is isomorphic to M.

We arrive now at one of the main results of this paper, the asymptotic freeness of all the (permuted) balanced marginals of a random multipartite quantum state. Theorem 4.13. In the asymptotic setting described above, the family of random matrices (W f,π ) f ∈( 1,n n/2 ), π∈S n/2 converges, as N → ∞, to a family of n n/2 (n/2)! free random variables, each having a Marčenko-Pastur distribution of parameter 1. At the level of moments, this reads lim

N →∞ EN -n(p+1)/2 Tr W f ,π = M∈M bal f ,π 1, ( 113
)
where f is a word of p (n/2)-subsets of 1, n , π is a word of p permutations and M bal f ,π := {M ∈ M p : M is planar and f , π are constant on the white vertices of M}.

Remark 4.14. The above result can be easily generalized to the case where the dimensions of the Hilbert spaces are not identical: one has to replace the "1" in the right hand side of (113) with the appropriate product of constants c i , where the dimensions of the Hilbert spaces scale as N i ∼ c i N , with N → ∞.

Proof. The statement about freeness follows easily from the moment formula, using the correspondence between planar maps and non-crossing partitions and the fact that a map M belongs to M bal f ,π iff its corresponding non-crossing partition α satisfies α ≤ (ker f ∧ ker π) (see also the end of the proof of Theorem 3.6).

We now show the moment formula (113). Note that we are actually in the setting of Theorem 4.3 and of Proposition 4.6, with the permutations σ being given by σ a = π -1 a+1 π a , ∀a ∈ 1, p . Applying now Proposition 4.6 (remember that all Hilbert space dimensions are N in our current setting), we get ⊆ M bal f ,π . First, from ∆ f ,σ (M) = 0, we get that f is constant on the white vertices of M, see (92). The fact that π is also constant on the white vertices of M follows form Lemma 4.11, while the planarity of M follows from Proposition 4.7.

For the other inclusion, consider a map M ∈ M bal f ,π . Since f is constant on the white vertices of M, we have ∆ f ,σ (M) = 0. Also, recall from the proof of Lemma 4.11 that, for any edge a and color i ∈ J a , the black vertex to which (a, i) is connected in M f ,σ has the following incident edges: {(c, π -1 c π a (i))} c∈ 1,p . In particular, using the fact that π is constant on the white vertices of M, all edges (b, i) incident to the same white vertex as (a, i) in M f ,σ are connected to the same black vertex, and thus M f ,σ has as many connected components as black vertices (i.e. k), so Σ(M f ,σ ) = 0. Finally, from Lemma 4.12, the map M f ,σ consists of k copies of M, so that if M is planar, M f ,σ is too. This concludes the proof.

Remark 4.15. The result above has a simple interpretation in the case of two marginals W f,id and W f ,id , for complementary sets f, f : f f = 1, n . Indeed, for all N , these two marginals have exactly the same eigenvalues, and their eigenvector unitary operators are independent (this follows from the fact that they act on non-overlapping tensor factors). Free independence is in this case a consequence of Voiculescu's classical result [START_REF] Voiculescu | Circular and semicircular systems and free product[END_REF]. The discussion above generalizes easily to an arbitrary number of non-overlapping marginals. ≤ n + (k + l)(p -2) + l + V • (f , σ) -µ.

The fact that µ ≥ 0 follows from the hypothesis n -2l -2k < 0, and from the bounds V (M) ≥ 2, g(M) ≥ 0, L(f , σ, M) ≥ 0.

If all the edges have the same set of colors, the minimum µ is achieved by the unique planar map in M p having two vertices (one black and one white). However we stress that in the general case, this map does not have a vanishing L(f , σ, M) ≥ 0, and therefore the minimization problem seems highly non-trivial.

Indeed, depending on the values taken by the coefficient of V (M) -2, two maps with different g, L, and V might have the same minimal value of 2lg + L + (2k + 2l -n)(V -2). We illustrate this with a very simple example in the case l = 1 (color A), k = 2 (colors B, C), and n is kept as a parameter. We take σ a = π -1 a+1 π a , π 1 = (12) (the transposition) and π 2 = (1)(2) (the identity), so that both σ 1 and σ 2 are transpositions. In that case,

g = 0 L = 0 V -2 = 1 g = 0 L = 1 V -2 = 0
• if n = 3, 2k + 2l -n = 3, and the only map in M > f ,σ is that with two vertices, • if n = 5, 6, 2k + 2l -n = 1, 0 and the only map in M > f ,σ is the tree, • if n = 4, 2k + 2l -n = 2 and the two maps belong to M > f ,σ .

4.3. The unbalanced asymptotical regime. We consider in this section the asymptotical regime where the dimensions of the Hilbert spaces that "move around" are fixed; in this respect, the setting here generalized the 4-partite situation considered in Section 3.4. In this case, we fix N J = m, and set ∀i ∈ 1, l , N i = N and ∀j ∈ r, n , N j ∼ c 1 n-r+1 N as N → ∞. Theorem 4.3 writes

E Tr σ W f = N pl+n-r+1
M∈Mp h(N ) V (M)-1 N -2lg(M)-(l+r-n-1)(V (M)-2) (115)

× m kp+V•(f ,σ)+(r-l-2k-1)(V (M)-1)-L(f ,σ,M) , where h is an arbitrary function whose limit at infinity is c. Again we have different behaviors depending on the factor l + r -n -1 = 2l + k -n in front of V (M) -2:

• If l + r -n -1 < 0, the large N limit is trivial, the only surviving map being the one with p + 1 vertices (1 black and p white). This situation was detailed in the beginning of Sec. 4.2. • If l + r -n -1 > 0, then l > 0, and the large N limit selects the only planar map with two vertices. The limit is also trivial.

Besides the two trivial cases above, the only remaining situation is l + r -n -1 = 0 (or, equivalently 2l + k = n) and we therefore assume it in the following. In this case, Theorem 4.3 amounts to "flipping" the two edges, as shown in the figure below (the half-edges incident to the white vertices remain the same, but the endpoints incident to the black vertices are exchanged 15 ). The flip clearly decreases the genus or the number of connected components. To see this, note that we can decompose the flip in edge-deletions and creations. The edges have (the same) two faces incident each. Delete one of them, and the other gets a single face visiting both sides. Delete the second edge: if it is a bridge the number of connected components increases (and Σ f ,σ goes from 1 to 0) and the genus is unchanged, so L decreases by two. If not, Σ f ,σ is unchanged, and the genus decreases by one, so L decreases by two, and we get the formula, after creating two new edges to complete the flip (creating the two new edges does not change the number of connected components nor the genus). We now assume that ∆ f ,σ = 1. The strategy is to pick some a such that σ a+1 contains B in its image (we can always find one, since ∆ f ,σ = 0), and replace σ a and σ a+1 so that B is a fixed point of σ a+1 . Then do the same for σ a-1 and σ a , and so on. At each step, one should be careful that the variation in the number of inversions is the same as L. In the end, B is a common fixed point of all the σ a , so we are in the ∆ f ,σ = 0 case treated above, and we can conclude. There are two cases, depending whether B is in the pre-image of σ a+1 or not:

(1) If B is in the pre-image of σ a+1 (and is not a fixed point of σ a+1 , in which case we go to the next step right away), we exchange it with the other color in the pre-image of σ a+1 , both in σ a and σ a+1 :

σ a σ a+1 c a → c a+1 → B c a → B → c a+2 c a → B → B c a → c a+1 → c a+2 . ( 121
)
There are two sub-cases: either the number of inversions decreases by two, in which case we verify that L decreases by two (6 possibilities), or the number of inversions is constant, in which case we verify that L is also constant (4 possibilities).

(2) If B is not in the pre-image of σ a+1 , we just exchange the pre-image of B with B, and leave the other color in the pre-image of σ a+1 untouched. Again, either the number of inversions decreases by two, in which case we verify that L decreases by two (12 possibilities), or the number of inversions is constant, in which case we verify that L is also constant (8 possibilities). This concludes the proof When trying to generalize the formula of the result above to V (M) > 2, we face the following difficulty: around a white vertex, there is at least one corner for which we do not have a σ (the white corner which does not face a black corner). On the example of Fig. 25, it is easily seen that the value of L for the full map is not given as the sum of the value L for the two submaps with one white vertex. This indicates that there is no factorization of L over the cycles of α. It is also 15 Indeed changing σa and σa+1 changes the cycles of Γ f ,σ while keeping α invariant 12), one has L(f , σ, M) = 2, while the value of L for both submaps with one white vertex is 2, so the sum over the one white vertex submaps of L does not match L(f , σ, M). However, if π 1 = (1)(2), π 2 = (1)(2), π 3 = (1)(2), π 4 = (1)(2), then L(f , σ, M) = 0 and the one white vertex submaps also have L = 0. This behavior indicates that there is no factorization over the cycles of α.

possible to display a formula for the case of maps whose edges are only of color CD. This formula shows that there is no factorization over the cycles of α. Thus it is not possible to write the free cumulants in the context of (scalar) free probability. It seems that the right framework needed to tackle this unbalanced scaling is the one of free probability with amalgamation over an algebra B that is contained in 1 A ⊗ M m (C) ⊗ M m (C) ⊃ B. In order to limit the length of this paper, we postpone the exploration and the exposition of such results to future work.

Figure 1 .

 1 Figure 1. The density of the Marčenko-Pastur distribution MP c for c = 1 (left) and c = 5 (right).

Figure 4 .

 4 Figure 4. A trace of p matrices can be represented as a vertex of valency p with a cyclic ordering of incident half-edges.

Figure 5 .

 5 Figure 5. Combinatorial maps for the example of Figure 3.

Figure 6 .

 6 Figure 6. Bijection between elements of M 0 p and labeled plane trees.

Figure 10 .

 10 Figure 10. Local duplication of white vertices to obtain the map M f .

  where γ B (resp. id C ) denotes the permutation γ (resp. id) acting on [p] B = {(1, B), . . . , (p, B)} (resp. on [p] C = {(1, C), . . . , (p, C)}), and δ C is a product of transpositions, S 2p δ C := a∈f -1 ({AC}) (a, B)(a, C) .

Figure 11 .

 11 Figure 11. Diagram showing that whenever f (a) = f (α -1 (a)), (a, B) and (a, C) belong to the same cycle of ε = δ C αδ C (γ B ⊕ id C ). Since the partition α Kr is noncrossing (blue partition), the element (a, C) (in red) must "escape" the interval [a, α(a)] through (a, B) (in blue).

Figure 13 .

 13 Figure 13. The density of the free multiplicative square of the Marčenko-Pastur distribution MP c for c = 1 (left) and c = 2 (right) versus Monte-Carlo simulations. The supports of the two probability measures are [0, 27/4] (left) and [(71 -17 √ 17)/8, (71 + 17 √ 17)/8] ≈ [0.11, 17.64] (right).

  Cat σ (r)Cat π Kr (s) Mob(σ, π), (58) where Cat σ is the multiplicative extension of Catalan numbers Cat σ (r) := b cycle of σ Cat i∈b r i , Kr denotes the Kreweras complementation (see [Kre72] or [NS06, Lecture 9]) and Mob is the Möbius function on the non-crossing partition lattice (see [NS06, Lecture 10]), defined for σ ≤ π by Mob(σ, π) := Mob(σ -1 π) = b cycle of σ -1 π (-1) |b|-1 Cat |b|-1 .

Figure 14 .

 14 Figure 14. Bijection between tricolored maps with one black vertex, and tricolored trees in T f .

Figure 18 .

 18 Figure 18. Above: Tutte duality. The initial map M on the left, its Tutte dual T (M) on the right. The middle map represents the intermediate steps of the construction. The labeling of the edges of T (M) is inherited from the map M. Below: Kreweras complementation map.

.

  and this equality does not hold in the planar case, In the limit N → ∞, thanks to Proposition 3.15, the equation (70) rewrites lim

Figure 19 .

 19 Figure 19. Plots of Monte-Carlo simulations of the eigenvalues of P = W 1/2 AB W AC W 1/2 AB (yellow histogram) versus the square of a Marčenko-Pastur distribution (red curve) and the free multiplicative square of the same Marčenko-Pastur distribution (blue curve). On the top row, we have N = 600, c = 1 and m = 1, 2, 5, while on the bottom row we have N = 600, c = 5, m = 1, 2, 5.

Figure 21 .

 21 Figure 21. Local duplication of white vertices.

Figure 22 .Figure 23 .

 2223 Figure22. The map M f ,σ for an example in the case (85). The white vertex in the original map M is tripled into 3 white vertices having half-edges of the same color attached. The order of the half-edges around the black vertex is given by the permutation Γ f ,σ .

Figure 24 .

 24 Figure 24. Diagrams for the permuted marginals W {1,2},(1)(2) and W {1,3},(12) .

  lim N →∞ EN -n(p+1)/2 Tr W f ,π = M∈M (0)free f ,σ 1,where we recall that (see (104)) M (0)free f ,σ = {M : L(f , σ, M) = 0}, the functional L being defined in (95). From Lemma 4.10, we get that for all maps M, V• (f , σ) = k. It is now enough to show that M (0)free f ,σ = M bal f ,π .Let us start with the inclusion M (0)free f ,σ

4.2. 2 .

 2 Case where less than half the colors are traced. Proposition 4.16. In the large N regime, if l + k > n 2 and with the previous notations,lim N →∞ 1 N n+(k+l)(p-2)+l+V•(f ,σ)-µ E Tr σ W f = M∈M > f ,σ c V (M)-1 , (114)where µ = min{2lg(M) + L(f , σ, M) + (2k + 2l -n)(V (M) -2) : M ∈ M p } ≥ 0 and M > f ,σ := {M ∈ M p : 2lg(M) + L(f , σ, M) + (2k + 2l -n)(V (M) -2) = µ}.Proof. It follows from (97), but let us detail the computations starting from Theorem 4.3. Under the usual balanced scaling, the exponent of N from Theorem 4.3 can be bounded as followsl(2 + p-V (M) -2g(M)) + (n -r + 1)(V (M) -1) + kp + V • (f , σ)+ + (r -l -2k -1)(V (M) -1) -L(f , σ, M) = (n -2l -2k)(V (M) -1) + (k + l)p + l -2lg(M) + V • (f , σ) -L(f , σ, M) = n + (k + l)(p -2) + l + V • (f , σ)--[2lg(M) + L(f , σ, M) + (2k + 2l -n)(V (M) -2)]

Figure 25 .

 25 Figure 25. For π 1 = (1)(2), π 2 = (12), π 3 = (1)(2), π 4 = (12), one has L(f , σ, M) = 2, while the value of L for both submaps with one white vertex is

  The diagram of the random variable Tr(W AB W AC ). Each tensor box has 4 decorations corresponding to H A (round shape), H B , H C (square shapes), and H D (diamond shape), displayed in this order from top to bottom.

	X	X *	X	X *
	Figure 9.

Note that by linking the newly added vertices to the corners around the black vertex instead of the white vertices, we have a duality between elements of M 0 p , see Def.

3.14.

See later sections.

This can actually be shown, however this lies out of the main scope of this paper.

The letter a stands for arêtes, which means edges in French.

In the following we omit color A from the support and image of the σa's as it is a common fixed point.

writes E Tr σ W f = N l(p+1) M∈Mp h(N ) V (M)-1 N -2lg(M) m kp+V•(f ,σ)+(r-l-2k-1)(V (M)-1)-L(f ,σ,M) , (116) where (r -l -2k -1)(V (M) -1) -L(f , σ, M) is a generalization of "alt" which we defined in Def. 51, Section 3 for n = 4. We recall that M (0) p is the subset of planar elements of M p . Proposition 4.17. In this regime, with the previous notations,

Note that if for a ∈ 1, p we assume V • (f , σ) = k, and we define Wf(a) = 1 N l m k W f (a) , as well as c = cm r-l-2k-1 , this rewrites as lim

An example for n = 5. We now focus on the quantity L(f , σ, M), in the following interesting case of "ABCDE": n = 5, l = 1 (color A), r = 5 (color E is always traced), and k = 2 (the edges carry two colors among {B, C, D}); the colors A and E correspond to spaces of dimension growing to infinity, while the colors B, C, D correspond to C m . As l + r -n -1 = 0, Prop. 4.17 and Eq. ( 118) apply. The study of the properties of L in this case shows that the corresponding cumulant functions do not factorize over the cycles of non-crossing partitions, but depend more subtly on these non-crossing partitions. In what follows, we study the special case of a map with two vertices (one black, one white), and we show, by the means of an example, that in the general case the factorization property does not hold.

To start, we recall the expression of L(f , σ, M)

Any permutation σ ∈ S 2 14 defines a partial permutation on l +1, r -1 = 2, 4 , which we complete canonically to have a permutation Φ(σ) of 2, 4 . Indeed, there is only one missing color in the support and image of each permutation σ a , so we just pair them (e.g. if

We recall that |σ| is the length of σ, i.e. its number of inversions.

Proposition 4.18. In the n = 5 example described above, if

Proof. To begin, note that the white vertex of M can be duplicated into 2 or 3 white vertices in the unfolded map M f ,σ , depending on whether the edges share the same set of colors or not; hence, ∆ f ,σ is either 0 or 1. We first assume that ∆ f ,σ = 0. This means that e.g. B is not in any of the color sets J a , the σ a are either the identity, either the transposition, on {C, D}, and that Φ(σ a )(B) = B for all a. We can assume that there are only transpositions in the corners, because when one has the identity permutation, one can collapse the two neighboring edges into a single one without changing the genus of M or M f ,σ , nor Σ(M f ,σ ). Here one can make a simple recurrence, which is quite similar to the one appearing in the proof of Proposition 3.5 in terms of maps: pick any a, and replace σ a and σ a+1 by two identities. The number of inversions decreases by two. In M f ,σ , this