Multiparticle entanglement criteria for nonsymmetric collective variances
Oliver Marty, Marcus Cramer, Giuseppe Vitagliano, Geza Toth, Martin B. Plenio
MMultiparticle entanglement criteria for nonsymmetric collective variances
O. Marty , M. Cramer , , G. Vitagliano , G. T´oth , , , and M.B. Plenio Institut f¨ur Theoretische Physik, Universit¨at Ulm, D-89081 Ulm, Germany Institut f¨ur Theoretische Physik, Leibniz Universit¨at Hannover, D-30167 Hannover, Germany Institute for Quantum Optics and Quantum Information (IQOQI),Austrian Academy of Sciences, A-1090 Vienna, Austria Department of Theoretical Physics, University of the Basque Country UPV/EHU, P.O. Box 644, E-48080 Bilbao, Spain IKERBASQUE, Basque Foundation for Science, E-48013 Bilbao, Spain Wigner Research Centre for Physics, Hungarian Academy of Sciences, P.O. Box 49, H-1525 Budapest, Hungary (Dated: August 24, 2017)We introduce a general scheme to detect various multiparticle entanglement structures from global non-permutationally invariant observables. In particular, we derive bounds on the variance of non-permutationallyinvariant and collective operators for the verification of k -party entanglement. For a family of observables re-lated to the spin structure factor, we give quantitative bounds on entanglement that are independent of the totalnumber of particles. We introduce highly non-symmetric states with genuine multipartite entanglement that isverifiable with the presented technique and discuss how they can be prepared with trapped ions exploiting thehigh degree of control in these systems. As a special case, our framework provides an alternative approachto obtain a tight relaxation of the entanglement criterion by Sørensen and Mølmer [Phys. Rev. Lett. 86, 4431(2001)] that is free from technical assumptions and allows to calculate the bounds with an improved scaling inthe detectable depth. PACS numbers: 03.65.Ud, 03.67.Ac
I. INTRODUCTION
One of the most fascinating challenges in quantum infor-mation science is to explore the prospects of quantum effectsto go beyond the capabilities of classical physics. An exam-ple is the concept of spin-squeezing that describes a collectiveproperty of an aggregation of spins [1, 2]. Originally, it hasbeen introduced to achieve performance gain using quantummetrology. The basic concept can be quantified in a multi-tude of spin-squeezing parameters [1–3]. The role of spin-squeezing in the context of quantum improved measurementscan be illustrated graphically, providing an intuitive picture ofthe spin-squeezing parameters [2, 3].On the other hand, a different central application of spinsqueezing parameters is the detection of many-body quan-tum correlations: strongly (anti-)correlated spins in squeezedstates exceeding the standard quantum limit are required tobe entangled as has been observed in [4–6]. Quantitativelythe degree of spin-squeezing is a measure of multiparty en-tanglement, i.e. in a strongly squeezed state the entanglementnecessarily spreads among a large number of spins. Entangle-ment criteria based on spin-squeezing parameters benefit fromthe fact that these parameters usually depend on simple andglobal observables only, in particular, typically on low-ordermoments of collective spin operators.This has two major implications. Experimentally, the ap-proach provides an accessible and robust way for entangle-ment detection that is free of any assumptions on the systemand may therefore be suitable for many different platforms.On the theoretical side, a criterion may be obtained from twomain ingredients: (i) local uncertainty relations due to thefew-body correlations involved and (ii) exploiting the permu-tation invariance of the observables. As a consequence, thisreduces the complexity of the task of determining entangle- ment criteria drastically and hence, for example, a completeset of inequalities useful for the detection of non-separabilityfor the first and second moments of the magnetization maybe given explicitly [7]. Yet, these simplifications also set thelimitations to the spin-squeezing criteria. Extending themis therefore desirable, in particular, to platforms with non-permutation invariant observables [9–11], to open them up toentanglement schemes that are established in the permutationinvariant setting.In this work, we focus on criteria which do not rely on per-mutation invariant observables. To this end, it is important tonote that the methodology of entanglement detection via spin-squeezing, in particular the application of local uncertaintyrelations, is also applicable to observable quantities other thancollective spin operators [8]. Here, we apply Lagrange-dualityto a specific constrained optimization problem to introduce ageneral scheme which allows for the detection of many-bodyentanglement via global observables. To give a concrete ex-ample, we will focus on Fourier-transformed spin operators.Such observables arise in scattering experiments and they areintrinsically non-permutation invariant. A prominent exampleis the static structure factor, which is accessible, e.g., by neu-tron scattering from magnetic materials. The structure factorhas been demonstrated useful for entanglement quantificationfor example in the vicinity of phased Dicke states [12, 13, 16]where usual spin-squeezing criteria are unable to confirm en-tanglement. The method presented hereinafter can be used todetect k -party and other forms of multipartite entanglementsuch as k -wide entanglement by means of the structure factor[21]. To this end, the route of Lagrange duality turns out to bea fruitful way of approaching the problem.Strikingly, this approach provides an alternative way toderive the multiparty entanglement criteria of Ref. [6] if itis combined with a numerically efficient method for lowerbounds to ground state energies introduced by Baumgratz and a r X i v : . [ qu a n t - ph ] A ug Plenio in [17]. With this, it is possible to close gaps in theproof of the entanglement bounds [6] and, moreover, under amild (and numerically testable) assumption that has also beenexploited in Ref. [6] we can calculate the bounds for any en-tanglement depth k without the need of increasing the Hilbertspace dimension of the underlying problem. This impacts ex-isting experimental and theoretical work that builds on Ref.[6], see e.g. [14, 15, 24, 25, 30–33], and paves the way forthe detection of larger, potentially macroscopic, numbers ofentangled particles.For the more general non-permutation symmetric observ-ables, in order to calculate the criteria explicitly, we consideran algorithm for global, non-convex eigenvalue optimizationwhich could be combined with matrix-product state meth-ods. The general case is of practical interest as, e.g., theymay be accessible in scattering experiments with neutrons oncrystalline magnetic compounds [34] or with X-Ray light oncold atoms [22, 23]. Finally, we construct states, that can beproven to be genuine multipartite entangled by our scheme,by demonstrating how they can experimentally be generatedwith trapped ions using the high degree of control over the in-teraction provided by these systems. These findings supportthe versatility and practical importance of our framework forthe field of controlled quantum systems. II. PRELIMINARIES
In order to obtain a detection scheme for k -party entangle-ment, Sørensen and Mølmer determined the minimal varianceof the collective spin operator of a many-body state as a func-tion of its magnetization in one of the orthogonal directions.To generalize these results we start by introducing the vari-ance of a (not necessarily Hermitian) operator ˆ O in a state ˆ (cid:37) as ∆ (cid:37) [ ˆ O ] := (cid:104) ˆ O † ˆ O (cid:105) ˆ (cid:37) − (cid:104) ˆ O † (cid:105) ˆ (cid:37) (cid:104) ˆ O (cid:105) ˆ (cid:37) . (1)Eq. (1) reduces to the usual definition of the variance if ˆ O isHermitian.The goal is now to find a function F C such that for states ˆ (cid:37) belonging to a certain class C of states (e.g., k -produciblestates) one has ∆ (cid:37) [ ˆ O ] ≥ F C ( (cid:104) ˆ M (cid:105) ˆ (cid:37) ) , (2)i.e. a lower bound to the variance in terms of an additional ob-servable ˆ M playing the role of the magnetization in [6]. Letus assume that we have access to ∆ (cid:37) [ ˆ O ] and (cid:104) ˆ M (cid:105) ˆ (cid:37) in an ex-periment. If the measurements happen to violate the aboveinequality then it is guaranteed that the state in the labora-tory is not in that class. E.g., if C corresponds to the set of k -producible states then such a violation shows that the stateis ( k + 1) -party entangled. We set out to determine F C fordifferent classes C . We start with considering N spin- S parti-cles and later, for concrete examples, focus on spin chains and S = 1 / . We will define classes of states as follows. Any state ˆ (cid:37) on N spins may be written as ˆ (cid:37) = (cid:88) n p n ˆ (cid:37) ( n )1 ⊗ · · · ⊗ ˆ (cid:37) ( n ) P n , p n ≥ , (cid:88) n p n = 1 , (3)where each ˆ (cid:37) ( n ) p corresponds to the state on a subset Z ( n ) p ofthe N spins and P n denotes the number of factors in the n ’thsummand. One may now define classes of states by restrictingthe subsets Z ( n ) p : E.g., if one demands that all |Z ( n ) p | = 1 then this defines the fully separable states. If one demands P n = k then all such states are k -separable. If one restrictsthe Z ( n ) p to contain at most k spins then this defines the set of k -producible states. III. MAIN OBSERVATION
Consider now a certain class of states on N spins, i.e.,all density matrices ˆ (cid:37) as in Eq. (3) with Z ( n ) p ∈ C , where C defines the class under consideration. Furthermore, let ˆ O = (cid:80) Ni =1 ˆ O i and ˆ M = (cid:80) Ni =1 ˆ M i with ˆ O i and ˆ M i actingonly on the i ’th spin but potentially different operators at each i , i.e., we do not demand ˆ O nor ˆ M to be permutation invari-ant. Our main observation is that for any such operators oneobtains (via Lagrange duality and the variational characteriza-tion of the variance [18] generalized to non-Hermitian opera-tor, see Appendix A for details) a lower bound as in Eq. (2)with F C ( m ) = sup λ ∈ R (cid:18) λm + N min Z∈C G Z ( λ ) |Z| (cid:19) , (4)where G Z ( λ ) = inf s ∈ C λ min (cid:104) ( ˆ O Z − s ) † ( ˆ O Z − s ) − λ ˆ M Z (cid:105) . (5)Here, λ min [ · ] denotes the smallest eigenvalue of the Hermitianmatrix in brackets and ˆ O Z = (cid:80) i ∈Z ˆ O i and similarly for ˆ M .We note that G Z is concave and F C is convex. Furthermore,the above holds for any collection of spins, such that, e.g. D -dimensional lattices are included. If the variance of ˆ O andthe mean value of ˆ M are experimentally accessible and vio-late the inequality in Eq. (2) then, without making any furtherassumptions, one can conclude that the state in the labora-tory is not in the class C . How strong the bound can be vio-lated depends on the observables. So far, we have introduceda framework for the detection of various multiparticle entan-glement structures by global measurements without makingany assumption on the underlying system. As the violationof these criteria necessitates a sufficiently small variance Eq.(1) it may be seen as a generalization of the spin-squeezingphenomenon to, both, arbitrary observables and more generalforms of entanglement.Notably, we only require that the operators are the sumof single site operators, so that, e.g., ˆ O = (cid:80) Ni =1 f i ˆ σ iz with f i ∈ C fits into our framework. Additionally, a bound to asum of variances in terms of the expectation values of multi-ple observables can directly be incorporated into equation Eq.(4), if these operators have the above local form [35].Now, how hard is it to actually compute F C ? First of all,to determine G Z , (a) one needs to be able to find the smallesteigenvalue of a potentially very large matrix: a priori, the di-mensions of the involved matrices are exponentially large in |Z| . If, e.g., the goal is to detect k -party entanglement then thedimensions of the involved matrices are exponentially large in k . Secondly, in order to obtain F C , (b) one needs to deter-mine G Z for all Z ∈ C . Considering the example of k -party entanglement again, one needs to compute G Z for all subsets Z containing at most k − spins as we allow fornon-permutation-invariant operators. If the involved oper-ators are permutation invariant this complexity is dramati-cally reduced: It is then sufficient to determine G Z for Z = { } , { , } , . . . , { , . . . , k − } .Finally, (c) the function to be minimized over s in Eq. (5)may exhibit local minima and thus calls for global non-convexoptimization.We remark that (a) may be addressed using efficientmethods such as the density-matrix renormalization group(DMRG) exploiting the local form of the observables. Impor-tantly, to guarantee a lower bound one has to carefully monitorconvergence. On the other hand, it is also possible to utilize ascheme based on a semi definite program (SDP) that providesa lower bound to λ min ( s ) [17]. Notably, for the permutation-ally invariant case where ˆ O and ˆ M are given by collectivespin- / operators in two orthogonal directions, we find a re-laxation of the SDP that calculates a lower bound to the small-est eigenvalue where the size of the configuration |Z| entersthe optimization soley as a parameter and present an applica-tion below. The technical details of this method are shifted toAppendix E.The third step, (c), can be tackled utilizing an algorithm in-troduced in Ref. [19] that is based on quadratic support func-tions, which are determined by (i) the value of the eigenvaluefunction, (ii) its derivative for specific values of s , and (iii) anestimate of the curvature that is given as an input to the algo-rithm. To reliably obtain a global optimum, the estimate of thecurvature is required to be a lower bound on the second deriva-tive of the eigenvalue function in the entire parameter range.Here, we use the algorithm heuristically, decreasing the esti-mated curvature until we do not observe any change in the re-sult, see Appendix C for details. Note, that in order to provideresults for large numbers of spins in a chain, one may combinethis algorithm with matrix-product states (MPSs) and opera-tors (MPOs) by reformulating steps (i) and (ii) in terms ofMPSs and MPOs. As mentioned above, step (i) is a simpleground state search as can be carried out using DMRG. Oncethe eigenvalue function has been evaluated using DMRG, thederivative may be obtained from the calculated optimal state,see Appendix C. As a general algorithm for eigenvalue opti-mization of matrix-valued functions, it might be a useful toolalso for other application in quantum science.On the other hand, when we use in step (a) the scheme ofRef. [17] to obtain a lower bound to the lowest eigenvalue, theminimization over s is carried out by computing the functionin the entire parameter range, supported by the standard opti-mization toolbox of MATLAB. For the permutation invariant collective spin observables we emphasize again that the sizeof the configuration |Z| enters the optimization as a parameteronly. IV. WITNESSING MULTIPARTY ENTANGLEMENT
A state is k -producible if it can be decomposed as in Eq.(3) with each ˆ (cid:37) ( n ) p corresponding to a state of at most k spins.Denoting by [ N ] = { , . . . , N } the set of all spins, we hencehave that k -producible states fulfil Eq. (2) with F C as in Eq. (4)and C = (cid:8) Z ⊂ [ N ] (cid:12)(cid:12) |Z| ≤ k (cid:9) . (6)Operationally, a pure k -producible state can be prepared froma fully separable state via an interaction that acts on non-overlapping sets of spins separately and the number of spinsin each set is upper bounded by k . Mixed k -producible statesare just mixtures of the states of the above type. States whichare not k -producible are ( k + 1) -party entangled. Experimen-tally, in different setups the presence of k -party entanglementhas been verified, see e.g. [14, 22–33]. k -party entanglementprovides a multiparticle entanglement hierarchy that can beverified for a large system even if only a subset of spins canbe accessed: the number of parties that may be confirmed tobe entangled on the subset also gives a lower bound to the en-tanglement of full system [20]. In contrast, k -partite entangle-ment cannot solely be specified by conditions on each subset Z ( n ) p separately, but needs knowledge of the full system thatone would like to characterize. We remark that the versatilityof Eq. (4) opens the possibility to find criteria for other en-tanglement structures. For example, note that k -producibilityis insensitive to the spatial distribution of entanglement: Sup-pose the N spins are arranged on a chain and are genuinely -party entangled (i.e., they are not fully separable). One mightthus want to be able to distinguish between whether the firsttwo spins are entangled or the first and last spin on the chainare entangled as it might be much less challenging to preparethe former case than the latter. This fine-grained form of mul-tiparticle entanglement is captured by the notion of k -wide en-tanglement introduced in [21] and easily incorporated in ourframework by defining C as the set of configurations wherethe spins are at most a distance k apart. V. NON-PERMUTATIONALLY SYMMETRIC EXTREMESPIN SQUEEZING
To give concrete examples, we now focus on spins arrangedon a chain and consider the operators ˆ O = 1 √ N ˆ S z ( q ) = 1 √ N N (cid:88) j =1 e i qj ˆ S ( j ) z (7)and ˆ M = ˆ S x (0) /N =: ˆ S x /N . Here, ˆ S ( j ) z denotes a spin- S operator along the z direction acting on the j ’th spin and F C ~ h S x i / N h S x i / N FIG. 1. Detecting k -party entanglement: Lower bounds ˜ F C to thefunctions F C in Eq. (8) for the three coinciding cases q = 0 , π/ , and π , with S = 1 / , and k = 1 , , , , , for any N ob-tained from the SDP described in Appendix E. To compare, we alsoshow the bounds from Ref. [6] for k even (light blue, dashed). Underthe assumption (2) as described in the main text the bounds can becalculated for large k with numerical effort independent of k . Fordemonstration, we show the bound for k = 2 · (inset). If a mea-surement lies below a curve corresponding to k in the plot then thestate is ( k + 1) -party entangled. q ∈ [0 , π ] and we choose ˆ M to be the magnetization per spinin the direction x . As a consequence of our main observationit follows that for every state ˆ (cid:37) in the class C we have ∆ (cid:37) [ ˆ S z ( q )] N ≥ F C (cid:32) (cid:104) ˆ S x (cid:105) ˆ (cid:37) N (cid:33) , (8)with F C as in Eq. (4).Inequality (8) includes the k -party bounds of Sørensen andMølmer as a special case with q = 0 and the correspondingset C . The inequality of Sørensen and Mølmer is maximallyviolated by the ground state of the so-called spin-squeezingHamiltonian, i.e. the one axis twisting plus external field, thathas been shown to be an enhancement over the simple one-axis twisting. Experimentally, the bound has been able to con-firm multiparticle entanglement in various setups [24, 25, 30–33].For this important special case, the SDP method presentedin Appendix E for spin-1/2 gives a reliable lower bound to F C that scales linearly with the size of the set C and, hence, lin-early with k . Clearly, these bounds can then be used to obtaincriteria for arbitrary spin- S particles as well. Notably, besidesthe scalability, our approach avoids all assumptions the proofof the bounds in [6] has been relying on and, hence, makestechnically subsequent publications [14, 15, 24, 25, 30–33]rigorous that refer to the original work.More specifically, (1) convexity of the bounds in [6] hasbeen one of the requirements of the proof and which is verifiednumerically can either be investigated numerically or follows as a result of assumption (3) discussed below. From numer-ical inspections for small k one may also infer that (2) theoptimal configuration in Eq. (4) is Z = { , . . . , k − } and,(3) for |Z| even, the infimum in Eq. (5) (with the above men-tioned observables) is achieved for (cid:104) ˆ S z (cid:105) ˆ (cid:37) = 0 . The latter as-sumption transforms the variance-minimization into a simpleground state search that can be solved efficiently, in particu-lar, if (4) the ground state is assumed to lie in the symmetricsubspace of dimension k + 1 . In contrast, the case |Z| odd re-mains numerically more costly [6, 33] and is usually omitted.(1) can be also considered the direct consequence of (3) andthe fact that the set of points corresponding to physical statesin the ( (cid:104) ˆ S x (cid:105) , (cid:104) ˆ S z (cid:105) ) -space is convex.Importantly, our method does not need these assumptionsand there is no technical difference between an even and oddnumber of spins. We can thus, indeed, consider the optimiza-tion over all configurations in Eq. (4) and, moreover, effi-ciently determine bounds for k even and odd. We observenumerically very good agreement with the bounds obtainedby previous methods and that they improve with increasing k ,see Fig. 1.As noted before, not only for this special case but for any q , the observables that appear in Eq. (8) are of practical im-portance since, e.g., the generalized variance may be acces-sible in scattering experiments [22, 23]. Then, depending onthe value of q , symmetries may simplify the different stepsthat are required to numerical compute F C significantly. Onone hand, the operators may exhibit an efficient parametriza-tion, see Appendix B for a discussion of symmetries of theoperators under consideration. On the other hand, in gen-eral, in order to obtain a criterion for k -party entanglement,the minimization over all subsets Z ⊂ [ N ] of cardinality atmost k needs to be considered. Here, the presence of symme-tries, that may also be present for q (cid:54) = 0 , reduces the numberof ways to select k out of N spins that may lead to distinctbounds. For example, taking translational symmetries into ac-count, the number of inequivalent subsets may be counted (seeAppendix B) and determined numerically [36]. To provide anexample, for q = π/ , S = 1 / , and k up to we com-pute the lower bound ˜ F C to F C for k -producible states (i.e., C as in Eq. (6)) using the SDP approach [17] and AppendixE, see Fig. 1. Note that the given bounds for k -local statesare valid for any number of spins N , again a consequence ofsymmetries, see Appendix B. VI. ENGINEERING k -PARTY ENTANGLED STATES BY AQUANTUM QUENCH IN ION TRAPS Modern experimental platforms such as trapped ions allowfor the implementation of quantum systems with a wide rangeof tunable interactions. This has raised the interest in controland study of spin systems with artificial interactions that, forexample, can result in exotic quantum phases and novel quan-tum states [37–40] where the quantification of entanglementcan help to characterize those quantum effects [41–44]. Wefind numerically and demonstrate below that states violatinginequality Eq. (8) can be generated under a time-evolution. In t [ J -1 ] ¢ [ S z ( q )]/ N FIG. 2. Detection of k -party entanglement after a quench under theHamiltonian in Eq. (9) with couplings as in Eq. (10). Shown as afunction of time are ∆ [ ˆ S z (2 π/ /N (red), and the lower bounds F C evaluated at the instantaneous magnetization (cid:104) ˆ S x (cid:105) /N (black) for N = 15 spins. At t ≈ . J − we observe the maximum of 15-party entanglement. particular, one may consider a quantum quench under an IsingHamiltonian with a transverse field ˆ H = (cid:88) i
In conclusion, we have introduced a method to derive crite-ria for the detection of various many-body entanglement struc-tures, with emphasis on k -party entanglement. Other entan-glement structures such as k -partite or k -wide [21] entangle-ment are immediately covered by our scheme if an additionaloptimization over specific spin configurations is taken into ac-count. The criteria make no assumptions on the state, requireto measure a few global observables only and are applicableto any number of total spins. In contrast to previous works,the observables do not have to be permutation invariant. In-stead, even if there is no symmetry identified, we show howto compute the multiparty entanglement criteria, where thealgorithmic approach we take can be further extended usingDMRG which may allow for the investigation of systems ofmany spins. We leave this exploration to future work. Wefind that for the case of permutation invariant observables, ourapproach enables to derive the bounds without any technicalassumption. As an application of the method we provide anexperimental protocol to test the entanglement criteria underrealistic conditions with trapped ions. VIII. ACKNOWLEDGEMENTS
We acknowledge the support of the ERC synergy grantBioQ, the EU project QUCHIP, the EU (ERC Starting Grant258647/GEDENTQOPT, CHIST-ERA QUASAR, COST Ac-tion CA15220), the Spanish Ministry of Economy, Indus-try and Competitiveness and the European Regional Devel-opment Fund FEDER through Grant No. FIS2015-67161-P (MINECO/FEDER), the Basque Government (Project No.IT986-16), the UPV/EHU program UFI 11/55 and the Aus-trian Science Fund (FWF) through the START project Y879-N27. This work was supported by the DFG through SFB 1227(DQ-mat) and the RTG 1991, the ERC grants QFTCMPS andSIQS, and the cluster of excellence EXC201 Quantum Engi-neering and Space-Time Research. The authors acknowledgesupport by the state of Baden-W¨urttemberg through bwHPC. [1] D. Wineland, et al., Spin squeezing and reduced quantum noisein spectroscopy. Phys. Rev. A 46, R6797 (1992).[2] M. Kitagawa, and M. Ueda, Squeezed Spin States, Phys. Rev.A 47, 5138 (1993).[3] J. Ma, X. Wang, C.P. Sun, and F. Nori, Quantum Spin Squeez-ing, Phys. Rep. 509, 89 (2011).[4] J. Korbicz, J. Cirac, M. Lewenstein, Spin Squeezing Inequali-ties and Entanglement of N Qubit States, Phys. Rev. Lett. 95,120502 (2005).[5] G. T´oth, C. Knapp, O. G¨uhne, H.J. Briegel, Spin Squeezing andEntanglement, Phys. Rev. A 79, 042334 (2009).[6] A.S. Sørensen, and K. Mølmer, Entanglement and ExtremeSpin Squeezing, Phys. Rev. Lett. 86, 4431 (2001).[7] G. Vitagliano, I. Apellaniz, I.L. Egusquiza, and G. T´oth, SpinSqueezing and Entanglement for an Arbitrary Spin, Phys. Rev.A 89, 032307 (2014).[8] G. Vitagliano, et al., Spin Squeezing Inequalities for ArbitrarySpin, Phys. Rev. Lett. 107, 240502 (2011).[9] R. Blatt, and C.F. Roos, Quantum simulations with trappedions, Nat. Phys. 8, 277 (2012).[10] I. Bloch, J. Dalibard, and S. Nascimb´ene, Quantum simulationswith ultracold quantum gases, Nat. Phys. 8, 267 (2012).[11] M.H. Devoret, and R.J. Schoelkopf, Superconducting Circuitsfor Quantum Information: An Outlook, Science 339, 1169(2013).[12] M. Cramer, H. Wunderlich, and M.B. Plenio, Measuring En-tanglement in Condensed Matter Systems, Phys. Rev. Lett. 106,020401 (2011).[13] O. Marty, et al., Quantifying entanglement with scattering ex-periments, Phys. Rev. B 89, 125117 (2014).[14] B. L¨ucke, J. Peise, G. Vitagliano, J. Arlt, L. Santos, G. T´oth,and C. Klempt, Detecting Multiparticle Entanglement of DickeStates, Phys. Rev. Lett. 112, 155304 (2014).[15] G. Vitagliano, et al., Entanglement and extreme spin squeezingof unpolarized states, New J. Phys. 19, 013027 (2017).[16] P. Krammer, H. Kampermann, D. Bruß, R.A. Bertlmann, L.C.Kwek, and C. Macchiavello, Multipartite Entanglement Detec-tion via Structure Factors, Phys. Rev. Lett. 103, 100502 (2009).[17] T. Baumgratz, and M.B. Plenio, Lower bounds for ground statesof condensed matter systems, New J. Phys. 14, (2012).[18] L. Dammeier, et al., Uncertainty relations for angular momen-tum, New J. Phys. 17, 093046 (2015); see also Apellaniz et. al,Optimal witnessing of the quantum Fisher information with fewmeasurements, Phys. Rev. A 95, 032330 (2017).[19] E. Mengi, et al., Numerical Optimization of Eigenvalues of Her-mitian Matrix Functions, SIAM 35, 2 (2014).[20] O. G¨uhne, G. T´oth, and H.J. Briegel, Multipartite entanglementin spin chains, New J. Phys. 7, 229 (2005).[21] S. W¨olk, and O. G¨uhne, Characterizing the Width of Entangle-ment, New J. Phys. 18, 123024 (2016).[22] R.A. Hart et al., Observation of antiferromagnetic correlationsin the Hubbard model with ultracold atoms, Nature 519, 211(2015).[23] A. Mazurenko, et al., Experimental realization of a long-rangeantiferromagnet in the Hubbard model with ultracold atoms,Nature 545, 462 (2017).[24] M.A. Riedel, et al., Atom-chip-based generation of entangle-ment for quantum metrology, Nature 464, 1170 (2010).[25] C. Gross, et al., Nonlinear atom interferometer surpasses clas-sical precision limit, Nature 464, 1165 (2010). [26] T. Monz, et al., 14-Qubit Entanglement: Creation and Coher-ence, Phys. Rev. Lett. 106, 130506 (2011).[27] J.T. Barreiro, et al., Demonstration of genuine multipartite en-tanglement with device-independent witnesses, Nat. Phys. 9,559 (2013).[28] F. Haas, et al., Entangled states of more than 40 atoms in anoptical fiber cavity, Science 344, 180 (2014).[29] R. McConnell, et al., Entanglement with negative Wigner func-tion of almost 3,000 atoms heralded by one photon, Nature 519,439 (2015).[30] O. Hosten, Measurement noise 100 times lower than thequantum-projection limit using entangled atoms, Nature 529,505 (2016).[31] K.C. Cox, Deterministic Squeezed States with Collective Mea-surements and Feedback, Phys. Rev. Lett. 116, 093602 (2016).[32] N.J. Engels, et al., Bell Correlations in Spin-Squeezed States of500000 Atoms, Phys. Rev. Lett. 118, 140401 (2017).[33] L. Dellantonio, et al., Multi-partite entanglement detection withnon symmetric probing, Phys. Rev. A 95, 040301(R) (2017).[34] J. Jensen, and A.R. Mackintosh, Rare Earth Magnetism: Struc-tures and Exciations, Clarendon Press, Oxford (1991).[35] G. Vitagliano, et al., Entanglement and extreme planar spinsqueezing, arXiv:1705.09090 (2017).[36] J. Sawada, Generating bracelets in constant amortized time,SIAM J. Comput., 31(1), 259–268 (2001); S. Karima, J.Sawada, Z. Alamgir, S.M. Husnine, Generating bracelets withfixed content, Theoret. Comp. Sc. 475 103, (2013).[37] P. Hauke, L. Bonnes, M. Heyl and W. Lechner, Probing En-tanglement in Adiabatic Quantum Optimization with TrappedIons, Front. Phys. 3:21, 1 (2015).[38] P. Hauke, D. Marcos, M. Dalmonte, and P. Zoller, QuantumSimulation of a Lattice Schwinger Model in a Chain of TrappedIons, Phys. Rev. X 3, 041018 (2013).[39] G.-D. Lin, C. Monroe, and L.-M. Duan, Sharp Phase Transi-tions in a Small Frustrated Network of Trapped Ion Spins, Phys.Rev. Lett. 106, 230402 (2011).[40] A. Chiuri, G. Vallone, N. Bruno, C. Macchiavello, D. Bruß, andP. Mataloni, Hyperentangled Mixed Phased Dicke States: Opti-cal Design and Detection, Phys. Rev. Lett. 105, 250501 (2010).[41] O. Marty, M. Cramer, and M.B. Plenio, Practical EntanglementEstimation for Spin-System Quantum Simulators Phys. Rev.Lett. 116, 105301 (2016).[42] M. Cramer, et al., Spatial entanglement of bosons in opticallattices, Nat. Comm. 4, 2161 (2013).[43] D.A. Abanin and E. Demler, Measuring Entanglement Entropyof a Generic Many-Body System with a Quantum Switch Phys.Rev. Lett. 109, 020504 (2012); A.J. Daley, et al., MeasuringEntanglement Growth in Quench Dynamics of Bosons in anOptical Lattice, Phys. Rev. Lett. 109, 020505 (2012).[44] R. Islam, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).[45] L. Novo, T. Moroder, and O. G¨uhne, Genuine multiparticle en-tanglement of permutationally invariant states, Phys. Rev. A 88,012305 (2013).[46] J.H. Redfield, The Theory of Group-Reduced Distributions,Amer. J. Math. 49, 3 (1927).[47] G. P´olya, Kombinatorische Anzahlbestimmungen f¨ur Gruppen,Graphen und chemische Verbindungen Acta Math. 68, 145(1937).[48] K. Kim, et al., Entanglement and Tunable Spin-Spin Couplingsbetween Trapped Ions Using Multiple Transverse Modes, Phys.
Rev. Lett. 103, 120502 (2009).[49] M. Grant, and S. Boyd, CVX: Matlab software for disciplinedconvex programming, version 2.0 beta. http://cvxr.com/cvx,September 2013. [50] J. L¨ofberg, YALMIP : A toolbox for modeling and optimizationin MATLAB, in Proc. IEEE International Symposium on Com-puter Aided Control Systems Design, Taipei, Taiwan, 2004.
Appendix A: Main observation
For any state (cid:37) and any operator O we have ∆ (cid:37) [ O ] = (cid:104) O † O (cid:105) (cid:37) − (cid:104) O † (cid:105) (cid:37) (cid:104) O (cid:105) (cid:37) = (cid:104) O † O (cid:105) (cid:37) − (cid:104) O (cid:105) ∗ (cid:37) (cid:104) O (cid:105) (cid:37) + inf s ∈ C | s − (cid:104) O (cid:105) (cid:37) | = inf s ∈ C (cid:2) (cid:104) O † O (cid:105) (cid:37) + | s | − s ∗ (cid:104) O (cid:105) (cid:37) − s (cid:104) O (cid:105) ∗ (cid:37) (cid:3) = inf s ∈ C (cid:104) ( O − s ) † ( O − s ) (cid:105) (cid:37) (A1)and (as the second line shows) the minimum is attained at s = (cid:104) O (cid:105) (cid:37) . Using this variational form of the generalized variance, wefind that for any state (cid:37) , any operator O , any λ ∈ R , and any Hermitian operator M , ∆ (cid:37) [ O ] = λ (cid:104) M (cid:105) (cid:37) + ∆ (cid:37) [ O ] − λ (cid:104) M (cid:105) (cid:37) = λ (cid:104) M (cid:105) (cid:37) + inf s ∈ C (cid:104) (cid:2) ( O − s ) † ( O − s ) − λM (cid:3) (cid:105) (cid:37) ≥ λ (cid:104) M (cid:105) (cid:37) + inf s ∈ C λ min (cid:2) ( O − s ) † ( O − s ) − λM (cid:3) =: λ (cid:104) M (cid:105) (cid:37) + G [ N ] ( λ ) , (A2)where λ min [ · ] denotes the smallest eigenvalue of the Hermitian matrix in brackets. Note that λ min (cid:2) ( O − s ) † ( O − s ) − λM (cid:3) is concave in λ such that G [ N ] ( λ ) is also concave in λ .Suppose now that (cid:37) is of the product form (cid:37) ⊗ · · · ⊗ (cid:37) P . (A3)This divides the N spins into sets of spins Z p ⊂ { , . . . , N } , p = 1 , . . . , P , which form a partition (cid:83) p Z p = { , . . . , N } and Z p denotes the set of spins that (cid:37) p acts on. If we further assume that O = (cid:80) Ni =1 O i and M = (cid:80) Ni =1 M i then (we use the shorthandnotation O Z = (cid:80) i ∈Z O i and similarly for M ), we find for states as in Eq. (A3) ∆ (cid:37) [ O ] = P (cid:88) p =1 ∆ (cid:37) p [ O Z p ] ≥ P (cid:88) p =1 (cid:18) λ (cid:104) M Z p (cid:105) (cid:37) p + |Z p | G Z p ( λ ) |Z p | (cid:19) = λ (cid:104) M (cid:105) (cid:37) + P (cid:88) p =1 |Z p | G Z p ( λ ) |Z p | . (A4)Now suppose that all the Z p are elements of some set C . Then for all states (cid:37) as in Eq. (A3) ∆ (cid:37) [ O ] ≥ λ (cid:104) M (cid:105) (cid:37) + P (cid:88) p =1 |Z p | min Z∈C G Z ( λ ) |Z| = λ (cid:104) M (cid:105) (cid:37) + N min Z∈C G Z ( λ ) |Z| . (A5)Finally, we may extend the above to states that are convex combinations of states as in Eq. (A3), i.e., for states of the form (cid:37) = (cid:88) n p n (cid:37) ( n )1 ⊗ · · · ⊗ (cid:37) ( n ) P n =: (cid:88) n p n (cid:37) n , p n ≥ , (cid:88) n p n = 1 , (A6)by the Cauchy–Schwarz inequality: ∆ (cid:37) [ O ] = (cid:88) n p n (cid:104) O † O (cid:105) (cid:37) n − (cid:88) n,m p n p m (cid:104) O † (cid:105) (cid:37) n (cid:104) O (cid:105) (cid:37) m ≥ (cid:88) n p n (cid:104) O † O (cid:105) (cid:37) n − (cid:88) n p n |(cid:104) O (cid:105) (cid:37) n | ≥ (cid:88) n p n ∆ (cid:37) n [ O ] ≥ λ (cid:104) M (cid:105) (cid:37) + N min Z∈C G Z ( λ ) |Z| . (A7)As this holds for all λ ∈ R , we may take the supremum to arrive at ∆ (cid:37) [ O ] ≥ F C ( (cid:104) M (cid:105) (cid:37) ) , F C ( m ) = sup λ ∈ R (cid:18) λm + N min Z∈C G Z ( λ ) |Z| (cid:19) , (A8)where we recall that G Z is concave and note that F C convex. Appendix B: Symmetries for q = 2 π/z In this section, we show how to exploit the symmetries of the operator Eq. (7) for q = 2 π/z and z integer in order to reducethe numerical effort to derive the bound F C in Eq. (8). By the periodicity of the phases, for some set Z ⊂ [ N ] , we may write ˆ O Z = 1 √ N (cid:88) j ∈Z e i qj ˆ S ( j ) z =: 1 √ N ζ (cid:88) n =1 e i qn (cid:88) j ∈Z n ˆ S ( j ) z , (B1)where, since F C in Eq. (8) (see also the definitions Eq. (4) and (5)) is invariant under local spin flips in z -direction, we mayconsider ζ = z for z odd and ζ = z/ for z even, and define Z n := (cid:8) j ∈ Z (cid:12)(cid:12) (e i qj = e i qn ) ∨ (e i qj = − e i qn ) (cid:9) ⊂ Z .Therefore, the bounds for z = 1 coincides with the bounds for z = 2 , since in both cases ζ = 1 . For z = 4 , i.e., ζ = 2 , ˆ O Z = 1 √ N i (cid:88) j ∈Z ˆ S ( j ) z + (cid:88) j ∈Z ˆ S ( j ) z . (B2)With this, we find that G Z ( λ ) = inf s ∈ C λ min (cid:104) ( ˆ O Z ( j ) − s ) † ( ˆ O Z ( j ) − s ) − λ ˆ M Z ( j ) (cid:105) , = inf s ∈ C λ min √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:61) ( s ) † √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:61) ( s ) − λN (cid:88) j ∈Z ˆ S ( j ) x + √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:60) ( s ) † √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:60) ( s ) − λN (cid:88) j ∈Z ˆ S ( j ) x = inf (cid:61) ( s ) ∈ R λ min √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:61) ( s ) † √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:61) ( s ) − λN (cid:88) j ∈Z ˆ S ( j ) x + inf (cid:60) ( s ) ∈ R λ min √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:60) ( s ) † √ N (cid:88) j ∈Z ˆ S ( j ) z − (cid:60) ( s ) − λN (cid:88) j ∈Z ˆ S ( j ) x . (B3)That is, the case z = 4 can also be expressed in terms of the permutation invariant case.Now, note that for all n the operators ˆ O Z and ˆ M Z = (cid:80) j ∈Z ˆ S ( j ) x are invariant under permutations of the sites in Z n . Anyoperator with these symmetry properties may be represented by an operator ˆ X with a decomposition as (see e.g. Ref. [45]) ˆ X = (cid:77) j =( j ,...,j ζ ) (cid:88) m ˆ X [1] m ( j ) ⊗ · · · ⊗ ˆ X [ ζ ] m ( j ζ ) =: (cid:77) j ˆ X ( j ) , (B4)where j n = j min n , j min n + 1 , . . . , |Z n | / with j min n = 0 ( j min n = 1 / if |Z n | is even (odd) and where ˆ X [ n ] m ( j n ) are of dimension (2 j n + 1) × (2 j n + 1) . Now, every eigenvector of ˆ X belongs to one of the blocks in the decomposition Eq. (B4). To calculate G Z ( λ ) we may therefore minimize each block separately, and obtain G Z ( λ ) = min j inf s ∈ C λ min (cid:104) ( ˆ O Z ( j ) − s ) † ( ˆ O Z ( j ) − s ) − λ ˆ M Z ( j ) (cid:105) , (B5)where ˆ O Z ( j ) = (cid:80) ζn =1 e i qn ˆ S z ( j n ) and M Z ( j ) = (cid:80) ζn =1 ˆ S x ( j n ) , with ˆ S α ( j ) the α -component spin- j operator.Notably, for every configuration of spins Z , G Z depends the cardinality of the subsets |Z n | only. Moreover, the varianceof ˆ O Z in Eq. (B4) only depends on e i q ( m − n ) , for m, n = 1 , . . . , ζ . Therefore we may consider soley configurations withcardinalities |Z n | inequivalent under cyclic shifts and reflection of the index n , i.e. inequivalent under n (cid:55)→ n + l for l ∈ N , and n (cid:55)→ ζ − n + 1 . Only configurations that cannot be obtained from one another by these operations will result in different bounds.Now, in order to derive the bound for k -party entanglement we need to consider all inequivalent configurations Z withcardinality at most k . This then gives the size of the set C to obtain bounds for multiparty entanglement for the observables ofconsideration. In general, this may be done numerically. For the special case where N is divisible by ζ we count the number ofthese configurations for for fixed k using P´olya’s enumeration theorem (PET) [46, 47]. We set out to count the number of waysone can assign a number ≤ k ( n ) ≤ k max , where k max := N/ζ , to every phase labelled by n = 1 , . . . , ζ with the constraint (cid:80) ζn =1 k ( n ) = N . For a particular configuration or ‘coloring’ k : { , . . . , ζ } → { , . . . , k max } , k ( n ) determines the numberof spins with phase factor n . As already noted, since only quantities of the form e i q ( m − n ) are relevant, combinations that areequal up to a cyclic shifts and reflection will result in the same bound. Mathematically, two colorings k and k (cid:48) are equivalent ifthere is a permutation τ of the set { , . . . , ζ } that belongs to the dihedral group D ζ and is such that k (cid:48) = k ◦ τ − . We introducethe generating function, a polynomial in k max variables F D ζ ( r , . . . , r k max ) = (cid:80) p ,...,p k max f D ζ ( p , . . . , p max ) r p · · · r p k max k max where f D ζ is the number of orbits, i.e. distinct configuartions, under D ζ with fixed content. Hence, in order to count all orbitsfor a constant number of sites, we are interested in the coefficient sum n k,N,ζ = (cid:88) (cid:80) np n = N f D ζ ( p , . . . , p k max ) . (B6)By the PET F D ζ ( r , . . . , r p max ) = Z D ζ (cid:32) k max (cid:88) n =0 r n , . . . , k max (cid:88) n =0 r ζn (cid:33) , (B7)where Z D ζ is the so-called cycle index polynomial of the dihedral group, here for two colors, given by Z D ζ ( t , t , . . . ) = ζ (cid:16)(cid:80) d | ζ ϕ ( d ) t ζ/dd + ζt t ( ζ − / (cid:17) , ζ odd, ζ (cid:16)(cid:80) d | z ϕ ( d ) t ζ/dd + ζ t t ( ζ − / + ζ t ζ/ (cid:17) , ζ even, (B8)where the sum runs over all divisors d of ζ and ϕ denotes Euler’s totient function. Hence, we may use the cycle index to calculateEq. (B6). We find n k,N,ζ = 12 ζ (cid:88) d | gcd ( k,ζ ) ϕ ( d ) (cid:88) p ∈ S ( d ) (cid:18) ζ/d p (cid:19) + (cid:40) (cid:80) k max n =0 (cid:80) p ∈ S ( n ) (cid:0) ( ζ − / p (cid:1) , ζ odd, (cid:16)(cid:80) p ∈ S (2) (cid:0) ζ/ p (cid:1) + (cid:80) ( p , ˜ p ) ∈ S (cid:0) p (cid:1)(cid:0) ( ζ − / p (cid:1)(cid:17) , z even, (B9)0 N k
80 20 200 0 060100 60200 n k , N , z N k
400 20 200 0
201 202 153 15104 105 50 0 ³ =3 ³ =4 ³ = N n k , N , z l og ( n ) k , N , z N k
FIG. 3. The function n k,N,ζ Eq. (B9) which gives the number of ways choosing k spins out of N that are inequivalent under the dihedralgroup D z . This gives an upper bound on the number of configuration one has to optimize over in order to find a bound for k -party bound fromobservables of the form Eq. (7) for ζ different phases. where S ( d ) = (cid:40) ( p , . . . , p k max ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k max (cid:88) n =0 p n = ζ/d, k max (cid:88) n =1 np n = k/d (cid:41) ,S ( m ) = (cid:40) ( p , . . . , p k max ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k max (cid:88) n =0 p n = ( ζ − / , k max (cid:88) n =1 np n = k − m (cid:41) ,S = (cid:40) ( p , . . . , p k max , ˜ p , . . . , ˜ p k max ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k max (cid:88) n =0 p n = 2 , k max (cid:88) n =0 ˜ p n = ( ζ − / , k max (cid:88) n =1 n ( p n + 2˜ p n ) = k (cid:41) . (B10) Appendix C: About the algorithm
Here, we describe an algorithm to solve numerically the optimization problem of Eq. (5). More specifically, we need to solvea global eigenvalue minimization, i.e. to minimize the lowest eigenvalue, λ min ( x ) = λ min ( A ( x )) of a matrix-valued function A ( x ) = ( O − s ) † ( O − s ) + λM, (C1)over a box B ⊂ R given by the conditions x ∈ [ λ min [( O † + O ) / , λ max [( O † + O ) / and x ∈ [ λ min [i( O † − O ) / , λ max [i( O † − O ) / and s = x + i x . Generally speaking, the difficulty of the optimization problem comes fromits nonconvexity. The algorithm described in [19] adresses this challenge by introducing quadratic support functions that pro-vide a lower bound to the eigenvalue function λ min ( x ) . The determination of the support functions relies on a global lower bound γ to λ min [ ∇ λ min ( x )] which requires the analyticity of the eigenvalue function. For any x ∈ B where Λ( x ) is non-degenerate,a support function is given by q ( x ) = λ min ( x ) + ∇ λ min ( x ) T ( x − x ) + γ (cid:107) x − x (cid:107) . (C2)Thus, to determine the support function we need to evaluate λ min ( x ) and the gradient λ min ( x ) , which is given by ( ∇ λ min ( x )) j = (cid:104) Ψ | ∂ A ( x ) ∂x j | Ψ (cid:105) , (C3)where | Ψ (cid:105) denotes the eigenstate of A to the lowest eigenvalue λ min ( x ) .For large systems, we can use DMRG to determine λ min ( x ) as well as | Ψ (cid:105) and calculate the gradient Eq. (C3) exploiting thefact that expectation values of matrix-product operators with matrix-product states can be determined efficiently. To study thesecond derivative of λ min ( x ) in more details, we assume that λ min ( x ) is non-degenerate for all x inside the parameter rangedefined above. The Hessian of λ min ( x ) is then given by (see section 3.2.3 of [19]) ( ∇ λ min ( x )) i,j = (cid:104) Ψ | ∂ A ∂x i ∂x j | Ψ (cid:105) − (cid:88) k> λ k ( x ) − λ ( x ) Re (cid:20) (cid:104) Ψ | ∂ A ∂x i | Ψ k (cid:105)(cid:104) Ψ k | ∂ A ∂x j | Ψ (cid:105) (cid:21) , (C4)1 -1 -14 -12 -10 -8 -6 -4 -14 -12 -10 -8 -6 -4 j F Z ° = - F Z ° = j ¸ h S x i / N j G Z ° = - G Z ° = j / j Z j FIG. 4. As an example we consider the error between the bound Eq. (8) F C to detect ( k + 1) -party entanglement with k = 15 and q = 2 π/ obtained for γ = 8 and γ = 10 . For a large magnetic field which corresponds to an expectation value of (cid:104) S x (cid:105) /N close to 0.5 the problemis convex, as may be seen from direct inspection whereas for smaller external fields the problem becomes nonconvex (not shown). In theexample, the algorithm stops either when the gap between the support function Eq. (C2) and the function value becomes smaller than − orwhen a maximum of 1500 iteration is reached. where λ k and | Ψ k (cid:105) denote the k th smallest eigenvalue and corresponding eigenvector, respectively, of A . Moreover, ∂ A ∂x i ∂x j = 2 δ i,j ,∂ A ∂x = − ( O † + O ) + 2 x ,∂ A ∂x = − i( O † − O ) + 2 x . (C5)In particular for O = S α ( q ) and M = S β one obtains ∂ A ∂x = 2 − N (cid:88) j =1 cos( qj ) S jα ,∂ A ∂x = 2 + N (cid:88) j =1 sin( qj ) S jα . (C6)A lower bound to the minimal eigenvalue of ∇ Λ in terms of the spectral gap may be given as λ min [ ∇ λ min ( x )] = 2 − λ max (cid:34)(cid:88) k> λ k ( x ) − λ ( x ) Re (cid:20) (cid:104) Ψ | ∂ A ∂x i | Ψ k (cid:105)(cid:104) Ψ k | ∂ A ∂x j | Ψ (cid:105) (cid:21)(cid:35) i,j ≥ (cid:18) − λ max [Re[ C ]] λ ( x ) − λ ( x ) (cid:19) , (C7)where C i,j = (cid:104) Ψ | ∂ A ∂x i ∂ A ∂x j | Ψ (cid:105) − (cid:104) Ψ | ∂ A ∂x i | Ψ (cid:105)(cid:104) Ψ | ∂ A ∂x j | Ψ (cid:105) . (C8)Hence, whenever A has a non-degenerate ground state for all x ∈ B , there is a global lower bound to λ min [ ∇ λ min ] . In thenumerical examples shown in Fig. 2 we increase γ until we do not observe any change in the bound. as shown in Fig. 4. Appendix D: Engineering the couplings in ion traps
For completeness, we start this section by summarizing the derivation of the effective spin couplings as they can be generatedwith trapped ions. We consider the ions to be confined in a linear trap with the interactions generated by one bichromatic laser2field for each of the two transversal directions α = x, y at frequencies ω atom ± µ α , respectively. Here, ω atom denotes thelevel splitting of the internal two-level system used to encode the spin, e.g. the hyperfine clock states of an Ytterbium ion. Thedesired laser field can be achieved by two Raman beams per direction with corresponding frequency differences [48]. The basicinteraction of the i ’th ion with a laser field at frequency ω and wave-vector k is given by (assume Ω ≥ ) H = (cid:126) Ω cos( kδx + ϕ + ωt ) σ ix , (D1)where δx denotes the deviation of the ion from its equilibrium position. Therefore the resulting interaction for the aboveconsidered laser fields with the ion chain is described by Hamiltonian H = (cid:126) (cid:88) α = x,y N (cid:88) i =1 Ω αi cos( δk α δx i + µ α t ) σ ix , (D2)where ( δk α ) ξ = δ ξ,α δk is the wave-vector difference of the Raman beams in the direction α and δx i denotes the deviation ofthe i ’th ion from its equilibrium position. We may write δk ( δx i ) α = (cid:80) Nn =1 η αi,n ( a α,n + a † α,n ) with the Lamb-Dicke parameter η αi,n = δk ( b α,n ) i (cid:112) (cid:126) / M ω αn , where M denotes the mass of the ions. Thereby, b α,n denotes the eigenvector of the n ’theigenmode in direction α and a α,n ( a † α,n ) the corresponding annihilation (creation) operators. Within Lamb-Dicke regime where η αi,n is small such that the condition | δk α δx i | (cid:28) holds, and under the rotating-wave approximation justified by the condition ω atom (cid:29) µ α (cid:29) Ω αi , for α = x, y and i = 1 , . . . , N , one finds a state dependent spin-spin interaction mediated by the transversaleigenmodes of the trap. If furthermore | ω αn − µ α | (cid:29) η αi,n Ω αi , also called the “slow” regime, excitation of the vibrational modesare only virtually excited and one obtains a spin Hamiltonian of the form (cid:80) Ni,j =1 J i,j σ ix σ jx with [48] J i,j = ( (cid:126) δk ) M (cid:88) α = x,y Ω αi Ω αj N (cid:88) n =1 ( b α,n ) i ( b α,n ) j µ α − ( ω αn ) . (D3)Next, we outline how Ising couplings of the form J i,j ∝ cos q ( i − j ) , with q = 2 π/z and z integer, may be designed withtrapped ions. Since for all i, j we have cos q ( i − j ) = cos( qi + ϕ ) cos( qj + ϕ ) + sin( qi + ϕ ) sin( qj + ϕ ) for any ϕ , and hencewe may write J ∝ λ c v c v Tc + λ s v s v Ts , (D4)where ( v c ) i = cos( qi + ϕ N ) and ( v s ) j = sin( qj + ϕ N ) where ϕ N = π (1 − ( N + 1) /z ) . Note that v c and v c are orthogonaland therefore proportional to the eigenvectors of J . On the other hand, the form of the matrix ( J i,j ) = J defined by the effectivecouplings Eq. (D3) is mathematically equivalent to J = ( (cid:126) δk ) M (cid:88) α,n β α,n β Tα,n (D5)where ( β α,n ) i := Ω αi ( b α,n ) i / (cid:112) µ α − ( ω αn ) for α = x, y and i, n = 1 , . . . , N . Suppose that we nearly resonantly excite twotransversal modes, one in each of two directions x and y , which we denote by m x and m y , respectively. The coupling matrix isapproximately described by the matrix J = ( (cid:126) δk ) M (cid:16) β x,m x β Tx,m x + β y,m y β Ty,m y (cid:17) and is hence of rank two. In order to mimickthe couplings Eq. (D4) for certain values of q we may choose the Rabi frequencies Ω αn in order to fulfil β x,m x ∝ v c , β y,m y ∝ v s . (D6)For concreteness we assume that the two transversal modes and frequencies are equal, i.e. b x,n = b y,n =: b n and ω xn = ω yn =: ω n for n = 1 , . . . , N . We choose q in order to obtain condition Eq. (D6) for the two transversal modes that correspond to thesecond and third highest frequencies, i.e. for m x = N − and m y = N − (recall that the center of mass mode has the highest3frequency for the transversal modes). In the example shown in the main text we choose N = 15 and z = 16 and find v c = 12 (cid:112) √ √ (cid:112) − √ − (cid:112) − √ −√ − (cid:112) √ − − (cid:112) √ −√ − (cid:112) − √ (cid:112) − √ √ (cid:112) √ and v s = 12 (cid:112) − √ √ (cid:112) √ (cid:112) √ √ (cid:112) − √ − (cid:112) − √ −√ − (cid:112) √ − − (cid:112) √ −√ − (cid:112) − √ (D7)Two modes that resemble these vectors, i.e. that have entries with the same sign, are given by b ∝ . . − . − . − . − . − . − . − − − . . . and b ∝ . . . . . . − . − . − . − . − . − . − . (D8)We can thus choose the Rabi frequencies appropriately such that conditition Eq. (D6) is fulfilled, which results in the Hamiltonian H ∝ (cid:88) i,j =1 cos (cid:16) π i − j ) (cid:17) σ iz σ jz + B (cid:88) i =1 σ ix . (D9) Appendix E: Lower bound to the ground state energy
In the following we describe how to obtain a lower bound using semidefinite programming as described in [17] for thepermutation invariant case. The ground state energy of a Hamiltonian H may be expressed as the minimization of tr[ H(cid:37) ] overdensity matrices (cid:37) ≥ with tr[ (cid:37) ] = 1 . Now, consider a set of operators { X n } Mn =1 with the property that H = (cid:88) m,n H mn X † m X n (E1)for an M × M matrix H . Note that, for any density matrix (cid:37) the matrix X with entries X mn = tr[ X † m X n (cid:37) ] is positivesemidefinite. Hence, a lower bound to the ground state energy of H is given by min X ≥ tr[ HX ] where the optimization is overpositive semidefinite matrices X that may follow additional constraints imposed by relations among the operators { X † m X n } m,n .In the present problem, when we consider spin-1/2 particles, the Hamiltonian is given by4 H = 14 (cid:88) m,n cos ( q ( m − n )) σ mx σ nx − (cid:88) n ( s cos qn + t sin qn ) σ nx + s + t − B (cid:88) n σ nz . (E2)In the following we concentrate on the case q = 0 . The Hamiltonian can be decomposed into a direction sum according thedecomposition of the Hilbert space into irreducible representations of SU (2) . With J x , J y , and J z denoting the spin- J operators,we may thus consider H J = J x − sJ x + s − B J z . (E3)Furthermore we define the set { X n } n to be { , J x , J y , J z , J x J x , J x J y , . . . , J z J z } , i.e. it consist of first and second momentsof the the spin operators, and the identity. The operators fulfil the commutator relations [ J α , J β ] = i (cid:15) αβγ J γ , where α, β, γ ∈{ x, y, z } and (cid:15) αβγ denotes the Levi-Civita symbol, and the relation J x + J y + J z = J ( J + 1) . As mentioned before, we usethese relations to linearly constrain the matrix X defined above, and may denote, generally, the set of all matrices that fulfil theseconstraints by B . Then the SDP min tr( H J X ) X ≥ , X ∈ B , (E4)where H J = − s − B O × − s s ... B O × · · · O × , (E5)provides a lower bound the ground state of H J for a fixed value of s . For a fixed number of spins, we have to take into account all J in the decomposition of the Hilbert space, i.e. J = N , N − , N − , . . . in order to obtain a lower bound for the HamiltonianEq. (E2). Here, O m × n denotes the matrix of dimension m × n with all zero entries. The problem Eq. (E4) can be solved usingtools from convex optimization [49, 50]. We use a separate optimization to find the minimum over s . Note that the dimension ofmatrix X does not increase with J .We can also use directly the relations among the Pauli matrices without decomposing the Hilbert space into irreducible repre-sentations. This will also result in an optimization where the involved matrices are of fixed dimension. This has the advantage,that for a fixed number of spins, we, indeed, only have to take into account a single optimization. To define the set { X n } n wechoose (in the following order) the operators σ lα , and σ mβ σ nγ , with α = x, y, z , ( β, γ ) = ( x, x ) , ( x, y ) , ( x, z ) , ( y, y ) , ( y, z ) and ≤ l, m (cid:54) = n ≤ N . With this choice we may write the coefficient matrix of the Hamiltonian H as H = I N − i B i B − i s O N × N ( N − i s s O N ( N − × N O N ( N − × N ( N − , (E6)where I d is defined as the matrix of dimension d × d with all entries equal to 1. The moment matrix X can be parametrized as X = X x,x X x,y X x,z X x X x,xx X x,xy X x,xz X x,yy X x,yz X † x,y X y,y X y,z X y X y,xx X y,xy X y,xz X y,yy X y,yz X † x,z X † y,z X z,z X z X z,xx X z,xy X z,xz X z,yy X z,yz X † x X † y X † z X , X † xx X † xy X † xz X † yy X † yz X † x,xx X † y,xx X † z,xx X xx X xx,xx X xx,xy X xx,xz X xx,yy X xx,yz X † x,xy X † y,xy X † z,xy X xy X † xx,xy X xy,xy X xy,xz X xy,yy X xy,yz X † x,xz X † y,xz X † z,xz X xz X † xx,xz X † xy,xz X xz,xz X xz,yy X xz,yz X † x,yy X † y,yy X † z,yy X yy X † xx,yy X † xy,yy X † xz,yy X yy,yy X yy,yz X † x,yz X † y,yz X † z,yz X yz X † xx,yz X † xy,yz X † xz,yz X † yy,yz X yz,yz , (E7)5where the block matrices are of the form X α,β ∈ C N × N , X α ∈ C N , X α,βγ ∈ C N × N ( N − and X αβ,γδ ∈ C N ( N − × N ( N − and every block represents a moment matrix with entries X α,β = (cid:104) ˆ σ α ˆ σ β (cid:105) , X α,βγ = (cid:104) ˆ σ α ˆ σ β ˆ σ γ (cid:105) and X αβ,γδ = (cid:104) ˆ σ α ˆ σ β ˆ σ γ ˆ σ δ (cid:105) ,respectively. The blocks are linked with each other through the algebraic relations σ α σ β = δ α,β + (cid:80) γ = x,y,z i (cid:15) αβγ σ γ that wetranslate into constraints on X . For example, ( X x,xx ) l,mn = σ lx σ mx σ nx and, also, ( X xy,xz ) l,m,n,m = i σ lx σ mx σ nx , for l (cid:54) = n and l (cid:54) = m . Therefore, we require that (cid:80) l (cid:54) = n,l (cid:54) = m ( X x,xx ) l,mn + i( X xy,xz ) l,m,n,m = 0 . We define the matrix J d := I d − d × d .Furthermore, let A ( n ) and B ( n ) , n = 1 , . . . , , be the matrices with entries, respectively, given by ( A (1) ) l,mn = δ l,m , ( A (2) ) l,mn = δ l,n , ( A (3) ) l,mn = 1 − δ l,m − δ l,n , (E8)and ( B (1) ) kl,mn = δ k,m δ l,n , ( B (2) ) kl,mn = δ k,m (1 − δ l,n ) , ( B (3) ) kl,mn = (1 − δ k,m ) δ l,n , ( B (4) ) kl,mn = δ k,n δ l,m , ( B (5) ) kl,mn = δ k,n (1 − δ l,m ) , ( B (6) ) kl,mn = (1 − δ k,n ) δ l,m , (E9)and B (7) = I N ( N − − (cid:80) n =1 B ( n ) . Next, we introduce the matrix W with the same block structure as Eq. (E7) and blocksof the form W α,β = w (1) α,β + w (2) α,β J , W α,βγ = (cid:80) n =1 w ( n ) α,βγ A ( n ) and W αβ,γδ = (cid:80) n =1 w ( n ) αβ,γδ B ( n ) , respectively, and werequire that w (2) yz,xz = w (7) xz,xz = w (7) yz,yz = w (7) yz,xy = w (7) xy,xy = w (7) xx,xx = w (7) yy,yy , w (3) x,xx + i w (3) xy,xz , w (7) yz,xy + w (7) yy,xz , w (7) xy,xy + 2 w (7) xx,yy , − i w (2) yz,xx − i w (5) yz,xx + w (3) z,xz − i w (6) xy,xz , w (3) y,yz + w (3) z,yy + i w (2) xy,yy + i w (5) xy,yy + i w (6) yz,xz , w (3) z,xx − i w (3) xy,xx − i w (6) xy,xx + w (3) x,xz + i w (5) yz,xz , w (3) xyy + w (3) y,xy + i w (3) yz,xx + i w (6) yz,xx − i w (3) yz,yy − i w (6) yz,yy − i w (5) xy,xz , w (2) x,z − i w (1) z,yz + i w (2) x,xy + w (1) yz,xy − i w (1) y,xx − i w (2) y,xx + ( n − w (5) yz,xy , w x + i w (1) y,z + ( N − (cid:16) w (2) y,xy + w (1) x,xx + w (2) x,xx − w (1) yz,yy − i w (4) yz,yy + w (2) z,xz + i w (1) xyxz (cid:17) , w z + i w (1) x,y + ( N − (cid:16) w (1) y,yz − i w (1) xy,xx − i w (4) xy,xx + i w (1) xy,yy + i w (4) xy,yy + w (1) x,xz − i w (1) yz,xz (cid:17) , w (2) z,z + w (4) xy,xy + 2i (cid:16) w (1) x,yz − w (1) xx,yy − w (1) xx,yy − w (4) xx,yy − i w (1) y,xz (cid:17) + ( N − (cid:16) w (2) yz,yz + w (2) xz,xz (cid:17) , w (2) y,y + w (4) xz,xz − (cid:16) w (2) x,yz + i w (1) z,xy + w yy (cid:17) + ( N − (cid:16) w (3) yz,yz + w (2) xy,xy + w (2) yy,yy + w (3) yy,yy + w (5) yy,yy + w (6) yy,yy (cid:17) , w (2) x,x + w (4) yz,yz − (cid:16) w xx + i w (2) z,xy + i w (2) y,xz (cid:17) + ( N − (cid:16) w (3) xy,xy + w (2) xx,xx + w (3) xx,xx + w (5) xx,xx + w (6) xx,xx + w (3) xz,xz (cid:17) , w (2) x,z − i w (1) z,yz + i w (2) x,xy + w (1) yz,xy − i w (1) y,xx − i w (2) y,xx + w xz + w (2) xx,xz + w (6) xx,xz + w (1) yy,xz + w (4) yy,xz + ( n − w (5) yz,xy . (E10)By the Pauli algebraic relations, the conditions Eq. (E10) guarantee that tr[ WX ] = N (cid:80) w (1) α,α + 1 + N ( N − (cid:80) w (1) αβ,αβ ,where the last sum runs over pairs ( α, β ) = ( x, x ) , ( x, y ) , ( x, z ) , ( y, y ) , ( y, z ) . Therefore, for all X tr[ HX ] ≥ max (cid:110) N (cid:88) w (1) α,α + 1 + N ( N − (cid:88) w (1) αβ,αβ : H ≥ W (cid:111) , (E11)where the maximum is over w ( n ) α,β , w ( n ) α,βγ and w ( n ) αβ,γδ .Now, we set out to show that the eigenvalues of H − W can be determined efficiently, i.e. the problem reduces to determinethe eigenvalues of a × matrix. We start by explicitly constructing a basis. It will be useful that we can represent any vector6 v ∈ C N ( N − with entries v n , n = 1 , . . . , N ( N − equivalently by the matrix v = v N · · · v ( N − +1 v v ( N − +2 ... ... . . . ... v N − v N − · · · . (E12)In the following (cid:8) e dj (cid:9) dj =1 will denote the standard basis of C d . With this, we introduce E dj,k = e dj ( e dk ) T . Now for any pair ( k, l ) we define T j,k = (cid:18) d − E dj,j E dj,k E dk,j d − E dk,k (cid:19) (E13)and T j,k = d − E dj,j E dj,k E dk,j d − E dk,k , (E14)for N even and odd, respectively. Now, let v ( n ) j,k , n = 1 , , be the vectors with matrix representation given by ( d = (cid:98) N/ (cid:99) ) v (1) = 1 √N J N , v (2) j,k = 1 √N T j,k (cid:18) J d −J d (cid:19) T j,k , v (3) j,k = 1 √N T j,k (cid:18) I d −I d (cid:19) T j,k , (E15)if N is even, where N = N ( N − , N = N ( N/ − and N = N / , and v (1) = 1 √N J N , v (2) j,k = 1 √N T j,k J d d Td − Td − d −J d T j,k , v (3) j,k = 1 √N T j,k d I d − Td Td −I d − d T j,k (E16)if N is odd, where N = N ( N − , N = ( N − N − / and N = N ( N − / . For N odd we additionally define v (2)1 , (cid:98) N (cid:99) = 1 √N T , (cid:98) N (cid:99) J d d Td − Td − d −J d T , (cid:98) N (cid:99) , v (2)1 , (cid:98) N (cid:99) = 1 √N T ,N/ d I d − Td Td −I d − d T , (cid:98) N (cid:99) , (E17)where T , (cid:98) N (cid:99) = d − E d , e d e dT d . (E18)Furthermore, let w (4) j,k ;1 = c k c Tl + c l c Tk + s k s Tl + s l s Tk , k ≤ l, k, l = 1 , . . . , (cid:22) N (cid:23) ,w (4) j,k ;2 = c k c Tl + c l c Tk − s k s Tl − s l s Tk , k ≤ l, k, l = 1 , . . . , (cid:22) N − (cid:23) ,w (4) j,k ;3 = c k s Tl + s k c Tl + c l s Tk + s l c Tk , k ≤ l, k, l = 1 , . . . , (cid:22) N − (cid:23) ,w (4) j,k ;4 = c k s Tl + s l c Tk − c l s Tk − s k c Tl , k < l, k, l = 1 , . . . , (cid:22) N (cid:23) , (E19)and w (5) j,k ;1 = c k c Tl − c l c Tk ,w (5) j,k ;2 = c k s Tl − s l c Tk ,w (5) j,k ;3 = s k s Tl − s l s Tk . (E20)7By Lemma 3 below, we can form linear combinations of, respectively, v ( n ) j,k for n = 1 , , and w ( n ) j,k ; n for n = 4 , so that thecorresponding vectors form an orthonormal basis { ˜ v ( n ) j } j,n of C N ( N − . Additionally, u j := ( A (1) + A (2) )˜ v (2) j / (cid:107) ( A (1) + A (2) )˜ v (2) j (cid:107) defines a basis of C N . We find u Tj J u j = ( N − δ ,j , u Tj A (1) ˜ v ( n ) j = (cid:112) N − δ n, δ j,k , u Tj A (2) ˜ v ( n ) j = − (cid:112) N − δ n, δ j,k u Tj = √ N δ j, , ˜ v ( n ) Tj = √ N − δ n, ˜ v ( n ) j B (2) ˜ v ( m ) k = ( N − δ j,k δ n, m, n = 1 , , v ( n ) j B (2) ˜ v ( m ) k ) m =2 , n =2 , = δ j,k (cid:18) α − β − β γ (cid:19) , (˜ v ( n ) j B (3) ˜ v ( m ) k ) m =2 , n =2 , = δ j,k (cid:18) α ββ γ (cid:19) , (˜ v ( n ) j B (5) ˜ v ( m ) k ) m =2 , n =2 , = δ j,k (cid:18) − α − ββ γ (cid:19) , (˜ v ( n ) j B (6) ˜ v ( m ) k ) m =2 , n =2 , = δ j,k (cid:18) − α β − β γ (cid:19) , (E21)Since all matrix elements in Eq. (E21) are proportional to δ j,k we can form a basis of C N +1+5 N ( N − with H − W blockdiagonal and with all blocks equal. Each block is of size × .We split the proof of Lemma 3 by considering first two preparatory Lemmas. Lemma 1 (i) There are linear combinations ˜ v (4) l = (cid:80) j,k,m c ljkm w (4) j,k ; m for l = 1 , . . . , N − , such that the vectors ˜ v (4) j defined via Eq. (E12) are orthonormal and span a ( N − N − / − -dimensional subspace. Furthermore, (ii) the vectorscorresponding to w (5) j,k ; m span a ( N − N − / -dimensional subspace. Proof. (i) The diagonal entries of w (4) j,k ; m for m = 1 , . . . , are given by diag( w (4) j,k ;1 ) l = ( c j − k ) l , diag( w (4) j,k ;2 ) l = ( c j + k ) l , diag( w (4) j,k ;3 ) l = ( s j − k ) l , diag( w (4) j,k ;4 ) l = ( s j + k ) l . (E22)Since {{ c j } j , { s j } j } are mutually orthogonal, we may consider linear combinations of w (4) j,k ; m that have equal diagonals. Count-ing reveals that these matrices may have one of N distinct diagonals. It thus follows by the fact that w (4) j,k ; m are mutually orthog-onal in the Hilbert-Schmidt inner product, that there are ( N − N − / − linear independent linear combinations withzeros on the diagonal. (ii) The matrices w (5) j,k ; m are mutually orthogonal with zeros on the diagonal. Hence, counting shows thatthey span a subspace of dimension ( N − N − / . Lemma 2
There are coefficients c ljk such that ˜ v ( n ) l = (cid:80) j,k c ljk v ( n ) jk define two sets of orthonormal vectors { ˜ v ( n ) j } j , for n = 2 , ,via the correspondence in Eq. (E12) , that, respectively, span ( N − -dimensional subspaces. Proof.
In order to show that v (2) j,k and v (3) j,k defined through the Eqs. (E15), (E16) and (E17) span a ( N − -dimensionalsubspace we calculate the Gramian matrices of the two sets of vectors. Defining for any NG = 1 d (( d − I d ⊗ I d + 2( d ⊗ I d + I d ⊗ d )) , (E23)with d = (cid:98) N/ (cid:99) , we find for N even that ( G ) jk,lm = (cid:104) v (2) j,k , v (2) l,m (cid:105) = (cid:104) v (3) j,k , v (3) l,m (cid:105) . An explicit calculation of the eigendecompo-sition of G reveals that it is of rank d − N − . For N odd we find that the Gramian matrices for both sets of vectors againcoincide and are given by (in the following the first row and column correspond to v (2)((3))1 ,d ) G = (cid:18) g T g G (cid:19) , (E24)where g = 2( N − (cid:18) N − d ⊗ d + e d ⊗ d (cid:19) (E25)8and G given in Eq. (E23). Since G has rank N − , it follows that G is of rank N − . The two sets, for n = 2 and n = 3 , haveidentical Gramian matrices and hence the statement follows. Lemma 3
The set { ˜ v ( n ) j } j,n is an orthonormal basis of C N ( N − . Proof.
For n = 1 , . . . , , ˜ v ( n ) j are eigenvectors of B (2) + B (3) to the eigenvalue λ n , where λ = 2( N − , λ = N − , λ = N − and λ , = − . In addition, these vectors are eigenvectors of B (4) to the eigenvalues τ n , with τ , , = 1 and τ , = − .Thus ˜ v ( n ) j mutually orthogonal for all j, nj, n