---
abstract: |
    Given a random quantum state of multiple (distinguishable or indistinguishable) particles, we provide an algorithm, rooted in symplectic geometry, to compute the joint probability distribution of the eigenvalues of its one-body reduced density matrices, and hence some associated physical invariants of the state.

    As a corollary, by taking the support of this probability distribution, which is a convex polytope, we recover a complete solution to the one-body quantum marginal problem, i.e., the problem of characterizing the one-body reduced density matrices that arise from some multi-particle quantum state. In the fermionic instance of the problem, which is known as the one-body $N$-representability problem, the famous Pauli principle amounts to one linear inequality in the description of the convex polytope.

    We obtain the probability distribution by reducing to computing the corresponding distribution of diagonal entries (i.e., to the quantitative version of a classical marginal problem), which is then determined algorithmically. This reduction applies more generally to symplectic geometry, relating invariant measures for a compact Lie group action to that for the maximal torus action; we state and prove our results in this more general symplectic setting. Our approach is in striking contrast to the existing solution to the computation of the supporting polytope by Klyachko and by Berenstein and Sjamaar, which made crucial use of non-Abelian features.

    In algebraic geometry, Duistermaat–Heckman measures correspond to the asymptotic distribution of multiplicities of irreducible representations in the associated coordinate ring. In the case of the one-body quantum marginal problem, these multiplicities include bounded height Kronecker and plethysm coefficients. A quantized version of the Abelianization procedure provides an efficient algorithm for their computation.
address:
- 'Matthias Christandl, Institute for Theoretical Physics, ETH Zürich, Wolfgang–Pauli–Strasse 27, CH-8093 Zürich, Switzerland'
- 'Brent Doran, Department of Mathematics, ETH Zürich, Rämistrasse 101, CH-8092 Zürich, Switzerland'
- 'Stavros Kousidis, Institute for Theoretical Physics, ETH Zürich, Wolfgang–Pauli–Strasse 27, CH-8093 Zürich, SwitzerlandInstitute of Physics, University of Freiburg, Rheinstrasse 10, 79104 Freiburg, Germany'
- 'Michael Walter, Institute for Theoretical Physics, ETH Zürich, Wolfgang–Pauli–Strasse 27, CH-8093 Zürich, Switzerland'
author:
- Matthias Christandl
- Brent Doran
- Stavros Kousidis
- Michael Walter
bibliography:
- 'dhmeasure.bib'
title: Eigenvalue Distributions of Reduced Density Matrices
---

Introduction
============

The pure state of a quantum system is described by a vector in a complex Hilbert space, or more precisely by a point in the corresponding projective space. Herein we consider the finite dimensional case, for instance a spin system. Since the Hilbert space for multiple particles is given by the tensor product of the Hilbert spaces of the individual particles, its dimension grows exponentially with the number of particles. This exponential behavior is therefore the key obstruction to classical modeling of quantum systems. The observation is as old as quantum theory itself, and physicists ever since have tried to find ways around it. In that spirit, our paper presents an effective method to extract those physical features that “only depend on the one-body eigenvalues” associated to randomly-chosen quantum states of any fixed number of particles. Our methods are rooted in geometry, and so further a nascent dialog between algebraic and symplectic geometry on the one hand, and the theory of quantum computation and quantum information on the other, with new results for each subject. We therefore begin with some context, and make an effort in the body of the paper to build a correspondence in terminologies.

Typically, to address the aforementioned exponential complexity, physicists make use of a very simple yet powerful observation: Important properties such as energy and entropy often do not depend on the whole wavefunction but rather on only a small part, namely the *reduced density matrix*, or quantum marginal, of a few particles. For instance, the binding energy of a molecule is given as a minimization over two-electron reduced density matrices arising from $N$-electron wavefunctions. Mathematically, the reduced density matrix is given as the contraction of (or trace over) the indices of the projection operator onto the wavefunction over the remaining particles (see for the precise definition). The problem of characterizing the set of possible reduced density matrices, known as the *quantum marginal problem* in quantum information theory and as the $N$-representability problem in quantum chemistry [@ruskai69; @colemanyukalov00], has therefore been considered one of the most fundamental problems in quantum theory [@stillinger95]. The general problem is computationally intractable, even on a quantum computer; more precisely it is QMA-complete and NP-hard [@liu06; @liuchristandlverstraete07]. However, the characterization of the one-body reduced density matrices of a pure global quantum state [@klyachko04; @daftuarhayden04; @klyachko06] admits a very elegant mathematical interpretation: it amounts to a description of the possible eigenvalues of the reduced density matrices; and this, just as for Horn’s problem, requires the computation of moment polytopes for coadjoint orbits of unitary groups, as first observed in [@christandlmitchison06; @daftuarhayden04; @klyachko04]. An interesting consequence is that defining linear inequalities for the polytope can have physical interpretation: for instance, the Pauli principle is simply one linear inequality bounding the polytope of a fermionic system (see [@coleman63; @borlanddennis72; @ruskai07; @klyachkoaltunbulak08; @klyachko09] and [@higuchi03; @higuchisudberyszulc03; @bravyi04] for further examples).

Random states play a fundamental role in physics. In classical *statistical mechanics*, the canonical state (or, “canonical ensemble”) is the marginal probability distribution (Gibbs measure, or Boltzmann distribution) of states of the system arising from a uniformly random configuration of system and bath, subject to an energy constraint and fixed particle numbers for both system and bath [@huang90]. In quantum statistical mechanics, the canonical state of the system is the reduced density matrix of the uniform state on a subspace (encoding the energy constraint) of system and bath. In fact it can be shown that the canonical state almost always approximates the reduced density matrix of a random pure state in this subspace, as for large systems a concentration of measure occurs [@popescushortwinter06; @lloyd06; @goldsteinlebowitztumulkaetal06]. Earlier, and also out of thermodynamic considerations, Lloyd and Pagels in a seminal work computed the distribution of eigenvalues of the reduced density matrix of the system when system and bath are in a random pure state [@lloydpagels88], i.e., in the simplest case of two large “particles”, one representing the system and one the bath. In this paper we are concerned with computing exact eigenvalue distributions for any number of particles, rather than with asymptotic results or with simple model systems of two particles. More precisely, we are concerned with the following question:

> *What is the probability distribution of the eigenvalues of the one-body reduced density matrices of a pure many-particle quantum state drawn at random from the unitarily invariant distribution?*

In this article, we answer this question completely by describing an explicit algorithm to compute the joint eigenvalue distribution of the reduced density matrices for an arbitrary number of particles of any statistics (distinguishable, Bose, Fermi). As a special case we easily recover the Lloyd–Pagels result for two distinguishable particles. As a corollary, our work naturally leads to a solution of the one-body quantum marginal problem in terms of a finite union of polyhedral chambers, thereby providing a complementary perspective on the work of Klyachko and of Berenstein–Sjamaar, and more recently of Ressayre, who each instead provided procedures based on geometric invariant theory (in particular the Hilbert–Mumford criterion) to list characterizing linear inequalities [@berensteinsjamaar00; @klyachko04; @ressayre10] for the moment polytope.

We emphasize that the computed eigenvalue distributions can be directly used to infer distributions of Rényi and von Neumann entropies, both of which play a fundamental role in statistical physics and quantum information theory and are functions of the eigenvalues only. In particular, one can recover the average entropy of a subsystem [@lubkin78; @page93], which featured in an analysis of the black hole entropy paradox [@page94; @HaydenPreskill]. For applications to the study of quantum entanglement, see .

From a mathematical perspective, the eigenvalue distributions that we compute are *Duistermaat–Heckman measures*, which are defined using the push-forward of the Liouville measure on a symplectic manifold along the moment map (see for the precise definition which we use in this article) [@heckman82; @guilleminsternberg82b; @guilleminsternberg84; @guilleminlermansternberg88; @guilleminlermansternberg96; @guilleminprato90]. The support of such a measure is a moment polytope, which in our physical context is the solution to the one-body quantum marginal problem. We will work and establish our results in this more general symplectic setting, in similar spirit to the existing solution to the one-body quantum marginal problem [@berensteinsjamaar00; @klyachko04]. We remark that, in particular, the one-body quantum marginal problem subsumes the well-known Horn’s problem, which asks for the possible eigenvalues of the sum $A + B$ of two Hermitian matrices $A$, $B$ with fixed eigenvalues [@helmkerosenthal95; @klyachko98; @knutsontao99; @fulton00; @knutsontao01; @knutsontaowoodward03]. The corresponding eigenvalue distribution for randomly-chosen matrices $A$ and $B$ has also been studied in the literature, and our methods allow us to recover the main results of [@dooleyrepkawildberger93] (); see also [@frumkingoldberger06] for a more concrete approach.

However, one striking difference between existing approaches to the one-body quantum marginal problem and the corollary to our approach is that the subtleties associated with the non-Abelian nature of their solutions can be completely bypassed. For instance, their reliance upon cohomology of Schubert cycles, and the interplay of different Weyl groups and sub-tori that feature because of repeated use of the Hilbert–Mumford criterion, may be seen as incidental. The problem, even in the full generality of the symplectic setting, is at heart an Abelian one whose essential combinatorics is encoded in the maximal torus action together with taking finitely many explicit derivatives. Remarkably, this is more than philosophy, as it has real import for computation.

The first mathematical contribution of this paper is an effective technique for providing explicit Duistermaat–Heckman measures under rather weaker assumptions than appear in the literature (). The fact that this subsumes our main question above, and in particular the one-body quantum marginal problem, while the existing literature does not, is the crucial second point of this work (). A third feature is our statement of the Abelianization procedure via the derivative principle for invariant measures (); though we have not seen this principle so formulated in the literature, there certainly are predecessors and the result should be equivalent to one of Harish-Chandra [@harishchandra57]. That a “quantized” version of our algorithm can be used to compute (efficiently, unlike existing algorithms, see ) multiplicities for the branching problem of representation theory may be seen as the fourth mathematical consequence of our approach.

The basic strategy we follow to address our main question above is a sequence of reductions. Firstly, our general quantum problem is replaced by an equivalent but more tractable one by “purification” of the quantum state (). This reduced problem is seen to satisfy a weak non-degeneracy assumption, so that the image of the moment map does not lie entirely on a wall in the relevant Weyl chamber. Under this assumption we can reduce via a *derivative principle* to the Duistermaat–Heckman measure for the maximal torus action (), which we evaluate by a “single-summand” algorithm along the lines of Boysal–Vergne [@boysalvergne09] ().

This derivative principle holds for more general $K$-invariant measures on the dual of the Lie algebra, $\mathfrak k^*$: Every invariant measure can be reconstructed from its projection onto the dual of the Lie algebra of the maximal torus by taking partial derivatives in the direction of negative roots (). Again, we remark that such a reduction is not possible on the level of the supports of the Duistermaat–Heckman measures, i.e., on the moment polytopes. We also remark that in the case of the one-body quantum marginal problem, the Duistermaat–Heckman measure for the maximal torus action is the joint distribution of the diagonal entries of the one-body reduced density matrices (as opposed to their eigenvalues). Amusingly, this distribution can be viewed as the solution to the quantitative version of a classical marginal problem ().

The single-summand algorithm alluded to above is an effective method for computing Duistermaat–Heckman measures for torus actions on projective spaces (). Its name stems from the fact that it amounts to evaluating a single summand of the kind that occurs in the well-known Heckman formula of Guillemin–Lerman–Sternberg [@guilleminlermansternberg88], which expresses the measure as an alternating sum of iterated convolutions of Heaviside measures. We also describe an algorithm based on this latter formula () that can in particular be applied to projections of coadjoint orbits (), and hence to the setting of Berenstein–Sjamaar [@berensteinsjamaar00].

Whenever a Hamiltonian group action can be quantized in a certain technical sense, Duistermaat–Heckman measures have an interpretation as the asymptotic limit of associated representation-theoretic quantities [@heckman82; @guilleminsternberg82b; @sjamaar95; @meinrenken96; @meinrenkensjamaar99; @vergne98]. In the second part of this article ( and ), we thus study the representation theory connected to the one-body quantum marginal problem; here, the relevant multiplicities include the Kronecker coefficients, which play a major role in the representation theory of the unitary and symmetric groups [@fulton97], as well as in Mulmuley and Sohoni’s geometric complexity theory approach to the P vs. NP problem in computer science [@mulmuleysohoni01; @mulmuleysohoni08; @mulmuley07; @burgisserlandsbergmaniveletal11]. It has been observed that the existence of a pure tripartite quantum state with given marginal eigenvalue spectra is equivalent to the asymptotic non-vanishing of an associated sequence of Kronecker coefficients [@christandlmitchison06; @klyachko04; @christandlharrowmitchison07], see also [@daftuarhayden04; @knutson09; @burgisserchristandlikenmeyer11; @burgisserchristandlikenmeyer11b]. For a similar connection in the context of Horn’s problem and Littlewood–Richardson coefficients, see [@lidskii82; @knutson00; @christandl08].

Our main results in this context are quantized versions of our earlier theorems: Kronecker coefficients can be computed by applying finite difference operators to weight multiplicities which are related to the classical marginal problem: Instead of measuring the volume of a polytope, one has to count the number of lattice points in the polytope. This can be computed efficiently using Barvinok’s algorithm [@barvinok94], and so leads to an efficient algorithm for computing Kronecker coefficients for Young diagrams with bounded height. Again, we shall establish the results in greater generality and recover the version for Kronecker coefficients as a special case.

Notation and Conventions {#notation}
------------------------

Throughout this article, $K$ will denote a compact, connected Lie group with maximal torus $T \subseteq K$, rank $r = \dim T$, Weyl group $W$, respective Lie algebras $\mathfrak k$ and $\mathfrak t$, and integral lattice $\Lambda = \ker {\left.\exp\vphantom{\big|}\right|_{\mathfrak t}}$ [@cartersegalmacdonald95; @kirillov08]. We write $\pi_{K,T} \colon \mathfrak k^* \rightarrow \mathfrak t^*$ for the projection dual to the inclusion $\mathfrak t \subseteq \mathfrak k$. We will think of weights as elements of the dual lattice, $\Lambda^* = \operatorname{Hom}_{{\mathbb Z}}(\Lambda, {\mathbb Z}) \subseteq \mathfrak t^*$, and identify a character $\chi \colon T \rightarrow \operatorname{U}(1)$ with the weight $d\chi/{2 \pi i} \in \mathfrak t^*$. We denote by $d\lambda$ the Lebesgue measure on $\mathfrak t^*$ that is normalized in such a way that any fundamental domain of the weight lattice $\Lambda^*$ has unit measure. Let us also choose a positive Weyl chamber $\mathfrak t^*_+ \subseteq \mathfrak t^*$; this determines a set of positive roots $\{ \alpha > 0 \} = \{ \alpha_1, \ldots, \alpha_R \} \subseteq \mathfrak t^*$. The set of negative roots is by definition $\{ -\alpha : \alpha > 0 \}$. We write $\mathfrak t^*_{>0}$ for the interior of the positive Weyl chamber, which contains the strictly dominant weights. We will often identify $\mathfrak k$ and its dual $\mathfrak k^*$, as well as $\mathfrak t$ and $\mathfrak t^*$, via some fixed $K$-invariant inner product $\braket{-,-}$ on $\mathfrak k$.

For the special unitary group $\operatorname{SU}(d)$, which we always take to be the group of unitary $d\times{}d$-matrices with unit determinant, we use the maximal torus consisting of diagonal matrices, on which the Weyl group $S_d$ acts by permuting diagonal entries. The Lie algebra ${\mathfrak{su}}(d)$ consists of anti-Hermitian matrices with trace zero, and our choice of invariant inner product is $\braket{X, Y} = -{\mathrm{tr}\left(X Y\right)}$. Using it to identify $\mathfrak t$ and $\mathfrak t^*$, a positive Weyl chamber $\mathfrak t^*_+$ is given by the set of diagonal matrices $\lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_d)$ with purely imaginary entries, summing to zero and arranged in such a way that $i \lambda_1 \geq \ldots \geq i \lambda_d$. This corresponds to choosing the positive roots $\alpha_{j,k}(\lambda) = i(\lambda_j - \lambda_k)$ with $j < k$. The points in the interior $\mathfrak t^*_{>0}$ are those $\lambda \in \mathfrak t^*_+$ with all distinct eigenvalues, i.e., $i \lambda_1 > \ldots > i \lambda_d$.

Let $M$ be a compact, connected Hamiltonian $K$-manifold of dimension $2n$, with symplectic form $\omega_M$ and a choice of moment map $\Phi_K \colon M \rightarrow \mathfrak k^*$ [@cannasdasilva08; @guilleminsternberg84]. The intersection $\Delta_K(M) = \Phi_K(M) \cap \mathfrak t^*_+$ of its image with the positive Weyl chamber is a compact convex polytope, called the *moment polytope* or *Kirwan polytope* [@guilleminsternberg82; @kirwan84b; @guilleminsjamaar05]. If $M$ is a coadjoint $K$-orbit $\mathcal O_\lambda$ through some $\lambda \in \mathfrak t^*_+$, it will always be equipped with the Kirillov–Kostant–Souriau symplectic form and the moment map induced by the inclusion $\mathcal O_\lambda \subseteq \mathfrak k^*$. Evidently, $\Delta_K(\mathcal O_\lambda) = \{\lambda\}$.

Throughout this article, we will always impose the following non-degeneracy condition:

\[main assumption\] The moment polytope $\Delta_K(M)$ has non-empty intersection with the interior of the positive Weyl chamber, $\mathfrak t^*_{>0}$.

In view of [@lermanmeinrenkentolmanetal98 Lemma 3.9] and well-known facts about compact Lie group actions, this assumption in fact implies the following: The set $\Phi_K^{-1}(K \cdot \mathfrak t^*_{>0})$ is an open, dense subset of $M$ whose complement has Liouville measure zero. We show in that does *not* restrict the applicability of our techniques to the problem of computing eigenvalue distributions of reduced density matrices.

The *Duistermaat–Heckman measure* $\operatorname{DH}^K_M$ is then defined as follows [@duistermaatheckman82]: Push forward the Liouville measure $\mu_M = \omega_M^n / ((2 \pi)^n n!)$ on $M$ along the moment map $\Phi_K$, compose with the push-forward along the quotient map $\tau_K \colon \mathfrak k^* \rightarrow \mathfrak t^*_+$ which sends all points in a coadjoint orbit $\mathcal O_\lambda$ to $\lambda$ in the positive Weyl chamber; then divide the resulting measure by the polynomial $p_K(\lambda) = \prod_{\alpha > 0} {\langle \lambda, \alpha \rangle} / {\langle \rho, \alpha \rangle}$, where $\rho$ is half the sum of positive roots. That is, $$\label{definition duistermaat-heckman measure}
  \operatorname{DH}^K_M = \frac 1 {p_K} (\tau_K)_* (\Phi_K)_* (\mu_M).$$ ensures that $\operatorname{DH}^K_M$ is a locally finite measure on the interior of the positive Weyl chamber. Its support is equal to the moment polytope. Note that $p_K(\lambda)$ is equal to the Liouville volume of a maximal-dimensional coadjoint orbit $\mathcal O_\lambda$ [@berlinegetzlervergne03 Proposition 7.26]. Therefore, the above is a natural definition to use in our context: It is normalized so that the Duistermaat–Heckman measure associated with the action of $K$ on a generic coadjoint orbit $\mathcal O_\lambda$ is a probability distribution concentrated at the point $\lambda$.

If $H$ is another compact, connected Lie group, with Lie algebra $\mathfrak h$, acting on $M$ via a group homomorphism ${\varphi}\colon H \rightarrow K$, then this action is also Hamiltonian, with a moment map given by the composition $$\label{restriction to subgroups}
  \Phi_H = (d{\varphi})^* \circ \Phi_K \colon M \rightarrow \mathfrak k^* \rightarrow \mathfrak h^*.$$ This in turn determines a Duistermaat–Heckman measure $\operatorname{DH}^H_M$. In particular, we can associate moment maps and Duistermaat–Heckman measures with all closed subgroups of $K$. In the case of the maximal torus $T \subseteq K$, we shall call $\Phi_T$ the *Abelian moment map* and $\operatorname{DH}^T_M$ the *Abelian Duistermaat–Heckman measure*, in order to distinguish them from the *non-Abelian moment map* $\Phi_K$ and the *non-Abelian Duistermaat–Heckman measure* $\operatorname{DH}^K_M$, respectively. Explicitly, $$\label{definition duistermaat-heckman measure abelian}
  \operatorname{DH}^T_M = (\Phi_T)_* (\mu_M) = (\pi_{K,T})_* (\Phi_K)_* (\mu_M).$$

Throughout this paper we shall assume for simplicity that the Abelian moment polytope $\Delta_T(M)$ is of maximal dimension. This can always be arranged for by replacing $T$ by the quotient $T / \bigcap_{m \in M} T_m$, where $T_m$ denotes the $T$-stabilizer of a point $m \in M$. If $T$ is the maximal torus of a semisimple Lie group, it follows already from that $\Delta_T(M)$ is of maximal dimension. As a consequence, $T$ acts locally freely on a dense, open subset whose complement has Liouville measure zero [@duistermaatheckman82 Lemma 3.1]. In particular, generic points in $M$ are regular for the Abelian moment map. Therefore the Abelian Duistermaat–Heckman measure is absolutely continuous with respect to Lebesgue measure on $\mathfrak t^*$. By the Duistermaat–Heckman Theorem, $\operatorname{DH}^T_M$ in fact has a polynomial density function of degree at most $n-r$ on each connected component of the set of regular values [@duistermaatheckman82 Corollary 3.3]. We shall call these components the *regular chambers*, and one can show that, except for the unbounded one, every such chamber is an open convex polytope. If the closures of two regular chambers have a common boundary of maximal dimension (i.e., of codimension one) then we shall say that the two chambers are *adjacent* and call the common boundary a *singular wall*.

All Hilbert spaces which we consider in this article are complex and finite-dimensional. We write $P_\psi$ for the orthogonal projection onto a one-dimensional subspace ${\mathbb C}\psi$, and ${\lVertX\rVert}_2 = \sqrt{{\mathrm{tr}\left(X^*X\right)}} = \sum_j s_j^2$ for the Hilbert–Schmidt norm of an operator $X$ with singular values $(s_j)$. We use $\braket{-,-}$ to denote inner products as well as the pairing between measures (or more general distributions) and test functions. We write $\delta_p$ for the *Dirac measure* at $p$, i.e., the probability measure concentrated at the point $p$, and $H_\omega$ for the *Heaviside measure* which is defined by $\langle H_\omega, f \rangle = \int_0^\infty f(t \omega) dt$. We sometimes use the letter ${\mathbf P}$ for probability distributions.

Throughout the paper when we speak of the quantum marginal problem we always refer to its one-body version as described in .

Finally, we offer a word of caution for people acquainted with the theory of geometric quantization [@guilleminsternberg77; @woodhouse92]: Our quantum states do not arise via some quantization procedure from a classical symplectic phase space. In contrast, herein, as detailed in below, the spaces of quantum states themselves are Hamiltonian manifolds. Probability distributions can be realized as quantum states of a special form, and the passage from quantum to classical is related to passing from a non-Abelian group to its maximal torus (see for precise statements). The “semiclassical limit” well-known in geometric quantization does not have an analogous physical meaning in our setting; its significance is solely to connect the symplectic geometry with representation theory (see ).

Density Matrices and Purification {#math phys dictionary}
=================================

The applicability of symplectic geometry to the quantum marginal problem relies on the close relation between the Lie algebra of $\operatorname{SU}(d)$ and the density matrices of quantum mechanics, and on the fact that restricting to certain subgroups has the physical meaning of passing to reduced density matrices, which describe the quantum state of subsystems. In this section we will describe this relationship in some detail (, ), and show how one can reduce the general problem of computing joint eigenvalue distributions of reduced density matrices to the case of globally pure quantum states, that is, to the Duistermaat–Heckman measure for a projective space (). We briefly discuss how probability distributions and the classical marginal problem are embedded in our setup () and describe some immediate physical applications ().

Density Matrices {#density matrices}
----------------

A *density matrix* is a positive-semidefinite Hermitian operator $\rho$ of trace one acting on a finite-dimensional Hilbert space $\mathcal H$. We will often choose coordinates and think of $\rho$ as a matrix. If $\rho$ is the orthogonal projection onto a one-dimensional subspace then we say that $\rho$ is a *pure state*; otherwise, it is a *mixed state*. An *observable* is an arbitrary Hermitian operator acting on $\mathcal H$.

Density matrices on $\mathcal H$ describe the state of a quantum system modeled by the Hilbert space $\mathcal H$: According to the postulates of quantum mechanics, the expectation value of an observable $O$ is given by the pairing ${\mathrm{tr}\left(O \rho\right)} \in {\mathbb R}$. Of course, $\rho$ is characterized by these expectation values, even if we use anti-Hermitian observables instead and restrict to trace zero (since the trace of $\rho$ is fixed). That is, we have an injection $$\label{density matrices to functionals}
  \rho \mapsto \left( X \mapsto i {\mathrm{tr}\left(X \rho\right)} \right) \in {\mathfrak{su}}(\mathcal H)^*$$ which extends to an isomorphism between the affine space of trace-one Hermitian operators on $\mathcal H$ and ${\mathfrak{su}}(\mathcal H)^*$. This isomorphism is $\operatorname{SU}(\mathcal H)$-equivariant; its inverse sends a coadjoint orbit $\mathcal O_\lambda$ to the set of Hermitian operators with eigenvalues $$\label{spectra to positive weyl chamber}
  \hat\lambda_j = \frac 1 {\dim \mathcal H} + i \lambda_j
  \quad
  (j=1,\ldots,\dim \mathcal H),$$ where the $\lambda_j$ are the eigenvalues of $\lambda$ (eigenvalues are repeated according to their multiplicity). Let us choose coordinates $\mathcal H \cong {\mathbb C}^d$ and identify $\operatorname{SU}(\mathcal H) \cong \operatorname{SU}(d)$ accordingly. Then the eigenvalues $\lambda_i$ are just the diagonal entries of the matrix $\lambda \in \mathfrak t^*_+$ labeling the coadjoint orbit (see for our conventions). It follows that defines a bijection between the positive Weyl chamber $\mathfrak t^*_+$ and the set of eigenvalue spectra of trace-one Hermitian operators, which we think of as elements of the set $\{ \hat\lambda \in \mathbf R^d : \hat\lambda_1 \geq \ldots \geq \hat\lambda_d, \sum_j \hat\lambda_j = 1 \}$.

Note that the set of pure states is identified with the coadjoint orbit through the highest weight of the defining representation of $\operatorname{SU}(\mathcal H)$, that is, with projective space ${\mathbb P}(\mathcal H)$: The density matrix corresponding to a point $[\psi] \in {\mathbb P}(\mathcal H)$ is simply the orthogonal projection $P_\psi$ onto ${\mathbb C}\psi$.

Moreover, the Liouville measure on a coadjoint orbit $\mathcal O_\lambda$ is identified via with the unique $\operatorname{SU}(\mathcal H)$-invariant measure on the set of Hermitian matrices with spectrum $\hat\lambda$, normalized to total volume $$\label{volume coadjoint su orbit}
  p_{\operatorname{SU}(\mathcal H)}(\lambda) =
  \prod_{j < k} \frac{i (\lambda_j - \lambda_k)}{k - j} =
  \prod_{j < k} \frac{(\hat\lambda_j - \hat\lambda_k)}{k - j}.$$ The non-Abelian Duistermaat–Heckman measure of $\mathcal O_\lambda$ corresponds to the Dirac measure at $\hat\lambda$, while the Abelian Duistermaat–Heckman measure corresponds to the distribution of the diagonal entries of a density matrix with spectrum $\hat\lambda$ chosen according to the invariant measure.

\[bloch ball\] The Lie algebra ${\mathfrak{su}}(2)$ is three-dimensional, generated by $-\frac i 2$ times the Pauli matrices $\sigma_x, \sigma_y, \sigma_z$. Functionals in its dual ${\mathfrak{su}}(2)^*$ can be identified with points in ${\mathbb R}^3$ by evaluating them at these generators. In this picture, the inverse of associates to a vector $\vec r = (x,y,z)$ the Hermitian matrix $\rho(\vec r) = \frac 1 2 ( {\mathbf 1}+ \vec r \cdot \vec\sigma )$, where $\vec\sigma$ is the Pauli vector $(\sigma_x, \sigma_y, \sigma_z)$. The $z$-axis is identified with the Lie algebra of the maximal torus, ${\mathbb R}\sigma_z$, its positive half-axis with our choice of positive Weyl chamber $\mathfrak t^*_+$, and $(0,0,2)$ with the corresponding positive root $\alpha > 0$. The coadjoint action of elements in $\operatorname{SU}(2)$ amounts to rotating the *Bloch vector* $\vec r$ via the two-fold covering map $\operatorname{SU}(2) \rightarrow \operatorname{SO}(3)$. Therefore, coadjoint orbits are spheres, commonly called *Bloch spheres* in quantum mechanics. They can be labeled by their radius $r$, that is, by their intersection with the positive half of the $z$-axis. Points on such a sphere correspond to Hermitian matrices with eigenvalue spectrum $(\frac {1+r} 2, \frac {1-r} 2)$. Note that $\rho(\vec r)$ is positive-semidefinite (i.e., a density matrix) if and only if $\vec r$ is contained in the unit ball of ${\mathbb R}^3$.

The non-Abelian moment map is just the inclusion map of a Bloch sphere into ${\mathbb R}^3$. Hence its composition with the quotient map $\tau_{\operatorname{SU}(2)}$ sends all points in a Bloch sphere of radius $r > 0$ to $r$, while the Abelian moment map projects all points onto the $z$-axis. The Liouville measure is equal to the usual round measure, normalized to total volume $r$. Therefore, the non-Abelian Duistermaat–Heckman measure of a Bloch sphere with radius $r$ is equal to the Dirac measure $\delta_r$, while the Abelian Duistermaat–Heckman measure is obtained by pushing forward the Liouville measure onto the $z$-axis (see ). As already observed by Archimedes, any two zones of the same height on a sphere have the same area. Hence this latter measure is proportional to Lebesgue measure on the interval $[-r,r]$. An analogous statement holds for arbitrary projective spaces ().

![A Bloch sphere, its height function and the induced measure: a generic coadjoint $\operatorname{SU}(2)$-orbit, its Abelian moment map and the Abelian Duistermaat–Heckman measure.[]{data-label="bloch ball figure"}](blochsphere.png){width="6cm"}

Observe that the components $\braket{X, \Phi_{\operatorname{SU}(\mathcal H)}} \colon \mathcal O_\lambda \rightarrow {\mathbb R}$ of the moment map send a quantum state to the expectation value of the corresponding observable $-i X$. Without loss of generality, we may assume that $X$ generates a one-dimensional torus and that $X$ has one-dimensional eigenspaces. Then $\braket{X, \Phi_{\operatorname{SU}(\mathcal H)}}$ is just the moment map for the action of the torus generated by $X$ and its distribution can be computed immediately by using the Abelian Heckman formula (). This gives a short and conceptual proof of the formula derived recently in [@venutizanardi12].

Reduced Density Matrices {#reduced density matrices}
------------------------

Composite quantum systems are modeled by the tensor product of the Hilbert spaces describing their constituents. It is useful to think of these subsystems as individual particles, although they can be of more general nature; for instance, the subsystems can describe different degrees of freedom such as position and spin. Depending on whether the particles are in principle distinguishable or indistinguishable, we distinguish two basic classes of composite systems, which are of fundamentally different nature.

If the quantum system is composed of $N$ *distinguishable particles*, its global quantum state is described by a density matrix on the tensor product $\mathcal H = \mathcal H_1 \otimes \ldots \otimes \mathcal H_N$, where the $\mathcal H_j$ are the Hilbert spaces describing the individual particles. Quantum mechanics also tells us that observables $O_j$ acting on a single subsystem $\mathcal H_j$ correspond to tensor product observables ${\mathbf 1}_{\mathcal H_1 \otimes \ldots \otimes \mathcal H_{j-1}} \otimes O_j \otimes {\mathbf 1}_{\mathcal H_{j+1} \otimes \ldots \otimes \mathcal H_N}$, which act by the identity on all other subsystems. By non-degeneracy of the inner product, there exists a unique density matrix $\rho_j$ on $\mathcal H_j$ such that $$\label{definition reduced density matrix}
  {\mathrm{tr}\left(({\mathbf 1}_{\mathcal H_1 \otimes \ldots \otimes \mathcal H_{j-1}} \otimes O_j \otimes {\mathbf 1}_{\mathcal H_{j+1} \otimes \ldots \otimes \mathcal H_N}) \rho\right)} = {\mathrm{tr}\left(O_j \rho_j\right)}$$ for all observables $O_j$. It describes the quantum state of the $j$-th subsystem.

The density matrix $\rho_j$ is called the *(one-body) reduced density matrix* or *quantum marginal* for the $j$-th particle of the quantum system.

Note that we can embed $\operatorname{SU}(\mathcal H_j)$ into $\operatorname{SU}(\mathcal H)$ by $U_j \mapsto {\mathbf 1}_{\mathcal H_1 \otimes \ldots \otimes \mathcal H_{j-1}} \otimes U_j \otimes {\mathbf 1}_{\mathcal H_{j+1} \otimes \ldots \otimes \mathcal H_N}$. This induces an embedding on the level of Lie algebras. The dual projection ${\mathfrak{su}}(\mathcal H)^* \rightarrow {\mathfrak{su}}(\mathcal H_j)^*$, given by restricting functionals to the subalgebra, is identified by with the map $\rho \mapsto \rho_j$. Similarly, the group homomorphism from the Cartesian product $\operatorname{SU}(\mathcal H_1) \otimes \ldots \otimes \operatorname{SU}(\mathcal H_N)$ to $\operatorname{SU}(\mathcal H)$ given by $(U_1,\ldots,U_N) \mapsto U_1 \otimes \ldots \otimes U_N$ induces the map $\rho \mapsto (\rho_1,\ldots,\rho_N)$ sending a density matrix to the tuple of all its one-body reduced density matrices.

The *(one-body) quantum marginal problem for distinguishable particles* asks for the possible tuples of one-body reduced density matrices $(\rho_1, \ldots, \rho_N)$ of an arbitrary density matrix $\rho$ with fixed spectrum, or, equivalently, for the possible tuples of their eigenvalues. By the above discussion, this is precisely equivalent to determining the moment polytope $\Delta_K(M)$ associated with the Hamiltonian action of the subgroup $K = \operatorname{SU}(\mathcal H_1) \otimes \ldots \otimes \operatorname{SU}(\mathcal H_N)$ on a coadjoint orbit $M = \mathcal O_{\tilde\lambda}$ for $\operatorname{SU}(\mathcal H)$, with moment map as defined in . The quantum marginal problem for globally pure states is the special case where $M = {\mathbb P}(\mathcal H)$. Moreover, the *joint eigenvalue distribution of reduced density matrices* we set out to compute in this article corresponds to the non-Abelian Duistermaat–Heckman measure $\operatorname{DH}^K_M$ as defined in : Up to the identification $\hat\lambda \mapsto \lambda$ between spectra of trace-one Hermitian operators and the positive Weyl chamber as defined in , it is given by $$\label{eigenvalue distribution dist}
  (\tau_K)_* (\Phi_K)_*\left(\frac {\mu_M} {\operatorname{vol}M}\right) =
  \frac 1 {\operatorname{vol}M} p_K \operatorname{DH}^K_M.$$ Note that we divide by the Liouville volume of $M$, which is just $p_{\operatorname{SU}(\mathcal H)}(\tilde\lambda)$, to obtain a probability measure.

Similarly, the joint distribution of the diagonal entries of the reduced density matrices corresponds to the Abelian Duistermaat–Heckman measure $\operatorname{DH}^T_M$.

If the quantum system is composed of *indistinguishable particles*, each particle is of course modeled by the same Hilbert space $\mathcal H_1$. The global state of the system is described by a density matrix supported on an irreducible sub-representation $\mathcal H \subseteq \mathcal H_1^{\otimes N}$, namely $\mathcal H = \operatorname{Sym}^N(\mathcal H_1)$ for *bosons* and ${\Lambda}^N(\mathcal H_1)$ for *fermions* (but we can in principle also consider other irreducible sub-representations which correspond to more exotic statistics). Note that since every such density matrix commutes with permutations, all the one-body reduced density matrices are equal.

Note that we can let single-particle observables $O$ act more intrinsically by the symmetric expression $\frac 1 N (O \otimes {\mathbf 1}_{\mathcal H_1^{\otimes (N-1)}} + \ldots + {\mathbf 1}_{\mathcal H_1^{\otimes (N-1)}} \otimes O)$, without changing their expectation values. Up to a factor $N$, this corresponds to the embedding of Lie algebras induced by the diagonal map $U \mapsto U \otimes \ldots \otimes U$, which is of course precisely the action of $\operatorname{SU}(\mathcal H_1)$ on the representation $\mathcal H$. This embedding therefore induces the map $\rho \mapsto \sum_j \rho_j = N \rho_1$ in the same way as described above. It follows that the *(one-body) quantum marginal problem for indistinguishable particles* amounts to determining $\frac 1 N \Delta_K(M)$ for the induced action of $K = \operatorname{SU}(\mathcal H_1)$ on a coadjoint orbit $M = \mathcal O_{\tilde\lambda}$ of $\operatorname{SU}(\mathcal H)$, and that, up to the identification $\hat\lambda \mapsto \lambda$, the *eigenvalue distribution of the reduced density matrix* is given by $$\label{eigenvalue distribution indist}
  \kappa_* (\tau_K)_* (\Phi_K)_*\left(\frac {\mu_M} {\operatorname{vol}M}\right) =
  \frac 1 {\operatorname{vol}M} p_K \kappa_*(\operatorname{DH}^K_M),$$ where the linear map $\kappa(\lambda) = \frac \lambda N$ counteracts the factor $N$ in the moment map.

                          Setting                   Hilbert space $\mathcal H$                                                 Group $K$
  ------------------------------- -------------------------------------------------------------- ----------------------------------------------------------------------
    $N$ distinguishable particles  ${\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N}$   $\operatorname{SU}(d_1) \times \ldots \times \operatorname{SU}(d_N)$
                       $N$ bosons             $\operatorname{Sym}^N({\mathbb C}^d)$                                      $\operatorname{SU}(d)$
                     $N$ fermions                  ${\Lambda}^N({\mathbb C}^d)$                                          $\operatorname{SU}(d)$

  : The quantum marginal problem is modeled by the action of the group $K$ on a coadjoint $\operatorname{SU}(\mathcal H)$-orbit $M = \mathcal O_{\tilde\lambda}$.[]{data-label="QMP summary table"}

In we have summarized the groups and spaces relevant for the quantum marginal problems of main physical interest, and we will focus on these in the remainder of this article. One can also combine both cases, e.g., to describe a quantum system composed of two different sorts of indistinguishable particles (as happens for the purified double of a bosonic or fermionic quantum marginal problem as defined in below), or a number of indistinguishable particles each of which have multiple internal degrees of freedom. In the latter case, arbitrary irreducible representations of the special unitary group can appear if one restricts to the reduced density matrices corresponding to only some of the degrees of freedom (see, e.g., [@klyachkoaltunbulak08]).

Purification
------------

Let $\mathcal H$ be an arbitrary finite-dimensional Hilbert space. It is well-known that every density matrix $\rho$ on $\mathcal H$ is the reduced density matrix of a pure state in $[\psi] \in {\mathbb P}(\mathcal H \otimes \mathcal H)$, called a *purification* of the quantum state $\rho$. Indeed, if $\rho = \sum_i p_i P_{v_i}$ is the spectral decomposition of $\rho$ then we can simply choose $\psi = \sum \sqrt{p_i} v_i \otimes v_i$. In this sense, the global state of a quantum system can always be described by a pure state; reduced density matrices occur only in the description of the states of its subsystems. This motivates the following definition:

For any unitary $K$-representation $\mathcal H$, we define the *purified double* to be the Hamiltonian $K \times \tilde K$-manifold ${\mathbb P}(\mathcal H \otimes \mathcal H)$, where $\tilde K = \operatorname{SU}(\mathcal H)$, equipped with the moment map constructed in the usual way by embedding into $\mathfrak u(\mathcal H \otimes \mathcal H)^*$ and “restricting” the functionals to elements in $\mathfrak k \oplus \mathfrak {\tilde k}$.

Observe that if $\mathcal H$ is one of the representations of modeling a setup of the quantum marginal problem, then the purified double corresponds to the pure-state quantum marginal problem where one has adjoined a single distinguishable particle modeled by $\mathcal H$.

The purification $[\psi] \in {\mathbb P}(\mathcal H \otimes \mathcal H)$ of a quantum state $\rho$ on $\mathcal H$ is unique up to a unitary acting on the second copy of $\mathcal H$. Evidently, such operations do not change the reduced density matrix $\rho = (P_\psi)_1$ and they leave the eigenvalue spectrum of $(P_\psi)_2$ invariant. In particular, the eigenvalue spectra of the reduced density matrices $(P_\psi)_1$ and $(P_\psi)_2$ are always equal. This implies that we can reduce the quantum marginal problem to the case of globally pure states, both for distinguishable and indistinguishable particles:

\[purification polytope\] Let $\mathcal O_{\tilde\lambda}$ be a coadjoint orbit of $\tilde K = \operatorname{SU}(\mathcal H)$, with $\tilde\lambda \in \mathfrak {\tilde t}^*_+$ corresponding to the eigenvalue spectrum of a density operator. Then $\lambda \in \Delta_K(\mathcal O_{\tilde\lambda})$ if and only if $(\lambda,\tilde\lambda) \in \Delta_{K \times \tilde K}({\mathbb P}(\mathcal H \otimes \mathcal H))$.

In other words, $$\Delta_{K \times \tilde K}({\mathbb P}(\mathcal H \otimes \mathcal H)) = \bigcup_{\tilde\lambda \in \tilde\Delta} \Delta_K(\mathcal O_{\tilde\lambda}) \times \{\tilde\lambda\},$$ where $\tilde\Delta = \{ \tilde\lambda \in \mathfrak {\tilde t}^*_+ : \widehat{\tilde\lambda}_j \geq 0 \}$ is the convex subset of the positive Weyl chamber corresponding to the eigenvalue spectra of density operators. We can similarly reduce the problem of determining the joint eigenvalue distribution to the case of globally pure states: For this, let us define probability measures $$\begin{aligned}
  \label{measures to be purified}
  {\mathbf P}&= (\tau_{K \times \tilde K})_* (\Phi_{K \times \tilde K})_* \left( \frac {\mu_{{\mathbb P}(\mathcal H \otimes \mathcal H)}} {\operatorname{vol}{\mathbb P}(\mathcal H \otimes \mathcal H)} \right),\\
  {\mathbf P}_{\tilde\lambda} &= (\tau_K)_* (\Phi_K)_* \left( \frac {\mu_{\mathcal O_{\tilde\lambda}}} {\operatorname{vol}O_{\tilde\lambda}} \right),
\end{aligned}$$ where $\operatorname{vol}$ denotes the Liouville volume.

\[purification measure\] The measures in are related by $$\langle {\mathbf P}, f \rangle
  = \frac 1 Z \int_{\tilde\Delta} d\tilde\lambda ~
    p^2_{\tilde K}(\tilde\lambda) ~
    \langle {\mathbf P}_{\tilde\lambda}, f(-,\tilde\lambda) \rangle,$$ for all test functions $f \in C_b(\mathfrak t^*_+ \oplus \mathfrak {\tilde t^*}_+)$, where $d\tilde\lambda$ is Lebesgue measure on $\mathfrak {\tilde t}^*_+$ and $Z$ a suitable normalization constant.

Each of the one-body reduced density matrices of a Liouville-distributed bipartite pure state in ${\mathbb P}(\mathcal H \otimes \mathcal H)$ is distributed according to the Hilbert–Schmidt measure restricted to the set of density matrices. In particular, its eigenvalues are distributed according to the well-known formula of [@lloydpagels88; @zyczkowskisommers01], so that $$(\tau_{\tilde K})_* (\Phi_{\tilde K})_* \left( \frac {\mu_{{\mathbb P}(\mathcal H \otimes \mathcal H)}} {\operatorname{vol}{\mathbb P}(\mathcal H \otimes \mathcal H)} \right)
    = \frac 1 Z \, p^2_{\tilde K}(\tilde\lambda) \, {\mathbf 1}_{\tilde\Delta}(\tilde\lambda) \, d\tilde\lambda,$$ with ${\mathbf 1}_{\tilde\Delta}$ the indicator function of $\tilde\Delta$ and $Z$ a suitable normalization constant. We have just seen that both reduced density matrices necessarily have equal eigenvalue spectrum. This implies that $$\langle (\Phi_{\tilde K \times \tilde K})_* \left( \frac {\mu_{{\mathbb P}(\mathcal H \otimes \mathcal H)}} {\operatorname{vol}{\mathbb P}(\mathcal H \otimes \mathcal H)} \right), g \rangle
    = \frac 1 Z \int_{\tilde\Delta} d\tilde\lambda \int_{\mathcal O_{\tilde\lambda} \times \mathcal O_{\tilde\lambda}} g.$$ See for an independent derivation using the techniques of this paper. It follows that $$\begin{aligned}
    \langle {\mathbf P}, f \rangle
    &= \langle (\tau_{K \times \tilde K})_* (\Phi_{K \times \tilde K})_* \left( \frac {\mu_{{\mathbb P}(\mathcal H \otimes \mathcal H)}} {\operatorname{vol}{\mathbb P}(\mathcal H \otimes \mathcal H)} \right), f \rangle \\
    &= \langle (\tau_{K} \Phi_{K} \times \tau_{\tilde K})_* (\Phi_{\tilde K \times \tilde K})_* \left( \frac {\mu_{{\mathbb P}(\mathcal H \otimes \mathcal H)}} {\operatorname{vol}{\mathbb P}(\mathcal H \otimes \mathcal H)} \right), f \rangle\\
    &= \frac 1 Z \int_{\tilde\Delta} d\tilde\lambda \int_{\mathcal O_{\tilde\lambda} \times \mathcal O_{\tilde\lambda}} (\tau_{K} \Phi_{K} \times \tau_{\tilde K})^* \big( f \big)\\
    &= \frac 1 Z \int_{\tilde\Delta} d\tilde\lambda ~ p_{\tilde K}(\tilde\lambda) ~ \int_{\mathcal O_{\tilde\lambda}} (\tau_{K} \Phi_{K})^* \left( f(-,\tilde\lambda) \right)\\
    &= \frac 1 Z \int_{\tilde\Delta} d\tilde\lambda ~ p^2_{\tilde K}(\tilde\lambda) ~ \langle {\mathbf P}_{\tilde\lambda}, f(-,\tilde\lambda) \rangle.
    \qedhere
  \end{aligned}$$

Note that ${\mathbf P}_{\tilde\lambda}$, and in particular the eigenvalue distributions and , vary continuously with the global spectrum $\tilde\lambda$. therefore implies that we can reconstruct them from the eigenvalue distribution for the purified double by taking limits. If the latter distribution has a continuous Lebesgue density, as will often be the case, then we can simply restrict this density function to the global spectrum $\tilde\lambda$ of interest.

We will now show that is always satisfied when working with the purified double. In quantum-mechanical terms, we have to show that there exists a global pure state $[\psi] \in {\mathbb P}(\mathcal H \otimes \mathcal H)$ such that the eigenvalue spectra of all the reduced density matrices are non-degenerate (with respect to the quantum marginal problem where we have added a single distinguishable particle with Hilbert space $\mathcal H$).

For distinguishable particles, where $\mathcal H \cong {\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N}$, this follows from the following general criterion, since the purified double is constructed by adding an additional Hilbert space of dimension $d_{N+1} = \dim \mathcal H = d_1 \cdots d_N$:

\[main assumption qmp\] Let $N \geq 1$ and $d_1 \leq \ldots \leq d_N \leq d_{N+1}$. Then there exists a global pure state in ${\mathbb P}({\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N} \otimes {\mathbb C}^{d_{N+1}})$ whose one-body reduced density matrices have non-degenerate eigenvalue spectra if and only if $$d_{N+1} \leq \left( \prod_{i=1}^N d_i \right) + 1.$$

The condition is clearly necessary, since it follows from the singular value decomposition that at most $\prod_{i=1}^N d_i$ eigenvalues of $\rho_{N+1}$ can be non-zero.

For sufficiency, let us construct a state with the desired property: For this, we consider the standard tensor product basis vectors $e_{i_1} \otimes \ldots \otimes e_{i_N}$ of ${\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N}$, labelled by integers $i_j \in \{1,\ldots,d_j\}$, $j=1,\ldots,N$. We choose a subset of $d_{N+1}-1$ many such basis vectors $e_{i_1(k)} \otimes \ldots \otimes e_{i_N(k)}$ in such a way that, for each subsystem $j=1,\ldots,N$, at least $d_j - 1$ of the $d_j$ integers occur. This clearly is possible by our assumptions. Finally, we set $$\psi = \sum_{k=1}^{d_{N+1}-1} 2^{-k} e_{i_1(k)} \otimes \ldots \otimes e_{i_N(k)} \otimes e_k.$$ Then $[\psi]$ is a pure state such that all of its one-body reduced density matrices have non-degenerate eigenvalue spectrum.

For bosons and fermions, we will use the following lemma to show that is satisfied:

\[boson fermion lemma\] The convex hull of the weights of $\operatorname{Sym}^N({\mathbb C}^d)$ has maximal dimension. The same is true for ${\Lambda}^N({\mathbb C}^d)$ if $N<d$.

In the following we write $\omega_j$ for the weight corresponding to the character $\operatorname{diag}(t_1, \ldots, t_d) \mapsto t_j$, with $j=1, \ldots, d$.

\(1) $\operatorname{Sym}^N({\mathbb C}^d)$: The vectors $e_j \otimes \ldots \otimes e_j$ are weight vectors of weight $N \omega_j$, with $j=1,\ldots,d$. Clearly, the convex hull of these weights already has maximal dimension.

\(2) ${\Lambda}^N({\mathbb C}^d)$, where $1 \leq N < d$: It is well-known that the weights are given by $\sum_{j \in J} \omega_j$ for all $N$-element subsets $J \subseteq \{1,\ldots,d\}$. (The corresponding weight vectors are the well-known occupation number basis vectors for fermions.) Fix any such weight, say, the one corresponding to $J = \{1,\ldots,N\}$. The difference vectors between this weight and the weights obtained by replacing a single element of $J$ are proportional to the positive roots $\alpha_{j,k}(\lambda) = i(\lambda_j - \lambda_k)$ with $j \in J$ and $k \in \{1,\ldots,d\} \setminus J$. There are at least $d-1$ such roots, and they form a basis of $\mathfrak t^*$.

implies that is also satisfied for the purified double of the bosonic and fermionic quantum marginal problems: Indeed, it clearly suffices to show that in each case there exists a density operator $\rho$ on $\mathcal H$ such that both $\rho$ and its one-body reduced density matrix $\rho_1$ have non-degenerate eigenvalue spectrum (then any purification of $\rho$ has the desired properties). By , there exists a convex combination of weights $\sum_k p_k \omega_k \in \mathfrak t^*_{>0}$. By perturbing slightly, we can arrange for the weights $(p_k)$ to be mutually disjoint. Choose corresponding (orthogonal) weight vectors $v_k \in \mathcal H$ and consider the density matrix $\rho = \sum_k p_k P_{v_k}$. Clearly, both $\rho$ and its one-body reduced density matrix $\rho_1$have non-degenerate eigenvalue spectrum.

To summarize, we have shown that the problem of computing the joint eigenvalue distribution of reduced density matrices is equivalent to the computation of Duistermaat–Heckman measures associated with certain Hamiltonian group actions (cf. ). Moreover, by passing to the purified double, we can always reduce to the case where $M = {\mathbb P}(\mathcal H)$ is a projective space satisfying .

Probability Distributions {#classical}
-------------------------

Under the identification , elements of the dual of the Lie algebra of the maximal torus correspond to diagonal density matrices. These are precisely the diagonal matrices with non-negative entries summing to one, and can therefore be interpreted as *probability distributions* of a random variable $Z$ with values in the orthonormal basis $(e_i)$ we have chosen. This interpretation is in agreement with quantum mechanics: If we perform an actual measurement of a density matrix $\rho$ with respect to this orthonormal basis then the probability of getting outcome $e_i$ is given precisely by the diagonal element ${\mathrm{tr}\left(P_{e_i} \rho\right)} = \braket{e_i, \rho \, e_i} = \rho_{i,i}$.

Note that the moment map for the action of the maximal torus $\tilde T \subseteq \operatorname{SU}(\mathcal H)$ on the projective space ${\mathbb P}(\mathcal H)$ corresponds to sending a pure state $[\psi]$ onto its diagonal. As we vary $[\psi]$ over all pure states in ${\mathbb P}(\mathcal H)$, the diagonal entries attain all possible probability distributions. In other words, the Abelian moment polytope $\Delta_{\tilde T}({\mathbb P}(\mathcal H))$ is just the simplex $\tilde\Delta$ defined in . The corresponding Duistermaat–Heckman measure is equal to a suitably normalized Lebesgue measure on $\tilde\Delta$ (this is a special case of below).

Now consider as in the case of $N$ distinguishable particles. Choose orthonormal bases to identify $\mathcal H_k \cong {\mathbb C}^{d_k}$, and therefore $\mathcal H \cong {\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N}$ using the tensor product basis. Note that we can interprete diagonal density matrices $\rho$ on $\mathcal H$ as the joint probability distribution of a tuple of random variables $(Z_1,\ldots,Z_N)$, where each $Z_k$ takes values in the standard basis of the corresponding ${\mathbb C}^{d_k}$, by setting $${\mathbf P}(Z_1 = e_{i_1}, \ldots, Z_N = e_{i_N}) =
  {\mathrm{tr}\left(P_{e_{i_1} \otimes \ldots \otimes e_{i_N}} \rho\right)}.$$ The marginal distributions of the random variables $Z_k$ in the sense of probability theory are then given by $$\begin{aligned}
  {\mathbf P}(Z_k = e_{i_k})
  &= \sum_{i_1, \ldots, \check{i_k}, \ldots, i_N} {\mathrm{tr}\left(P_{e_{i_1} \otimes \ldots \otimes e_{i_N}} \rho\right)}\\
  &= {\mathrm{tr}\left( \left( {\mathbf 1}_{{\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_{k-1}}} \otimes P_{e_{i_k}} \otimes {\mathbf 1}_{{\mathbb C}^{d_{k+1}} \otimes \ldots \otimes {\mathbb C}^{d_N}} \right) \rho\right)}
   = {\mathrm{tr}\left(P_{e_{i_k}} \rho_k\right)},\end{aligned}$$ where for the second identity we have used that $\rho$ is a diagonal matrix. That is, the marginal distributions of the $Z_k$ are precisely described by the reduced density matrices $\rho_k$ (i.e., by the quantum marginals), which are also diagonal if $\rho$ is diagonal.

Accordingly, the moment polytope $\Delta_T({\mathbb P}(\mathcal H))$ for the action of the maximal torus $T \subseteq \operatorname{SU}(d_1) \times \ldots \operatorname{SU}(d_N)$ on the set of pure states describes the tuples of marginal probability distributions that arise from joint distributions of the $(Z_1,\ldots,Z_N)$. This *(univariate) classical marginal problem* is of course trivial, since there are no constraints on the joint distribution. However, its quantitative version, which corresponds to computing the Abelian Duistermaat–Heckman measure $\operatorname{DH}^T_{{\mathbb P}(\mathcal H)}$, is interesting and not at all trivial to solve. In fact, the problem of computing joint eigenvalue distributions of reduced density matrices, which we set out to solve in this article, can be reduced to the computation of $\operatorname{DH}^T_{{\mathbb P}(\mathcal H)}$. This reduction, or rather the generalization which we describe in below, is at the core of the algorithms presented in .

Physical Applications {#physical applications}
---------------------

As indicated in the introduction, the eigenvalue distributions and have direct applications to quantum physics. In quantum statistical mechanics, among others, one typically studies bipartite setups $\mathcal H = \mathcal H_S \otimes \mathcal H_E$ composed of a system $S$ and an environment (or bath) $E$. Randomly-chosen pure states give rise to a distribution of reduced density matrices $\rho_S$, whose properties vary with the size of the environment. Physical motivations have lead to the computation of the corresponding eigenvalue distribution [@lloydpagels88], which we can easily re-derive using the techniques of this paper (). Note that many basic physical quantities are functions of the eigenvalues, such as the *von Neumann entropy* $$H(S) = H(\rho_S) = -{\mathrm{tr}\left(\rho_S \log \rho_S\right)} = \sum_j -\hat\lambda_j \log \hat\lambda_j,$$ where $(\hat\lambda_j)$ are the eigenvalues of $\rho_S$, or more general Rényi entropies and purities (cf. ). The average von Neumann entropy of a subsystem [@lubkin78; @page93] in particular has featured in the analysis of the black hole entropy paradox [@HaydenPreskill].

We can also consider other coadjoint orbits such as Grassmannians: Here, the density matrix corresponding to a $d$-dimensional subspace $\mathcal H' \subseteq \mathcal H_A \otimes \mathcal H_E$ is the normalized projection operator $\rho = {\mathbf 1}_{\mathcal H'}/d$, and the reduced density matrix $\rho_A$ is interpreted as a *canonical state* in the sense of statistical mechanics [@popescushortwinter06; @lloyd06; @goldsteinlebowitztumulkaetal06]. The probability distributions we compute can therefore be used to analyze the typical behavior of canonical states.

The tripartite case, in itself already interesting from the perspective of the quantum marginal problem, is also highly relevant to applications: It corresponds to the situation where $S$ itself is composed of two particles $A$ and $B$, so that $\mathcal H = \mathcal H_A \otimes \mathcal H_B \otimes \mathcal H_E$. In the study of quantum entanglement, remarkable recent progress has been made by analyzing the entanglement properties of the two-body reduced density matrix $\rho_{AB}$ of a randomly-chosen pure state in large dimensions, where the concentration of measure phenomenon occurs [@HaydenRandomizing; @haydenleungwinter06; @aubrunszarekye11b; @aubrunszarekye11; @collinsnechitaye11]. In particular, a negative resolution of the additivity conjecture of quantum information theory [@shor-additivity] has recently been obtained by related methods [@hastings-additivity; @aubrun-hastings]. The joint eigenvalue distribution of the reduced density matrices in particular determines *quantum conditional entropies* and *quantum mutual informations*, that is, the quantities $$\begin{aligned}
  H(A|B) &= H(AB) - H(B) = H(E) - H(B),\\
  I(A:B) &= H(A) + H(B) - H(AB) = H(A) + H(B) - H(E),\end{aligned}$$ since the eigenvalue spectra of $\rho_{AB}$ and $\rho_E$ are equal (cf. ). They have immediate applications to entanglement theory; for example, the quantum mutual information provides an upper bound on the amount of entanglement that can be distilled from a quantum state [@christandlwinter04].

In all these applications, most known results are for large Hilbert spaces, since the techniques employed rely on asymptotic features such as measure concentration. Our algorithms require no such assumption. In particular, they are well-suited for low-dimensional systems, which previously remained inaccessible.

Derivative Principle for Invariant Measures {#derivative principle}
===========================================

In this section we will describe a fundamental property of $K$-invariant measures on $\mathfrak k^*$ that are concentrated on the union of the maximal-dimensional coadjoint orbits (that is, on $K \cdot \mathfrak t^*_{>0}$). Every such invariant measure can be reconstructed from its projection onto $\mathfrak t^*$ by taking partial derivatives in the direction of negative roots (). In particular, this implies that the non-Abelian Duistermaat–Heckman measure $\operatorname{DH}^K_M$ can be recovered from the Abelian Duistermaat–Heckman measure $\operatorname{DH}^T_M$ ().

For the invariant probability measure supported on a single coadjoint orbit of maximal dimension, this follows from a well-known formula of Harish-Chandra, as was already observed by Heckman:

\[ableitungsformel coadjoint orbit\] Let $\mathcal O_\lambda$ be a coadjoint orbit through $\lambda \in \mathfrak t^*_{>0}$. Then, $${\left.\left( \prod_{\alpha > 0} \partial_{-\alpha} \right) \operatorname{DH}^T_{\mathcal O_\lambda}\vphantom{\big|}\right|_{\mathfrak t^*_{>0}}} = \delta_\lambda,$$ where the partial derivatives and the restriction are in the sense of distributions.

Harish-Chandra’s formula for the Fourier transform of a coadjoint orbit states that $$\langle \operatorname{DH}^T_{\mathcal O_\lambda}, e^{i \langle -, X \rangle} \rangle =
    \sum_{w \in W} (-1)^{l(w)} e^{i \langle w\lambda, X \rangle} \prod_{\alpha > 0} \frac 1 {i \langle \alpha, X \rangle}$$ for every $X \in \mathfrak t$ which is not orthogonal to a root (see [@berlinegetzlervergne03 Corollary 7.25] for a recent account). Here, $l(w)$ is the length of the Weyl group element $w$. This implies that the Abelian Duistermaat–Heckman measure is given by the following alternating sum of convolutions $$\label{harish chandra}
    \operatorname{DH}^T_{\mathcal O_\lambda} = \sum_{w \in W} (-1)^{l(w)} \delta_{w \lambda} \star H_{-\alpha_1} \star \ldots \star H_{-\alpha_R}.$$ Recall that $H_\omega$ is the Heaviside measure defined in by $\langle H_\omega, f \rangle = \int_0^\infty f(t \omega) dt$. By the fundamental theorem of calculus, we have $\partial_\omega H_\omega = \delta_0$. Therefore, $$\label{harish chandra derivative}
    \left( \prod_{\alpha > 0} \partial_{-\alpha} \right) \operatorname{DH}^T_{\mathcal O_\lambda} = \sum (-1)^{l(w)} \delta_{w \lambda},$$ and the assertion follows if we restrict to the interior of the positive Weyl chamber.

Every Bloch sphere of radius $r > 0$ is a coadjoint orbit of maximal dimension (cf. ). We have seen that $\operatorname{DH}^T_{\mathcal O_r}$ is equal to $\frac 1 2 {\mathbf 1}_{[-r,r]}(z) dz$, where $dz$ is Lebesgue measure on the $z$-axis. In agreement with , we observe that $${\left.\partial_\alpha \operatorname{DH}^T_{O_r}\vphantom{\big|}\right|_{\mathfrak t^*_{>0}}} =
    {\left.\partial_z {\mathbf 1}_{[-r,r]}(z) dz\vphantom{\big|}\right|_{{\mathbb R}_{> 0}}} =
    \delta_r.$$

\[main theorem measures\] Let $\nu$ be a $K$-invariant Radon measure on $\mathfrak k^*$ satisfying $\nu(\mathfrak k^* \setminus K \cdot \mathfrak t^*_{>0}) = 0$. Then, $${\left.\left(\prod_{\alpha > 0} \partial_{-\alpha} \right) (\pi_{K,T})_*(\nu)\vphantom{\big|}\right|_{\mathfrak t^*_{>0}}} =
    {\left.\frac 1 {p_K} (\tau_K)_*(\nu)\vphantom{\big|}\right|_{\mathfrak t^*_{>0}}},$$ where the partial derivatives and the restriction are in the sense of distributions.

Let $f \in C_c^\infty(\mathfrak t^*_{>0})$ be a test function, which we extend by zero to all of $\mathfrak t^*$, and set $g := (\pi_{K,T})^*\left(\left(\prod_{\alpha > 0} \partial_\alpha \right) f\right)$. By definition and assumption, respectively, $$\langle {\left. \left(\prod_{\alpha > 0} \partial_{-\alpha} \right) (\pi_{K,T})_*(\nu)\vphantom{\big|}\right|_{\mathfrak t^*_{>0}}}, f \rangle =
    \langle \nu, g \rangle =
    \langle {\nu}\big|_{K \cdot \mathfrak t^*_{>0}}, g \rangle.$$ Since $\nu$ is a $K$-invariant measure, we can use Fubini’s theorem to replace $g$ by its $K$-average. On each maximal-dimensional coadjoint orbit $\mathcal O_\lambda \subseteq K \cdot \mathfrak t^*_{>0}$, this average is given by $$\frac 1 {\operatorname{vol}\mathcal O_\lambda} \langle \mu_{\mathcal O_\lambda}, g \rangle =
    \frac 1 {p_K(\lambda)} \langle \mu_{\mathcal O_\lambda}, (\pi_{K,T})^*\left(\left(\prod_{\alpha > 0} \partial_\alpha \right) f\right) \rangle =
    \frac 1 {p_K(\lambda)} \langle \left(\prod_{\alpha > 0} \partial_{-\alpha} \right) \operatorname{DH}^T_{\mathcal O_\lambda}, f \rangle,$$ which by is precisely equal to ${f(\lambda)}/{p_K(\lambda)}$. In other words, the averaged function is on $K \cdot \mathfrak t^*_{>0}$ equal to the pullback $(\tau_K)^* \left( f/{p_K} \right)$. We conclude that $$    \langle {\nu}\big|_{K \cdot \mathfrak t^*_{>0}}, g \rangle =
    \langle {\nu}\big|_{K \cdot \mathfrak t^*_{>0}}, (\tau_K)^* \left( \frac f {p_K} \right) \rangle =
        \langle {\left.\frac 1 {p_K} (\tau_K)_* (\nu)\vphantom{\big|}\right|_{\mathfrak t^*_{>0}}}, f \rangle.
    \qedhere$$

\[main theorem\] The Duistermaat–Heckman measures as defined in are related by $${\left.\left(\prod_{\alpha > 0} \partial_{-\alpha} \right) \operatorname{DH}^T_M\vphantom{\big|}\right|_{\mathfrak t^*_{>0}}} =
    {\operatorname{DH}^K_M}\bigg|_{\mathfrak t^*_{>0}}.$$

guarantees that we can apply to the push-forward of the Liouville measure along the non-Abelian moment map $\Phi_K$.

This is the *derivative principle* alluded to in the title of this section. As we shall see in the following, it is a powerful tool for lifting results about the Duistermaat–Heckman measure for torus actions to general compact Lie group actions.

According to [@woodward05 §3.5], was already known to Paradan and also follows from a different result of Harish-Chandra. In we will describe another way to establish it by using the connection between Duistermaat–Heckman measures in algebraic geometry and multiplicities in group representations.

Note that completely determines the measure $\nu$ from its projection onto $\mathfrak t^*$, since $\nu$ is by assumption concentrated on the union of the coadjoint orbits of maximal dimension. Similarly, the non-Abelian Duistermaat–Heckman measure $\operatorname{DH}^K_M$ can be fully reconstructed from $\operatorname{DH}^T_M$ by using .

We stress that it is oftentimes not necessary to explicitely compute the non-Abelian Duistermaat-Heckman measure. Indeed, is of course by definition equivalent to $$\langle \operatorname{DH}^K_M, f \rangle = \langle \operatorname{DH}^T_M, \left(\prod_{\alpha > 0} \partial_\alpha \right) f \rangle$$ for all $f \in C_c^\infty(\mathfrak t^*_{>0})$, so that we can reduce the computation of averages over $\operatorname{DH}^K_M$ directly to integrations with respect to the Abelian Duistermaat–Heckman measure (cf. proof of ).

\[finite union of regular chambers\] It follows from and the discussion in that, on each (open) regular chamber, the non-Abelian Duistermaat–Heckman measure also has a polynomial density, namely the partial derivative in the directions of the negative roots of the density of the Abelian measure. However, there could still be non-zero measure on the singular walls separating the regular chambers. If we would like to exclude this then we need to understand the smoothness properties of the Abelian density function in the vicinity of singular walls, or, equivalently, the nature of the term by which the polynomial density changes when crossing a singular wall. If this jump term vanishes to order at least $R$ on the wall, then the Abelian density function is at least $R$-times weakly differentiable in the vicinity of the wall, and therefore the non-Abelian Duistermaat–Heckman density is also absolutely continuous there. This vanishing condition can be checked explicitly for each singular wall using the jump formula described in .

In case the vanishing condition is satisfied, the non-Abelian moment polytope $\Delta_K(M)$ is equal to the closure of a finite union of regular chambers for the Abelian moment map: Indeed, on each regular chamber the density polynomial is either equal to zero, or it is non-zero on an open, dense subset.

We cannot resist giving an easy application of to the generalized Horn’s problem: Here, the goal is to describe the sum of two coadjoint orbits $\mathcal O_\lambda + \mathcal O_\mu$ (Horn’s original problem referred to coadjoint orbits of the special unitary group). In other words, one considers the diagonal action of $K$ on $\mathcal O_\lambda \times \mathcal O_\mu$, which is Hamiltonian with moment map $(X,Y) \mapsto X+Y$, and one would like to describe the associated moment polytope or Duistermaat–Heckman measure.

\[main theorem horn\] Let $\lambda \in \mathfrak t^*_{>0}$ and $\mu \in \mathfrak t^*_+$. Then, $$\operatorname{DH}^K_{\mathcal O_\lambda \times \mathcal O_\mu} =
    \sum_{w \in W} (-1)^{l(w)} \delta_{w \lambda} \star \operatorname{DH}^T_{\mathcal O_\mu},$$ where $l(w)$ is the length of the Weyl group element $w$.

Clearly, since $T$ is Abelian, $$\operatorname{DH}^T_{\mathcal O_\lambda \times \mathcal O_\mu} =
    \operatorname{DH}^T_{\mathcal O_\lambda} \star \operatorname{DH}^T_{\mathcal O_\mu}.$$ Since $\lambda + \mu \in \mathfrak t^*_{>0}$, is satisfied. Therefore, is applicable, and the assertion follows together with , $$\operatorname{DH}^K_{\mathcal O_\lambda \times \mathcal O_\mu} =
    \left( \prod_{\alpha > 0} \partial_{-\alpha} \right) \operatorname{DH}^T_{\mathcal O_\lambda} \star \operatorname{DH}^T_{\mathcal O_\mu} =
    \sum_{w \in W} (-1)^{l(w)} \delta_{w \lambda} \star \operatorname{DH}^T_{\mathcal O_\mu}.
    \qedhere$$

The general case where both $\lambda$ and $\mu$ are contained in the boundary of the positive Weyl chamber can be treated as in [@dooleyrepkawildberger93] by taking limits. Of course we can also expand $\operatorname{DH}^T_{\mathcal O_\mu}$ as an alternating sum of convolutions by using or its version for lower-dimensional coadjoint orbits [@berlinegetzlervergne03 Theorem 7.24].

Algorithms for Duistermaat–Heckman Measures {#algorithms}
===========================================

In this section we present two algorithms for computing Duistermaat–Heckman measures. Both algorithms are based on the derivative principle from , in that they first compute the Abelian measure and then take partial derivatives according to .

The first algorithm, the *Heckman algorithm*, is based on the Heckman formula by Guillemin, Lerman and Sternberg, which expresses the Abelian measure as an alternating sum of iterated convolutions of Heaviside measures. The density function of each such convolution is piecewise polynomial and can be evaluated inductively using recent work of Boysal and Vergne. While very useful for computing low-dimensional examples, the resulting algorithm is rather inefficient due to the large number of summands.

Our second algorithm, the *single-summand algorithm*, is based on another formula for the Abelian Duistermaat–Heckman measure in the case where $M$ is the projective space of an arbitrary finite-dimensional representation. It turns out that this formula is equivalent to evaluating a single iterated convolution of the above form (hence the name of the algorithm). It can therefore be computed in a similar way, but much more efficiently. Since by passing to the purified double the quantum marginal problem can always be reduced to the case where $M$ is a projective space (), this solves the problem of computing eigenvalue distributions of reduced density matrices in complete generality.

Heckman Algorithm {#heckman}
-----------------

Before stating the Heckman formula by Guillemin, Lerman and Sternberg, let us recall the following renormalization process as described in [@guilleminlermansternberg88]:

Suppose that there are only finitely many fixed points of the action of the maximal torus $T$ on $M$. For each such fixed point $p \in M^T$, consider the induced representation of $T$ on the tangent space $T_p M$. The weights of this representation are called isotropy weights and we can always choose a vector $\gamma \in \mathfrak t^*$ which is non-orthogonal to all isotropy weights (for all tangent spaces). The process of multiplying by $-1$ those isotropy weights that have negative inner product with $\gamma$ is then called *renormalization*, and the resulting weights are called renormalized weights. See for a discussion of the case where $M$ is a projective space and for examples.

\[abelian heckman\] Suppose that there are only finitely many torus fixed points $p \in M^T$. Denote by $n_p$ the number of isotropy weights in $T_p M$ that are multiplied by $-1$ during renormalization and by $\hat\omega_{p,1}, \ldots, \hat\omega_{p,n}$ the resulting renormalized weights. Then, $$\operatorname{DH}^T_M = \sum_{p \in M^T} (-1)^{n_p} \delta_{\Phi_T(p)} \star H_{\hat\omega_{p,1}} \star \ldots \star H_{\hat\omega_{p,n}},$$ with $H_{\hat\omega}$ the Heaviside measure defined by $\langle H_{\hat\omega}, f \rangle = \int_0^\infty dt f(\hat\omega t)$.

In other words, the stationary phase approximation for the Fourier transform of an Abelian Duistermaat–Heckman measure is exact. This generalizes the Harish-Chandra formula for coadjoint orbits , which we used to establish .

Observe that each summand of the Heckman formula can be written as the push-forward of the standard Lebesgue measure $dt$ on ${\mathbb R}^n_{\geq 0}$ along a linear map of the form $\hat P \colon (t_k) \mapsto \sum_{k=1}^n t_k \hat\omega_k$, translated by $\Phi_T(p)$, since $$\label{convolution translation}
  H_{\hat\omega_1} \star \ldots \star H_{\hat\omega_n} =
  {\hat P}_*(H_{e_1} \star \ldots \star H_{e_n}) =
  {\hat P}_*({\left.dt\vphantom{\big|}\right|_{{\mathbb R}^n_{\geq 0}}}).$$ In a recent paper [@boysalvergne09], Boysal and Vergne have analyzed general push-forward measures of this form under the assumption that the vectors $\hat\omega_k$ span a proper convex cone (i.e., a convex cone of maximal dimension that does not contain any straight line). This ensures that the measure is locally finite and absolutely continuous with respect to Lebesgue measure on $\mathfrak t^*$. This assumption is certainly satisfied for the renormalized isotropy weights occurring in the Heckman formula (by the very definition of renormalization and our assumption that the Abelian moment polytope has maximal dimension).

Let us briefly review their results: It is well-known that the push-forward measure has a piecewise homogeneous polynomial density function of degree $n-r$. Here, the *chambers* are the connected components of the complement of the cones spanned by at most $r-1$ of the weights $(\hat\omega_k)$. Except for the unbounded chamber, they are open convex cones. *Walls* are by definition the convex cones spanned by $r-1$ linearly independent weights.[^1] Similarly to , if the common boundary of the closure of two chambers is of maximal dimension then this common boundary is a wall; moreover, every wall arises in this way. Note that the union of the walls is precisely the complement of the union of the chambers.[^2]

Let $\hat\Delta_\pm$ be two adjacent chambers which are separated by a wall $\hat W$, and choose a normal vector $\hat\xi \in \mathfrak t^*$ pointing from $\hat\Delta_-$ to $\hat\Delta_+$. Order the weights such that precisely $\hat\omega_1, \ldots, \hat\omega_m$ lie on the linear hyperplane spanned by $\hat W$. In the following, we shall freely identify differential forms and the measures induced by them. Denote by $d\hat w$ the Lebesgue measure on the hyperplanes parallel to $\hat W$, normalized in such a way that $$\label{wall quotient measure cone}
  d\lambda = d\hat w \wedge d\hat\xi,$$ where $d\hat\xi$ is the pullback of the standard volume form of ${\mathbb R}$ along the coordinate function $\langle -, \hat\xi\rangle$. Denote by $\hat f_\pm$ the homogeneous polynomials describing the density function $\hat f$ on $\hat\Delta_\pm$. Finally, consider the push-forward of Lebesgue measure on ${\mathbb R}^m_{\geq 0}$ along the linear map $\hat P_{\hat W} \colon (u_k) \mapsto \sum_{k=1}^m u_k \hat\omega_k$. Its density with respect to $d\hat w$ is given by a single homogeneous polynomial on the wall $\hat W$, since $\hat W$ is always contained in the closure of a chamber for $\hat P_{\hat W}$. Denote by $\hat f_{\hat W}$ any polynomial function extending it to all of $\mathfrak t^*$. Then the result of Boysal and Vergne is the following [@boysalvergne09 Theorem 1.1]: The jump of the density function across the wall is given by $$\label{boysal vergne}
    \hat f_+(\hat\lambda) - \hat f_-(\hat\lambda)
  = {\left.\operatorname{Res}\vphantom{\big|}\right|_{z=0}} \left(
      \hat f_{\hat W}(\partial_{\hat x})
      \frac
        {e^{\langle \hat\lambda, {\hat x} + z\hat\xi \rangle}}
        {\prod_{k=m}^n \langle \hat\omega_k, {\hat x} + z \hat\xi \rangle}
    \right)_{{\hat x}=0},$$ where ${\left.\operatorname{Res}\vphantom{\big|}\right|_{z=0}} g = a_{-1}$ is the residue of a formal Laurent series $g = \sum_k a_k z^k$. (The residue appears as part of an inversion formula for the Laplace transform.)

In the case where only a minimal number of weights lie on the linear hyperplane spanned by $\hat W$ ($m=r-1$), the wall polynomial $\hat f_{\hat W}$ can be chosen as a constant, since the corresponding push-forward map is merely a change of coordinates:

\[minimal wall jump cones\] Suppose that precisely $r-1$ weights $\hat\omega_1, \ldots, \hat\omega_{r-1}$ lie on $\operatorname{span}{\hat W}$. Then, $$\hat f_{\hat W}^{-1} \equiv |d\lambda\left(\hat\omega_1, \ldots, \hat\omega_{r-1}, \frac {\hat\xi} {{\lVert\hat\xi\rVert}^2}\right)|.$$

Since the map $\hat P_{\hat W} \colon {\mathbb R}^{r-1} \rightarrow \operatorname{span}\hat W, (u_k) \mapsto \sum_{k=1}^{r-1} u_k \hat\omega_k$ along which we push forward is a linear isomorphism, the polynomial $\hat f_{\hat W}$ can be chosen as the constant of proportionality between the push-forward of Lebesgue measure on ${\mathbb R}^{r-1}$ and the measure $dw$. We can compute its value by comparing the volume of the parallelotope spanned by the $(\hat\omega_k)$ with respect to the two measure. For the former measure, this is of course one, while for the latter it follows from that $$d\lambda(\hat\omega_1, \ldots, \hat\omega_{r-1}, \hat\xi) =
    dw(\hat\omega_1, \ldots, \hat\omega_{r-1}) \,{\lVert\hat\xi\rVert}^2.
    \qedhere$$

This immediately gives rise to the following inductive algorithm:

\[boysal vergne algorithm\] The following algorithm computes the piecewise polynomial density of the push-forward of Lebesgue measure on ${\mathbb R}^n_{\geq 0}$ along $(t_k) \mapsto \sum_{k=1}^n t_k \hat\omega_k$:

1.  Start with the unbounded chamber, where $\hat f \equiv 0$.

2.  Iteratively jump over walls $\hat W$ separating the current chamber with an adjacent chamber:

    1.  Denote by $\hat\omega_1, \ldots, \hat\omega_m$ the weights which lie on the hyperplane through $\hat W$.

    2.  If the wall is minimal ($m=r-1$), compute $\hat f_{\hat W}$ via .

    3.  Otherwise, recursively apply to compute the piecewise polynomial density of the push-forward of Lebesgue measure on ${\mathbb R}^m_{\geq 0}$ along $(u_k) \mapsto \sum_{k=1}^m u_k \hat\omega_k$.[^3] On $\hat W$ itself, it is given by a single homogeneous polynomial. Choose any polynomial extension $\hat f_{\hat W}$ to all of $\mathfrak t^*$.

    4.  Compute the density on the adjacent chamber using .

By combining with the Heckman formula, we arrive at the following algorithm for computing Duistermaat–Heckman measures. We shall call it the *(Abelian) Heckman algorithm*.

\[abelian heckman algorithm\] Under the assumptions and using the notation of , the following algorithm computes the piecewise polynomial density function of the Abelian Duistermaat–Heckman measure:

1.  Compute the density of each of the $|M^T|$ iterated convolutions $\delta_{\Phi_K(p)} \star H_{\hat\omega_{p,1}} \star \ldots \star H_{\hat\omega_{p,n}}$ using .

2.  Form their alternating sum according to .

The non-Abelian Duistermaat–Heckman measure can then be computed via . By passing to its support, we can also determine the non-Abelian moment polytope (cf. ).

The algorithm as we have stated it assumes that the fixed-point data is part of the input. Let us describe it in the situations we are interested in:

\[fixed-point data projective space\] Consider the projective space $M = {\mathbb P}(V)$ associated with an arbitrary finite-dimensional, unitary $K$-representation $V$. Torus fixed points in $M$ correspond to weight vectors in $V$. Therefore, $M^T$ is finite if and only if all the weight spaces of $V$ are one-dimensional. If this is the case, let $V = \bigoplus_{k=0}^n {\mathbb C}v_k$ be the weight-space decomposition, with $v_k$ weight vectors of pairwise distinct weight $\omega_k$, so that the torus fixed points are precisely the points $[v_0], \ldots, [v_n] \in M$. Then, *before renormalization*, the isotropy weights in $T_{[v_k]} M$ are given by the vectors $\omega_l - \omega_k$ for $l \neq k$.

Note that the representations associated with the pure-state quantum marginal problems displayed in indeed have one-dimensional weight spaces, so that is directly applicable: This is obvious for ${\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N}$ and can also be verified for $\operatorname{Sym}^N({\mathbb C}^d)$ and ${\Lambda}^N({\mathbb C}^d)$ (e.g., by observing that any single-row or single-column semistandard tableaux is already determined by its weight vector). However, other irreducible representations of $\operatorname{SU}(d)$, which correspond to indistinguishable particles of more exotic statistics, typically have weight spaces of dimension larger than one [@fulton97].

\[heckman for coadjoint orbit reductions\] Consider more generally the action of $T$ on a coadjoint $\tilde K$-orbit $M = \mathcal O_{\tilde\lambda}$ induced by a group homomorphism $\varphi \colon T \rightarrow \tilde T \subseteq \tilde K$. Even though this action might have infinitely many fixed points, there is an obvious way to write down an alternating sum formula for $\operatorname{DH}^T_{\mathcal O_{\tilde\lambda}}$: Note that it follows directly from that $$\operatorname{DH}^T_{\mathcal O_{\tilde\lambda}} =
    \pi_* \operatorname{DH}^{\tilde T}_{\mathcal O_{\tilde\lambda}},$$ where $\pi = (d\varphi)^*$ is the dual map $\mathfrak {\tilde t}^* \rightarrow \mathfrak t^*$. Therefore, we can simply take the Abelian Heckman formula for the $\tilde T$-action (which is always applicable since the fixed point set of $\tilde T$ is the Weyl orbit of $\tilde\lambda$, hence finite), and push forward each summand along $\pi$. In the case of a maximal-dimensional coadjoint orbit and for a suitable choice of renormalization direction, the result is just the push-forward of the Harish-Chandra formula , $$\label{heckman for maximal dimensional coadjoint orbit reductions}
    \operatorname{DH}^T_{\mathcal O_{\tilde\lambda}} = \sum_{\tilde w \in \tilde W} (-1)^{l(\tilde w)} \delta_{\pi(\tilde w \tilde\lambda)} \star H_{-\pi(\tilde \alpha_1)} \star \ldots \star H_{-\pi(\tilde \alpha_{\tilde R})},$$ with $\tilde\alpha_1, \ldots, \tilde\alpha_{\tilde R}$ the positive roots of $\tilde K$. The formula for lower-dimensional coadjoint orbits can be obtained by using [@berlinegetzlervergne03 Theorem 7.24] instead of .

In particular, this approach allows the computation of the Abelian Duistermaat–Heckman measure for arbitrary setups of the quantum marginal problem by an obvious variant of .

While and the variant described in are useful for computing low-dimensional examples, any approach relying on the Heckman formula has the major problem that the number of summands in the Heckman formula is typically very large (e.g., it is exponential in the number of distinguishable particles or fermions). Moreover, even though the Boysal–Vergne algorithm computes the density of a single summand chamber-by-chamber, this is less straightforward for the alternating sum, where all summands have to be evaluated in parallel. In below we will therefore derive an algorithm which does not suffer from these problems.

There is also a non-Abelian Heckman formula due to Guillemin and Prato [@guilleminprato90] (which suffers from the same problems). It can be deduced directly from the Abelian one by applying the derivative principle:

\[non-abelian heckman\] Suppose that there are only finitely many torus fixed points $p \in M^T$ and that in each tangent space $T_p M$ each positive root $\alpha > 0$ or its negative occurs as an isotropy weight. Denote by $n_p$ the number of isotropy weights in $T_p M$ that are multiplied by $-1$ during renormalization. For each positive root $\alpha > 0$ and in each $T_p M$, remove either $\alpha$ or $-\alpha$ from the list of renormalized isotropy weights. Denote the remaining weights by $\hat\omega_{p,1}, \ldots, \hat\omega_{p,n-R}$, and let $k_p$ be the number of negative roots that have been removed. Then, $$\operatorname{DH}^K_M = \sum_{p \in M^T} (-1)^{n_p+k_p} \delta_{\Phi_K(p)} \star H_{\hat\omega_{p,1}} \star \ldots \star H_{\hat\omega_{p,n-R}} \Big|_{\mathfrak t^*_+}.$$ In particular, the second assumption is satisfied when the moment map $\Phi_K$ sends each torus fixed points to the interior of a Weyl chamber.

Since $\partial_{\hat\omega} H_{\pm{\hat\omega}} = \pm\delta_0$ (cf. the proof of ), the asserted formula follows at once by combining with .

Only the final remark needs elaboration: As observed by Guillemin and Prato, the assumption that $\Phi_K(p) \in W \cdot \mathfrak t^*_{>0}$ implies that the $K$-stabilizer at each fixed point $p$ is precisely $T$, so that the infinitesimal action of $K$ generates a copy of $\mathfrak k / \mathfrak t$ inside the tangent space $T_p M$. Therefore, at any fixed point $p$, each positive root $\alpha > 0$ or its negative occurs as an isotropy weight.

This gives rise to an obvious non-Abelian variant of :

\[non-abelian heckman algorithm\] Under the assumptions and using the notation of , the following algorithm computes the non-Abelian Duistermaat–Heckman measure:

1.  Compute the $|M^T|$ iterated convolutions $\delta_{\Phi_K(p)} \star H_{\hat\omega_{p,1}} \star \ldots \star H_{\hat\omega_{p,n-R}}$ using (see ).

2.  Form their alternating sum according to .

By passing to its support, we can also determine the non-Abelian moment polytope (cf. ).

\[subtle remark\] There is a slight subtlety involved with the formulation of step (1) of : In case the renormalized isotropy weights $\hat\omega_{p,1}, \ldots, \hat\omega_{p,n-R}$ in some $T_p M$ do not span all of $\mathfrak t^*$, the corresponding iterated convolution is of course not absolutely continuous with respect to $d\lambda$, and cannot be applied directly (see, e.g., the first proof of ). Instead, we need to replace $\mathfrak t^*$ by the span of the $\hat\omega_{p,k}$ and apply accordingly.

In we will use both the Abelian and the non-Abelian version of the Heckman algorithm to compute the eigenvalue distribution of the reduced density matrices of a random pure state of two qubits () and of $N$ bosonic qubits (), as well as of random mixed states of two qubits ().

Single-Summand Algorithm for Projective Space {#projective space}
---------------------------------------------

We will now derive explicit formulas for the Duistermaat–Heckman measure associated with a projective space, $M = {\mathbb P}(V)$, where $V$ is a $(n+1)$-dimensional unitary representation of $K$, and where $M$ is equipped with the Fubini–Study symplectic form $\omega_\text{FS}$, normalized in such a way that its Liouville measure is equal to $\frac 1 {n!}$. The $K$-action is Hamiltonian, and a canonical moment map is given by [@kirwan84] $$\label{projective space non-abelian moment map}
  \Phi_K \colon
  {\mathbb P}(V) \rightarrow \mathfrak k^*, \quad
  [v] \mapsto \left( X \mapsto \frac 1 i  \frac {\langle v, X v \rangle} {\langle v, v \rangle} \right).$$

We start by decomposing the representation $V$ into one-dimensional weight spaces, $V = \bigoplus_{k=0}^n {\mathbb C}v_k$, where $v_k$ is a weight vector of weight $\omega_k$ (repetitions allowed). In the corresponding homogeneous coordinates, the Abelian moment map has the following simple form, $$\label{projective space abelian moment map}
  \Phi_T \colon {\mathbb P}(V) \rightarrow \mathfrak t^*, \quad
  [z_0 : \ldots : z_n] \mapsto \frac {\sum_{k=0}^n |z_k|^2 \omega_k} {\sum_{k=0}^n |z_k|^2},$$ and it is straightforward to see that the Abelian Duistermaat–Heckman measure can be written as the push-forward of Lebesgue measure on the standard simplex along a linear map:

\[projective space abelian via standard simplex\] We have $$\operatorname{DH}^T_{{\mathbb P}(V)} = P_*({\left.dp\vphantom{\big|}\right|_{\Delta_n}}).$$ Here, $P$ is the linear map ${\mathbb R}^{n+1} \rightarrow \mathfrak t^*, (t_k) \mapsto \sum_k t_k \omega_k$, and $dp$ is Lebesgue measure on the affine hyperplane ${\mathbf H}:= \{ (t_k) : \sum_k t_k = 1 \} \subseteq {\mathbb R}^{n+1}$, normalized in such a way that the standard simplex $\Delta_n :=   \{ (p_k) : p_k \geq 0, \sum_{k=0}^n p_k = 1 \}$ has measure $\frac 1 {n!}$.

The Fubini-Study measure is the push-forward of the usual round measure on the unit sphere $S^{2n+1} \cong \{ (z_0,\ldots,z_n) : {\lvertz_0\rvert}^2 + \ldots + {\lvertz_n\rvert}^2 = 1 \} \subseteq V$ along the quotient map $(z_0,\ldots,z_n) \mapsto [z_0:\ldots:z_n]$, normalized to total volume $\frac 1 {n!}$. On the other hand, the round measure on the unit sphere also induces Lebesgue measure on the standard simplex by pushing forward along the map $(z_0,\ldots,z_n) \mapsto (|z_0|^2,\ldots,|z_n|^2)$, since $dx dy = d(r^2) d\theta$ in polar coordinates. The claim follows from comparing with $P \colon (t_k) \mapsto \sum_{k=0}^n t_k \omega_k$.

can also be established by applying the Heckman formula as described in .

Denote by $dp/d\lambda$ a differential form corresponding to Lebesgue measure on the affine subspaces $P^{-1}(\lambda) \cap {\mathbf H}$, normalized in such a way that $$\label{quotient form polytope}
  dp = dp/d\lambda \wedge P^*(d\lambda)$$ when restricted to the affine hyperplane ${\mathbf H}$.

\[density in polytope picture\] The density function $f \colon \mathfrak t^* \rightarrow [0,\infty)$ of the Abelian Duistermaat–Heckman measure is given by $$f(\lambda) = \operatorname{vol}\, \{ p_k \geq 0 : \sum_{k=0}^n p_k \omega_k = \lambda, \sum_{k=0}^n p_k = 1 \},$$ where the volume is measured with respect to the measure induced by $dp/d\lambda$ on $P^{-1}(\lambda) \cap {\mathbf H}$.

For all test functions $g \in C_b(\mathfrak t^*)$, we have $$\langle \operatorname{DH}^T_{{\mathbb P}(V)}, g \rangle =
    \int_{\Delta_n} dp \, g(P(p)) =
    \int_{\mathfrak t^*} d\lambda \left( \int_{P^{-1}(\lambda) \cap \Delta_n} dp/d\lambda \right) g(\lambda),$$ by using and Fubini’s theorem for the fibration ${\left.P\vphantom{\big|}\right|_{{\mathbf H}}}$ [@guilleminsternberg77 pp. 307].

That is, the Abelian Duistermaat–Heckman density measures the volume of a family of convex polytopes parametrized by $\mathfrak t^*$. This is also true for the density of the iterated convolutions studied in (see below). There are exact numerical schemes that can be used to compute the polynomial density functions on each regular chamber which have already been implemented in software packages, e.g., the parametric extension of Barvinok’s algorithm [@barvinok93] described in [@verdoolaegeseghirbeylsetal07; @verdoolaegebruynooghe08]. We will not pursue this route any further. However, in we will show that its “quantized” counterpart gives rise to an efficient way of computing the corresponding representation-theoretic quantities (in particular, the Kronecker coefficients).

In the following, we will instead describe a combinatorial algorithm based on the same principles as our Heckman algorithm. Before doing so, let us determine explicitly the regular chambers for the Abelian moment map, i.e., the connected components of the set of regular values of $\Phi_T$, each on which the measure is given by a polynomial. For this, we define the *support* of a point $p = [v] \in {\mathbb P}(V)$ as the set of weights which contribute to the weight-space decomposition of $v$, $$\operatorname{supp}p := \{ \omega_k : z_k \neq 0, p = [z_0 : \ldots : z_n] \}.$$ The significance of this definition is that the support of a point already fully determines whether it is regular or singular:

\[regular points\] Let $p \in P(V)$. Then $p$ is a regular point of the Abelian moment map if and only if $$\operatorname{span}\{ \omega - \omega' : \omega, \omega' \in \operatorname{supp}p \} = \mathfrak t^*.$$

It follows readily from the definition of the moment map that a point $p$ is regular if and only if $\mathfrak t_p$, the Lie algebra of its stabilizer, is trivial [@guilleminsternberg82 Lemma 2.1]. But $\mathfrak t_p$ is already determined by the support of $p$: $$\mathfrak t_p = \{ X \in \mathfrak t: \omega(X) = \omega'(X) \quad \forall \omega, \omega' \in \operatorname{supp}p \}$$ This is the annihilator of the linear span in the statement of the lemma.

We arrive at the following characterization of the set of singular values of the Abelian moment map:

\[singular values\] The set of singular values of $\Phi_T$ is the union of all convex hulls of subsets containing (at most) $r$ weights, $$\bigcup_{\#I = r} \operatorname{conv}\{ \omega_k : k \in I \} =
    \bigcup_{\#I \leq r} \operatorname{conv}\{ \omega_k : k \in I \}.$$

It is clear from and that the convex hull of any subset of weights of cardinality at most $r$ consists of singular values. The converse follows from Carathéodory’s theorem.

In particular, the singular walls are convex hulls of $r$ weights in general position. From this description we can easily determine the regular chambers. Observe again that there is a single unbounded regular chamber.

We will now use the result of Boysal and Vergne described in to derive intrinsic formulas for the jumps of the Duistermaat–Heckman density when crossing a singular wall. Recall that the measures they consider are push-forwards of Lebesgue measure on the convex cone ${\mathbb R}^{n+1}_{\geq 0}$ rather than of Lebesgue measure on the standard simplex $\Delta_n$, which is of course the intersection of ${\mathbb R}^{n+1}_{\geq 0}$ with the affine hyperplane ${\mathbf H}= \{ (t_k) : \sum_{k=0}^n t_k = 1 \}$. It is however straightforward to translate between both pictures: In order to avoid confusion, we shall use the same convention as in that hatted quantities correspond to the Boysal–Vergne picture. Let us consider the “extended” weights $\hat\omega_k := (\omega_k,1) \in \mathfrak t^* \oplus {\mathbb R}$ ($k=0,\ldots,n$) together with the corresponding linear map $$\hat P \colon {\mathbb R}^{n+1} \rightarrow \mathfrak t^* \oplus {\mathbb R}, \quad (t_k) \mapsto \sum_{k=0}^n t_k \hat\omega_k = (P(t_0, \ldots, t_n), \sum_{k=0}^n t_k).$$ Denote by $dt$ standard Lebesgue measure on ${\mathbb R}^{n+1}$ and equip $\mathfrak t^* \oplus {\mathbb R}$ with the measure $d\hat\lambda = d\lambda ds$, where $ds$ is standard Lebesgue measure on ${\mathbb R}$. Choose a differential form $dt/d\hat\lambda$ inducing Lebesgue measure on the fibers of $\hat P$, normalized in such a way that $$\label{quotient form cone}
  dt = dt/d\hat\lambda \wedge {\hat P}^*(d\hat\lambda) = dt/d\hat\lambda \wedge P^*(d\lambda) \wedge (dt_0 + \ldots + dt_N).$$ Then one can establish just as in the proof of the following formula for the density function of the push-forward of Lebesgue measure on ${\mathbb R}^{n+1}_{\geq 0}$ along $\hat P$ with respect to $d\hat\lambda = d\lambda ds$, $$\label{single summand density}
  \hat f(\lambda, s) = \operatorname{vol}\, \{ t_k \geq 0 : \sum_{k=0}^n t_k \omega_k = \lambda, \sum_{k=0}^n t_k = s \},$$ where the volume is measured with respect to $dt/d\hat\lambda$. But comparing and and noting that $dt = dp \wedge (dt_0 + \ldots + dt_N)$ on ${\mathbf H}$, we see that in fact $dt/d\hat\lambda$ and $dp/d\lambda$ induce the same measure on the fibers $P^{-1}(\lambda) = \hat P^{-1}(\lambda, 1)$, so that $$\label{density transfer eqn}
  \operatorname{DH}^T_{{\mathbb P}(V)} =
  f(\lambda) \, d\lambda =
  \hat f(\lambda, 1) \, d\lambda.$$ This shows that we can work equivalently in the convex cone picture of Boysal and Vergne.[^4]

We shall now describe the jump formula. Let $W$ be a singular wall separating regular chambers $\Delta_\pm \subseteq \mathfrak t^*$, and choose a normal vector $\xi \in \mathfrak t^*$ pointing from $\Delta_-$ to $\Delta_+$. Order the weights such that precisely $\omega_0, \ldots, \omega_{m-1}$ lie on $W$. Denote by $dw$ Lebesgue measure on the hyperplanes parallel to $W$, normalized in such a way that $$\label{wall quotient measure polytope}
  d\lambda = dw \wedge d\xi,$$ where $d\xi$ is the pullback of the standard volume form of ${\mathbb R}$ along the coordinate function $\langle -, \xi\rangle$. Denote by $f_\pm$ the polynomials describing the density function $f$ on the regular chambers $\Delta_\pm$. Finally, consider the Duistermaat–Heckman measure for the action of $T$ on the projective space over $V_W = \bigoplus_{k=0}^{m-1} {\mathbb C}v_k$, the direct sum of the weight spaces corresponding to the weights which lie on the hyperplane through $W$. Its density with respect to $dw$ is given by a single polynomial on the singular wall $W$, since $W$ is always contained in the closure of a regular chamber for ${\mathbb P}(V_W)$. Choose any polynomial function $f_W$ extending it to all of $\mathfrak t^*$.

\[wall jump abelian\] The jump of the Abelian Duistermaat–Heckman density across the singular wall is given by $$f_+(\lambda) -f_-(\lambda)
  = {\left.\operatorname{Res}\vphantom{\big|}\right|_{z=0}} \left(
      \hat f_{\hat W}(\partial_x, \partial_y)
      \frac
        {e^{z \langle \lambda - \omega_0, \xi \rangle + \langle \lambda, x \rangle + y}}
        {\prod_{k=m}^n z \langle \omega_k - \omega_0, \xi \rangle + \langle \omega_k, x \rangle + y}
    \right)_{x=0, y=0}.$$ Here, $\hat f_{\hat W}(\lambda,s) = s^{m-r} f_W(\frac \lambda s)$ is the homogeneous “extension” of $f_W$ to $\mathfrak t^* \oplus {\mathbb R}$.

The convex cones $\hat\Delta_\pm$ through $\Delta_\pm \times \{1\}$ are chambers in the sense of Boysal and Vergne. They are separated by a wall $\hat W$, namely the convex cone through $W \times \{1\}$. Note that $\hat\xi = (\xi,-\braket{\omega_0,\xi})$ is a normal vector to $\hat W$. Denote by $\hat f_\pm$ the homogeneous polynomials describing the density function of the push-forward of Lebesgue measure on ${\mathbb R}^{n+1}_{\geq 0}$ along $\hat P$. It is clear that $d\hat w = ds \wedge dw$ induces Lebesgue measure on $\hat W$ and that it is normalized in such a way that $d\hat\lambda = d\hat w \wedge d\hat\xi$. By and the jump formula of Boysal and Vergne, we have $$\begin{aligned}
    f_+(\lambda) - f_-(\lambda)
  = \hat f_+(\lambda, 1) - \hat f_-(\lambda, 1)
  = {\left.\operatorname{Res}\vphantom{\big|}\right|_{z=0}} \left(
      \hat f_{\hat W}(\partial_{\hat x})
      \frac
        {e^{z \langle (\lambda, 1), \hat x + z \hat \xi \rangle}}
        {\prod_{k=m}^n \langle \hat\omega_k, \hat x + z \hat\xi \rangle}
    \right)_{\hat x=0}.
  \end{aligned}$$ The polynomial $\hat f_{\hat W}$ as defined above agrees with its original definition in , since it is a homogeneous polynomial and can thus be reconstructed from $f_W$, which by is its restriction to the slice $\mathfrak t^* \times \{1\}$, by the formula given above. Writing $\hat x = (x,y) \in \mathfrak t^* \oplus {\mathbb R}$ and expanding the hatted quantities, we arrive at the assertion.

As in , the case where only a minimal number of weights lie on the affine hyperplane through $W$ is particularly simple to evaluate:

\[minimal wall jump constant polytope\] Suppose that precisely $r$ weights $\omega_0, \ldots, \omega_{r-1}$ lie on the affine hyperplane through $W$. Then, $$\hat f_{\hat W}^{-1} \equiv f_W^{-1} \equiv |d\lambda\left(\omega_1-\omega_0, \ldots, \omega_{r-1}-\omega_0, \frac{\xi}{{\lVert\xi\rVert}^2}\right)|.$$

We argue as in the proof of : In view of and the minimality assumption, the map $(q_k) \mapsto \sum_{k=0}^{r-1} q_k \omega_k$ along which we push forward is an isomorphism, and $f_W$ is equal to the constant of proportionality between the push-forward of Lebesgue measure on ${\mathbf H}$ (normalized in such a way that the standard simplex has measure $\tfrac 1 {d!}$) and the measure $dw$. We can compute this constant by comparing the volume of the parallelotope spanned by the $(\omega_k)$: For the former measure this constant is one (by its very normalization), while for the latter it follows from that $$d\lambda(\omega_1-\omega_0, \ldots, \omega_{r-1}-\omega_0, \xi) = d\omega(\omega_1-\omega_0, \ldots, \omega_{r-1}-\omega_0) \,{\lVert\xi\rVert}^2. \qedhere$$

These results give rise to the following inductive algorithm for computing the Abelian and non-Abelian Duistermaat–Heckman measure of a projective space. We will call it the *single-summand algorithm*, since in view of it amounts to computing a push-forward measure that is equivalent to a single summand of the Abelian Heckman formula (cf. ).

\[projective space algorithm\] The following algorithm computes the piecewise polynomial density function of the Abelian Duistermaat–Heckman measure of the projective space ${\mathbb P}(V)$:

1.  Start with the unbounded regular chamber, where $f \equiv 0$.

2.  Iteratively jump over singular walls $W$ separating the current regular chamber with an adjacent regular chamber:

    1.  Denote by $\omega_0, \ldots, \omega_{m-1}$ the weights which lie on the hyperplane through $W$.

    2.  If the wall is minimal ($m=r$), compute $f_W$ via .

    3.  Otherwise, recursively apply to compute the piecewise polynomial density of the Abelian Duistermaat–Heckman measure of ${\mathbb P}(V_W)$, where $V_W = \bigoplus_{k=0}^{m-1} {\mathbb C}v_k$ is the direct sum of the weight spaces for the weights in (a).[^5] On $W$ itself, it is given by a single polynomial. Choose any polynomial extension $f_W$ to all of $\mathfrak t^*$.

    4.  Compute the density on the adjacent chamber using .

The non-Abelian Duistermaat–Heckman measure can then be computed via . By passing to its support, we can also determine the non-Abelian moment polytope (cf. ).

\[complete solution\] In view of and and by passing to the purified double (), solves the problem of computing the eigenvalue distribution of reduced density matrices in complete generality.

We conclude this section by explicitly stating the Abelian and non-Abelian jump formula for the case where only a minimal number of weights lie on the affine hyperplane through the wall. They will be used later for computing examples.

\[minimal wall jump abelian\] Suppose that precisely $r$ weights $\omega_0, \ldots, \omega_{r-1}$ lie on the affine hyperplane through the singular wall $W$. Then the jump of the Abelian Duistermaat–Heckman density across the wall is given by $$f_+(\lambda) - f_-(\lambda) =
    f_W
    \left( \prod_{k=r}^n \langle \omega_k - \omega_0, \xi \rangle \right)^{-1}
    \frac{\langle \lambda - \omega_0, \xi \rangle^{n-r}}{(n-r)!},$$ where $f_W$ is the constant from .

This follows immediately from by pulling out the constant $f_W$, setting $x = y = 0$ and evaluating the residue at $z=0$.

The non-Abelian formula follows directly by applying :

\[minimal wall jump non-abelian\] Suppose that precisely $r$ weights $\omega_0, \ldots, \omega_{r-1}$ lie on the affine hyperplane through the singular wall $W$, and that $n-r \geq R$, so that the non-Abelian Duistermaat–Heckman measure of ${\mathbb P}(V)$ is absolutely continuous in the vicinity of $W$. Denote by $f^K_\pm$ the polynomials describing its density on the regular chambers. Then the jump across the wall is given by $$\begin{aligned}
    f^K_+(\lambda) - f^K_-(\lambda)
    = f_W
    \left( \prod_{k=r}^n \langle \omega_k - \omega_0, \xi \rangle \right)^{-1}
    \left( \prod_{\alpha > 0}  - \langle \alpha, \xi \rangle \right)
    \frac{\langle \lambda - \omega_0, \xi \rangle^{n-r-R}}{(n-r-R)!},
  \end{aligned}$$ where $f_W$ is the constant from .

has already been established in [@guilleminlermansternberg88], where the authors also envisaged an algorithm similar to our Heckman algorithm. They did however not have a general jump formula such as at their avail. Instead, they had to resort to an inexact formula which in general only holds in highest order (in the distance to the wall).

Examples
========

In this section we illustrate our algorithms by computing some eigenvalue distributions of reduced density matrices. The global quantum states will always be chosen according to one of the invariant probability measures described in . Many of our examples will involve *qubits*, i.e., quantum systems modeled by two-dimensional Hilbert spaces, so that the algorithms can be nicely visualized. But of course our algorithms can be used to determine the eigenvalue distributions for arbitrary instances of the quantum marginal problem (see ).

Pure States of Multiple Qubits {#pure states qubits}
------------------------------

We start by considering pure states of $N$ qubits, where $K = \operatorname{SU}(2)^N$ acts on $M = {\mathbb P}(({\mathbb C}^2)^{\otimes N})$ by tensor products (cf. ). It will be convenient to identify $\mathfrak t^* \cong {\mathbb R}^N$ in such a way that the positive Weyl chamber corresponds to the cone ${\mathbb R}^N_{\geq 0}$ and the fundamental weights to the standard basis vectors $e_j = (\delta_{j,k})$ ($k=1,\ldots,N$). That is, if $\lambda = (\lambda^{(j)}) \in \mathfrak t^*$ then we will by slight abuse of notation identify $\lambda^{(j)}$ with the scalar $i(\lambda^{(j)}_1 - \lambda^{(j)}_2) =$ $2 i \lambda^{(j)}_1$. It follows that $d\lambda$ is simply the usual Lebesgue measure on ${\mathbb R}^N$, that the symplectic volume polynomial is given by $p_K(\lambda) = \lambda^{(1)} \cdots \lambda^{(N)}$, and that the positive roots are $2 e_1, \ldots, 2 e_N$ (cf.  and ). Moreover, amounts to assigning to a point $(\lambda^{(j)}) \in {\mathbb R}^N$ the tuple $(\rho_1,\ldots,\rho_N)$ of diagonal density matrices acting on ${\mathbb C}^2$, where $\rho_j$ has maximal eigenvalue $\hat\lambda^{(j)}_{\max} = \frac 1 2 + i \lambda^{(j)}_1 = \tfrac {1+\lambda^{(j)}} 2$.

We first discuss in detail the toy example of $N=2$ qubits, demonstrating both the non-Abelian Heckman algorithm and the single-summand algorithm.

$$\begin{tikzpicture}
      \draw[help lines] (-1.3,-1.3) grid (1.3,1.3);
      \draw[->] (-1.6,0) -- (1.6,0) node[right] {\tiny $\lambda^{(1)}$};
      \draw[->] (0,-1.6) -- (0,1.6) node[left] {\tiny $\lambda^{(2)}$};
      \draw[fill=gray,opacity=0.3] (1,1) -- (-1,1) -- (-1,-1) -- (1,-1) -- (1,1);
      \draw[very thick] (1,1) -- (-1,1);
      \draw[very thick] (1,-1) -- (-1,-1);
      \draw[very thick] (1,-1) -- (1,1);
      \draw[very thick] (-1,-1) -- (-1,1);
      \draw[very thick] (1,1) -- (-1,-1);
      \draw[very thick] (-1,1) -- (1,-1);
      \draw(1,1) node[below right] {\tiny $(1,1)$};
      \draw(1,-1) node[below right] {\tiny $(1,-1)$};
      \draw(-1,1) node[below left] {\tiny $(-1,1)$};
      \draw(-1,-1) node[below left] {\tiny $(-1,-1)$};
    \end{tikzpicture}$$

\[two qubits non-abelian\] The non-Abelian Duistermaat–Heckman measure for the action of $\operatorname{SU}(2) \times \operatorname{SU}(2)$ on ${\mathbb P}({\mathbb C}^2 \otimes {\mathbb C}^2)$ is given by $$\langle \operatorname{DH}^{\operatorname{SU}(2) \times \operatorname{SU}(2)}_{{\mathbb P}({\mathbb C}^2 \otimes {\mathbb C}^2)}, f \rangle =
    \frac 1 2 \int_0^1 f(t,t) dt,$$ i.e., by a one-dimensional Lebesgue measure supported on the diagonal between the origin and $(1,1)$.

The four fixed points of the action correspond to the standard basis vectors $e_j \otimes e_k$ ($j,k=1,2$), which are weight vectors of weight $(\pm 1, \pm 1)$ using the conventions fixed above (the vertices of the grey rectangle in ). Let us choose the direction $\gamma = (-2,-1)$ for renormalization. After removal of the positive and negative roots, $(\pm 2,0)$ and $(0,\pm 2)$, only a single renormalized isotropy weight remains at each fixed point (cf. ). Therefore, shows that the non-Abelian Duistermaat–Heckman measure is given by the restriction to the positive Weyl chamber of $$\delta_{(1,1)} \star H_{(-2,-2)}
    - \delta_{(1,-1)} \star H_{(-2,2)}
    + \delta_{(-1,1)} \star H_{(-2,2)}
    - \delta_{(-1,-1)} \star H_{(-2,-2)}.$$ Only the first summand contributes to the positive Weyl chamber, and its restriction is given precisely by the formula displayed above (cf. ).

Note that the non-Abelian wall jump formula () is not directly applicable, since $n-r \not\geq R$. Indeed, as we have seen above, the non-Abelian measure does not have a Lebesgue density, since it is concentrated on the diagonal.

Therefore, we will follow , which uses the Abelian wall jump formula, and afterwards takes partial derivatives in direction of the negative roots according to : The decomposition of $\mathfrak t^*$ into regular chambers is indicated in . We start in the unbounded chamber, where the density is equal to the zero polynomial and cross the horizontal singular wall at the top. Evaluating the Abelian jump formula (; say, with $\omega_0 = (1,1)$ and $\xi = (0,-1)$), we find that the density polynomial on the upper regular chamber is equal to $\frac 1 8 (1 - \lambda^{(2)})$. Next, we cross the diagonal singular wall separating the upper and the right-hand side regular chamber. Using the Abelian jump formula once again, we see that the density polynomial changes by $\frac 1 8 (\lambda^{(2)} - \lambda^{(1)})$. Therefore, the Abelian Duistermaat–Heckman measure has the following piecewise polynomial density on the positive Weyl chamber: $$\frac 1 8 \left( 1 - \max(\lambda^{(1)}, \lambda^{(2)}) \right)$$ Taking partial derivatives in the direction of the negative roots, $(-2,0)$ and $(0,-2)$, we arrive at the measure asserted above.

\[two qubits marginals\] The joint distribution ${\mathbf P}_{\operatorname{eig}}$ of the maximal eigenvalues of the reduced density matrices of a randomly-chosen pure quantum state of two qubits is given by $$\langle {\mathbf P}_{\operatorname{eig}}, f \rangle = 24 \int_{\frac 1 2}^1 f(s, s) \left( s - \frac 1 2 \right)^2 ds,$$ for all test functions $f(\hat\lambda^{(1)}_{\max}, \hat\lambda^{(2)}_{\max})$.

According to , multiply the non-Abelian Duistermaat–Heckman measure by the symplectic volume polynomial $p_{\operatorname{SU}(2) \times \operatorname{SU}(2)}(\lambda) = \lambda^{(1)} \lambda^{(2)}$, divide by $\frac 1 {3!}$, the volume of ${\mathbb P}({\mathbb C}^2 \otimes {\mathbb C}^2)$. Finally, push forward along $(\lambda^{(j)}) \mapsto (\hat\lambda^{(j)}_{\max} = \tfrac {1+\lambda^{(j)}} 2)$.

This eigenvalue distribution is in fact known more generally for bipartite pure states chosen at random [@lloydpagels88; @zyczkowskisommers01]. We will later show how to compute its generalization using the techniques of this paper ().

For higher tensor powers, evaluating the Heckman formula quickly becomes unwieldy. However, it can still be used to compute the Duistermaat–Heckman measure locally:

\[n qubits local\] The non-Abelian Duistermaat–Heckman measure for the action of $\operatorname{SU}(2)^N$ on ${\mathbb P}(({\mathbb C}^2)^{\otimes N})$ is on the closures of the regular chambers that contain the vertex $(1,\ldots,1)$ given by the convolution product $$\delta_{(1,\ldots,1)} \star H_{\omega_1} \star \ldots \star H_{\omega_{2^N-N-1}},$$ where $\{ \omega_k \}$ is the set of weights of the form $(-2,\ldots,-2,0,\ldots,0)$ (at least two non-zero entries) as well as their $S_N$-permutations.

If we renormalize with respect to the direction $\gamma \approx (-1,\ldots,-1)$ then just as in the first proof of only a single summand in the non-Abelian Heckman formula contributes in the vicinity of the vertex $(1,\ldots,1)$ and, moreover, this summand is of the above form: Indeed, the weights $\{ \omega_k \}$ are precisely the isotropy weights with the negative roots removed (cf. ). Since the density function of $\operatorname{DH}^K_M$ is polynomial on each regular chamber adjacent to the vertex, we can extend the local formula to their closures.

It is in fact easy to see that the domain of validity of this formula is the intersection of the half-space $$\left\{ \lambda : \sum_{j=1}^N \lambda^{(j)} \geq N-2 \right\}$$ with the positive Weyl chamber (the regular chambers not adjacent to $(1,\ldots,1)$ lie in the complement of this half-space).

gives a local description of the non-Abelian moment polytope, namely by the cone based at $(1,\ldots,1)$ and spanned by the rays with direction vectors $\{\omega_k\}$. By convexity, its intersection with the positive Weyl chamber is an outer approximation to the moment polytope.

Let us specialize to the case $N=3$: Here, precisely the rays with the direction vectors $(-2,-2,0)$, $(-2,0,-2)$ and $(0,-2,-2)$ are extremal. Their intersection with the positive Weyl chamber has to be contained in the non-Abelian moment polytope: Otherwise, there would be additional vertices in the interior of the positive Weyl chamber — but only $(1,1,1)$ is the image of a torus fixed point. Since also the origin is contained in the moment polytope (the Greenberger–Horne–Zeilinger state, $[\psi] = [e_1 \otimes e_1 \otimes e_1 + e_2 \otimes e_2 \otimes e_2]$, is a preimage of the origin [@greenbergerhornezeilinger89]), we conclude that the convex hull $$\operatorname{conv}\{ (0,0,0), (1,0,0), (0,1,0), (0,0,1), (1,1,1) \}$$ is an inner approximation to the moment polytope. Both approximations are in fact equal and therefore describe the moment polytope precisely (). Inequalities characterizing the moment polytope for $N$ qubits have been determined in [@higuchisudberyszulc03].

\(O) at (0,0,0); (W) at (0.3333,0.3333,0.3333); (W1) at (1,0,0); (W2) at (0,1,0); (W3) at (0,0,1); (SEP) at (1,1,1);

(-0.3,0,0) – (1.3,0,0); (0,-0.3,0) – (0,1.3,0); (0,0,-0.3) – (0,0,1.3); (1,-0.15,0) – (1,1.15,0); (-0.15,1,0) – (1.15,1,0); (1,1,-0.15) – (1,1,1.15);

(-0.15,1,1) – (1.15,1,1); (0,-0.15,1) – (0,1.15,1); (0,1,-0.15) – (0,1,1.15);

\(W) – (W1); (W) – (W2); (W) – (W3); (W) – (O); (W) – (W1); (W) – (W2); (W) – (W3); (W) – (SEP); (O) – (W1) – (W2) – cycle; (SEP) – (W1) – (W2) – cycle; (SEP) – (W2) – (W3) – cycle; (O) – (W3) – (W1) – cycle; (SEP) – (W3) – (W1) – cycle;

(W1) node\[below\][$(1,0,0)$]{}; (W2) node\[left\][$(0,1,0)$]{}; (W3) node\[left\][$(0,0,1)$]{}; (SEP) node\[anchor=west\][$(1,1,1)$]{}; (O) node\[below left\][$0$]{};

\[three qubits non-abelian density\] The non-Abelian Duistermaat–Heckman measure for the action of $\operatorname{SU}(2)^3$ on ${\mathbb P}({\mathbb C}^2 \otimes {\mathbb C}^2 \otimes {\mathbb C}^2)$ has the piecewise linear Lebesgue density $$\begin{cases}
      \frac{1}{16} \min \lambda^{(j)} & \text{in the lower pyramid},\\
      \frac{1}{32} \left( 1 - \sum_{j=1}^3 \lambda^{(j)} + 2 \min \lambda^{(j)} \right) & \text{in the upper pyramid},\\
      0 & \text{otherwise}
    \end{cases}$$ (compare ).

By , the non-Abelian Duistermaat–Heckman measure is on the closures of the regular chambers containing $(1,1,1)$ given by the convolution $$\delta_{(1,1,1)} \star H_{(-2,-2,-2)} \star H_{(-2,-2,0)} \star H_{(-2,0,-2)} \star H_{(0,-2,-2)}.$$ Using we can readily compute its density: $$\begin{aligned}
     &\int_0^\infty dt_1 \cdots \int_0^\infty dt_4 \, \delta\big({\left(\begin{smallmatrix}1\\1\\1\end{smallmatrix}\right)} + t_1 {\left(\begin{smallmatrix}-2\\-2\\-2\end{smallmatrix}\right)} + t_2 {\left(\begin{smallmatrix}-2\\-2\\0\end{smallmatrix}\right)} + t_3 {\left(\begin{smallmatrix}-2\\0\\-2\end{smallmatrix}\right)} + t_4 {\left(\begin{smallmatrix}0\\-2\\-2\end{smallmatrix}\right)} - \lambda\big)\\
     = &\frac{1}{32} \int_0^\infty ds_1 \int_{-\infty}^\infty ds_2 \cdots \int_{-\infty}^\infty ds_4 \, {\mathbf 1}_C(s_2,s_3,s_4) \, \delta\big((1-s_1) {\left(\begin{smallmatrix}1\\1\\1\end{smallmatrix}\right)} - {\left(\begin{smallmatrix}s_2\\s_3\\s_4\end{smallmatrix}\right)} - \lambda\big)\\
     = &\frac{1}{32} \int_0^\infty ds_1 \, {\mathbf 1}_C\big((1-s_1) {\left(\begin{smallmatrix}1\\1\\1\end{smallmatrix}\right)} - \lambda\big)\\
          = &\frac{1}{32} \max \{ 0, 1 - \sum_{j=1}^3 \lambda^{(j)} + 2 \min(\lambda^{(j)}) \},
  \end{aligned}$$ where ${\mathbf 1}_C$ is the indicator function of the cone spanned by ${\left(\begin{smallmatrix}1\\1\\0\end{smallmatrix}\right)}$, ${\left(\begin{smallmatrix}1\\0\\1\end{smallmatrix}\right)}$, and ${\left(\begin{smallmatrix}0\\1\\1\end{smallmatrix}\right)}$, i.e., $${\mathbf 1}_C(a,b,c) =
    \begin{cases}
      1 &\text{ if } a+b \geq c \text{ and } a+c \geq b \text{ and } b+c \geq a,\\
      0 &\text{ otherwise}.
    \end{cases}$$ We have therefore established the claimed density on the complement of the lower pyramid.

According to , the jump across the hyperplane separating the upper and the lower pyramid is given by $$\frac 1 4 \left(- \frac 1 {64}\right) 8 \left(1 - \sum_{j=1}^3 \lambda^{(j)} \right) = \frac 1 {32} \left(\sum_{j=1}^3 \lambda^{(j)} - 1 \right)
                                                                                    $$ (the left-hand side terms are ordered just like in the jump formula). This is precisely the difference between the densities on the upper and lower pyramid as asserted in the statement of the proposition.

It is straightforward to deduce from this the eigenvalue distribution (cf. the proof of ):

\[three qubits marginals\] The joint distribution of the maximal eigenvalues of the reduced density matrices of a randomly-chosen pure quantum state of three qubits has Lebesgue density $$8! \left( \prod_{j=1}^3 \hat\lambda^{(j)}_{\max} - \frac 1 2 \right)
    \begin{cases}
      \min \hat\lambda^{(j)}_{\max} - \frac 1 2 & \text{if } \sum_{j=1}^3 \hat\lambda^{(j)}_{\max} \leq 2,\\
      \max \left\{ 0, \frac 1 2 \left( 1 - \sum_{j=1}^3 \hat\lambda^{(j)}_{\max} \right) + \min \hat\lambda^{(j)}_{\max} \right\} & \text{if } \sum_{j=1}^3 \hat\lambda^{(j)}_{\max} \geq 2,
    \end{cases}$$ on the space of maximal eigenvalues $(\hat\lambda^{(j)}_{\max}) \in [\frac 1 2,1]^3$.

Our use of the local convolution formula () and of the non-Abelian wall jump formula () were merely convenient shortcuts: It is clear that we could have completely algorithmically computed the measure by following .

Mixed States of Two Qubits {#bravyi example}
--------------------------

We will now use the non-Abelian Heckman algorithm to treat the case of random two-qubit states with fixed, non-degenerate global eigenvalue spectrum. That is, we consider the action of $K = \operatorname{SU}(2) \times \operatorname{SU}(2)$ on a coadjoint $\operatorname{SU}(4)$-orbit through a point $\tilde\lambda$ contained in the interior of the positive Weyl chamber. Recall that the Weyl group of $\operatorname{SU}(4)$ is the symmetric group $S_4$, with $(-1)^{l(\tilde w)}$ equal to the signum of a permutation $\tilde w \in S_4$. By , $$\label{bravyi abelian}
  \operatorname{DH}^{T}_{\mathcal O_{\tilde\lambda}} = \sum_{\tilde w \in S_4} \operatorname{sign}(\tilde w) \, \delta_{\pi(\tilde w \tilde\lambda)} \star H_{-\pi(\tilde\alpha_1)} \star \ldots \star H_{-\pi(\tilde\alpha_6)},$$ where $\tilde\alpha_1, \ldots, \tilde\alpha_6$ are the positive roots of $\operatorname{SU}(4)$ (see for our conventions), and where $\pi$ is the restriction map $\mathfrak {\tilde t}^* \rightarrow \mathfrak t^*$, with $\mathfrak {\tilde t}^*$ the dual of the Lie algebra of the maximal torus of $\operatorname{SU}(4)$. With respect to our identification $\mathfrak t^* \cong {\mathbb R}^2$ fixed in , the map $\pi$ is given by $$\label{bravyi projection}
  \pi \colon
  \mathfrak {\tilde t}^* \rightarrow {\mathbb R}^2, \quad
  (\tilde\lambda_1,\ldots,\tilde\lambda_4) \mapsto
    2i(\tilde\lambda_1+\tilde\lambda_2, \tilde\lambda_1+\tilde\lambda_3).$$ One computes readily that the $-\pi(\tilde\alpha_k)$ are precisely the weights $(-2,2)$, $(-2,0)$ (twice), $(-2,-2)$ and $(0,-2)$ (twice). In particular, the two negative roots of $\operatorname{SU}(2) \times \operatorname{SU}(2)$ are contained in this list (each of them is in fact contained twice). By applying we arrive at the following formula:

\[bravyi non-abelian\] The non-Abelian Duistermaat–Heckman measure for the action of $\operatorname{SU}(2) \times \operatorname{SU}(2)$ on a coadjoint $\operatorname{SU}(4)$-orbit $\mathcal O_{\tilde\lambda}$ with $\tilde\lambda \in \mathfrak {\tilde t}^*_{>0}$ is given by $$\operatorname{DH}^{\operatorname{SU}(2) \times \operatorname{SU}(2)}_{\mathcal O_{\tilde\lambda}} =
    {\left.
    \left( \sum_{\tilde w \in S_4} \operatorname{sign}(\tilde w) \, \delta_{\pi(\tilde w \tilde\lambda)} \right) \star H_{(-2,2)} \star H_{(-2,0)} \star H_{(-2,-2)} \star H_{(0,-2)}
    \vphantom{\big|}\right|_{\mathfrak t^*_+}}.$$

Following , we evaluate the right-hand side iterated convolution using . The result is the following:

\[bravyi convolution\] The measure $H_{(-2,2)} \star H_{(-2,0)} \star H_{(-2,-2)} \star H_{(0,-2)}$ has Lebesgue density $$f(\lambda^{(1)}, \lambda^{(2)}) =
    \begin{cases}
      0 & \text{in chamber 0},\\
      \frac 1 {64} \left( \lambda^{(1)} + \lambda^{(2)} \right)^2 & \text{in chamber 1},\\
      \frac 1 {64} \left( \left( \lambda^{(1)} \right)^2 + 2 \lambda^{(1)} \lambda^{(2)} - \left( \lambda^{(2)} \right)^2 \right) & \text{in chamber 2},\\
      \frac 1 {32} \left( \lambda^{(1)} \right)^2 & \text{in chamber 3}.
    \end{cases}$$ See for the labelling of the chambers and an illustration of the density.

(-2.3,-2.3) grid (1.3,2.3); (-2.6,0) – (1.6,0) node\[right\] [$\lambda^{(1)}$]{}; (0,-2.6) – (0,2.6) node\[above\] [$\lambda^{(2)}$]{}; (0,0) node\[below right\] [$0$]{}; (1,0) node\[below\] [$1$]{}; (0,1) node\[left\] [$1$]{}; (0,0) – (-2.3,2.3); (0,0) – (-2.3,0); (0,0) – (-2.3,-2.3); (0,0) – (0,-2.3); (0,0) – (-2.3,2.3) – (-2.3,-2.3) – (0,-2.3) – (0,0); (-0.5,1.5) circle (0.15) node [$0$]{}; (-1.5,0.5) circle (0.15) node [$1$]{}; (-1.5,-0.5) circle (0.15) node [$2$]{}; (-0.5,-1.5) circle (0.15) node [$3$]{}; (0,-3) node [ ]{};

The density of the non-Abelian Duistermaat–Heckman measure is thus given by the restriction to the positive Weyl chamber of an alternating sum of $24$ copies of the density described in , one copy attached to each of the points $\pi(\tilde w \tilde\lambda)$. In view of the geometry of the support of the latter density, it is clear that in fact only summands for points in the right halfplane $\{\lambda^{(1)} > 0\}$ contribute (i.e., at most half of the points). Using , one finds that the points $\pi(\tilde w \tilde\lambda)$ are the six points whose coordinates are equal to any two out of the three values $c_1 = |2i(\tilde\lambda_1+\tilde\lambda_4)|$, $c_2 = 2i(\tilde\lambda_1+\tilde\lambda_3)$, or $c_3 = 2i(\tilde\lambda_1+\tilde\lambda_2)$ (without repetitions), as well as their Weyl conjugates. See for illustration.

(-0.3,-0.3) grid (5.3,5.3); (-0.6,0) – (5.6,0) node\[right\] [$\lambda^{(1)}$]{}; (0,-0.6) – (0,5.6) node\[above\] [$\lambda^{(2)}$]{}; (0,0) node\[below left\] [$0$]{}; (5,0) node\[below\] [$1$]{}; (0,5) node\[left\] [$1$]{}; (1,0) node\[below\] [$c_1$]{}; (2,0) node\[below\] [$c_2$]{}; (4,0) node\[below\] [$c_3$]{}; (0,1) node\[left\] [$c_1$]{}; (0,2) node\[left\] [$c_2$]{}; (0,4) node\[left\] [$c_3$]{}; (4,1) – (3,0) – (0,0) – (0,3) – (1,4) – (2,4) – (4,2) – (4,1); (4,1) – (3,0) – (0,0) – (0,3) – (1,4) – (2,4) – (4,2) – (4,1);

(4,2) circle (0.15) node [$+$]{}; (4,1) circle (0.15) node [$-$]{}; (2,1) circle (0.15) node [$+$]{}; (2,4) circle (0.15) node [$-$]{}; (1,4) circle (0.15) node [$+$]{}; (1,2) circle (0.15) node [$-$]{};

Moreover, we can deduce that the non-Abelian moment polytope has the form described in the figure. To do so, we simply need to check in each regular chamber whether the density polynomial vanishes. By doing so and describing the resulting polytope in terms of inequalities, we recover a well-known result by Bravyi [@bravyi04]:

The non-Abelian moment polytope for the action of $\operatorname{SU}(2) \times \operatorname{SU}(2)$ on a generic coadjoint $\operatorname{SU}(4)$-orbit is given by $$\Delta_{\operatorname{SU}(2) \times \operatorname{SU}(2)}(\mathcal O_{\tilde\lambda}) =
    \{
      (\lambda^{(1)},\lambda^{(2)}) :
      0 \leq \lambda^{(1)}, \lambda^{(2)} \leq c_3,
      \lambda^{(1)} + \lambda^{(2)} \leq c_2 + c_3,
      {\lvert\lambda^{(1)} - \lambda^{(2)}\rvert} \leq c_3 - c_1
    \}.$$

In the limit where the global state becomes pure, the moment polytope converges to the diagonal between the origin and $(1,1)$. This is in agreement with . One can similarly recover the eigenvalue distribution of the reduced density matrices of a random pure state of two qubits by taking a corresponding limit.

In view of , the distributions computed in can be assembled to give the joint eigenvalue distribution of the reduced density matrices of a randomly-chosen pure state in ${\mathbb P}({\mathbb C}^2 \otimes {\mathbb C}^2 \otimes {\mathbb C}^4)$.

Pure States of Bosonic Qubits {#pure states bosonic qubits}
-----------------------------

We now turn to random pure states of $N$ bosonic qubits, where $K = \operatorname{SU}(2)$ and $M = {\mathbb P}(\operatorname{Sym}^N({\mathbb C}^2))$. We will use the Abelian Heckman algorithm:

\[sym N qubits abelian\] The Abelian Duistermaat–Heckman measure for ${\mathbb P}(\operatorname{Sym}^N({\mathbb C}^2))$ has Lebesgue density $$\frac 1 {2^N (N-1)! N!} \sum_{k=-N,-N+2,\ldots,N} (-1)^{\frac{N+k}{2}} {\binom{N}{\frac{N+k}{2}}} (\lambda - k)_+^{N-1}.$$ Here, we set $(\lambda-k)^{N-1}_+ = (\lambda-k)^{N-1}$ for $\lambda \geq k$ and $0$ otherwise.

Equivalently, $N!$ times the Abelian Duistermaat–Heckman measure is equal to the probability distribution of the sum of $N$ independent random variables that are uniformly distributed on the interval $[-1,1]$ [@feller71 §I.9, Theorem 1a].

The weights of $\operatorname{Sym}^N({\mathbb C}^2)$ are $\{ -N, -N+2, \ldots, N \}$; let us write $v_k$ for a weight vector of weight $k$. The associated projective space has precisely $N+1$ torus fixed points. At any such fixed point $[v_k]$, the isotropy weights are given by $$\{ (l-k) : l = -N, -N+2, \ldots, \check{k}, \ldots, N \},$$ and we will denote them by $\hat\omega_{k,1,}, \ldots, \hat\omega_{k,N}$ (cf. ). Observe that precisely $n_k = \frac{N+k}{2}$ of them are negative with respect to the renormalization direction $\gamma = +1$. By , the corresponding summand of the Heckman formula is equal to the push-forward of Lebesgue measure on ${\mathbb R}^N_{\geq 0}$ along $$P \colon {\mathbb R}^N_{\geq 0} \rightarrow \mathfrak u_1^*, \quad
    (s_1, \ldots, s_N) \mapsto \sum_{i=1}^N s_i {\lvert\hat\omega_{k,i}\rvert} + k.$$ We first compute its cumulative distribution function: $$\begin{aligned}
    P_*(ds)\left((-\infty,(k+\lambda)]\right)
    &= ds(\{ s_1, \ldots, s_N \geq 0 : \sum_{i=1}^N s_i {\lvert\hat\omega_{k,i}\rvert} \leq \lambda \})\\
    &= ds(\{ s_1, \ldots, s_N \geq 0 : \sum_{i=1}^N s_i \leq 1 \}) \, \frac 1 {\prod_i {\lvert\hat\omega_{k,i}\rvert}} \, \lambda_+^N\\
        &= \frac 1 {N!} \frac 1 {2^N (\frac{N+k}{2})! (\frac{N-k}{2})!} \lambda_+^N\\
    &= \frac 1 {2^N N! N!} {\binom{N}{\frac{N+k}{2}}} \lambda_+^N.
  \end{aligned}$$ The density is then given by the derivative, $$f_k(\lambda) = \frac 1 {2^N (N-1)! N!} {\binom{N}{\frac{N+k}{2}}} (\lambda-k)_+^{N-1},$$ and by forming the alternating sum of these terms we arrive at the formula displayed above.

In we give an alternative proof of using representation theory and combinatorics. It follows directly from the derivative principle that the non-Abelian Duistermaat–Heckman measure is given by the following formula:

\[sym N qubits non-abelian\] The non-Abelian Duistermaat–Heckman measure for the action of $\operatorname{SU}(2)$ on ${\mathbb P}(\operatorname{Sym}^N({\mathbb C}^2))$ with $N \geq 2$ has Lebesgue density $$\frac 1 {2^{N-1} (N-2)! N!} \sum_{k=-N,-N+2,\ldots,N} (-1)^{\frac{N+k}{2}+1} {\binom{N}{\frac{N+k}{2}}} (\lambda-k)_+^{N-2}$$ on $[0,\infty)$.

Again, it is clear how to translate the above into the eigenvalue distribution of the one-body reduced density matrix by using . See for an illustration in the case of $N=2$ bosonic qubits.

\[sym two qubits figure\] $$\begin{aligned}
    \begin{tikzpicture}
      \draw[help lines] (-2.3,-0.3) grid (2.3,1.3);
      \draw[->] (-2.6,0) -- (2.6,0) node[right] {\tiny $\lambda$};
      \draw[->] (0,-0.6) -- (0,1.6);
      \draw(0,0) node[below left] {\tiny $0$};
      \draw(1,0) node[below left] {\tiny $1$};
      \draw(2,0) node[below left] {\tiny $2$};
      \draw(-1.6,0.2) node[above] {\tiny $\frac {2+\lambda} 8$};
      \draw(1.6,0.2) node[above] {\tiny $\frac {2-\lambda} 8$};
      \draw[very thick,dashed] (-2,0) -- (0,1);
      \draw[very thick] (2,0) -- (0,1);
    \end{tikzpicture}
    & \quad \quad
    \begin{tikzpicture}
      \draw[help lines] (-0.3,-0.3) grid (2.3,1.3);
      \draw[->] (-0.6,0) -- (2.6,0) node[right] {\tiny $\lambda$};
      \draw[->] (0,-0.6) -- (0,1.6);
      \draw(0,0) node[below left] {\tiny $0$};
      \draw(1,0) node[below left] {\tiny $1$};
      \draw(2,0) node[below left] {\tiny $2$};
      \draw(2,1) node[right] {\tiny $\frac 1 4$};
      \draw[very thick] (0,1) -- (2,1);
    \end{tikzpicture}
    & \quad \quad
    \begin{tikzpicture}
      \draw[help lines] (-0.3,-0.3) grid (1.3,1.3);
      \draw[->] (-0.6,0) -- (1.6,0) node[right] {\tiny $\hat\lambda_{\max}$};
      \draw[->] (0,-0.6) -- (0,1.6);
      \draw(0,0) node[below left] {\tiny $0$};
      \draw(1,0) node[below left] {\tiny $1$};
      \draw[very thick] (0,0) -- (0.5,0) -- (1,1);
      \draw(1,0.7) node[right] {\tiny $8 (\hat\lambda_{\max} - \tfrac 1 2)$};
    \end{tikzpicture}
  \end{aligned}$$

As an application, let us compute the average value of the reduced purity of a randomly-chosen pure state of bosonic qubits. The reduced purity is by definition equal to $$\label{definition reduced purity}
  {\lVert\rho_1\rVert}_2^2 = \left( \hat\lambda_{\max}^2 + (1-\hat\lambda_{\max})^2 \right),$$ where $\rho_1$ denotes the one-body reduced density matrix and $\hat\lambda_{\max}$ its maximal eigenvalue.

\[purity example\] The average reduced purity of a randomly-chosen pure state of $N$ bosonic qubits is given by $$\frac 1 2 + \frac 1 {2 N}.$$

We will not use directly, but instead work with the Abelian Duistermaat–Heckman measure: Denote by ${\mathbf P}$ the probability distribution of the maximal eigenvalue of the one-body reduced density matrix. By using and $\hat\lambda_{\max} = \frac 1 2 + \frac \lambda {2 N}$, we find that the average reduced purity is given by $$\begin{aligned}
           &\int \left( \hat\lambda_{\max}^2 + \left( 1-\hat\lambda_{\max} \right)^2 \right) \; d{\mathbf P}(\hat\lambda_{\max})\\
    = &\frac 1 {4 N^2} \int \lambda \left( \left( N + \lambda \right)^2 + \left( N - \lambda \right)^2 \right) \, N! \, d\operatorname{DH}^{\operatorname{SU}(2)}_M(\lambda)\\
    = &\frac 1 {2 N^2} \int \left( N^2 \lambda + \lambda^3 \right) \, N! \, d\operatorname{DH}^{\operatorname{SU}(2)}_M(\lambda).
  \end{aligned}$$ By the derivative principle, , this is equal to $$\begin{aligned}
    &\frac 1 {N^2} \int_0^{\infty} \left( N^2 + 3 \lambda^2 \right) \, N! \, d\operatorname{DH}^{\operatorname{U}(1)}_M(\lambda)\\
    = &\frac 1 {2 N^2} \int \left( N^2 + 3 \lambda^2 \right) \, N! \, d\operatorname{DH}^{\operatorname{U}(1)}_M(\lambda).
  \end{aligned}$$ In we have seen that $N! \, \operatorname{DH}^{\operatorname{U}(1)}_{M}$ is the probability distribution of the sum of $N$ independent random variables that are uniformly distributed on the interval $[-1,1]$. Since the variance of any such random variable is $\frac 1 3$ and since variances of independent random variables are additive, the above is equal to $$\frac 1 {2N^2} \left( N^2 + 3 \frac N 3 \right)
    = \frac 1 2 + \frac 1 {2N}.
    \qedhere$$

In accordance with the concentration of measure phenomenon, $\rho_1 \rightarrow {\mathbf 1}/2$ in distribution as $N \rightarrow \infty$. We remark that our result matches [@mullerdahlstenvedral11 Theorem 34] if one works out the quantities left uncalculated therein. Note that our proof illustrates the power of the derivative principle: Instead of explicitly computing the eigenvalue distribution, we can reduce to the Abelian Duistermaat–Heckman measure by differentiating the quantity we are interested in.

Pure States of Bipartite Systems
--------------------------------

We conclude this series of examples by re-deriving the eigenvalue distribution of the reduced density matrices of a randomly-chosen pure state in the case of a general bipartite quantum system, corresponding to the action of $\operatorname{SU}(a) \times \operatorname{SU}(b)$ on $M = {\mathbb P}({\mathbb C}^{a} \otimes {\mathbb C}^{b})$. Instead of following one of the algorithms it will be most convenient to directly work with the formula given in .

Suppose that $b \geq a$. If $b > a+1$ then is not satisfied (): Indeed, it always follows from the singular value decomposition that $$\label{bipartite spectra}
  \operatorname{eig}\rho_B = (\operatorname{eig}\rho_A, 0, \ldots, 0),$$ so that in this case the non-Abelian moment polytope is contained in the boundary of the positive Weyl chamber. However, of course implies that for any choice of $b \geq a$ the joint eigenvalue distribution is already determined by the eigenvalue distribution of $\rho_A$. That is, it suffices to compute the Duistermaat–Heckman measure for the action of $K = \operatorname{SU}(a)$. Denote by $T$ the standard maximal torus of $\operatorname{SU}(a)$. As in , we identify points $\lambda \in \mathfrak t \cong \mathfrak t^*$ with diagonal density matrices $\hat\lambda = \frac 1 a + i \lambda$.

Clearly, the Abelian moment polytope consists of those $\lambda$ with $\hat\lambda \in \Delta_{a-1}$, and the non-Abelian moment polytope is its intersection with the positive Weyl chamber (cf. ).

\[lloyd pagels abelian\] On the Abelian moment polytope, the Duistermaat–Heckman measure for the action of the maximal torus of $\operatorname{SU}(a)$ is proportional to $$\prod_{j=1}^a \hat\lambda_j^{b-1} ~ d\lambda =
    \prod_{j=1}^a \left(\frac 1 a + i \lambda_j\right)^{b-1} ~ d\lambda.$$

Choose the weight-space decomposition of ${\mathbb C}^a \otimes {\mathbb C}^b$ given by the standard basis vectors $e_j \otimes e_k$ ($j=1,\ldots,a$ and $k=1,\ldots,b$). According to , the density of the Abelian Duistermaat–Heckman measure is given by $$f(\lambda) = \operatorname{vol}\, \{ p_{1,1}, \ldots, p_{a,b} \geq 0 : \sum_{k=1}^b p_{j,k} = \hat\lambda_j = \frac 1 a + i \lambda_j \quad (j=1,\ldots,a) \}$$ with respect to the volume measure $dp/d\lambda$ defined therein. Note that the right-hand side set is the Cartesian product of $a$ rescaled standard simplices. The measure factorizes accordingly, and it is easy to see that $$\begin{aligned}
    f(\lambda)
    = \prod_{j=1}^a \operatorname{vol}\, \{ p_1, \ldots, p_b \geq 0 : \sum_{k=1}^b p_k = \hat\lambda_j \}
    = \frac 1 Z \prod_{j=1}^a \hat\lambda_j^{b-1}
  \end{aligned}$$ for $\hat\lambda \in \Delta_{a-1}$, and zero otherwise, with $Z$ a suitable normalization constant.

\[lloyd pagels\] On the non-Abelian moment polytope, the Duistermaat–Heckman measure for the $\operatorname{SU}(a)$-action on ${\mathbb P}({\mathbb C}^a \otimes {\mathbb C}^b)$ is proportional to $$\prod_{j=1}^a \hat\lambda_j^{b-a}
    \prod_{j < k \leq a} (\hat\lambda_j - \hat\lambda_k) ~
    d\lambda.$$

According to the derivative principle, we have to apply $\prod_{j < k} i(\partial_{\lambda_k} - \partial_{\lambda_j}) = \prod_{j < k} ( \partial_{\hat\lambda_k} - \partial_{\hat\lambda_j} )$ to the Abelian Duistermaat–Heckman density as computed in .

This is a partial differential operator of order $\binom a 2$, therefore the resulting non-Abelian density polynomial has total degree at most $d_{\max} = a(b-1) - a(a-1)/2$. Since we differentiate each variable at most $a-1$ times, it is a multiple of the symmetric polynomial $\prod_{j=1}^a \hat\lambda_j^{b-a}$. On the other hand, the result is evidently antisymmetric, and therefore a multiple of the Vandermonde determinant $\prod_{j < k} (\hat\lambda_j - \hat\lambda_k)$. Since the total degrees add up to $d_{\max}$, this implies the assertion.

In view of , this result implies the following well-known formula [@lloydpagels88; @zyczkowskisommers01]:

The distribution of the eigenvalue spectrum $\hat\lambda = \operatorname{eig}\rho_1$ of a randomly-chosen bipartite pure state $\rho$ on ${\mathbb C}^a \otimes {\mathbb C}^b)$ has Lebesgue density proportional to $$\prod_{j=1}^a \hat\lambda_j^{b-a}
    \prod_{j < k \leq a} (\hat\lambda_j - \hat\lambda_k)^2$$ on the space of eigenvalue spectra $\{ \hat\lambda \in \Delta_{a-1} : \hat\lambda_1 \geq \ldots \geq \hat\lambda_a \}$.

It is also easy to deduce the corresponding formula for the action of $\operatorname{SU}(a) \times \operatorname{SU}(b)$:

\[lloyd pagels diagonal\] Denote by $\Delta = \{ \lambda \in \mathfrak t^*_+ : \hat\lambda \in \Delta_{a-1} \}$ the non-Abelian moment polytope for the $\operatorname{SU}(a)$-action. Then the push-forward of Liouville measure along the moment map for the $\operatorname{SU}(a) \times \operatorname{SU}(b)$-action is given by $$\left\langle (\Phi_{\operatorname{SU}(a) \times \operatorname{SU}(b)})_*(\mu_{{\mathbb P}({\mathbb C}^a \otimes {\mathbb C}^b)}), g \right\rangle
    = \frac 1 Z \int_{\Delta} d\lambda \int_{\mathcal O^{\operatorname{SU}(a)}_\lambda \times \mathcal O^{\operatorname{SU}(b)}_{((\hat\lambda,0,\ldots,0) - \frac {{\mathbf 1}} b)/i}} g,$$ where $Z$ is a suitable normalization constant.

By , each coadjoint orbit $\mathcal O_\lambda$ for $\operatorname{SU}(a)$ is paired with the coadjoint $\operatorname{SU}(b)$-orbit through $((\hat\lambda,0,\ldots,0) - \frac{{\mathbf 1}}{b})/i$. By its $\operatorname{SU}(a) \times \operatorname{SU}(b)$-invariance, on each such pair of coadjoint orbits the push-forward measure is just a multiple of the usual Liouville measure. The assertion follows by observing that the density in is at any point $\lambda$ proportional to the symplectic volume of the corresponding coadjoint $\operatorname{SU}(b)$-orbit.

Multiplicities of Representations {#multiplicities section}
=================================

All results discussed so far can be considered as asymptotic limits of corresponding statements in representation theory, at least if the Hamiltonian $K$-manifold $M$ can be linearized (“quantized”) in a certain technical sense. This is in the spirit of Kirillov’s orbit method and the theory of geometric quantization [@guilleminsternberg77; @guilleminsternberg84; @woodhouse92; @guilleminlermansternberg96; @kirillov99].

In particular, this is the case when $M$ is a $K$-invariant smooth irreducible complex projective subvariety of ${\mathbb P}(V)$ for a finite-dimensional unitary $K$-representation $V$. In this situation, the Fubini–Study form of ${\mathbb P}(V)$ restricts to a non-degenerate symplectic form on $M$, and the $K$-action is Hamiltonian with moment map the restriction of . We still assume that is satisfied.

Coadjoint orbits $\mathcal O_\lambda$ through dominant integral weights $\lambda \in \Lambda^*$ (and only these) can be realized in this setup [@kirillov99]: They are in a natural way projective subvarieties of ${\mathbb P}(V_\lambda)$, where $V_\lambda$ is the unitary $K$-representation with highest weight $\lambda$.

In particular, the quantum marginal problem can be analyzed in this framework: Coadjoint $\operatorname{SU}(d)$-orbits through integral highest weights correspond to Hermitian matrices with integral eigenvalue spectra, and it suffices to consider these, since we can always rescale and take limits, or simply pass to the purified double ().

In we will recall the limit alluded to above. We then proceed to describe the representation-theoretic analogue of the derivative principle: Multiplicities of irreducible $K$-representations can be computed from weight multiplicities by taking finite differences (). In the case of the projective space associated with a unitary $K$-representation, the relevant weight multiplicities are those for the symmetric powers of the representation. In , we give a concrete formula describing these weight multiplicities as the number of integer points in certain rational convex polytopes; we indicate that this is again amenable to algorithmic implementation. Finally, we show that in the limit we recover the corresponding statements of .

The Semi-Classical Limit {#semiclassicallimit subsection}
------------------------

Since the $K$-action on $M \subseteq P(V)$ originates from a linear action on $V$, each graded part of the homogeneous coordinate ring ${\mathbb C}[M]$ is naturally a finite-dimensional $K$-representation and can thus be decomposed into irreducible sub-representations, $${\mathbb C}[M]
    \cong \bigoplus_{k=0}^\infty \bigoplus_{\lambda \in \Lambda^* \cap \mathfrak t^*_+} V_\lambda \otimes \operatorname{Hom}_K(V_\lambda, {\mathbb C}[M]_k).$$ We shall encode their multiplicities, suitably re-scaled, in the following sequence of discrete measures, $$\label{discrete irrep measure}
  \mu^K_{M,k} := \frac 1 {k^{n - R}} \sum_{\lambda \in \Lambda^* \cap \mathfrak t^*_+} \dim \operatorname{Hom}_K(V^*_\lambda, {\mathbb C}[M]_k) \, \delta_{\lambda/k}.$$ The factor $k^R$ accommodates for the growth of the dimension of a generic irreducible representation in the coordinate ring, which has highest weight in $\mathfrak t^*_{>0}$ (we still assume that is in place).

It is well-known that in the *semi-classical limit* $k \rightarrow \infty$ this sequence of measures converges in distribution to the non-Abelian Duistermaat–Heckman measure [@heckman82; @guilleminsternberg82b; @sjamaar95; @meinrenken96; @meinrenkensjamaar99; @vergne98], $$\label{semiclassicallimit}
  \mu^K_{M,k} \rightarrow \operatorname{DH}^K_M.$$ In fact, one can show using the Hirzebruch–Riemann–Roch theorem that the piecewise polynomial density function of the Duistermaat–Heckman measure is at any regular point $\lambda \in \mathfrak t^*_{>0}$ for $\Phi_K$ (equivalently, for $\Phi_T$) given by $$\label{semiclassicallimit pointwise}
  \lim_{k \rightarrow \infty}
  \frac 1 D
  \frac 1 {k^{n-R-r}}
  \dim \operatorname{Hom}_K(V_{k \lambda^*}, {\mathbb C}[M]_k).$$ Here, $D$ is the number of elements in the generic stabilizer of points in $M$,[^6] and the right-hand side is understood as the limit over the subsequence with $k \lambda$ integral. The additional factor $k^r$ comes from the fact that we now consider the density with respect to the $r$-dimensional Lebesgue measure $d\lambda$.

It is well-known that the right-hand side multiplicity in is a *quasi-polynomial* in $k$, i.e., a polynomial whose coefficients are periodic functions of $k$. Observe that the existence of the limit implies that (a) the degree of this quasi-polynomial is at most $n-R-r$, and (b) if it is of maximal degree then its leading order coefficient is in fact a constant independent of $k$ (namely, the limit). We shall therefore call the limit the *maximal-order growth coefficient* of the quasi-polynomial.

The rational points of the moment polytope $\Delta_K(M)$ are precisely those of the form $\lambda / k$ with $V^*_\lambda \subseteq {\mathbb C}[M]_k$ [@brion87]. In other words, $\operatorname{supp}\operatorname{DH}^K_M \cap \mathfrak t^*_{{\mathbb Q}} = \bigcup_k \operatorname{supp}\mu^K_{M,k}$.

Restricting the action to the maximal torus $T \subseteq K$, we observe that the Abelian Duistermaat–Heckman measure $\operatorname{DH}^T_M$ captures the asymptotic distribution of weights in the homogeneous coordinate ring of $M$, i.e., the asymptotics of the character of ${\mathbb C}[M]_k$ as $k \rightarrow \infty$.

For strictly dominant and integral $\lambda \in \Lambda^* \cap \mathfrak t^*_{>0}$, the Borel–Weil theorem shows that the homogeneous coordinate ring of the coadjoint orbit $\mathcal O_\lambda \subseteq P(V_\lambda)$ is equal to $${\mathbb C}[\mathcal O_\lambda] = \bigoplus_{k=0}^\infty V^*_{k \lambda}.$$ Therefore, all the multiplicity measures $\mu_{\mathcal O_{\lambda},k}^K$ (and hence their limit) are equal to the Dirac measure at $\lambda$.

Multiplicities of Irreducible Representations via Finite Differences {#multiplicities of irreducibles}
--------------------------------------------------------------------

Multiplicities of weights and highest weights in finite-dimensional $K$-representations are related by iteratively taking (negative) finite differences in the directions of the positive roots. This can be seen as the “quantized” version of the derivative principle (). Its proof is in essence a rephrasing of the Weyl character formula, an idea which goes back at least to Steinberg [@steinberg61].

\[steinberg lemma\] Denote by $m_K$ and $m_T$ the highest weight and weight multiplicity function, respectively, of a finite-dimensional $K$-representation $V$. Then on the positive Weyl chamber we have $$m_K = {\left.\left(\prod_{\alpha > 0} - D_\alpha \right) m_T\vphantom{\big|}\right|_{\mathfrak t^*_+}},$$ where $(D_\alpha m)(\lambda) = m(\lambda + \alpha) - m(\lambda)$ is the finite-difference operator in direction $\alpha$. Note that any two of the operators $D_\alpha$ commute, so that their product is independent of the order of multiplication.

By linearity of the finite-difference operators it suffices to establish the lemma for a single irreducible representation $V_\lambda$ of highest weight $\lambda$. It will be convenient to work with the formal character $\mathrm{ch}(V_\lambda) = \sum_\mu m_T(\mu) \, e^\mu$ [@cartersegalmacdonald95; @knapp02]. By the Weyl character formula, $$\prod_{\alpha > 0} \left( 1 - e^{-\alpha} \right) \mathrm{ch}(V_\lambda) =
    e^{-\rho} \sum_{w \in W} (-1)^{l(w)} e^{w(\lambda + \rho)}$$ where $W$ is the Weyl group, $l(w)$ the length of a Weyl group element $w$, and $\rho$ half the sum of the positive roots.

Now observe that the left-hand side is the generating function of $\left( \prod_{\alpha > 0} - D_\alpha \right) m_T$, since taking finite differences corresponds to multiplying the generating function by $1 - e^{-\alpha}$. Up to terms corresponding to non-dominant weights, the right-hand side is equal to $e^\lambda$, which is the generating function of $m_K = \delta_{\lambda,-}$. The assertion follows from this.

Note that the Weyl character formula can be seen as the representation-theoretic analogue of the Harish-Chandra formula that was used to establish . In the semi-classical limit , the finite differences become infinitesimal and we recover an alternative proof of in the algebro-geometric setting. This argument can also be turned around to establish for general compact Lie groups $K$ from its Abelian version [@guilleminsternberg84 (34.8)] and .[^7]

Multiplicities for Projective Spaces {#multiplicities for projective space}
------------------------------------

As in , let $M = {\mathbb P}(V)$ be the complex projective space for a unitary $K$-representation $V$. Its homogeneous coordinate ring is equal to the symmetric algebra, $${\mathbb C}[{\mathbb P}(V)] = \operatorname{Sym}(V) = \bigoplus_{k=0}^\infty \operatorname{Sym}^k(V).$$

Choose a weight-space decomposition $V = \bigoplus_{k=0}^n {\mathbb C}v_k$, and identify $V \cong {\mathbb C}^{n+1}$ and $\operatorname{U}(V) \cong \operatorname{U}(n+1)$ accordingly. Observe that the maximal torus $T \subseteq K$ acts via the standard maximal torus of $\operatorname{U}(n+1)$, that is, the set of unitary diagonal matrices, which we denote by $\tilde T$.

Each symmetric tensor power $\operatorname{Sym}^k({\mathbb C}^{n+1})$ is an irreducible representation of $\operatorname{U}(n+1)$. Its weight spaces are all one-dimensional, and the weights that occur are precisely the $\lambda = \operatorname{diag}(\lambda_0, \ldots, \lambda_n)$ with $i \lambda_j \in {\mathbb Z}$, $\lambda_j \geq 0$, and $\sum_j \lambda_j = k$ [@fulton97]. Clearly, we can identify this set of weights with the integral points in $k \Delta_n$, where $\Delta_n$ is the $n$-dimensional standard simplex in ${\mathbb R}^{n+1}$. In the language of Young diagrams, these are the weight vectors corresponding to semistandard Young tableaux of shape $(k)$ with entries in $\{0,\ldots,n\}$.

To determine the weight multiplicities with respect to $T \subseteq K$, we have to “restrict” each weight to $\mathfrak t$. This corresponds precisely to applying the map $P \colon {\mathbb R}^{n+1} \rightarrow \mathfrak t^*, (t_k) \mapsto \sum_k t_k \omega_k$ introduced in . Therefore, the multiplicity in $\operatorname{Sym}^k(V)$ of a weight $\lambda \in \Lambda^*$ is given by counting integral points in a rational convex polytope parametrized by $k$ and $\lambda$: $$\label{weight multiplicities for projective space}
  m_{T,k}(\lambda) = \#\left(\Delta(\lambda, k) \cap {\mathbb Z}^{n+1}\right),$$ where $$\Delta(\lambda, k) = \Big\{ (t_j) \in {\mathbb R}^{n+1} : t_j \geq 0, \sum_{j=0}^n t_j \omega_j = \lambda, \sum_{j=0}^n t_j = k \Big\}.$$

\[future remark\] Such *vector partition functions* can be evaluated efficiently using Barvinok’s algorithm if the group $K$ and the ambient dimension $\dim V = n + 1$ is fixed [@barvinok93; @barvinok94; @barvinokpommersheim99], namely in time $O(\mathrm{poly}(\log k))$. In fact, $m_{T,k}$ is a piecewise quasi-polynomial function in both $\lambda$ and $k$, and there are parametric generalizations of Barvinok’s algorithm for computing these quasi-polynomials [@verdoolaegeseghirbeylsetal07; @verdoolaegebruynooghe08]. Since we can compute multiplicities of irreducible $K$-representations by taking finite differences of weight multiplicities in the direction of positive roots (), this can also be done efficiently if $K$ is fixed. We will report on a generalization of this technique to the general branching problem for compact connected Lie groups in a forthcoming article [@christandldoranwalter12].

There is also a jump formula by Boysal and Vergne [@boysalvergne09], which as in can be used to inductively compute the quasi-polynomials chamber by chamber.

We now turn to the semi-classical limit. As $k \rightarrow \infty$, it is clear that $$\mu^{\tilde T}_{{\mathbb P}(V), k} = \frac 1 {k^n} \sum_{\lambda \in \Delta_n \cap \frac 1 k {\mathbb Z}^{n+1}} \delta_\lambda$$ converges to Lebesgue measure on the standard simplex $\Delta_n$, normalized to total volume $$\lim_{k \rightarrow \infty} \frac 1 {k^n} \dim \operatorname{Sym}^k(V) = \lim_{k \rightarrow \infty} \frac 1 {k^n} {\binom{n+k}{n}} = \frac 1 {n!}.$$ Therefore, $\mu^T_{{\mathbb P}(V),k}$ converges to the push-forward of Lebesgue measure on $\Delta_n$ along the map $P$. By the semi-classical limit , this is of course equivalent to the assertion of . Moreover, note that the quantity $$m_{T,k}(k \lambda) =
  \#\left(\Delta(k \lambda, k) \cap {\mathbb Z}^{n+1}\right) =
  \#\left(\Delta(\lambda, 1) \cap \tfrac 1 k {\mathbb Z}^{n+1}\right)$$ is the Ehrhart quasi-polynomial associated to rational polytope $\Delta(\lambda, 1)$ [@beckrobins09]. It is intuitively clear that its growth in $k$ should be related to the volume of this polytope. Indeed, $$m_{T,k}(k \lambda) = k^{n-r} \operatorname{vol}\Delta(\lambda, 1) + O(k^{n-r-1}),$$ where $\operatorname{vol}$ is the $(n-r)$-dimensional volume with respect to the measure $dt/d\hat\lambda$ defined in [@beckrobins09 Exercise 3.29]. Observe that this agrees with and : The maximal-order growth coefficient is a constant equal to the Abelian Duistermaat–Heckman density at point $\lambda$.

Kronecker and Plethysm Coefficients {#kronecker section}
===================================

In this section, we describe the representation theory of the quantum marginal problem in more detail. For distinguishable particles, the relevant multiplicities can be expressed in terms of decomposing tensor products of irreducible representations of the symmetric group (). In particular, the joint eigenvalue distribution of the reduced density matrices of a tripartite pure state is determined by the asymptotics of the Kronecker coefficients (see ). We emphasize that by specializing the method described in we get a novel algorithm for computing Kronecker coefficients which is efficient for Young diagrams of bounded height. Indistinguishable particles correspond to certain plethysm coefficients and we conclude by illustrating this connection ().

Kronecker Coefficients {#kronecker coefficients}
----------------------

Recall that for $N$ distinguishable particles we have to consider the action of $K = \operatorname{SU}(d_1) \times \ldots \times \operatorname{SU}(d_N)$ on a coadjoint $\operatorname{SU}(d_1 \cdots d_N)$-orbit $M = \mathcal O_{\tilde\lambda}$, where we now assume that $\tilde\lambda$ is an integral weight in $\mathfrak t^*_{>0}$. The multiplicity measures $\mu^K_{M,k}$ are determined by the decomposition of the homogeneous coordinate ring $${\mathbb C}[M] = \bigoplus_{k=0}^\infty V^{d_1 \cdots d_N}_{k \tilde\lambda}$$ into $K$-isotypical components (the superscript labels the corresponding $\operatorname{SU}$).

We can express this equivalently using the representation theory of the symmetric group $S_m$. Recall that by Schur–Weyl duality the diagonal action of $\operatorname{SU}(d)$ and the permutation action of $S_m$ on $({\mathbb C}^d)^{\otimes m}$ generate each other’s commutant, so that $$\label{schur weyl decomposition}
  ({\mathbb C}^d)^{\otimes m} \cong \bigoplus_{\mu} V^d_\mu \otimes [\mu].$$ Here, the sum runs over all Young diagrams $\mu=(\mu_1,\ldots,\mu_d)$ with ${\lvert\mu\rvert} := \sum_j \mu_j = m$ boxes and at most $d$ rows, $V^d_\mu$ is the irreducible representation of $\operatorname{SU}(d)$ with highest weight $X \mapsto i \sum_j X_j \mu_j$, and $[\mu]$ is the corresponding irreducible representation of $S_m$ (see [@fulton97] for details). We shall freely identify Young diagrams and the corresponding highest weights.

In particular, we can realize the irreducible representation $V^{d_1 \cdots d_N}_{k \tilde\lambda}$ in $({\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N})^{\otimes {\lvertk \tilde\lambda\rvert}}$. Comparing the Schur–Weyl decomposition for the full Hilbert space with the tensor product of the decompositions for the individual subsystems, we find that $$V^{d_1 \cdots d_N}_{k \tilde\lambda}
  \cong
  \bigoplus_{\lambda_1,\ldots,\lambda_N}
  V^{d_1}_{\lambda_1} \otimes \ldots \otimes V^{d_N}_{\lambda_N} \otimes
  \operatorname{Hom}_{S_{{\lvertk\tilde\lambda\rvert}}}([k \tilde\lambda], [\lambda_1] \otimes \ldots \otimes [\lambda_N]),$$ where the sum runs over the Young diagrams $\lambda_i$ with ${\lvertk\tilde\lambda\rvert}$ boxes and at most $d_i$ rows. Therefore, $$\begin{aligned}
  \label{generalized kroneckers}
  \mu^K_{M,k}
  = &\frac 1 {k^{n-R}}
  \sum_{\lambda_1,\ldots,\lambda_N} \dim \operatorname{Hom}_{S_{{\lvertk\tilde\lambda\rvert}}}([k \tilde\lambda], [\lambda_1]^* \otimes \ldots \otimes [\lambda_N]^*) ~ \delta_{(\lambda_1/k,\ldots,\lambda_N/k)}\\
  = &\frac 1 {k^{n-R}}
  \sum_{\lambda_1,\ldots,\lambda_N} \dim \operatorname{Hom}_{S_{{\lvertk\tilde\lambda\rvert}}}([k \tilde\lambda], [\lambda_1] \otimes \ldots \otimes [\lambda_N]) ~ \delta_{(\lambda_1/k,\ldots,\lambda_N/k)},
\end{aligned}$$ where the latter identity holds due to the self-duality of the representations of the symmetric group. In particular, the rational points of the non-Abelian moment polytope $\Delta_K(\mathcal O_{\tilde\lambda})$ are precisely $$\begin{aligned}
  \label{rational points QMP moment polytope}
    \bigcup_k \left\{ (\lambda_1/k,\ldots,\lambda_N/k) : \begin{array}{l}
    [k\tilde\lambda] \subseteq [\lambda_1] \otimes \ldots \otimes [\lambda_N], \text{ where the}\\
    \text{$\lambda_i$ have ${\lvertk\tilde\lambda\rvert}$ boxes and at most $d_i$ rows}
  \end{array} \right\}.\end{aligned}$$ See for an illustration of the multiplicity measures corresponding to the mixed-state quantum marginal problem for two qubits discussed in .

![(a) and (b) Illustration of the multiplicity measures $\mu^K_{M,k}$ for the mixed-state quantum marginal problem of two qubits with global spectrum $(4/7,2/7,1/7,0)$ and $k=28, 56$, which have been computed by the algorithm described in . (c) Their semi-classical limit, i.e., the corresponding Duistermaat–Heckman measure as computed in .[]{data-label="hedgehog figure"}](irrep_mults_4210_4.png "fig:"){width="4cm"} ![(a) and (b) Illustration of the multiplicity measures $\mu^K_{M,k}$ for the mixed-state quantum marginal problem of two qubits with global spectrum $(4/7,2/7,1/7,0)$ and $k=28, 56$, which have been computed by the algorithm described in . (c) Their semi-classical limit, i.e., the corresponding Duistermaat–Heckman measure as computed in .[]{data-label="hedgehog figure"}](irrep_mults_4210_8.png "fig:"){width="4cm"} ![(a) and (b) Illustration of the multiplicity measures $\mu^K_{M,k}$ for the mixed-state quantum marginal problem of two qubits with global spectrum $(4/7,2/7,1/7,0)$ and $k=28, 56$, which have been computed by the algorithm described in . (c) Their semi-classical limit, i.e., the corresponding Duistermaat–Heckman measure as computed in .[]{data-label="hedgehog figure"}](bravyi_polytope_4210.png "fig:"){width="4.8cm"}

We can write the multiplicities in in the following symmetric form: $$\dim \operatorname{Hom}_{S_{{\lvertk\tilde\lambda\rvert}}}([k \tilde\lambda], [\lambda_1]^* \otimes \ldots \otimes [\lambda_N]^*) =
  \dim \left(
    [\lambda_1]^* \otimes \ldots \otimes [\lambda_N]^* \otimes [k \tilde\lambda]^*
  \right)^{S_{{\lvertk\tilde\lambda\rvert}}}.$$ Observe that the right-hand side is a multiplicity for the pure-state quantum marginal problem for ${\mathbb C}^{d_1} \otimes \ldots \otimes {\mathbb C}^{d_N} \otimes {\mathbb C}^{d_1 \cdots d_N}$. Indeed, the homogeneous coordinate ring of a projective space is just the symmetric algebra (), whose graded parts correspond to the trivial representations of the symmetric groups. This is the representation-theoretic perspective on purification (cf. , in particular ).

For the tripartite pure-state quantum marginal problem (equivalently, the mixed-state bipartite quantum marginal problem), the relevant multiplicities are the well-known *Kronecker coefficients* of the symmetric group, $$g_{\lambda,\mu,\nu} = \dim \left( [\lambda] \otimes [\mu] \otimes [\nu] \right)^{S_k}.$$ They are the symmetric group analogue of the Littlewood–Richardson coefficients of the unitary group (in fact, the latter can be considered as a special case) but much harder to compute in general, since there is no combinatorial description like the Littlewood–Richardson rule.

The corresponding characterization of the non-Abelian moment polytope has already been observed in [@christandlmitchison06; @klyachko04; @christandlharrowmitchison07], as well as in [@daftuarhayden04] for the projection onto two of the subsystems. The semi-classical limit refines this characterization: Not only can one read off the existence of quantum states with given marginal eigenvalue spectra from the asymptotic non-vanishing of the corresponding Kronecker coefficients $g_{k\lambda,k\mu,k\nu}$, but their growth also encodes the probability of finding these eigenvalue spectra when the global state is chosen according to the invariant probability measure. Explicitly, states that $$\label{kroneckerlimit}
  \frac 1 {k^p} \sum_{\lambda,\mu,\nu} g_{\lambda,\mu,\nu} \, \delta_{\lambda/k, \mu/k, \nu/k}
  \rightarrow \operatorname{DH}^{{\mathbb P}({\mathbb C}^a \otimes {\mathbb C}^b \otimes {\mathbb C}^c)}_{\operatorname{SU}(a) \times \operatorname{SU}(b) \times \operatorname{SU}(c)},$$ where $p = n - R = {abc - 1 - a(a-1)/2 + b(b-1)/2 + c(c-1)/2}$, and where the sum runs over all Young diagrams $\lambda, \mu, \nu$ with $k$ boxes and at most $a$, $b$ and $c$ rows, respectively.

The method described in in particular provides a novel algorithm for computing the Kronecker coefficients which is efficient for Young diagrams of bounded height: Using the finite-difference formula of , we can reduce to the computation of a bounded number of weight multiplicities , which using Barvinok’s algorithm can be evaluated in polynomial time in the input size, i.e., in time $O(\mathrm{poly}(\log k))$, where $k$ is the number of boxes of the Young diagrams. As mentioned in , we will elaborate on this algorithm in a forthcoming article [@christandldoranwalter12].

Plethysm Coefficients {#plethysms}
---------------------

While the quantum marginal problem for distinguishable particles is connected to (generalized) Kronecker coefficients, it is for indistinguishable particles related to certain plethysm coefficients. Indeed, if $M = {\mathbb P}(V_\lambda)$ for an irreducible $\operatorname{SU}(d)$-representation $V_\lambda$ then its coordinate ring consists of the *plethysms* $${\mathbb C}[M] = \bigoplus_{k=0}^\infty \operatorname{Sym}^k(V_\lambda).$$ See e.g. [@macdonald95] for more information on plethysms, which are in general defined as the composition of Schur functors. In particular, the bosonic and fermionic pure-state marginal problem are related to the asymptotics of $\operatorname{Sym}^k(\operatorname{Sym}^N({\mathbb C}^d))$ and $\operatorname{Sym}^k({\Lambda}^N({\mathbb C}^d))$, respectively, as $k \rightarrow \infty$.

Let us illustrate this by giving an alternative derivation of the Duistermaat–Heckman measures for $N$ bosonic qubits (cf. ). We will explicitly compute the asymptotic weight multiplicity distribution of the plethysm $\operatorname{Sym}^k(\operatorname{Sym}^N({\mathbb C}^2))$ as $k \rightarrow \infty$, and then apply the derivative principle. The main combinatorial tool we shall employ are the $q$-binomial coefficients $${\genfrac{[}{]}{0pt}{}{n}{k}}_q = \frac{[n]_q !}{[k]_q ! [n-k]_q !}
  .$$ Recall that these are defined in terms of the $q$-integers $[n]_q = \frac{1-q^n}{1-q}$ and $q$-factorials $[n]_q ! = [n]_q [n-1]_q\ldots [1]_q$. We start with the following description of the character of $\operatorname{Sym}^k(\operatorname{Sym}^N({\mathbb C}^2))$ in terms of $q$-binomial coefficients:

\[springer character qbinomial\] Let $k,N \in {\mathbb N}$ and $q=e^{\omega_1}$. Then, $$\mathrm{ch} \left( \operatorname{Sym}^k \left( \operatorname{Sym}^N \left( {\mathbb C}^2 \right) \right) \right) = {\genfrac{[}{]}{0pt}{}{k+N}{N}}_{q^2} q^{-kN}.$$

\[asy eq q-stuff\] As functions on the open unit disk $\{ q \in {\mathbb C}: |q|<1 \}$ one has for fixed $N \in {\mathbb N}$ and $k \rightarrow \infty$ the following asymptotic equivalence $$[N]_q ! {\genfrac{[}{]}{0pt}{}{k+N}{N}}_q \sim [k]_q^N
    .$$

Following [@kaccheung02 (9.1)], for any fixed $c \in {\mathbb N}$ and $|q|<1$ one has $$\lim_{k \rightarrow \infty} \frac{1-q^{k+c}}{1-q} = \frac{1}{1-q}.$$ By applying this identity both to the numerator and the denominator, $$\lim_{k \rightarrow \infty} \frac{[k+c]_q}{[k]_q} = 1
    .$$ Hence $$\lim_{k \rightarrow \infty} \frac{ [N]_q !{\genfrac{[}{]}{0pt}{}{k+N}{N}}_q}{[k]_q^N} =
    \lim_{k \rightarrow \infty} \frac{[k + N]_q [k+N-1]_q \ldots [k+1]_q}{[k]_q^N} =
    1
    .
    \qedhere$$

The following corollary is an easy application of Osgood’s theorem [@osgood01] to (see e.g. [@beardonminda03]).

\[osgood\] Fix $N$, and define $f_k(q) = [N]_q ! {\genfrac{[}{]}{0pt}{}{k+N}{N}}_q$ and $g_k(q) = [k]_q^N$. Then, the sequences $f_k$ and $g_k$ converge, as $k \rightarrow \infty$, pointwise to the same holomorphic function on some open dense subset in the open unit disk $\{ q \in {\mathbb C}: |q|<1 \}$. In particular, the limits are equal as power series.

We can use this result to extract asymptotic multiplicity information.

\[asymptotic uniformness\] The discrete measures $\mu_{(N),k}^{\operatorname{U}(1)} := \mu_{\operatorname{Sym}^N({\mathbb C}^2),k}^{\operatorname{U}(1)}$ as defined in tend in the limit $k \rightarrow \infty$ to $\frac 1 {N!}$ times the probability distribution of the sum of $N$ independent random variables uniformly distributed on $[-1,1]$.

By , $\mu_{(N),k}^{\operatorname{U}(1)}$ is a finite measure with generating function $$\int q^x \, d\mu_{(N),k}^{\operatorname{U}(1)}(x \, \omega_1) =  \frac 1 {k^N} {\genfrac{[}{]}{0pt}{}{k+N}{N}}_{q^{2/k}} q^{-N}.$$ Let $\nu_k$, $\omega_k$ be finite measures with generating functions $[N]_{q^{2/k}}!$ and $$\frac 1 {k^N} [k]^N_{q^{2/k}} q^{-N} = \left( \frac{q^{-1} + q^{-1+2/k} + \ldots + q^{1-2/k}}{k} \right)^N,$$ respectively. Obviously, $\omega_k$ is asymptotically distributed like the sum of $N$ independent random variables uniformly distributed on the interval $[-1,1]$, and by so is $\mu^{\operatorname{U}(1)}_{(N),k} \star \nu_k$. Since $\nu_k \rightarrow N! \, \delta_0$ as $k \rightarrow \infty$, this implies our assertion.

By using the semi-classical limit , we conclude once again that the Abelian Duistermaat–Heckman measure is given by the formula that was established in . The non-Abelian Duistermaat–Heckman measure is obtained as in by applying the derivative principle.

For the plethysms $\operatorname{Sym}^k(\operatorname{Sym}^2({\mathbb C}^2))$ we can also illustrate the semi-classical limit for the $\operatorname{SU}(2)$-action, since the decomposition into irreducible $\operatorname{SU}(2)$-representations is well-known [@macdonald95 §1.5, Example 6 (a)]: $$\operatorname{Sym}^k(\operatorname{Sym}^2 ({\mathbb C}^2)) \cong
    \operatorname{Sym}^{2k}({\mathbb C}^2) \oplus
    \operatorname{Sym}^{2k-4}({\mathbb C}^2) \oplus
    \ldots \oplus
    \operatorname{Sym}^{2 / 0}({\mathbb C}^2)$$ The last summand is $\operatorname{Sym}^2 ({\mathbb C}^2)$ for odd $k$, and $\operatorname{Sym}^0 ({\mathbb C}^2)$ for even $k$. Therefore, the discrete measures as defined in are given by $$        \mu^{\operatorname{SU}(2)}_{{\mathbb P}(\operatorname{Sym}^2({\mathbb C}^2)),k} = \sum_{l=2k, 2k-4, \ldots, 2 \vert 0} \frac 1 {k^{2-1}} \delta_{\frac l k}.$$ In the limit $k \rightarrow \infty$, they converge to the non-Abelian Duistermaat–Heckman measure as computed in , $$\operatorname{DH}^{\operatorname{SU}(2)}_{{\mathbb P}(\operatorname{Sym}^2({\mathbb C}^2))} = - \tfrac 1 4 \left( (x+2)^0_+ - 2 x^0_+ + (x-2)^0_+ \right) dx =
    \tfrac 1 4 {\mathbf 1}_{[0,2)}(x) dx.$$ See for an illustration.

The description of the character of $\operatorname{Sym}^k(\operatorname{Sym}^N({\mathbb C}^2))$ via $q$-binomial coefficients has the additional advantage that one is able to compute all higher cumulants and moments of the associated distribution for any fixed $k$ and $N$ (see [@MR2019639]). This is due to a method by Panny [@MR845446].

Acknowledgements {#acknowledgements .unnumbered}
================

We would like to thank Alonso Botero, Emmanuel Briand, Peter Bürgisser, Benoît Collins, David Gross, Christian Ikenmeyer, Markus P. Müller, Mercedes Rosas, and Volkher Scholz for helpful discussions. We thank Graeme Mitchison for joint initial discussions on the topic of asymptotics of Kronecker coefficients. The second author would like to express his particular gratitude to Frances Kirwan for many fruitful discussions regarding moment maps and invariant theory and their many uses.

This work is supported by the Swiss National Science Foundation (grants PP00P2–128455 and 200021\_138071), the German Science Foundation (SFB/TR12, and grants CH 843/1–1 and CH 843/2–1), the National Center of Competence in Research ‘Quantum Science and Technology’, and the Excellence Initiative of the German Federal and State Governments through the Junior Research Group Program within the Institutional Strategy ZUK 43.

[^1]: \[non-compact footnote\]In fact, the $p$-th summand of is precisely the Duistermaat–Heckman measure corresponding to the isotropy representation of $T$ on the symplectic vector space $T_p M$, which is of course a non-compact symplectic manifold and, strictly speaking, does not fit into our setup. The decomposition of $\mathfrak t^*$ into regular chambers for the moment map of $M$ is refined by the common refinement of the chamber decompositions for the $T_p M$ (cf. ).

[^2]: This is our reason for choosing a different definition for walls than the one used in [@boysalvergne09]. There, walls were defined as linear hyperplanes spanned by $r-1$ linearly independent vectors.

[^3]: This density of course only depends on the hyperplane through $\hat W$, and can therefore re-used for all other walls that span the same hyperplane.

[^4]: The push-forward of Lebesgue measure on ${\mathbb R}^{n+1}_{\geq 0}$ along $\hat P$ can also be understood as the Duistermaat–Heckman measure associated with the Hamiltonian $T \times \operatorname{U}(1)$-action on the complex vector space $V$, where $\operatorname{U}(1)$ acts by scalar multiplication (cf. ).

[^5]: This density of course only depends on the hyperplane through $W$, and can therefore be re-used for all other singular walls that lie on the same hyperplane.

[^6]: Quotienting out a discrete subgroup leaves the Duistermaat–Heckman measure invariant, but changes the weight lattice, and therefore the normalization of the Lebesgue measure $d\lambda$.

[^7]: We thank Allen Knutson for pointing this out, as well as for sketching a self-contained proof of .
