The quantum marginal problem for symmetric states: applications to variational optimization, nonlocality and self-testing

In this paper, we present a method to solve the quantum marginal problem for symmetric $d$-level systems. The method is built upon an efficient semi-definite program that determines the compatibility conditions of an $m$-body reduced density with a global $n$-body density matrix supported on the symmetric space. We illustrate the applicability of the method in central quantum information problems with several exemplary case studies. Namely, (i) a fast variational ansatz to optimize local Hamiltonians over symmetric states, (ii) a method to optimize symmetric, few-body Bell operators over symmetric states and (iii) a set of sufficient conditions to determine which symmetric states cannot be self-tested from few-body observables. As a by-product of our findings, we also provide a generic, analytical correspondence between arbitrary superpositions of $n$-qubit Dicke states and translationally-invariant diagonal matrix product states of bond dimension $n$.


I. INTRODUCTION
The quantum marginal problem (QMP) is ubiquitous not only in modern physics, but also in modern chemistry, where it is usually referred to as the n-representability problem [1]. The QMP can be stated as determining whether a set of reduced density matrices (RDMs) are compatible with a global wavefunction. The QMP arises naturally when computing physically important quantities such as the energy of a system or its entropy, as they often only depend on few particles. As an illustrative example, let us imagine one is interested in computing the ground energy of a k-local Hamiltonian H = i H i , where each H i acts nontrivially on at most k particles. The solution to this problem is ψ|H|ψ , where |ψ is an eigenvector of H with lowest corresponding eigenvalue. Unfortunately, the amount of computational resources to describe |ψ grows, in general, exponentially with the system size, rendering this approach impractical. Alternatively, one can exploit the fact that H is a sum of much simpler terms, and compute instead ψ|H|ψ = i Tr[H i ρ i ], where each ρ i is the reduced density matrix of |ψ ψ| on the particles H i acts upon. The latter formulation, however, only appears to circumvent the exponential cost of describing |ψ . As a matter of fact, it actually comes at the cost of knowing the compatibility conditions of {ρ i } with a global state |ψ ψ|.
Despite this apparent simplification, the QMP has challenged the physics and chemistry communities since the 60s and every nontrivial advance has already supposed a milestone in the field [2]. The QMP is strongly believed to be very hard, even for a quantum computer: It is complete for the complexity class Quantum Merlin-Arthur (QMA) [3] which, roughly speaking, is the analogous of NP for a quantum computer. This may be not so surprising, since the k-local Hamiltonian problem itself is QMA-complete [4], even for k = 2 [5] or for quantum systems on a 1D geometry [6], and existing quantum algorithms take typically exponential time to solve it [7][8][9][10]. In spite of these intractability results, tremendous progress has been achieved over the years in the QMP, even before the advent of the quantum information processing era [11][12][13].
A great deal of the aforementioned progress has centered onto the one-body RDM problem; i.e., determining if a set of onebody RDMs is compatible with a global (pure) state. Klyachko showed that in this case it is sufficient to characterize the QMP compatibility conditions solely from the eigenvalues of the one-body RDMs [2,14]. For instance, the case of RDMs of a bipartite system is completely solved in terms of linear inequalities on the spectra [14], a rather surprising fact, since it is not obvious that compatible spectra form convex sets, let alone polytopes. This formalism has allowed for an elegant and mathematically tractable charcterization, from which further results have stemmed [15,16], for instance, in the context of witnessing genuinely multipartite entanglement from one-body RDMs [13].
The few-body QMP is encompassed with substantial additional challenges. First of all, the relevant information cannot be extracted solely from the eigenvalues of the RDMs, as it was the case on the one-body case. In the one-body case, the supports of different RDMs were necessarily disjoint; therefore, the action of local unitaries did not affect global compatibility and, in consequence, only the spectrum of the one-body RDMs was relevant. In the case that one considers few-body RDMs, their supports may intersect. Therefore, these must, at least, have the same reduced density matrix on the intersection of their supports.
Despite the additional requirements, some progress has been made [17]: for instance, almost all four-partite pure states are determined by their two-body marginals [18]. The QMP has also been extensively studied under the bosonic and fermionic formalism [2,12,[19][20][21]: There, one uses the assumption of the global pure state being either fully symmetric or antisymmetric and obtains conditions on the RDMs on a given subset of parties. These RDMs are all equal due to the symmetry of the global state.
In this work we consider the following problem: Given a reduced density matrix σ of m qudits, i.e., acting on the Hilbert space (C d ) ⊗m , is it compatible with a global reduced density matrix ρ acting on the symmetric space of n qudits Sym(C d ) ⊗n ? Note that we denote the symmetric space as the subspace which is spanned by the Dicke states defined in Section II. We present the solution to this problem by analytically writing the compatibility conditions between σ and ρ and we show that they can be efficiently determined numerically as a feasibility semidefinite program (SdP) (Section III).
The core results of our work can be summarized as follows: • We give the analytical conditions for any m-qudit RDM to be compatible with a larger n-qudit symmetric state • We show how these compatibility conditions can be efficienly solved (polynomially in n, with degree d − 1) via a SDP, thus enabling fast optimization over symmetric states Our work has implications in several aspects of quantum information processing, as we show in subsequent sections with exemplary case studies. We show in Section III A how this provides a computationally undemanding variational approach to the ground state energy of any local Hamiltonian. We benchmark our results with different physical models, such as the Lipkin-Meshkov-Glick from nuclear physics [22][23][24], an Ising chain with power-law interactions, and various XXZ chains with transverse and longitudinal magnetic fields (Section IV A). We further showcase how our method leads to a natural tool to optimize symmetric, few-body, Bell inequalities [25][26][27] and we apply our method to show the ground state of some XXZ one-dimensional model with 128 particles contains Bell correlations (Section IV B). In order to benchmark our results in large system sizes, we have also developed an analytical correspondence between pure symmetric states of n qubits and translationally invariant diagonal matrix product states of bond dimension n (Section IV C) which may be of independent interest. Another exemplary case study that stems from our method may have implications on the self-testing of symmetric states from few-body correlators. Our method allows us to find sufficient conditions to certify which symmetric states cannot be self-tested from their marginals, by analyzing when the compatibility conditions do not lead to a unique solution (Section IV D). We discuss the implications of our results and discuss further research directions in Section V.

II. PRELIMINARIES
Symmetric states constitute one of the most prominent classes of quantum states [28]. These are linear combinations of the so-called Dicke states, which arise naturally from the superradiance effect [29]. Dicke states and symmetric states have been successfully prepared in the laboratory in a plethora of systems, ranging from photons [30] to ultracold atoms [31,32]. Their entanglement properties have been extensively studied [33][34][35]. Furthermore, Device-Independent (DI) self-testing protocols exist for Dicke states [36,37], and symmetric states provide an advantage in quantum metrology [38]. Symmetric states are simple to describe, as their permutational invariance allows one to circumvent the exponential growth of the Hilbert space representability problem, therefore being confined to a subspace of the multipartite Hilbert space whose dimension scales only polinomially with n, with degree d − 1: In the case of n qubits, symmetric states are those spanned by the Dicke states, denoted |D n where S n denotes the symmetric group (the group of permutations of n elements), and τ is a permutation acting on the different local Hilbert spaces. For instance, the half-filled Dicke state of n = 4 qubits has two excitations k = 2, but they are delocalized among the different subsystems, in an equally-weighted coherent superposition of the same phase: The number of different terms in a qubit Dicke state of k excitations is given by the combinatorial expression n k and there are n + 1 Dicke states for d = 2.
In the general case of qudits, now one needs to specify how many |1 excitations there are, how many |2 excitations, etc. Hence, it is a natural choice to index Dicke states by partitions of n. We denote λ n as a partition of n in d elements, and it consists of a vector of d non-negative integers that sum d: λ = (λ 0 , . . . , λ d−1 ), where λ i ∈ Z ≥0 counts the number of excitations in state |i , so that i∈[d] λ i = n. We will omit mentioning d whenever it is clear from the context. There are n+d−1 d−1 such partitions and a qudit Dicke state is denoted |D λ , where The number of different terms in Equation (3) is given by the multinomial combinatorial expression For the sake of a more compact notation, in the rest of this manuscript we will denote the qudit Dicke state |D λ simply as |λ , since D is void of meaning.

III. COMPATIBILITY CONDITIONS WITH A GLOBAL SYMMETRIC STATE
Here we outline the fundamental building block of our work. We derive a set of necessary and sufficient conditions for compatibility of an m-qudit (symmetric) density matrix with a global n-qudit symmetric density matrix.
Let us consider a symmetric quantum system ρ of n qudits. Let us denote the components of ρ as ρ λ µ , where λ, µ n. Thus, Our first goal is to compute the partial trace of ρ of n − m subsystems, both in the computational and in the Dicke basis. Note that since ρ is symmetric, we can take the partial trace over any n − m subsystems, as any choice will give the same result.
Before stating Theorem 1, we introduce the following notations. We define [d] := {0, . . . , d − 1}, and we shall use indices with an overhead arrow (e.g. i) to denote matrix elements in the computational basis, as opposed to indices in boldface (e.g. λ) that denote matrix elements in the Dicke basis. We recall once more that such boldface indeces consist in partitions of n, and they are associated to the Dicke basis elements through Eq. (3).
Theorem 1 Let ρ be a symmetric state of n qudits, and m ≤ n. Let i, j ∈ [d] m and α, β m, so that we define σ := Tr n−m (ρ); i.e., Then, we have the coefficients in the computational basis and the coefficients in the Dicke basis a. Proof of Theorem 1: We begin by noting that, in the computational basis, the partial trace of n − m subsystems of an n-qudit density matrix ρ has the following expression: where the operator | denotes index concatenation, and the indices have been properly rearranged so that the n − m traced out parties are the last ones. Hence, our goal is to express the symmetric state in the computational basis in order to apply Equation (9) and then go back to the symmetric space. Since Dicke states are enumerated by the partitions of m it will be useful to define the following function: In words, w k ( i) counts how many coordinates of i are equal to k. It is then natural to define w( i) := (w 0 ( i), . . . , w d−1 ( i)). Note that, by construction, w( i) n for every i ∈ [d] n . The weight counting function w is useful in representing the Dicke state |λ in the computational basis, which we denote |D λ : where δ is the Kronecker Delta function, which is 1 if, and only if, its argument is the zero vector; 0 otherwise. Thanks to Equation (11) we can now define the inclusion operator Π : Sym(C d ) ⊗n → (C d ) ⊗n onto the symmetric space We note that Π † Π = 1 Sym(C d ) ⊗n and that Π † : (C d ) ⊗n Sym(C d ) ⊗n is the projector onto the symmetric space.
Let now ρ be a density matrix on Sym(C d ), whose components are labeled ρ λ µ in the Dicke basis and ρ i j in the computational basis. The relation between them is given by Now we are ready to trace out n − m parties of ρ, e.g. the last. Note that w( a| b) = w( a) + w( b). Hence, using Equation (9) we obtain σ := Tr n−m (ρ) as yielding Equation (7). Finally, Equation (8) is obtained via the transformation (σ α β ) = Π † (σ i j )Π. From Equation (8) we can now obtain a set of necessary and sufficient conditions for a reduced density matrix σ of m qudits to be compatible with a global (possibly mixed) symmetric state of n qudits, which can be expressed as the following semidefinite program (SDP): where we have defined the coefficients a λ,α µ,β from Equation (8); namely, If we think of a λ,α µ,β as the entries of a matrix A α β indexed by λ and µ, then we can express the SDP Equation (15) in canonical form as where ·, · denotes the Hilbert-Schmidt inner product. If the SDP Equation (17) is feasible, then σ admits an extension to a symmetric Dicke state of n qudits ρ which is precisely the solution of that SDP. In Section IV D we discuss about the uniqueness of that solution.

A. Variational ansatz
We can now easily modify the SDP Equation (17) to optimize any linear functional H on ρ, while maintaining compatibility over a given marginal state σ. This is done by considering the SDP The most interesting case arises when such a functional can be expressed as a sum of terms with support on, at most, m qudits. This includes many cases of physical interest, such as Hamiltonians or Bell operators composed of, at most, m-body interactions/correlators. In this case, let us denote H = i H i . Then, H, ρ can be expressed as a linear combination of terms of the form H i , σ , namely: Note that in Equation (19), both ρ and σ are treated as positive-semidefinite variables. The positive-semidefiniteness of σ is automatically implied by that of ρ. In fact, σ can be completely removed from Equation (19) and embedded into the objective function; however we keep it in this form for clarity of exposition. The form of Equation (19) is thus useful to optimize a functional H that depends only on the marginal information contained in the reduced states, while keeping compatibility with a global symmetric state. Recall that the size of ρ depends polinomially on n, with a degree d − 1, so that this procedure is efficient for systems of qudits of large n and fixed d.

IV. SOME APPLICATIONS
The aim of this section is to illustrate several applications in various, apparently uncorrelated, problems in quantum information, all of which have deep roots in the QMP.
In Section IV A we apply Equation (19) to benchmark our method as a variational ansatz to find a fast, upper bound, to the ground state energy and, in some cases, to well approximate the ground state of several paradigmatic Hamiltonians. In Section IV B we adapt our method to optimize Bell functionals composed of symmetric, few-body observables. In Section IV C we provide a method to approximate any n-qubit Dicke state with a translationally-invariant (TI) diagonal matrix product state (MPS) of bond dimension n. Finally, in Section IV D we show how our method can be used to show which symmetric states cannot be self-tested from few-body marginals.

A. Benchmarking the variational ansatz
Here we consider some Hamiltonians of exemplary spin models, in order to benchmark the performance of the variational ansatz presented in Section III A. Furthermore, we provide in Section IV A 7 the runtime and estimated resources consumed by our method, compared to DMRG as a benchmark.
The variational ansatz approximates the ground state and energy of an m-local Hamiltonian by an m-qudit RDM (denoted m-RDM for short) compatible with a global many-body symmetric state. We expect the variational ansatz to provide a good approximation when the coupling interactions are similar between all pairs and to provide exact results when the ground state lies in the symmetric space.

Lipkin-Meshkov-Glick model
Let us start by considering a spin model for which our variational method recovers the ground state exactly. We picked as an example the Lipkin-Meshkov-Glick (LMG) model [22][23][24], which involves long-range interactions that result in ground states that are symmetric under any permutation of the particles. The LMG model was originally proposed in nuclear physics to describe phase transitions in nuclei. However, nowadays it also serves to describe e.g. two-mode Bose-Einstein condensates experiments, since it captures the physics of interacting bosons in a double-well trapping potential. Furthermore, in its isotropic version, the ground states of the LMG model are pure Dicke states, which have been shown to display Bell correlations [25,27,39]. The phase transitions of the general model are also well understood [40]. Exact solutions are known [41][42][43], which we use here to benchmark our method.
In particular, we consider the following LMG Hamiltonian which describes a set of n spin−1/2 particles with anisotropic long-range interactions under an external transverse magnetic field h: where σ (i) k denotes the Pauli matrix in position i and direction k, λ > 0 correspond to ferromagnetic (FM) interactions, λ < 0 correspond to anti-ferromagnetic (AFM) interactions, and γ marks the anisotropy in the coupling terms, with γ = 1 being the isotropic case.
By adapting the SDP optimization problem in Equation (19), we can use the 2-RDM compatibility constraints to approximate its ground state by means of SDP. In particular, the optimization problem to find the best approximation within the variational ansatz can be reduced to the following SDP: where we emphasize thatH : In Figure 1 we show that the ground state energy of the model we considered is faithfully recovered using our variational method.
We note that this method unlocks the possibility to have access to an efficient description of the n-qubit state (given in the symmetric basis Equation (8)), since it is a matrix of size (m + 1) × (m + 1). This, in turn, allows to obtain any associated m-RDM ∀1 ≤ m ≤ n. This enables us to study the method against extensive quantities such as entropy, since it is now easy to obtain the m-block size entanglement entropy S m,n = − m i=0 p i log 2 p i by finding the eigenvalues p i of the m-RDM [40]. In Figure 1 we have used such a procedure to obtain the half-system entanglement entropy. The method reproduces the features of the LMG model phase diagram, as expected.  [44]) the values coincide. Center: Always for n = 128, we use the RDM compatibility constraints to obtain the half-system RDM from the ground state found by the variational method. This enables us to compute the entanglement entropy, and to characterize the phase diagram of the model [40]. Right: Half-system entanglement entropy for γ = 0 and different values of n obtained from the variational method using the RDM compatibility constraints in the symmetric basis. For γ = 1, one can appreciate the anomaly in h = 1 as we increase n, which becomes a critical point in the asymptotic limit.

Ising chain with variable-range interactions under a transverse field
As a second example, we consider the Ising model with variable-range interactions in a transverse field. The tuneable interaction range allows us to explore how our variational method performs as the range of interactions decreases from the infinite-range case (equivalent to the LMG model) to the nearest-neighbour case. In particular, we consider the following Hamiltonian for an Ising chain with decaying power-law interactions: where J ij = |i − j| −α , the paramether α tunes the range of interactions, and θ > 0 (θ < 0) results in AFM (FM) interactions. We note that in the limit α → 0 all pairs interact with the same strength [45], whereas in the other extreme α → ∞ we have interactions only between nearest neighbours. The phase diagram for this model has been extensively characterised [46][47][48]. In particular, for α > 0 the model exhibits three phases: an ordered ferromagnetic phase for −π/2 ≤ θ < θ − c (α); a disordered paramagnetic phase for θ − c (α) < θ < θ + c (α); and an ordered anti-ferromagnetic phase for θ + c (α) < θ ≤ π/2. Notably, such a model has also been shown to display Bell correlations at the critical points for the ferromagnetic couplings [49].
In order to construct the variational method for Equation (22), we consider the SDP in Equation (21) with the effective In Figure 2 we compare the ground states obtained using our variational method (VM) with those obtained from exact diagonalisation (ED), in terms of relative energy and fidelity. To compare the ground state energies we look at their ratio, E VM 0 /E ED 0 . For this, a few comments are in order: First of all, note that the ground state energy is negative and sufficiently far from zero to constitute a good approximation of 1 − δ, where δ is the relative error. Second, we have chosen the ratio as a figure of merit, instead of the relative error, as it gives a better visual comparison with the fidelity also plotted in Figure 2. To compute the ground state fidelity we use the definition F (ρ ED , ρ VM ) := Tr √ ρ ED ρ VM √ ρ ED 2 . For this model we expect that the analytical solutions for the long-range interaction regime are well approximated by the one for the LMG model, where our variational method yields exact results. However, in the transition from α 1 to α 0 the quality of the approximation is a priori not so clear. In Figure 2 it can be appreciated that for α > 0 the method fails to capture the anti-ferromagnetic phase θ + c (α) < θ, eventually yielding fidelities close to zero. On the other hand, for values θ < θ + c (α) the method provides a good approximation, even though for sufficiently large α the method is less accurate near the critical point θ − c . In Figure 3 we compare half-system entanglement entropies obtained from the VM and ED. As expected, we observe a discrepancy in the anti-ferromagnetic phase, because of the little overlap between the ground state and the symmetric space in  (22) with n = 10, compared to exact diagonalization. Left: energy ratio, Right: ground state fidelity. As expected, the case α = 0 is faithfully recovered by the variational method. However, we observe that as the value of α increases (the range of interaction decreases) the variational method fails to capture the anti-ferromagnetic phase for which the fidelity eventually drops to zero. The ground state energy and fidelity in the ferromagnetic and paramagnetic phases are well approximated, although for large values of α there is a little discrepancy near the critical points θ ≈ θ − c (α).
that region. Interestingly, we see from Figure 3 that in the ferromagnetic regime one can use the VM (without ED) to approximate its phase transition (between paramagnetic and ferromagnetic) for different number of particles, and extrapolate its asymptotic limit. Such an approximation naturally works better as the range of interactions increases. The VM has the potential to identify a phase transition in the model due to the following argument: One might expect that in many cases of interest, and for a finite number of particles, the sudden change in nature (e.g. symmetry) of the ground state is actually "smeared" rather smoothly around the critical point. Therefore, everywhere around this point the ground state may still have some overlap with the symmetric space, which is the one considered by our variational method. The numerical evidence presented in Figure 3 supports this conjecture. Although the paramagnetic to anti-ferromagnetic phase transition is not reproduced by the VM, one observes that the transition from ferromagnetic to paramagnetic is well approximated by looking at the discrepancies for values of α 1.5. Right: Half-system entanglement entropy for α = 2. Note that the transition from ferromagnetic to paramagnetic manifest itself in the vicinity of θ −0.1π for n = 10. From this observation we conjecture that our VM can be used to extrapolate some critical points where phase transitions occurs. We remark that the ED has not been used in this case and that the behavior observed arises from the VM alone.

Ising chain for nearest neighbours interactions
So far, in Section IV A 1 we focused on the extreme case of infinite-range interactions (the LMG), and in Section IV A 2 we explored how our variational method behaves as we decrease the range of interactions. Let us now investigate what happens in the other extreme case: an Ising model with nearest-neighbours interactions in a transverse field. The Hamiltonian we consider is:  (23) with n = 10. Left and center: comparison of the ground state energy and fidelity with respect to results from exact diagonalization. We observe that energies disagree mostly in the AFM case, due to its inherent asymmetry. The FM case is in general well approximated, with only some minor discrepancies near the values that hint at a critical point in the asymptotic limit. Right: The VM is used in the FM case (Jz = 1/2) in order to investigate the scaling of the half-system entanglement entropy, hinting at the existence of a critical point when extrapolating to the asymptotic limit.
where J z > 0 (J z < 0) corresponds to FM (AFM) coupling, and h tunes the transverse field strength. For Equation (23), the variational method is taken with the effective HamiltonianH : the previous sections, in Figure 4 we compare the VM with ED for low number of particles. As expected, the VM yields almost orthogonal solutions in the AFM region, while it provides fidelities close to unity in the FM region. Still, slight discrepancies arise in the FM case. In particular, one observes that in the vicinity of what is a critical point (in the asymptotic limit) the fidelity drops, nevertheless still providing a good upper bound to the ground state energy (see Figure 4). We remark that, by using the VM to determine the half-system entropy scaling in the FM case (J z = 1/2), we observe an anomaly at h ≈ 1 (see Figure 4), which signals the presence of critical point. Therefore, despite the discrepancy in fidelity, the FM critical point can still be well approximated by our VM.

XXZ model under a transverse field
In Section IV A 3 we have looked at a spin system in one dimension and with nearest-neighbour interactions. This model is actually solvable via Jordan-Wigner transformation [50], as it can be equivalently described as a system of free fermions (see e.g. [39]). Here we consider the validity of our variational ansatz beyond the free fermion scope, in an XXZ spin chain with an homogeneous magnetic field in the X direction: where we take J = 1 (J = −1) for the FM (AFM) couplings, ∆ marks the anisotropy (with ∆ = 1 corresponding to the isotropic case, the XXX model) and h tunes the transverse field strength. In this case, the effective Hamiltonian for the variational method In Figure 5 we compare the ground states obtained from our VM with those obtained from ED, in terms of relative energy and fidelity. For this case we observe that the VM provides a faithful approximation for values ∆ 1 (∆ −1) when considering FM (AFM) interactions. In particular, in such a regime the VM yields exact results except around a line which likely corresponds to critical points in the asymptotic limit [51,52].

Ferromagnetic XXZ with periodic boundary conditions
Let us now consider the following instance of an XXZ model: a periodic anisotropic ferromagnetic spin-1/2 chain, placed in an homogeneous magnetic field in the z direction. This model is described by the Hamiltonian where σ (n+1) = σ (1) , J x , J z ≥ 0 are the exchange coupling constants, and b tunes the strength of the external magnetic field. In Ref. [53] it was investigated the fidelity of preparation of Dicke states as ground states of Equation (25). In Table I we show the ground state distribution predicted with the perturbative results of [53]. Our variational ansatz comes as a natural tool to benchmark the fidelity of the preparation of Dicke states as the ground states of the Hamiltonian in Equation (25). In Figure 6 we see that, up to n = 10, the fidelity remains > 85%, which is consistent with the predictions of [53]. For this case, in the VM we take the effective HamiltonianH : Table I. Ground state for the model in Equation (25), as a function of the magnetic field b and the coupling parameter ∆J = Jx − Jz, according to the perturbative results presented in [53]. In Figure 6 we recover and strengthen the result.
The perturbative prediction from [53] splits the phase diagram among different regions, each having substantial overlap to a different Dicke state. However, there are already some discrepancies observed: In between these regions, the approximation with Dicke states does not have good overlap with the ground state (see Figure 1 in [53]), however it has good overlap with other basis Dicke states (see Figure 2 in [53]). Here, our method enables us to understand this discrepancy from a different perspective: in Figure 6 we see that the regions in which the perturbative approach fails actually correspond to a good overlap for the first excited state and a Dicke state. The low-end of the spectrum of the Hamiltonian considered has a good overlap with the symmetric space. However, it can happen that the VM chooses to approximate the first excited state instead of the ground state if the energy obtained becomes more favorable. This depends on both the overlap with the ground space and the energy gap of the Hamiltonian. Let us denote E 0 and E 1 the ground and first excited state energies, respectively. Let us also denote F 0 (F 1 ) the fidelity between the ground state (first excited state) and the symmetric space. The discontinuities may happen when F 0 E 0 = F 1 E 1 . Indeed, in Figure 6 we observe that, while the energy ratio is smooth, the fidelity may suddenly drop to zero or jump to almost one due to the above mentioned reason. For the larger n limit, the ground state of Equation (25) can be found exactly using matrix product states (MPS) and the density matrix renormalization group (DMRG) algorithm. In order to compute the fidelity of our variational solution with respect to the exact one, the most straightforward way is to contract the MPS representation of the ground state with the solution of the variational method. However, the latter is not given in a MPS form and it therefore has to be converted in a MPS form. In order to do so, we need to establish a correspondence between arbitrary superpositions of Dicke states and MPS. To the best of our knowledge, in general such a correspondence has not been established before, and this result may be of independent interest. Therefore, we have devoted a full section (Section IV C) to it.

Many-body spin-1 Hamiltonian with collective interactions
Finally, let us illustrate that our proposed variational method can be also easily applied to systems with local Hilbert space of dimension d > 2. To do so, we consider the three-orbital Lipkin-Meshkov-Glick Hamiltonian [54] (equivalently, the generalisation of the Lipkin Hamiltonian as proposed in [55]). Similarly to Section IV A 1, the variational ansatz is expected to recover exact results due to the long-range interactions resulting in ground states with permutation symmetry. The model is constructed by n identical but distinguishable three-level atoms, and it is also commonly used in nuclear shell models. It can also arise for three-level atoms collectively coupled to electromagnetic field modes of a cavity. Concretely, the model is described by the Hamiltonian: ij with τ ij = |i j| for i, j = {0, 1, 2}. In this case the effective Hamiltonian used in the SDP is In Figure 7 we show that the variational ansatz reproduces exactly the ground state energy, as expected. Furthermore, we use the compatibility conditions to obtain the half-system entanglement entropy, which is useful to provide insights about the phase diagram of the model.  [44]). Center: Half-system entanglement entropy obtained through the compatibility conditions Equation (8), capturing features of phase transitions. Right: Half-system entanglement entropies for different number of particles. The scalability can be used to extrapolate the peak anomalies to the asymptotic limit.

Benchmarking performance with existing methods
One of the most appealing features of the variational method proposed here is its ability to yield results for large system sizes with very modest time and memory requirements (and, consequently, energy consumption). Indeed, since the computations take place in the symmetric space, its dimension grows only polynomially with the system size. In particular, it is linear for qubits, quadratic for qutrits, etc. In the previous sections we have argued that the method yields results that capture traces of some quantities of physical interest. Therefore, it can be a good candidate to a first order exploration before trying more numericallyintensive results. To make this comparison quantitative, we here benchmark the computational requirements of the method with other existing techniques.
In this section we briefly comment on the time, memory and energy consumption devoted to the variational ansatz. The runtime of the variational method can be split in two steps: 1) to precompute the A matrices in Equation (17) for a fixed n, m and d and 2) to load and solve the SDP. The most expensive task, both in time and memory, is to compute the compatibility constraints. Hence, in order to agilitate the process, one would first preallocate and store the compatibility constraints for a fixed number of particles n, of local Hilbert space dimension d and with RDMs of size m. Then, once the compatibility constraints are preallocated, one can scan the phase diagram of the desired parametrized Hamiltonian model just by loading the compatibility constraints and proceeding to solve the corresponding SDP.
In Figure 8 we present a representative sample of the computing runtimes we have observed in order to prellocate the compatibility constraints of the 2-RDMs and half-system n/2-RDMs for different number of particles n and different local Hilbert space dimensions d. We have considered the 2-RDMs compatibility constraints case both in the computational basis (Equation (7)) and the symmetric basis (Equation (8)). Apart from requiring less memory, one observes a clear advantage in runtime when obtaining the compatibility constraints directly in the symmetric basis, as expected. Therefore, for the variational method it is desirable to project the effective Hamiltonian onto the symmetric subspace. For the n/2-RDMs case, we have considered only the symmetric representation which decreases the memory storage limitations in order to find, for instance, the half-system entanglement entropies. An additional comment is in order: For a constant value of d, note that the multinomial coefficients in Equation (7) or Equation (8) do not require a full expansion of the factorials, but there exist closed analytical formulas for them (see e.g. [27]). This has been taken into consideration in our calculations. Furthermore, it is desirable to apply such closed expressions, not only for speed, but more importantly for numerical stability issues (quotients of factorials of large numbers may give problems in floating-point arithmetic if these numbers are of the order of ∼ 100).
Number of Qudits

Runtimes [sec]
n/2-RDM, Symmetric basis Figure 8. Runtimes observed in order to prellocate the 2-RDMs and n/2-RDMs compatibility constraints in the computational and symmetric basis. Apart from the memory storage advantage, it is observed that the symmetric spaces offers a significant advantage also in runtimes. The computations have been carried out on a 64-bit operating system with 32GB ram and a 3.70GHz processor. No parallelization has been used, which can be easily be implemented, significantly speeding up the process. The runtimes might slightly vary at each run and are not meant to be taken as exact, but as an illustration of their order.
In Figure 9 we present some of the computing runtimes in order to load the constraints and solve the SDP for the Ising chain with decaying power-law interactions previously considered in Section IV A 2. We have considered the constraints and effective Hamiltonian in the symmetric basis, and in order to solve the SDP we have set the solver SDPT3 [56] to its maximal precision providing a numerical error up to O(10 −14 ) when the variational ansatz can reach the exact solution. We have carried out the comparison with the solution provided by DMRG. In order to find the fidelity between the DMRG solution and the VM solution, we have used the auxiliary results developed in Section IV C in order to represent the VM solution as a translationally invariant diagonal matrix product state.

B. Bell non-local correlations
The proposed variational ansatz is also convenient to investigate Bell non-local correlations in many-body systems. The permutational symmetry naturally synergizes with the so-called 2-body Permutation Invariant Bell Inequalities (PIBIs) presented in [25]. These type of Bell inequalities involve at most 2-body correlation functions, and some of them are violated by symmetric states [27]. One can now consider two approaches: On the one hand, to obtain the quantum state that gives the maximal violation of such equalities within the variational ansatz. On the other hand, to find quantum states that also have Bell correlations by using the VM to appoximate the ground state of a many-body Hamiltonian. In the latter case, the Hamiltonian considered need not Figure 9. Numerical results of the VM for the Ising Hamiltonian Equation (22) with n = 64, compared with DMRG. The DMRG algorithm follows [49,57,58] increasing the bond dimension up to 20. Left and Center: Energy ratio and fidelity of the VM solution with respect to the DMRG solution. In order to compute the fidelity we use the result in Section IV C to transform the symmetric basis representation of the VM solution into a matrix product state representation. Right: Runtimes comparison to achieve convergence with DMRG and VM method. The VM is significantly faster, making it a good candidate for a first rough exploration of large phase diagrams, and to upper bound ground state energies. Note that runtimes might slightly vary at each run. The total runtimes are 525.04s (1s preallocation) for the VM and 2854.25s for the DMRG. With a power consumption of 425W in our workstation, amounts to an environmental impact of around 18g vs 100g of CO2 into the atmosphere for the VM vs the DMRG methods (we have taken the 0.296 EU coefficient of kWh to kg of CO2 given by the European Environment Agency for 2016). We remark that the DMRG is an extremely optimized and efficient method that cannot be applied beyond 1 geometric dimension. In these cases, the benchmark with existing methods would be separated by even more orders of magnitude.
correspond to the Bell operator [39]. It is worth mentioning that the measurement settings might need to be optimized in order to increase the visibility of the Bell correlations.

Optimizing permutationally invariant two-body Bell inequalities
We first focus on two particular classes of 2-body PIBI. These inequalities satisfy the following condition for all correlations that can be described under local-realism assumptions (meaning that their violation signals the presence of non-local correlations): and (n mod 2)(n − 1)(nS 0 + S 1 ) + n 2 S 00 + nS 01 − S 11 ≥ n 2 (n + 2 + n mod 2), l are the one-and two-body symmetric correlators with M (i) k denoting the measurement in direction indexed by k = {0, 1} corresponding to particle i. The first inequality Equation (27) is particularly fitted to detect non-local correlations in superpositions of Dicke states, while the second inequality Equation (28) is tailored to detect non-local correlations in half-filled pure Dicke states [27,59].
In order to know if there exists a quantum state that violates such Bell inequalities, one still needs to find appropiate n pairs of measurement settings. This gives rise to the so-called Bell operator, a quantum observable in the n-partite Hilbert space whose expectation value with respect to a quantum state corresponds to the value of the Bell inequality for the chosen measurement settings. In general, finding such measurements consists in a very demanding non-convex optimization problem with no trivial solution. However, since the variational ansatz provides a global state, the complexity of the problem gets greatly reduced when we restrict the optimization to the case where all the measurement settings are the same for each party, in the same reference frame. This gives rise to a permutationally invariant Bell operator, whose extremal values within the symmetric space we be found using our variational method.
Contrary to previous approaches (see e.g. [27,60]), where one could use representation theory methods such as Schur-Weyl duality to project the permutationally invariant Bell operator onto the different symmetric blocks (see also [61]), here the VM circumvents this intermediate step: It is enough to consider the effective Hamiltonian where σ can be the 2-RDM obtained with the variational ansatz. We parametrize the measurements as M k := sin (θ k ) σ x + cos (θ k ) σ z , where k ∈ {0, 1} and σ x , σ z are the Pauli matrices, and use the VM to find the symmetric state minimizing the energy of the effective Hamiltonian. We note that this approach becomes particularly useful in the case of large d, since the number of blocks arising from the symmetry-adapted basis increases with d. We also remark that, since one can always apply a dual U ⊗n symmetry to both state and measurements without departing from the symmetric space, it is enough to optimize over the difference between measurement directions, e.g. θ 0 − θ 1 . Furthermore, since Equation (27) and Equation (28) have two inputs and two outputs per party, Jordan's lemma guarantees that using d = 2 is sufficient to find find its maximal violation.

Looking for Bell correlations in a direction specified by a Hamiltonian
Here we propose a two-step process to find Bell correlations in the ground state of Hamiltonians of physical interest (e.g. an XXZ chain). First, we use the VM to do a quick scan over the parameter space of a given Hamiltonian family, in order to find potential candidates whose ground state might display Bell correlations, see Figure 10. If we found Bell correlations with the VM, we would have obtained a symmetric state displying them, albeit we have no guarantee of the fidelity with the ground state of the model at this point. Then, we narrow down the search by computing the actual ground state with other methods, such as DMRG, also in Figure 10. As we show in Figure 10, we observe for the first time, to the best of our knowledge, that the ground state of such an XXZ chain violates the Bell inequality Equation (27), thus exhibiting Bell correlations, in the corresponding parameter regime. Bear in mind that in order to carry out the measurement optimization one can no longer consider a single 2-RDM as we posed in Equation (29), but one has to sum over all the different 2-RDMs σ ij obtained with the DMRG: Note that, since the RDMs are fixed in this case (by the exact solution), the measurement directions should not be restricted to the XZ plane in the Bloch sphere, but allowed to point in any direction. Furthermore, the measurement settings for each party do not need to coincide. Similarly, we have seen in Section IV A 5 that pure Dicke states provide a good approximation of the ground state for such ferromagnetic XXZ chain with periodic boundary conditions and longitudinal magnetic field b. As we have previously mentioned, inequality Equation (28) is tailored to half-filled Dicke states which happen to approximate the ground state in the range Therefore, in such region we expect to witness non-local correlations with Equation (28). Indeed, in Figure 11 we show the witnessed nonlocality where we have used the variational ansatz to approximate the ground state and optimize using Equation (29). We also observe that nonlocality detection goes beyond the specified region where the half-filled Dicke state approximates the ground state. As previously discussed in Section IV A 5, such extra range of nonlocality detection seems to arise from the variational ansatz approximating the first excited state instead of the ground state.

C. Generically expressing symmetric states as TI diagonal MPS
In this section we present an analytical method to generically represent any n-qubit Dicke state with a translationally-invariant (TI) matrix product state (MPS). More precisely, the goal is that, given a state of the form Some representations of important symmetric states have been known since their inception. For instance, the GHZ state (up to normalization) can be generated with D = 2 [62] using the following TI MPS: On the other hand, for the |W state, there exists no TI representation with bond dimension 2 [63]. However, it does admit (up to normalization) the following non-TI representation [63]: We observe that in Equation (34), all the coefficients of the MPS are either 1 or 0. Indeed, some MPS can be used to represent Boolean formula solutions [64]. More generally, the representability of quantum states with MPS of a particular form has deep connections with modern algebraic geometry [65][66][67]. For instance, the |W state can be arbitrarily well-approximated with a diagonal TI MPS of bond dimension D = 2: where The n-qubit |W state can be obtained as the limit of ε → 0 of the TI MPS given by Equation (35).
a. Proof: One simply notes that, since the MPS is diagonal, the element corresponding to the physical index (i 1 , . . . , i n ) with k := j i j is given by which amounts to 0 if k ≡ 0 mod 2 and ε (k−1)/(n−1) if k ≡ 1 mod 2. Hence, noting yields the result. Inspired by Lemma 1, we propose now a TI MPS of bond dimension n to approximate generically any superposition of Dicke states of the form Equation (31). We propose the following parameterization of A 0 and A 1 : For simplicity let us denote A 0 = y1 and k = j i j . It is then easy to see that Equation (32) leads to the following system of equations: Note that for states of the form of Equation (31) we only require n + 1 equations. Hence, it is natural to choose D = n. This motivates the following Lemma: Lemma 2 Consider the system of Equations where z 1 , . . . , z n ∈ C. The solutions of Equation (41) are the roots of the polynomial P (X) defined in Equation (A10).

Lemma 3
The system of equations that determines the coefficients y and x 1 . . . x n of the diagonal tensors A 0 and A 1 , used to represent the linear combination of Dicke states Equation (31), is given by The value of y is readily determined from the first equation, and the x's are found by finding the roots of the polynomial P (X) constructed from Lemma 2 using the remaining set of equations.
Note that generically, we will have n complex solutions, up to n! permutations. However, there is the possibility that some of the solutions lie at infinity in some pathological cases. Nevertheless, these cases form a zero-measure set which can in practice be avoided by adding an ε-perturbation to d k , in the same spirit as in Equation (36).
It is now clear that the bottleneck is solving Equation (41) in Lemma 2. We propose two approaches in order to do so: In this section, a variant of the Faddeev-Leverrier algorithm to solve Newton's identities in order to find the roots power-sum symmetric polynomials. In Appendix A, we propose a step-by-step computation of the solutions via Gröbner basis which could provide solutions for more general systems of equations.
We make two observations in order to solve Equation (41). The first one is that, if we consider a matrix A with eigenvalues {x 1 , . . . , x n }, then Equation (41) can be thought of as The second observation is that there exists a way to express the characteristic polynomial of a matrix A in terms of Tr[A k ]. Indeed, if P (X) = det (X1 − A) = (X − x 1 ) · · · (X − x n ) = n k=0 c k X k is the characteristic polynomial of A, then P (A) = 0, and taking the trace in both sides yields such an equation. The Faddeev-Leverrier algorithm provides an easy to compute form for the coefficients c k : They are given by The values of x are found by finding all the roots of the polynomial P (X) = n k=0 c k X k with c n = 1. In Appendix A we give an alternative approach to solve Equation (41) from a more algebraic point of view that gives more insight to the cominatorial structure underlying Equation (41).
Having an efficient way to represent a state of the form Equation (31), we can now use its MPS form to efficiently compute the fidelity of DMRG solutions, which are already given in the MPS formalism, for large n, thus being able to benchmark our method.
D. Determining which symmetric states cannot be self-tested from their marginals Self-testing is one of the most stringent protocols in the paradigm of device-independent quantum information processing. Self-testing consists in inferring, solely from the statistics of a Bell experiment, which quantum states and measurements are being used and performed, respectively [68][69][70]. Much of the existing work has been centered around the bipartite case [71,72], partly motivated by its more accessible physical implmentations [73,74], but also motivated by its more accessible theoretical analysis, exploiting in most cases properties of the maximally entangled state of two qudits [75][76][77]. In the multipartite case, the analysis becomes more complicated, although some ideas for the bipartite inspired some extensions [36]. In the multipartite case, symmetric states constitute a natural candidate to begin their study: for instance, the robust self-testing of the W state (|D 1 ) [78] inspired schemes to self-test Dicke states of the form |D k [36,37,79]. Nevertheless, these schemes use full-body correlators and require individual addressing, thus being less appealing from an experimental point of view. Therefore, some studies have been carried to find out whether self-testing is possible using only marginal information [80] (see also [81,82]): In [80], some efforts showed that the three-qubit states maximally violating some of the translationally invariant, two-body Bell inequalities from [83] could be self-tested using two-body correlators, thus giving a positive answer to this question.
Interestingly, the question of how much information from the statistics is needed (i.e., how many parties can one trace out) in order to self-test a quantum state is still open. In this section, we aim at showing how our method can be used to guarantee a negative answer to the previous question and to give evidence towards a positive answer as well, depending on the uniqueness of the solution to Equation (17). More precisely, we show how our method can be used to determine which symmetric states could potentially be self-tested from marginals and which symmetric states definitely could not, because their marginals do not have a unique (modulo local unitaries) extension in the symmetric space.
Let us first consider an n-qubit density matrix being a projector onto a n-qubit Dicke state |D k as in Equation (1). We shall denote it ρ n,k . Let σ = Tr 1 (ρ) be the resulting density matrix from tracing out a single particle. In virtue of Equation (7), we have To gain some intuition, let us first study under which conditions is it possible to show that the purification of σ is unique. Following the spirit of Lemma 5.2 of [84], we begin by considering a purification with an auxiliary system of the form where the Dicke states act on n − 1 qubits and It follows from elementary algebra that the (n − 1)-body RDM of |Φ Φ| is equal to σ if, and only if, Hence, in this case there exists a purification, it is unique, and it must be of the form Corollary 1 The (n − 1)−partite reduced state of |D k uniquely determines |D k in the symmetric space.
In Appendix B we show how the above example can be generalized to tracing out any number of parties. The uniqueness of the extension is in one-to-one correspondence to the uniqueness of a linear program (see Equation (B9)). We have numerically observed that such a solution is unique if we trace out up to n − 2 parties for a basis Dicke state. However, it is not a priori clear how generic the above property is. A more in-depth study suggests that generically, the uniqueness property depends on both the rank of the global density matrix and the number of parties traced out. In Figure 12 we provide numerical evidence that generically the uniqueness of the symmetric extension depends on the number of particles n, the number of parties in the RDMs m and the rank of the global density matrix ρ. We have followed the procedure below: 1. Generate a random symmetric state whose density matrix has a given rank 2. Use the compatibility conditions to obtain its RDM 3. Choose a random Hermitian matrix A and find a new global state compatible with the RDM by making use of the SDP in Equation (17), but with objective function − A, ρ and check its fidelity with the original global state. The matrix A forces the SDP to explore the feasible set in the direction given by A 4. We repeat step 3) a sufficient number of times (in our case, 100 times). If the fidelity remains always one up to numerical accuracy error, this is strong numerical evidence that the global state is unique. On the other hand, if the fidelity falls below one in some case, this indicates that too many parties have been traced out and the RDM would no longer be sufficient to self-test the original state, since it does not have a unique global extention of size n We note that, although mixed states in their generality cannot be self-tested, the fact that for some rank configurations and sizes of the RDM the extension to the symmetric state seems to be unique could open the door to a weaker form of self-testing, under the assumption that the global state is symmetric. Numerical results on the dependence between n, m, rank(ρ) for an m-RDM to have a unique symmetric extension ρ of n = 15 and n = 30 qubits respectively with a given rank(ρ). For each case we have carried out 100 trials forcing the SDP to explore the feasible set in a random direction A at each trial. The black squares correspond to the configurations for which the recovered global symmetric state has fidelity > 0.9999 with the original global state for 100% of the trials, thus providing evidence of having a unique symmetric extension. The numerical tolerance has been set to take into account the imprecision of the SDP solver. For the non-black squres, some of the trials have exhibited a fidelity < 0.9999. For those cases, we show the minimal fidelity obtained out of all the trials as a way to illustrate the tolerance. One clearly observes a certain correlation between size of the RDM m and rank(ρ), showing more chances to have a unique extension for low rank(ρ) by tracing out few particles.

V. CONCLUSIONS AND OUTLOOK
In the present work we have presented a study of the QMP restricted to symmetric states. We have provided analytical compatibility conditions for an m-qudit RDM σ to be compatible with an n-qudit global symmetric state ρ. We then use said compatibility conditions to answer the question of whether a given reduced density matrix σ is compatible with a global symmetric state ρ by turning it into a feasibility problem efficiently solvable via an SDP. Our results have implications in different fields. We have explored some of them in several case-studies: • We have developed a computationally efficient variational optimization method to upper bound the ground state energy of any local Hamiltonian. This method considers the resulting marginals to be compatible with a global symmetric state in order to carry out the optimization by means of SDP with the compatibility conditions as constraints. In order to benchmark the VM, we have considered several paradigmatic Hamiltonian spin models, that go from long-range to nearestneighbor interactions. In general, we observe that the VM provides a good upper bound for ferromagnetic and long-range interactions, yielding exact results in the infinite-range limit; while it misses to capture antiferromagnetic short-range interactions, where the ground state has poor overlap with the symmetric space. We have also used the compatibility conditions in order to obtain the half-system entanglement entropy in the symmetric space, which is a insightful quantity for many-body systems. Remarkably, we have observed that for some cases our variational method can also be used to approximately locate phase transitions. The reason for that lies that the change of properties of the ground state in a phase transition also maifest, to some extent, in the symmetric space projection and are therefore captured by our method. Another observed feature is that for some specific models the VM has recovered the first excited state, instead of the ground state, in some regions of the phase diagram. Finally, we have observed a significant speed advantage of the VM compared to a typical DMRG algorithm. The advantage of the VM lies on the low memory storage required and high speed, making the VM a suitable candidate for a first order exploration of large sets of parameters characterizing the phase diagram of spin Hamiltonians. While we have only considered qubits, qutrits and chain configurations, our VM is straightforwardly applicable to any qudit and lattices of arbitrary geometry and dimension. We leave open to implement and explore the VM in corresponding cases of interest.
• We have considered the VM in the context of Bell non-local correlations. In particular, we have explored its synergy with the so-called 2-body permutationally invariant Bell inequalities. The results in this context are two-fold: First, we have shown how the variational method comes as a natural tool to optimize a multipartite 2-body PIBI in order to find whether the inequality detects non-local correlations; Second, we have used the low computational cost of the variational method to look for non-local correlations in a spin-1/2 XXZ chain under a transverse field, narrowing the parameters to be considered and eventually leading to the detection of non-local correlations the ground state with n = 128 parties. We have also considered another spin-1/2 XXZ chain this time with periodic boundary conditions with longitudinal magnetic field, detecting non-local correlations on its ground state and first excited state of a specific phase. The tool we have here presented can be readily used in the context of Bell correlation depth [85] or DI entanglement depth certification [86,87] in the context of 2-body PIBIs.
• We have developed an analytical methodology to derive a translationally invariant diagonal matrix-product state representation of bond dimension n for pure symmetric states. This result is generic, and could be of independent interest. For our purposes, we have used it to transform the symmetric state solution obtained with the VM into a translationally invariant, diagonal, matrix product state, allowing us to check its fidelity with the DMRG solution.
• Finally, we have shown how the compatibility conditions can be used to determine which symmetric states ρ cannot be self-tested solely from their marginals. Remarkably, we present numerical evidence suggesting a correlation between the size of the global state n, its rank rank(ρ), how many particles m remain in the observed RDM and the uniqueness of a symmetric global state. This uniqueness property could open the way to a weaker form of self-testing, that usese the assumption that the global state is symmetric.
One may wish to impose the condition of the global state being pure; i.e., ρ being a rank-1 projector. Although this condition clearly breaks the convexity of the SDP program, one can combine the VM with an iterative projection onto the principal component of ρ [83,88] in order to converge to a rank-1 solution. Furthermore, one may also wish to explore the role of different symmetries in the SDP. Whether there exists a SDP invariant formulation of our problem [61,89] that could allow it to be formulated for other symmetry groups is unclear and we leave it for future research. In a following work, we shall investigate a variations of the VM in order to perform tomography/fidelity estimates with respect to a target symmetric state. This is of wide experimental relevance, as in the case of Bose-Einsten condensate where only partial information (e.g. not an informationally complete set of measurements) is available [26,90].
i.e., the roots of P correspond to the values of x i , up to a permutation.
Therefore, we only need to find the general form of P (X). The coefficients of P (X) are closely related to the partitions of n. Let us define the following: Definition 1 Let λ m denote a partition of m; i.e., λ = (λ µ1 1 , . . . λ µ k k ) where k i=1 µ i λ i = m and λ i > λ i+1 with λ i , µ i ∈ N.
We define the polynomial where We define by convention Q 0 := 1.
Note that λ m |ξ λ | = m! since ξ λ counts (with sign) the number of permutations of m elements of cycle type λ. In addition, we remark that the number of partitions p(m) of a given integer m scales as log p(m) ∼ C √ m, where C is a universal constant. This makes the sum Equation (8) prohibitive to evaluate already for modestly large values of m. However, as shown in Section IV C, it is possible to efficiently compute Q m (z) without splitting it into its different summands. Now that we know how to obtain the x that satisfy z in Equation (41), let us turn to the system of equations that actually arises from Equation (32). Note that Equation (41) does not take into consideration the z 0 term, but we incorporating the condition that A 0 ∝ 1. Then the system of equations of interest becomes Equation (42).
The system of equations Equation (41) is also known the the power sum ideal. Its reduced Groebner Basis is found as the elimination ideal of the power sums.

Corollary 3
The elimination ideal of the power sums gives the compatibility conditions on the weights d k of Equation (31) to be representable with a diagonal TI MPS of bond dimension D < n.
Indeed, let us consider n = 4 and D = 3. The elimination ideal of the power sums of three variables and degree four is 3 − z 4 ∩ K[z 1 , z 2 , z 3 , z 4 ] = q(z 1 , z 2 , z 3 , z 4 ) , (A11) where q(z 1 , z 2 , z 3 , z 4 ) = z 4 1 − 6z 2 1 z 2 + 3z 2 2 + 8z 1 z 3 − 6z 4 (A12) Note that the compatibility polynomial in Equation (A12) is precisely the same polynomial Q 4 (z) in the constant term of the univariate polynomial in Equation (A6). Hence, all the symmetric Dicke states for which the z obtained from their d k belongs to the elimination ideal of the power sums of D variables with degree n are representable as a diagonal TI MPS of the form A 0 ∝ 1 and A 1 = diag(x).
-If we trace out n − m parties, then we generalize the last two points: For α > k or α < k − (n − m) the condition on the diagonal is n−m p=0 ρ α+p α+p ξ α+p = 0, which implies ρ k+1+p k+1+p = ρ k−1−p k−1−p = 0 for p ≥ 0. This condition is nontrivial as long as the number of equations (m − 1) is greater than the number of nonzero left hand sides (n − m + 1), i.e., whenever m > n/2. Therefore, the condition ρ 0 implies that all the off-diagonal elements must be zero and therefore ρ k k = 1.
where we have defined for simplicity x p := ρ p p n k / n p . If m > n/2, there must necessarily be zeroes in the right hand side of Equation (B9).
The question about uniqueness of solutions of linear programs [91] and semidefinite programs [92,93] is an intensive field of research, due to its connection to rigidity theory. For instance, the general solution to the uniqueness of Equation (17) can be expressed via Theorem 2 [92] If ρ is a max-rank solution of Equation (17), and we write ρ = L † L, where L ∈ C r×n , then ρ is the unique solution of Equation (17) if, and only if, the kernel of the linear space spanned by L † A α β L is trivial.
Corollary 4 [92] If all the solutions to Equation (17) share the same rank, then the solution must be unique.