Exact eigenvalue order statistics for the reduced density matrix of a bipartite system

We consider the reduced density matrix $\rho_{A}^{(m)}$ of a bipartite system $AB$ of dimensionality $mn$ in a Gaussian ensemble of random, complex pure states of the composite system. For a given dimensionality $m$ of the subsystem $A$, the eigenvalues $\lambda_{1}^{(m)},\ldots, \lambda_{m}^{(m)}$ of $\rho_{A}^{(m)}$ are correlated random variables because their sum equals unity. The following quantities are known, among others: The joint probability density function (PDF) of the eigenvalues $\lambda_{1}^{(m)},\ldots, \lambda_{m}^{(m)}$ of $\rho_{A}^{(m)}$, the PDFs of the smallest eigenvalue $\lambda_{\rm min}^{(m)}$ and the largest eigenvalue $\lambda_{\rm max}^{(m)}$, and the family of average values $\langle \mathrm{Tr}\big(\rho_{A}^{(m)}\big)^{q}\rangle$ parametrised by $q$. Using values of $m$ running from $2$ to $6$ for definiteness, we show that these inputs suffice to identify and characterise the eigenvalue order statistics, i.e., to obtain explicit analytic expressions for the PDFs of each of the $m$ eigenvalues arranged in ascending order from the smallest to the largest one. When $m = n$ (respectively, $m<n$) these PDFs are polynomials of order $m^{2}-2$ (respectively, $mn - 2$) with support in specific sub-intervals of the unit interval, demarcated by appropriate unit step functions. Our exact results are fully corroborated by numerically generated histograms of the ordered set of eigenvalues corresponding to ensembles of over $10^{5}$ random complex pure states of the bipartite system. Finally, we present the general solution for arbitrary values of the subsystem dimensions $m$ and $n$, namely, formal exact expressions for the PDFs of every ordered eigenvalue.


December 2021
Abstract. We consider the reduced density matrix ρ we show that these inputs suffice to identify and characterise the eigenvalue order statistics, i.e., to obtain explicit analytic expressions for the PDFs of each of the m eigenvalues arranged in ascending order from the smallest to the largest one. When m = n (respectively, m < n) these PDFs are polynomials of order m 2 − 2 (respectively, mn − 2) with support in specific sub-intervals of the unit interval, demarcated by appropriate unit step functions. Our exact results are fully corroborated by numerically generated histograms of the ordered set of eigenvalues corresponding to ensembles of over 10 5 random complex pure states of the bipartite system. Finally, we present the general solution for arbitrary values of the subsystem dimensions m and n, namely, formal exact expressions for the PDFs of every ordered eigenvalue.

Introduction
A basic motif in the study of entanglement involves a bipartite system AB comprising subsystems A and B with Hilbert space dimensions m and n( m), respectively. Given a specified Gaussian ensemble of random pure states of the composite system, the task is to deduce the statistical properties of the reduced density matrix ρ (m) a . A fairly extensive literature exists in this regard (see Refs. [1]- [19] and references therein). The properties investigated include the joint probability density function (PDF) of the set of eigenvalues {λ [1,2], the leading large-m behaviour of the average values of the set of eigenvalues [9], the PDFs of the smallest and largest eigenvalues λ (m) min [3,9,11,15] and λ (m) max [17], their mean values and higher moments (λ a ) q where q is any positive integer [11,17,19], and the associated entropies that quantify the extent of entanglement of A and B such as the average subsystem von Neumann entropy − Tr(ρ as already stated. The natural extension of extreme value statistics is the order statistics of the eigenvalues, the task being to deduce the PDFs of the individual eigenvalues identified by their positions in an ordered sequence. This is the purpose of the present paper. The complications involved in the statistics of extreme values in the case of correlated random variables persist, of course, for order statistics as well. In order to avoid any confusion, we shall denote the eigenvalues λ The investigations cited in the opening paragraph above rely on the following basic result [2]. Let λ  (listed in no particular order). Their joint PDF is then given by Here the Dyson index β = 2 (respectively, 1) for a Gaussian ensemble of complex (respectively, real) pure states, and the normalizaton constant is As is well known in random matrix theory, the eigenvalues {λ (m) k } form a set of correlated random variables both because of the bunching effect arising from the requirement that their sum be equal to unity, as well as the level repulsion implied by the presence of the factor |λ k | β in the joint PDF of Eq. (2). In principle, the PDF p (m) k (x) of the k th eigenvalue in the ordered eigenvalue sequence can be found by multiplying the joint PDF in Eq. (2) by the product of step functions integrating over λ 1 , . . . , λ k−1 , λ k+1 , . . . , λ m , and normalizing the resulting function of x to unity. This approach, however, presents formidable technical problems, and is not feasible. We shall see that there is an alternative, simpler procedure to arrive at the result sought.
It is evident from Eqs. (2) and (4) that some simplification occurs when β = 2 (complex pure states), and we shall consider this case. While we shall finally present exact expressions for the PDFs of the ordered eigenvalues for arbitrary subsystem dimensions m and n, those expressions are somewhat complicated. The manner in which the structure of the solutions arises is best elucidated by explicit illustration for small values of m. Accordingly, we shall start with m = 2 and increase it step by step up to m = 6, to demonstrate how the PDFs build up. In the interests of clarity, we shall further set n = m for the most part in this demonstration, in order to take advantage of the simplification that ensues from the fact that α = 0 in this case.
The plan of the paper is as follows. In Section 2, we write down the ranges of the random variables {Λ (m) k }, followed by some properties of Mellin transforms that will be used in the sequel. We also list three specific, already known results that will be used to deduce the PDFs of the ordered eigenvalues. Next, in Section 3 we state (for the sake of completeness) the existing results in the trivial case m = n = 2. We then discuss in Section 4 the case m = n = 3, which is again special because there is just one eigenvalue in between the smallest and largest eigenvalues. In Sections 5 and 6, we illustrate our procedure to obtain the analytical expressions for PDFs corresponding to m = 4, 5, and 6 when m = n. We also show that our procedure works for the case m = n by considering the case m = 4, n = 5. Finally, in Section 7 we present the solution for the PDFs {p

Preliminaries
Integrating P (λ k } is a technically complicated task. It yields a so-called 'single-particle' PDF for a single eigenvalue [13]. But this procedure automatically averages over the location of the eigenvalue in the ordered set of eigenvalues. Thus, it leads, for instance, to the expected result that the average value λ (m) k is just 1/m. It is evident that this single-particle PDF is quite different from the individual PDFs corresponding to the ordered eigenvalues.
We start with some general properties of the ordered set of eigenvalues {Λ In particular, Λ (m) m−1 1/2. Moreover, it is evident that the largest eigenvalue cannot be smaller than 1/m, so that its range is given by The ranges in Eqs. (6) and (7) therefore specify the support I As we shall also be concerned with the higher moments of the eigenvalues Λ (m) k , it is helpful to recall very briefly some properties of Mellin transforms. The Mellin transform f (q) of a function f (x) (x ∈ [0, 1]) and its inverse are defined as where the Bromwich contour C runs from c − i∞ to c + i∞ to the right of all the singularities of f (q). The following easily established results will be used in the sequel: (a) If f (q) is a rational function of q with simple poles at the negative integers (b) Let r be a positive integer > 1. If f (q) is r −q times a rational function of q with simple poles at the negative integers q = −1, . . . , −ν, then f (x) is a polynomial in x of order ν − 1, multiplied by the unit step function Θ(1 − rx).
In the present context, we note that the PDF p (m) k (x) and the q th moment of Λ comprise a Mellin transform pair.
We turn now to the three known results that we require to deduce the PDFs of the ordered eigenvalues, setting β = 2 and m = n.
is given by [11] p (m) where Θ denotes the unit step function. The corresponding moments of the smallest eigenvalue are then given by a ) q is given by [19] Tr(ρ We note that this average is the Mellin transform of the sum of PDFs, m k=1 p (m) k (x). Inverting the transform will therefore yield that sum.
(iii) The third and most crucial ingredient is the PDF p m , for which an implicit formula has been derived in Ref. [17]. As pointed out therein, the determination of p  1 (x). The procedure [17] leading to the result sought may be summarized in the following sequence of steps. First, one defines the set of m 2 functions Ψ jl (s) = 1 0 e −su u j+l du (j, l = 0, 1, . . . , m − 1) (13) and evaluates the (m × m) determinant The inverse Laplace transform P(t) of P(s) is then obtained, and t is set equal to 1/x in the result. Then, the quantity is the cumulative distribution function of Λ (m) m , from which its PDF follows according to The Mellin transform of p We note, for future reference, that the counterparts of Eq. (12) and Eqs. (13)- (15) to general values of m and n ( m) are also available (Refs. [19] and [17], respectively). We proceed now to find the full set of PDFs {p 3. Two qubits: m = n = 2 A bipartite system of two qubits is a trivial case as far as order statistics are concerned, since there are only two eigenvalues Λ (2) 1 and Λ 1 . Setting m = 2 in Eq. (10), the PDF of the smaller eigenvalue is given, in this case, by Working out the steps outlined in Eqs. (13) to (16) with m set equal to 2, we obtain Apart from the step functions, these are polynomials of order m 2 − 2 = 2, with support in [0, 1/2] and [1/2, 1] respectively. The expression for p 2 (x) could have been written down from that for p (2) 1 (x) in this case: Since Λ In order to verify the expressions above and those to be obtained for p computed from an ensemble of random complex pure states of the composite system. We consider a randomly chosen pure state of the full system AB to be an mn-dimensional column vector, with the real and imaginary parts of each element of the vector drawn from a standard normal distribution. The state is then normalized. The moments (and hence the cumulants) of the numerically generated histograms, obtained with 1.001 × 10 5 random pure states, agree up to the third decimal place with those computed from the analytical expressions. We have also verified that this agreement improves on increasing the number of random pure states in the ensemble to 10 6 . These statements remain valid in all the cases to be considered in the sections that follow. Figure 1 shows that there is excellent agreement between the numerically generated histograms of Λ and (Λ 2 ) q = 3!q!(q 2 + q + 2 − 2 −q )/(q + 3)!
tallies with the corresponding expression for Tr (ρ 3 . There is, however, a simple strategy to find the PDF of Λ 2 in this instance. In brief, we first find (Λ 3 ) q , and use these results along with that for Tr (ρ 2 ) q . The inverse Mellin transform of the latter then yields the PDF p When m = n = 3, Eq. (10) gives for the PDF of the smallest eigenvalue the expression Exact eigenvalue order statistics for the reduced density matrix of a bipartite system 8 As before, working through the steps in Eqs. (13) to (16) with m = 3, we obtain The Mellin transforms of the two PDFs in Eqs. (22) and (23) yield, respectively, (Λ and (Λ 3 ) q = 8!q! 2 6 (q + 8)! 2 4 (q 4 + 2q 3 + 11q 2 + 10q + 12) On the other hand, setting m = 3 in Eq. (12) gives An important point that we note here for future reference is the following. After the ratio q!/(q + 8)! is simplified, the expression on the right-hand side of Eq. (26) is a rational function of q. There are no transcendental functions like r −q present in Tr(ρ Inverting the Mellin transform, we obtain for the PDF of the middle eigenvalue the explicit expression Figure 2 again shows that the three PDFs p  As the PDFs are essentially polynomials with compact support in ranges whose endpoints are rational numbers, all the moments of these PDFs (and hence their cumulants) are rational numbers. (This feature remains valid for all values of m and n.) Table 1 lists the values of the basic descriptors of the distributions concerned in terms of the corresponding cumulants κ r : the mean κ 1 , the variance κ 2 , the skewness κ 2 3 /κ 3 2 , and the excess of kurtosis κ 4 /κ 2 2 . k (x) for every one of these eigenvalues. The case m = 4 serves as the simplest illustration of this method. As before, we start with the PDF of the smallest eigenvalue, obtained by setting m = n = 4 in Eq. (10). We have Next, we find the explicit expression for the PDF p where We observe that ensuring that p (4) 4 (x) vanishes identically for x < 1/4, (As we know, its support is [1/4, 1]). Next, setting m = 4 in Eq. (12), we get Tr(ρ (4) a ) q = 15!q! 36(q + 15)! (144 + 156q + 184q 2 + 57q 3 + 31q 4 + 3q 5 + q 6 ).(33) Once again, we note that the expression on the right-hand side of Eq. (33) is a rational function of q (after the ratio q!/(q + 15)! is simplified). Hence, by the property (a) of Mellin transforms noted in Section 2, the step functions Θ(1 − 4x), Θ(1 − 3x) and Θ(1 − 2x) cannot appear in its inverse Mellin transform 4 k=1 p (4) k (x). The coefficients of these step functions must therefore vanish identically when the individual PDFs are added up. Note also that p where c 1 and c 2 are constants. They are determined from the normalization (to unity) of p 2 (x) and p We observe from the foregoing (and from all the cases to be considered in the sequel) that the constants multiplying the coefficients A (m,n) j (x) for a given j in the different PDFs p The higher cumulants can also be calculated, and they are all rational numbers. We have also verified that the sum of the q th moments of the eigenvalues tallies with the known expression for Tr(ρ

Other cases
As further checks of the method used, we have carried out similar calculations to determine the PDFs of the ordered eigenvalues in the cases m = n = 5, 6 and 7, respectively. The algebraic expressions become considerably more lengthy as m increases. The expressions for the PDFs when m = n = 5 are given in Appendix A, and these expressions agree very well with the numerically generated histograms, as shown in Figure 4. As already pointed out, we find that the constants multiplying the coefficient functions A  The expressions obtained for the PDFs in the case m = n = 6 are also recorded in Appendix A. Once again, we have also verified that there is very good agreement between the analytical expressions for the PDFs and the numerically generated histograms. Similarly, the expressions for m = n = 7 are also precisely along expected lines, and will not be given here.
Finally, in order to show that our method works even when m = n, we have found the analytical expressions for the PDFs when m = 4 and n = 5. We must now take into account the fact that the index α = 1 in this case, and use the corresponding generalizations of Eqs. (12)- (16). The details are given in Appendix B. Once again, the plots of the calculated PDFs are in complete agreement with the numerical histograms, as shown in Figure 5. Table 2 lists the averages Λ

Solution for general m and n
We now proceed to the exact formal expression for the PDF p The procedure followed is the same as that for the case m = n. As already mentioned, the counterparts of Eqs. (12) [19] and Eqs. (13)- (15) [17] for the case n m are now required. The pattern in the structure of the PDFs found in the foregoing sections aids us considerably in deducing the structure for general m and n. We obtain, finally, where A  Let σ = (σ 0 , σ 1 , . . . , σ m−1 ) be a permutation of the sequence {0, 1, . . . , m − 1), with sgn σ = ±1 depending on whether σ is an even or odd permutation of the natural order, and let S denote the set of all permutations σ. Setting a = m − j where 1 j m, we have where and It is evident that the general solution for the PDF p (m,n) k (x), while exact and explicit, is algebraically quite involved. This fact further corroborates the usefulness of displaying in detail the results for several small values of m, as has been done in the foregoing.
To summarize: We have obtained the probability density functions of the eigenvalue order statistics Λ In all the cases considered, the analytic expressions obtained for the PDFs are in excellent agreement with the numerically generated histograms of the eigenvalues concerned. As further corroboration, we also find that, in every case, the Mellin transform of the sum of the q th moments of these PDFs matches the known expression [19] for Tr(ρ Based on the explicit analytic solutions in the cases m = 3, 4, 5, 6 and 7, we deduce the following general properties. When m = n, the PDF p    In terms of these polynomials, the PDFs {p