Certifying the quantum Fisher information from a given set of mean values: a semidefinite programming approach

We introduce a semidefinite programming algorithm to find the minimal quantum Fisher information compatible with an arbitrary dataset of mean values. This certification task allows one to quantify the resource content of a quantum system for metrology applications without complete knowledge of the quantum state. We implement the algorithm to study quantum spin ensembles. We first focus on Dicke states, where our findings challenge and complement previous results in the literature. We then investigate states generated during the one-axis twisting dynamics, where in particular we find that the metrological power of the so-called multi-headed cat states can be certified using simple collective spin observables, such as fourth-order moments for small systems, and parity measurements for arbitrary sys-tem sizes.


Introduction
Quantum states of multipartite systems may represent useful resources for specific quantum technology applications [14].It is therefore a central task for quantum information science to provide robust and reliable protocols to estimate the resource content of a quantum state ρ prepared in a given device.In several situations, the quan-Guillem Müller-Rigat: guillem.muller@icfo.euIrénée Frérot: irenee.frerot@lkb.upmc.frtum advantage is a genuine quantum many-body effect, such as in quantum computing [37], quantum simulation [3,17,29] or quantum metrology [34].Hence, in the context of quantum certification -namely, the task of estimating the resource content of a given experimental preparation -, a pervasive challenge is that of scalability.Given that quantum state tomography quickly becomes unfeasible for systems of increasing size, this challenge signifies that one must base the certification protocols on partial information, typically the knowledge of a few mean values of collective observables.A certification algorithm then represents a conservative approach, which quantifies the minimal resource content over all quantum states compatible with the available partial information.A key property of all quantum resource quantifiers is that they are convex functions f (ρ) of the quantum state; namely, a quantum resource [11], be it quantum coherence [43], thermodynamic resources [24], quantum entanglement [23], Einstein-Podolski-Rosen steering [46], Bell non-locality [10,40], etc., cannot be produced (nor amplified) by mere statistical mixing of different preparations.Convexity turns out to be a key technical property for practical certification algorithms, as imposing compatibility of the quantum state with a given set of mean values -namely, linear constraints of the form Tr( Da ρ) = ⟨ Da ⟩ for a = 1, . . ., K where Da are some observables forming the available data -, maintains the convexity of f .Therefore, the conservative estimate of the resource content of ρ compatible with the data can be cast altogether as a convex optimization algorithm, hence presenting a unique and well-defined optimum.
In this work, we focus on estimating the quantum Fisher information (QFI) associated to unitary transformations of ρ [7,9].The QFI quantifies the usefulness of the state for quantum metrology (namely, the task of measuring small parameters, such as magnetic fields or gravitational fields, beyond the capabilities of classical systems of the same size) [34].While previous works have achieved significant progress in this respect [5,18], no existing approach fully exploits the convexity of the QFI to achieve an unbiased and certified estimate of the QFI given a set of mean values.In the present paper, exploiting the link of the QFI with the Uhlmann fidelity, we develop a semi-definite programming (SDP) algorithm to solve this problem.SDP represents a well-studied special class of convex optimization problems, for which efficient numerical routines are available [41].We apply our new approach to several problems of relevance to quantum metrology, such as Dicke states [34,52].We then explore the metrological resource produced during the out-of-equilibrium dynamics of the celebrated one-axis-twisting dynamics (OAT) [26].In this case, we obtain the unexpected result that fewparticle multi-headed cat states and their metrological power can be fully assessed by low-order moments of collective spin operators.We also find that Heisenberg-scaling metrological power can be optimally exploited by parity measurements throughout the OAT evolution.We finally illustrate our approach in small quantum spin chains with nearest-neighbour interactions.

Semidefinite programming solution
Quantum Fisher information.-Inthe context of quantum interferometry, the QFI evaluates the sensitivity of the state ρ to unitary transformations of the form ρ(θ) = e −iθ Ĝ ρe +iθ Ĝ (linear encoding of the phase θ) [34].Here, Ĝ is the observable generating the unitary transformation.In the specific case of magnetometry, for instance, Ĝ is some component of the magnetic moment of the system, and θ would incorporate both the time of evolution and the strength of the magnetic field.Via the quantum Cramér-Rao bound [13,20,38], the QFI (F Q [ρ, Ĝ]) quantifies the ultimate sensitivity that can be achieved on the estimation of θ by making measurements on the state ρ(θ), namely [34]: Here, (∆θ) 2 est. is the variance of the estimator of the phase θ, and Ô is some observable whose mean value ⟨ Ô⟩ = Tr[ Ô ρ(θ)] changes with θ as ∂ θ ⟨ Ô⟩ = i⟨[ Ĝ, Ô]⟩, and (∆ Ô) 2 = ⟨ Ô2 ⟩ − ⟨ Ô⟩ 2 is the variance of Ô.We denote by Ξ O the sensitivity of the observable Ô in the phase estimation task.The Cramér-Rao bound [Eq.(1)] is valid for arbitrary observables Ô.
Statement of the certification problem.-Theproblem we aim at solving is to certify the minimal QFI compatible with a set of K expectation values ⟨ D⟩ = (1, ⟨ D1 ⟩, ⟨ D2 ⟩, ..., ⟨ DK ⟩) (thereafter called data): where ρ is a quantum state acting in a Hilbert space H of a given dimension (namely: a Hermitian, semidefinite positive, unit-trace operator; note that we implement the unit-trace constraint with D0 = I).As mentioned in the introduction, this problem is well-posed since the QFI is a convex function of the state and the linear constraints corresponding to the data maintain the convexity property.
SDP solution to the certification problem.-Toarrive at our practical algorithm, we use the fact that the fidelity can be computed as a semidefinite-program (SDP) [41,47], Furthermore, it is straightforward to impose compatibility of the state with the data in the SDP formulation, yielding the final SDP algorithm (for a given fixed and sufficiently small δθ): (5) The SDP algorithm delivers both the maximal fidelity and an optimal quantum state ρ * , whose QFI is evaluated either as ), or directly using the optimal state ρ * .In Appendix A.1, we further discuss the accuracy of this certificate for small yet finite δθ.
Alternative SDP algorithm.-Wepropose an alternative SDP formulation to solve our central problem Eq. (2).The starting point is the definition of the inverse QFI as the minimal variance of a locally unbiased observable Ô [21,54]: Here, the operator Ô can be decomposed as Ô = i θ(i) |i⟩ ⟨i|, where {|i⟩} is the measurement basis and θ(i) is the value of the estimator of θ associated with the i-th measurement outcome.The maximum of Eq. ( 6) over all states ρ compatible with a given set of expectation values ⟨ D⟩ is equal to the inverse of the solution to problem Eq. (2).The minimization over Ô, Eq. (6), can be reformulated as a maximization problem using the strong duality theorem for SDP.Finally, the resulting double maximization problem (over ρ and over variables dual to Ô) can be recasted as a single SDP-see Appendix A for the derivation and exact formulation of this SDP, and for the comparison with the fidelity based approach.
Exploiting the metrological usefulness.-Afterobtaning the minimal QFI for fixed arbitrary data D, one may ask how to operationally exploit the corresponding metrological resource.In other words, one aims at finding an explicit observable Ô with sensitivity Ξ −2 O as close as possible to the QFI.From the optimal state ρ * , such observable can be provided with the certificate.This privileged observable corresponds to the symmetric logarithmic derivative (SLD) [33] S which is derived from the relation Notice, however, that in practice the SLD observable might be impractical to measure, and the practical optimal observable must be found on a case-by-case basis [18].
Taking advantage of symmetries.-Whilebeing a convex-optimization algorithm, the SDP of Eq. (5) involves the complete density matrix ρ as a variable of the problem, and hence is hardly scalable beyond about N = 10 qubits.However, in several cases of practical interest, both the generator Ĝ and the observables D respect elementary symmetries, such as translation invariance, or even invariance under permutation of the elementary degrees of freedom.Such symmetries can be accommodated in the SDP by splitting the optimization into different symmetry sectors, which yields a more scalable algorithm.As further discussed in Appendix B, if both the observables D and the generator Ĝ obey a certain symmetry, it is sufficient to optimize over symmetric states of the form ρ = ⊕ J ρ(J) , where ρ(J) is a (non-normalized) quantum state acting on the symmetry sector J (with Tr[ρ (J) ] the probability to find the system in the symmetry sector J).In essence, taking advantage of the symmetries allows one to block-diagonalize the problem with respect to the symmetry sectors, whose dimensions are typically much smaller than the complete manybody Hilbert space, yielding an SDP algorithm with a more favorable scaling (see Appendix B for a specific example on permutation invariance).

Alternative approaches.-The central problem
Eq. (2) has appeared in different contexts throughout the literature.In Appendix C, we give a concise review of past efforts to solve it and discuss its comparison to the algorithms proposed here.
Beyond unitary evolutions.-Inthis paper, we focus on situations where the metrological protocol consists in a unitary evolution of the probe: ρ(θ) = e −iθ Ĝ ρe +iθ Ĝ.However, the SDP of Eq. ( 5) can be extended to a generic (namely: non-unitary) parameter encoding operation by defining ρ± (θ) = ρ(θ) ± δθ L[ρ(θ)], where L is the corresponding Lindbladian, This can be used to incorporate noise and other imperfections of the metrological task [28], or to model more general open-system dynamics.Notice that in this case, the QFI becomes θ-dependent.
In the remainder of the paper, we illustrate our approach by providing the certified metrological usefulness in several classes of states of experimental relevance.We focus on states of an ensemble of N quantum spin-1/2 particles.In Section 3, we consider Dicke states, which in contrast to spin-squeezed states have a vanishing mean spin.In Section 4, we then consider the family of states generated during the so-called one-axis twisting dynamics [26] using collective spin observables of increasing orders.

Metrological usefulness of unpolarized states
In this section, we focus on unpolarized states, namely states of vanishing mean spin: individual spin-1/2 operators.Such unpolarized states may offer a metrological advantage, the most prominent example being balanced Dicke states, such that ⟨ Ĵ2 ⟩ = (N/2)(N/2 + 1) (maximal total spin) and ⟨ Ĵ2 x ⟩ = 0 (eigenstate of the collective spin along x with eigenvalue 0) [22].
The balanced Dicke state displays an exquisite sensitivity to rotations around any axis in the yz plane, with QFI given by F Q [ρ, Ĵz ] = 4⟨ Ĵ2 z ⟩ = N (N/2 + 1).To exploit this sensitivity, a possible approach is to measure ⟨ Ĵ2 x ⟩(θ), allowing one to reach the sensitivity bound of the QFI [4,25,30,55].In this protocol, one needs to estimate the variance of Ĵ2 x , and hence to access fourth-order moments of the collective spin.Here, we are instead interested in quantifying the metrological usefulness of unpolarized states using only second moments of collective spins.Our motivation for doing so is that in an experiment, it is impossible to prepare the balanced Dicke state with unit fidelity, and it is natural to certify the QFI using only the knowledge of ⟨ Ĵ2 x ⟩, ⟨ Ĵ2 y ⟩, ⟨ Ĵ2 z ⟩ .

Case ⟨ Ĵ2
x ⟩ = 0.-In particular, for ⟨ Ĵ2 x ⟩ = 0, as we prove in Appendix D, the QFI with respect to the generator Ĵn ⊥ is equal to 4⟨ Ĵ2 n ⊥ ⟩, where n ⊥ is any direction in the yz plane.The constraint ⟨ Ĵ2 x ⟩ = 0 implies that any state compatible with it is an eigenstate of Ĵx and therefore it is rotationally invariant around the x axis.Consequently, Hence, such observation extends the property of the balanced Dicke state (which is a pure state of maximal total spin ⟨ Ĵ2 ⟩ = (N/2)(N/2 + 1)) to general mixed states with arbitrary total spin.

Case ⟨ Ĵ2
x ⟩ > 0.-We then investigate the effect of a finite but small nonzero value of ⟨ Ĵ2 x ⟩ on the QFI of Ĵz .In Appendix E, we study the minimal F Q ( Ĵz ) compatible with a given ⟨ Ĵ2 x ⟩ and a given total spin ⟨J 2 ⟩ ∼ N 2 .Our main conclusion is that for ⟨ Ĵ2 x ⟩ ≳ 1/4, the Heisenberg scaling F Q ( Ĵz ) ∼ N 2 is lost.This result suggests that certifying the metrological power of Dickelike states requires knowledge beyond second moments.
Questioning the Dicke squeezing parameter.-Herewe question the validity of a proposed parameter to estimate the metrological usefulness of unpolarized states, namely the so-called Dicke squeezing parameter introduced in Ref. [52]: This parameter was observed numerically to quantify the metrological usefulness of a specific class of states (namely: superpositions of Dicke states close to the balanced Dicke state), but to our knowledge, no proof of the general relevance of the Dicke squeezing parameter has been given in the literature.Note that, as per the previous observation, for ⟨ Ĵ2 x ⟩ = 0 one has that Our approach allows us to provide a rigorous lower bound to the QFI, as opposed to ξ D which only provides a conjectured lower bound based on numerical observations.In Figure 1, we compare (2ξ D ) −1 with the minimum QFI of Ĵz compatible with a reduced total spin of ⟨ Ĵ2 ⟩ = 0.7(N/2)(N/2 + 1) as a function of ⟨ Ĵ2 x ⟩. x ⟩ we obtain the minimum QFI compatible with ⟨ Ĵx ⟩ = 0, ⟨ Ĵ2 x ⟩, ⟨ Ĵ2 y + Ĵ2 z ⟩ = 0.7(N/2)(N/2 + 1) − ⟨ Ĵ2 x ⟩, and N = 16 (solid orange line).We compare it with the Dicke squeezing parameter [Eq.( 7)] (dashed blue line).
We find that the Dicke squeezing parameter [Eq.(7)] only captures the metrological usefulness for very small values of ⟨ Ĵ2 x ⟩, but in general it largely overestimates the QFI of the state.Therefore, the main outcome of this study is that this parameter cannot be used as a faithful estimate of the metrological usefulness of the system without further assumptions.Notice though that in Ref. [4], a proper squeezing parameter tailored to Dicke states has been derived, which constitutes a valid lower bound of the QFI, and allows one to strongly improve the estimation of the metrological usefulness of the state.However, this criterion involves fourth moments of the collective spin, while our goal here was to include data up to the second moments only.It would be interesting in a future work to investigate the tightness of the squeezing parameter of Ref. [4] by including fourth-order moments in our SDP algorithm.
The OAT protocol has been widely investigated both theoretically and experimentally to prepare metrologically useful entangled many-body states in atomic ensembles [34].Its realization in different regimes of structured systems such as optical lattices has also been explored in the context of cold-atom and Rydberg-atom quantum simulators [8,12,15,35,36,50].The OAT dynamics generates spin squeezing at short times t ≲ 1/ √ N , with strength ξ −2 R ∼ N 3/2 .At longer times, the mean spin ⟨ Ĵx ⟩ vanishes and quantum correlations are not simply reflected in spin squeezing.While the theoretical metrological usefulness continues to increase with time, it becomes more demanding to certify or exploit it, as the correlations are spread in higher moments of the collective spin fluctuations.In particular, at times t = π/q, the state is a so-called q-headed cat state (MHCS), namely an equal-weight superposition of q CSS in the xy plane (θ = π/2) at angles ϕ m = 2πm/q for m ∈ {0, 1, .., q−1} [1,42].This includes, for t = π/2, the conventional cat state (with two heads) marking half of the total period of the dynamics.Before discussing the metrological usefulness of MHCS, and more general states generated during the OAT dynamics, we briefly review the case of spin-squeezed states.
Spin-squeezed states.-Spin-squeezedstates represent the most widely studied class of manybody spin states useful for metrology, and are generated at short time in the OAT dynamics.Those states have a nonzero mean spin ⟨ Ĵ⟩ = ⟨ Ĵn ⟩ ̸ = 0 (with n the direction of the mean spin) and a reduced variance for a spin component transverse to the mean spin ⟨ Ĵ2 n ⊥ ⟩ < ⟨ Ĵn ⟩ 2 /N (with n ⊥ orthogonal to n) [49].For simplicity, let us assume that the mean spin is along x (as is the case in the OAT dynamics), and that the minimal variance orthogonal to x is along s (squeezing axis).The QFI for rotations around as = s × x (antisqueezed axis) is lower-bounded by the spin squeezing parameter [34]: It is natural to investigate the tightness of this bound, namely, to find the minimal F Q [ρ, Ĵas ] compatible with given ⟨ Ĵx ⟩ and ⟨ Ĵ2 s ⟩.This question has been investigated previously in Ref. [5] (see Appendix C for a discussion of the algorithm used in this paper), and the bound is known to be tight except for data points at the boundary of the feasible set (namely, the set of values for ⟨ Ĵx ⟩, ⟨ Ĵ2 s ⟩ for which there exists a quantum state compatible with them), where the bound is tight within ∼ 2% accuracy.
In the remainder of this section, we will use expectation values as produced during the OAT dynamics [Eq.( 8)] to derive, using our SDP algorithm, the sensitivity bound to collective spin rotations.We will compare our obtained bounds with well-known lower bounds of the QFI, such as spin squeezing [Eq.( 9)], or nonlinear generalizations thereof [18].
Second moments.-Westart by taking the first and second moments of collective spin observables as input data to our algorithm, namely ⟨ Ĵa ⟩, and Re{⟨ Ĵa Ĵb ⟩} for a, b ∈ {x, y, z}.Note that since ⟨ Ĵ2 ⟩ = (N/2)(N/2 + 1) is given explicitly by the data, it is sufficient to optimize over the totally symmetric sector whose dimension is N + 1, instead of the complete Hilbert space of dimension 2 N .Such expectation values contain also implicitly the spin squeezing bound of Eq. (9), where the squeezed component Ĵs is used to sense rotations generated by Ĝ = Ĵas .At short times, the orientations {s, as} are contained in the yz plane (orthogonal to the mean spin, which is along x) and depend on time.At longer times, the state wraps around the Bloch sphere and becomes unpolarized.In this case, we use as a generator Ĝ = Ĵy .Minimizing the QFI with respect to the generator Ĝ allows us to compare both our SDP bound and the spin-squeezing bound.The result is displayed in Figure 2.
The results presented in Figure 2 lead us to several important observations.First, there is a gap between our bound and the spin squeezing parameter, demonstrating that using the same second moments as used to build the optimal spin-squeezing parameter, one can in fact certify more metrological usefulness than predicted by the same squeezing parameter.As a matter of fact, adding to the data the second moment along the mean spin (⟨ Ĵ2 x ⟩) in addition to the mean spin itself (⟨ Ĵx ⟩) and to the squeezed second moment (⟨ Ĵ2 s ⟩) is already sufficient to open this gap between the spin squeezing parameter and the SDP bound, and gives the same QFI lower-bound as using the full correlation matrix of the collective spin.We note that an illustrative special case is at As it has a vanishing mean spin, its spin squeezing bound is zero.It features , which implies that it exhibits maximal QFI with respect to rotations generated by Ĵx , F Q [|GHZ⟩ , Ĵx ] = 4(∆ Ĵx ) 2 = N 2 .However, using only second moments it is not possible to certify this metrological usefulness, as they are always compatible with a state invariant under rotations around x such as a statistical mixture of two coherent spin states [GHZ] = (|π/2, 0⟩⟨π/2, 0| + |π/2, π⟩⟨π/2, π|)/2.Nonetheless, one can certify the metrological usefulness of rotations around y (or z), for which , Ĵy ] = N , which is equivalent to considering rotations of a coherent spin state polarized along ±x.Finally, we notice that more generally, at all times at which the mean spin is zero, our SDP algorithm certifies metrological usefulness for rotations around y with a QFI lower bound ∼ N (see Appendix F for details).
In order to further reveal the metrological usefulness of non-spin-squeezed states beyond the small QFI lower-bound reported in Figure 2, it is necessary to consider data involving higher moments of the collective spin.
Fourth moments.-Wethen add to our SDP algorithm the information on third and fourth moments: Re{⟨ Ĵa Ĵb Ĵc ⟩}, Re{⟨ Ĵa Ĵb Ĵc Ĵd ⟩} for a, b, c, d ∈ {x, y, z}, in addition to the data considered in Figure 2. We emphasize that in Ref. [18], a method is proposed to incorporate systematically such quantities in generalized squeezing bounds, quantifying the optimal sensitivity to collective spin rotations for all observables in the span of the data operators.Namely, if the data include mean values of the form {⟨ Ôa ⟩} a = ⟨ Ô⟩ and second moments of it {⟨ Ôa Ôb ⟩} a,b = ⟨ Ô ÔT ⟩}, then Ref. [18] allows one to build the optimal observable of the form O is maximal, as well as the optimal direction for collective spin rotations.Here, if one considers Ô up to second moments, using fourthorder moments in D is sufficient to derive the corresponding bound.The details of the approach of Ref. [18] can be found in Appendix C. Using the said approach, we find that at short times, the optimal generator is Ĝ = Ĵas (as was already the case using only up to second moments of the collective spin); then until t ≲ 0.15π the optimal generator is Ĝ = Ĵx ; then for t ≳ 0.27π the optimal generator is Ĝ = Ĵy .We compare the nonlinear spin squeezing parameter obtained using the method of Ref. [18] to the QFI lower-bound obtained by our SDP, using up to fourth-order moments, and for the same optimal generator.Notice that in general the optimal generator as given by Ref. [18] need not be the same as the one for our approach.The result of the comparison of both sensitivity bounds for N = 10 particles is shown in Figure 3.We also plot the QFI of the OAT state |Ψ(t)⟩ with respect to Ĝ. (Note that in Figure 3, we have added infinitesimal noise to the quantum state which modifies the data as ⟨ D⟩ → (1 − p)⟨ D⟩ + pTr( D)/(N + 1) with p = 10 −4 ; this has a negligible effect on the QFI of the state or the generalized squeezing bound, and avoids singular results in the SDP bound.) Figure 3: Evolution of the metrological bounds from fourth moments of the collective spin during OAT dynamics for N = 10 particles.From bottom to top: in blue (dashed) nonlinear spin squeezing parameter as introduced in Ref. [18], in orange (solid line), the SDP bound, based on the same data.The green dotted line is 4(⟨ Ĝ2 ⟩ − ⟨ Ĝ⟩ 2 ) and corresponds to the QFI of the OAT-evolved state |Ψ(t)⟩ of Eq. ( 8) with respect to the same generator Ĝ as used to compute the other bounds.
We observe consistently that our bound systematically improves over the non-linear squeezing parameter of Ref. [18] as soon as the non-gaussianity of the state becomes significant (that is, at times t/π ≳ 0.1).However, we notice that in contrast to the previous case considering only up to second-order observables, here the gap appears to close when increasing N (see Appendix F for a detailed discussion on the scaling).Furthermore, a striking fact is that close to times t = π/4, π/3, where MHCS are realized, our bound nearly saturates the full QFI of the state.This is remarkable, for MHCS are highly non-gaussian states and are rather characterized by high-order interference fringes stemming from the coherent superpositions of several well-separated coherent spin states, while our bound requires using spin fluctuations only up to order four.We observe that for the 4-MHCS at time t = π/4 and for N = 2, 4, 6, 8, 10, our bound actually reaches the maximal QFI, namely N (N + 1)/2.This result stems from the unexpected fact that by fixing the collective spin fluctuations up to order four, one actually fully specifies the state: there exists an operator Ŵ , built as a linear combination of the data, such that the state corresponds to its extremal eigenvalue: min ρ Tr(ρ Ŵ ) = ⟨4MHCS| Ŵ |4MHCS⟩ = λ min ( Ŵ ), where λ min denotes the minimum eigenvalue of Ŵ with corresponding eigenvector |4MHCS⟩.In Appendix G we outline a systematic procedure to find such operator (if it exists).However, we observe that this result is a finite-size effect, only valid in few-particle systems: indeed, as the number of spins is increased, the peaks for our lower-bound to QFI at times t = π/4, π/3 drop.For the 4-MHCS, the certified QFI drops below N (coherent-state limit) for N = 24.Regarding the 3-MHCS, the certification drops below N already for N = 10 (see Appendix F for further details).For larger particle numbers, higher order moments then become necessary to certify the metrological resource of highly non-gaussian states such as MHCS.
Parity measurements.-Asa last illustration of our algorithm, we consider the use of parity operators Πa = (−1) Ĵa+N/2 , which are paradigmatic operators to exploit the maximal metrological usefulness of the cat (GHZ) state generated at time t = π/2 in the OAT dynamics [27].In analogy with Eq. ( 9), one can construct the parity squeezing parameter: For the GHZ state, |GHZ⟩ = (|π/2, 0⟩ + e iΘ |π/2, π⟩)/ √ 2, the optimal axes are a = x and b contained in the yz plane, yielding Ξ −2 Π = N 2 , which is the maximal sensitivity allowed by quantum mechanics (the so-called Heisenberg limit).Such a result can be also obtained with the approach of Ref. [18] using as operators Ô = { Ĵx , Ĵy , Ĵz , Πx , Πy , Πz }, and including their expectation values and the correlation matrix ⟨ Ô ÔT ⟩.Using this approach, the optimal generator throughout the OAT dynamics is in the anti-squeezed direction Ĵas for times t ≲ 0.27π, and then Ĵx for longer times.As in the previous illustrations, our proposed SDP-based algorithm is used to find a stronger metrological sensitivity based on the same data and with respect to such optimal generator.In Figure 4 the corresponding bounds are compared.There, we depict the maximal QFI associated with collective spin rotations, that is, (four times) the largest eigenvalue of the covariance matrix C ab = Re{⟨ Ĵa Ĵb ⟩} − ⟨ Ĵa ⟩⟨ Ĵb ⟩ for a, b ∈ {x, y, z}.From bottom to top: in blue (dashed) non-linear squeezing parameter derived from Ref. [18] using Ô = { Ĵx , Ĵy , Ĵz , Πx , Πy , Πz } as measurement operators; orange curve (solid line)the corresponding SDP bound using the same data.In green (dotted) parity squeezing parameter of Eq. ( 10) with a = y (generator direction) and b = x (measurement direction).In red (dashed-dotted) is plotted the maximal QFI of the OAT-generated states with respect to collective spin rotations.
We observe that at the times where the generator is Ĵas , our bound captures the full QFI of the state, a strong improvement with respect to the prediction of Ref. [18].In particular, we verify that it is sufficient in our algorithm to input the data ⟨ Πx ⟩, Im{⟨ Πx Ĵas ⟩} to saturate the full QFI of the state.In addition, we have verified that this QFI can be saturated using simply the parity squeezing parameter introduced in Eq. (10), with directions a = as, b = x.This simple observation highlights the possibility to use parity operators to exploit the maximal sensitivity during the whole OAT dynamics before the generation of the GHZ state.
At longer times (t ≳ 0.27π), when the optimal generator as predicted by Ref. [18] is Ĵx , our bound only yields a modest improvement with respect to the non-linear squeezing parameter of Ref. [18].However, we notice our bound can be drastically improved up to times t ≲ 0.35π, again saturating the full QFI of the state, by just changing the generator from Ĵx to Ĵy .Also in this case, the bound is saturated by the parity squeezing parameter of Eq. (10), taking as orientations a = y, b = x (see the green dotted curve of Figure 4).Related findings on the relevance of parity measurements to saturate the QFI through-out the evolution were reported in previous works [19,32], pointing out the crucial role played by parity conservation during the dynamics.

Conclusions
Certifying the resource content of a multipartite quantum system for metrology applications is an important task for quantum technologies.In this work, we introduced semidefinite programming algorithms to find the minimal sensitivity to unitary transformation, as quantified by the quantum Fisher information, compatible with an arbitrary set of mean values.These algorithms are general and improve over previous approaches by fully exploiting the convexity of the problem.
We have shown that one can take advantage of symmetries to drastically reduce the numerical cost of the implementation.As illustrations, we have focused on relevant families of states of spin-1/2 ensembles.In the case of Dicke states, we have shown that a so-called Dicke squeezing parameter introduced in the past [52] overestimates the certifiable metrological usefulness.We have then focused on states generated in the one-axis-twisting dynamics [26], where we obtained several significant results.First, we showed that using only the first and second moments of the collective spin, one can certify more metrological usefulness than given by the Wineland spin squeezing parameter [49].Including up to fourth-order moments, we made the surprising observation that, for up to a dozen of spins, the metrological power of so-called multi-headed cat states can be almost optimally certified.Finally, we found that using the sensitivity of parity operators (which are N -body correlation functions), close to the maximal QFI of the state can be certified throughout the whole dynamical evolution.Our results are directly relevant to many experiments where collective spin fluctuations are sampled and could be used to certify the metrological resource of such systems.Future applications of our approach could include spatially-structured systems, such as those studied in quantum simulators, and Appendix H proposes a first numerical study in this direction.expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union, European Commission, European Climate, Infrastructure and Environment Executive Agency (CINEA), nor any other granting authority.Neither the European Union nor any granting authority can be held responsible for them.

A Two SDP algorithms
In this appendix, we detail the main technical considerations of two SDPs solving the central problem as formalized in Eq. (2).Firstly, for the algorithm based on maximizing fidelity as presented in the main text, we discuss the effect of the small parameter δθ.Secondly, we provide a comprehensive derivation of the second SDP, which solves exactly the problem of Eq. (2) by maximizing the variance of a locally unbiased estimator.Finally, we compare the performance of both algorithms in terms of accuracy and running computational times.

A.1 Finite δθ considerations of the main-text SDP
In order to illustrate the error introduced by a finite δθ in the SDP algorithm, we consider the following example: a system of N = 4 qubits, in which we bound the QFI with respect to the generator Ĝ = Ĵz , imposing compatibility with the arbitrarily-chosen data: ⟨ Ĵx + 0.2 Ĵz ⟩ = 0.6078, ⟨ Ĵ2 y ⟩ = 0.5795.As in the main text, we denote by Ĵa = 4 i=1 σa i /2 the collective spin along orientation a, with σa i the Pauli matrices (a ∈ {x, y, z}).The SDP is solved for a range of values for δθ.For each value of δθ, we compute the QFI via two methods: 1.By inverting the relation ), the minimal QFI can be estimated from the fidelity resulting from the SDP, F Q [ρ * δθ , Ĝ] = (1 − F max )/(δθ) 2 (up to second-order corrections in δθ).

By computing the QFI directly from the optimal ρ *
δθ which is itself obtained as an output of the SDP algorithm.With the spectral decomposition of the state ρ = p p Π(p) , where p are the eigenvalues and Πp the projectors to the respective eigenspaces, the QFI with respect to Ĝ can be obtained with the closed formula, By construction, this method always gives a QFI realizable by a physical quantum state.
We compare such ways to lower-bound the QFI in Figure 5, where we plot both fidelity and QFI as a function of δθ.From Figure 5 we observe that for smaller values of δθ than a certain finite δθ critical , the fidelity is a concave function of δθ.After convergence of the SDP, this leads by means of method 1 to underestimation of the QFI bound.On the other hand, using method 2, a QFI value is given that it is exactly realizable by ρ * δθ .Thus, this method provides an overestimation, generally an upper bound to the method 1.In the limit δθ → 0, in practice δθ ∼ 10 −2 , the two methods converge yielding a tight lower bound.Hereinafter, we will use method 2, because as shown in Figure 5b it is the tightest in the limit of small δθ.

A.2 Alternative SDP algorithm
In this section, we derive an alternative SDP algorithm computing the minimal QFI compatible with a set of mean values.This second algorithm does not involve the fidelity, and has the advantage that it is then no longer necessary to introduce the small parameter δθ as discussed in the previous section.Let us start with the following optimization problem whose solution is the inverse of the QFI [21,54] : The operator Ô can be interpreted as an observable representing the estimator for the unknown phase of the unitary transformation Û (θ) = e −iθ Ĝ.It can be decomposed as Ô = i θ(i) |i⟩ ⟨i|, where {|i⟩} is the measurement basis and θ(i) is the value of the estimator of θ associated with the i-th measurement outcome.The objective function of the problem is the variance of the observable Ô, while the constraints ensure local unbiasedness.The last constraint, Tr(ρ Ô) = 0, can be dropped without affecting the result of the minimization.Indeed, any Ô can be decomposed as Ô = α1 + Ô0 , where Tr(ρ Ô0 ) = 0, changing α does not affect the first constraint, and Tr(ρ Ô2 ) = Tr(ρ Ô2 0 ) + α 2 Tr(1) is clearly minimal when α = 0, which implies that Tr(ρ Ô) = 0. Let us reformulate as an SDP the former problem: As a consequence of the last constraint, matrices Ô and Ô2 are Hermitian, and Ô2 ⪰ Ô2 due to Schur's complement condition.Therefore, the minimum of the objective function is achieved when Ô2 = Ô2 (because ρ is positive-semidefinite), which means that problems Eq.( 12) and Eq. ( 13) are equivalent.
Our task is to minimize the QFI over ρ satisfying certain constraints, so the inverse of the QFI must be maximized and the problem we want to solve is: To simplify the problem, we are going to change it from "max-min" to "max-max" -to achieve this, we will use a strong duality theorem for SDPs, and change the minimization over Ô, Ô2 to a maximization.Let us start with rewriting the minimization problem of Eq. ( 13) into a canonical form: where Φ is a linear, Hermitian-preserving map (notice that X and B are not necessarily of the same dimension).The dual maximization problem, whose solution is the same as of the primal one, reads [48] max where Φ * is the adjoint map, satisfying Tr Φ( X) Ŷ = Tr XΦ * ( Ŷ ) (17) for all Hermitian X, Ŷ .Let d be the dimension of ρ, and Φ be a map from 2d × 2d to (2d + 1) × (2d + 1) Hermitian matrices, defined as The problem Eq. ( 13) is equivalent to Eq. ( 15) and Since X is constrained to be positive-semidefinite, it must also be Hermitian, which, together with equalities X01 = X † 01 and X11 = 1 (which follow from condition Φ( X) = B) implies that X has structure X = Ô2 Ô Ô 1 , where Ô, Ô2 are both Hermitian.Moreover, Tr X Ĉ = Tr ρ Ô2 , and from Φ( X) = B follows that Tr i[ρ, Ĝ] Ô = 1, which shows that the problem Eq. ( 13) is indeed equivalent to the canonical form Eq. ( 15) with parameters given by Eq. (18), Eq. ( 19), Eq. (20).Let us write down the Hermitian matrix Ŷ as where Ŷij are d×d blocks, and plug in this representation of Ŷ into Eq.( 17) with Φ defined by Eq. ( 18).
After expanding both sides of Eq. ( 17) and comparing them, we obtain Moreover, Tr( B Ŷ ) = y 00 + Tr( Ŷ22 ), so the SDP problem dual to Eq. ( 13) is y 00 is real, and Ŷ11 , Ŷ22 are Hermitian because Ŷ is Hermitian.We are now almost ready to express problem Eq. ( 14) as a single SDP, in which the maximization is performed also over ρ.The only problem is that variables y 00 and ρ are multiplied in the constraint of Eq. (23).To deal with this issue, let us make the following observation: when block matrix satisfies A 00 A 01 A 10 A 11 ⪰ 0, then x 2 A 00 xA 01 xA 10 A 11 ⪰ 0 for any real x.This is a direct consequence of Schur's complement condition for matrix positivity.This observation allows us to change the constraint from Eq. ( 23) to where M = y 00 Ŷ11 is an arbitrary Hermitian matrix, which will be our new independent variable (notice that neither Ŷ11 nor M enter into the problem's objective).Let us introduce a new matrix variable ρ′ = y 2 00 ρ satisfying Tr(ρ ′ ) = y 2 00 , and impose a condition which is equivalent to Tr(ρ ′ ) ≥ y 2 00 .When the objective function y 2 00 + Tr( Ŷ22 ) is maximized, then the above inequality becomes equality.This allows us to treat y 00 and ρ′ as independent variables.Importantly, constraint Eq. ( 24) is linear with ρ′ .Finally, the constraints Tr( D ρ) = ⟨ D⟩ are equivalent to Tr( D ρ′ ) = ⟨ D⟩Tr(ρ ′ ).All these observations allow us to write the expression for the inverse of minimal QFI, as a single SDP:

A.3 Comparison
In this section, we compare the efficiency and accuracy of the two previously-introduced SDP algorithms.As an example, we consider a system of N qubits in a totally symmetric state, i.e. ⟨ Ĵ2 ⟩ = (N/2)(N/2 + 1), where Ĵ2 is the total spin.In the data, we include a mean spin ⟨ Ĵx ⟩ = 0.8 and the fluctuations in an orthogonal direction ⟨ Ĵ2 y ⟩ = 0.4.We find the minimal QFI with respect to the generator Ĵz compatible with such expectation values via the two methods, varying N between N = 4 and N = 46 (for even values).For each N we record both the computing time and corresponding QFI lower-bound.The number of particles N is related to the dimension of the totally symmetric subspace as D = N + 1.The corresponding results are shown in Figure 6.
In Figure 6, we observe how the fidelity-based SDP is marginally more efficient than the variancebased SDP.This fact stems from the extra parameter and constraint existing in the latter approach.In panel b) we study the discrepancy between the two bounds.We find that due to the finite δθ the fidelity-based method slightly overestimates the variance bound, within a 10 −4 relative accuracy.

B Permutation invariance
In this section, we discuss the implementation of symmetries in order to reduce the computational complexity of our SDP algorithm.For the sake of concreteness, we consider a system of N qubits.A generic density matrix describing the system is fully specified by exponentially-many parameters.We will show that if the observables and the generator are invariant under permutations of the qubits (PI), it is sufficient to optimize over that is: polynomially-many) real parameters representing a state which is block-diagonal with respect to the symmetry sectors.We hence obtain a drastic scaling reduction.
We first define the local basis.The state of each spin i ∈ {1, 2, .., N } is fully described by the Pauli vector σ(i) = (σ , where σ(i) a applies Pauli matrix a in party i with trivial action in the remaining qubits.From these operators, typical PI observables that we can form are collective spin Ĵ = N i=1 S (i) (elementary), where S = σ/2, and moments, i.e., polynomials thereof.All those operators have in common that they commute with a unitary representation of any permutation of N parties, π ∈ S N , Ûπ .Such unitary is described by its action in the computational basis Sufficiency of optimizing over PI states.-Inthis paragraph we show that it is sufficient to optimize over PI states if both the observables in D and the generator Ĝ are PI.The proof consists in two steps.First, we notice that if there is a state compatible with the data, then there always exists also a PI state compatible with the data.Indeed, for any state ρ we can generate the PI state by applying the averaging map ρPI = (N !) −1 π∈S N Û † π ρ Ûπ .Both states ρ and ρPI will yield the same expectation values for PI observables.The second step is to establish a relation between the QFI of ρ and that of ρPI .Since ÛPI commute with the generator Ĝ, for any π ∈ ŜN .Finally, by the convexity of QFI from the averaging map, we see that concluding that in our optimization for the search of the minimum QFI, PI states are sufficient.
Block-diagonal representation of PI states.-Here,we derive an efficient representation of PI states.First, one can check that any elementary PI operator ÂPI commutes with the total spin Ĵ2 .Hence, it can be decomposed in the block-diagonal form, ÂPI = ⊕ Jmax J=J min Â(J) , where J label the degenerate eigenspaces of Ĵ2 with eigenvalue J(J + 1), where J min = 0, 1/2, for even, odd N and J max = N/2.All sectors J have coarse-grained dimension 2J + 1, which is the number of PI inequivalent ways to combine N spins-1/2 to form a collective state of total spin J. Now, each vector represents a set of highly degenerated orthogonal states that cannot be distinguished by any PI observable.With with the data with minimal QFI.
Loong Len et al.-In Ref. [28] a method was proposed to lower bound the classical Fisher information (FI), F , of an observable Ô from moments of it.Up to K-th moments, they defined: where the vector b has components b p = ∂ θ ⟨ Ôp ⟩ = ⟨i[ Ĝ, Ôp ]⟩, and the matrix A, A pq = ⟨ Ôp+q ⟩ for p, q ∈ {0, 1, .., K}.The sensitivity Ξ −2 O is recovered in the first order (K = 1).By construction, F (K) ≤ F (K+1) , and in the limit K → ∞, it converges to F .The FI is a lower bound of the QFI.Thus, F (K) constitutes also a proper lower bound of QFI.Such values can be compared in our approach with the minimal QFI compatible with the data ⟨ D⟩ = {⟨i[ Ĝ, Ôp ]⟩, ⟨ Ôp+q ⟩} p,q∈{0,1,..,K} .
Gessner et al.-In Ref. [18] they ask which generator Ĝ ∈ Span( Ĝ) and observable Ô ∈ Span( Ĝ ⊕ Ô) capture the maximal sensitivity Ξ −2 O , which reads: where P Ĝ projects to the subspace generated by Ĝ and C = ⟨ Ô ÔT ⟩ − ⟨ Ô⟩⟨ Ô⟩ T is the complex covariance matrix.The optimal generator is G = n G • Ĝ where n G is the maximal eigenvector of M and the observable is This approach is used in the main text to find the optimal generator for a set of expectation values during OAT dynamics.For example, if spin and its fluctuations are given, the method recovers the squeezing and anti-squeezing axes for n O and n G respectively.By taking Ô up to second moments, we can derive a nonlinear squeezing parameter beyond spin.Such bound is also compared with our SDP value in the main text with the same data, i.e. considering data ⟨ D⟩ up to the fourth moments of the collective spin.In general, this approach yields a valid but not optimal lower bound to the QFI given the data.
Tóth et al.-In Ref. [45] an SDP-based approach was proposed to find lower bounds on quantities defined as convex roofs of polynomials of operator expectation values.It can also find a lower bound of the QFI given few arbitrary expectation values.The starting point is to notice that, as proven in Ref.
[51], the QFI is the convex roof of the variance: Then, the variance of a pure state |ψ⟩ ∈ H can be linearized in the two-copy space, (∆ Ĝ) In this embedding, Eq. ( 32) is rewritten as a minimization of separable symmetric states σ ∈ SEP(H ⊗ H).The separability condition in general cannot be incorporated as a positive-semidefinite constraint.However, one can relax the previous restriction to states with positive partial trace (PPT) and formalize the following SDP: where Ŝ is the exchange operator (swapping the two subsystems) with action Ŝ |ψ⟩ ⊗ |ϕ⟩ = |ϕ⟩ ⊗ |ψ⟩, t the partial transpose and tr the partial trace with respect to a subsystem.Note that this method requires the optimization over a matrix of size D 2 (with D the dimension of the Hilbert space E), while in our proposed algorithm, the largest matrices are of size 2D.This fact compromises the scalability of the method.For example, with D = 11 and the data used in Figure 6, we verify that the present algorithm is ∼ 1300 times slower than our proposed solutions.As for generic Hilbert spaces H, SEP ⊊ PPT, the solution of the problem Eq. (33) will generally underestimate the physical bound.However, for most practical purposes, e.g. the setting used in Figure 6, the bound obtained here is compatible with the bound given by our methods.Furthermore, the approach explained here is tailored to linear encoding.Its extension to arbitrary maps is not direct as the algorithm is based on the convex-roof construction Eq.(32) which was only proven valid for unitary transformations (see in particular the last remark in Ref. [51]).Note that using the recent results of Ref. [44], relevant for calculating the so-called quantum Wasserstein distance of two states, one can actually remove the restriction to symmetric states in the optimization of Eq. (33).
Fadel et.al.-In Ref. [16], the authors report a method to bound the fidelity with respect to a target state σ from partial information.As the quantum fidelity is lower-bounded by the Hilbert-Schmidt inner product, F(ρ, σ) ≥ Tr(ρσ), the problem can be linearized to an SDP: In the N -partite scenario, when the data D are symmetrized few-body correlations (or equivalently, low-order moments of collective observables), the problem can be solved efficiently in N [2].Such low-order moments are sufficient to certify spin-squeezed states generated via one-axis twisting (OAT) dynamics with almost unit fidelity.Within this regime, the QFI of the state would coincide with our bound based on the same moments (see e.g. Figure 3 for t ≲ 0.15π).Finally, they also suggest to use the said approach to bound the QFI with respect to transformations generated by collective operators, Ĝ , where ζ = 8 in general or ζ = 6 if one of the states is pure [6].

C.2 Exact methods
Here, we describe the approaches leading to a tight bound, i.e., exactly reproducible by a quantum state.
Suggestion by Kołodyński.-Another possibility is to use the definition of QFI from SLD, Ŝ.
From Eq. (35) we realize that the optimization over the state ρ for a given Ŝ is an SDP.Conversely, for a fixed state ρ, finding Ŝ is a linear problem.Consequently, one could alternate between these two problems in a see-saw strategy to find the optimal global solution.In our work, we incorporate implicitly both problems as a single SDP.For completeness we shall introduce yet another way to compute the QFI which was derived in Ref. [31] from the error propagation formula (c.f.Eq. (1) main text).
The RHS of Eq. ( 36) is thus a valid lower bound of the QFI given the data D = { X2 , i[ Ĝ, X]}.Such approach has been exploited to compute, also with a see-saw strategy, the so-called QFI of a channel Λ, defined as the optimization problem max ρ F Q [Λ(ρ), Ĝ] over all proper density matrices ρ.
Apellaniz et al.-We compare our approach with the method introduced in Ref. [5].The solution of the central problem of the main text, Eq. (2), is there recovered with the Legendre transform formalism as a maximization over lower bounds, inf µ W µ (r) := β(r), where with λ min denoting the smallest eigenvalue.For any value of r, β(r) will lower bound the minimum QFI with respect to Ĝ over states ρ compatible with Tr( Dρ) = ⟨ D⟩, namely, the solution of the central problem Eq. (2) (main text).In other words, r parameterizes a family of valid lower bounds, the stronger of which will be the max r β(r) = QFI Ĝ.As verified in the same Ref.[5] one can prove that by construction β is tightly coinciding with the solution of the problem.While the maximization over r is convex, this is not the case for µ, which turns the problem nonconvex.Now, for each update of r in the maximization algorithm, one must search the global minimum on µ over the large spectral range on Ĝ, µ ∈ [λ min ( Ĝ), λ max ( Ĝ)].In Figure 7, we show for the same example described in the Appendix A the landscape which is needed to explore to find the minimum QFI.We check that the minimum, min µ W µ (r * ) = 0.2401, is compatible with our previous value introduced in Appendix A. In green (dotted), the same for a random (not optimal) setting r = (0.89, 0.08) whose global minimum constitutes generally a lower bound of the QFI.
We verify several local minima appearing forcing to sweep over all the range of µ.This fact contrasts with our proposal, which just requires considering an infinitesimal parameter δθ.Once more, such method is based on the QFI as a convex roof Eq. (32), and it is especially suited to evaluate the metrological usefulness in phase estimation tasks such as interferometry.Its generalization to arbitrary parameter-encoding protocols is not addressed in Ref. [5].

D Bound for Dicke states of arbitrary spin
In this appendix, we prove the observation in the main text regarding the QFI of the manifold of states with ⟨ Ĵ2 x ⟩ = 0. First, we notice that such condition implies that the state is an eigenstate of Ĵx with eigenvalue 0. That is, ρ = where Π(J) Jx=0 are projectors properly orthonormalized projecting onto subspaces generated by states with well-defined zero spin projection along x, and total spin J.

E Scaling of the QFI bound with N close to Dicke states
As detailed in the main text, the condition ⟨ Ĵ2 x ⟩ = 0 is quite stringent.Here, we use our approach to find the minimal QFI compatible with some finite ⟨ Ĵ2 x ⟩ and ⟨ Ĵ2 ⟩.In particular, we fix the total spin to its maximal value, ⟨ Ĵ2 ⟩ = J max (J max + 1) with J max = N/2, so that we can restrict the optimization to the totally symmetric subspace and reach easily N = 70 on a standard computer.This enables us to study the scaling of the QFI with N .For ⟨ Ĵ2 x ⟩ = 0, we verified the Heisenberg scaling as the QFI is In Figure 8, we apply a small perturbation ⟨ Ĵ2 x ⟩ = aN/4 with a ≪ 1 and study how the QFI bound is modified.

F Scaling of the QFI bound with N with respect to squeezing bounds
In this appendix, we analyze the scaling of the advantage offered by our method with respect to the analytical lower bound given by the squeezing parameter ensuring the OAT dynamics.In the main text, for N = 10 particles, we show a gap existing between such two bounds when up to second moments and fourth moments are considered.Here, we investigate its robustness as the number of parties is increased.Specifically, in Figure 9 we display the gap per spin between the two bounds for both cases for N up to 30.In Figure 9a, we observe that at intermediate times, the gap is proportional to the particle number, leading to a robust improvement with respect to the usual spin squeezing parameter.At times t ≳ 0.15π, when the state becomes unpolarized, spin squeezing vanishes while our method certifies metrological usefulness of QFI/N ≈ 0.32, which is finite but smaller than the coherent-state limit QFI/N = 1.
The gap when up to fourth-order moments are considered with respect to the quadratic spin squeezing parameter, introduced in Ref. [18], is shown in Figure 9b.Unlike the linear case, we observe at intermediate times the gap saturating to zero as N is increased, signaling no metrological advantage with respect to the squeezing bound, except in the vicinity of the MHCS at t = π/3, π/4.It allows us to assess the robustness of the MHCS detection as reported in the main text.On the one hand, the usefulness of the 3-MHCS is already lost for N = 14.On the other, the detection of the 4-MHCS is stronger and does not drop below N (coherent limit) until N = 24.For larger system sizes, higher moments become necessary to detect the metrological usefulness of such states.

G Derivation of the data-driven witness Ŵ
As noticed in the main text, using up to fourth-order moments in the 4-MHCS for N ≤ 10 (even), the very same state is recovered with unit fidelity as the output of our SDP approach.This suggests that there is only a single state, namely the 4-MHCS, compatible with the data.
Here, in order to verify this fact, we outline a generic procedure to construct a witness operator Ŵ which is based on a certain data D (e.g. up to fourth order moments) and whose unique ground state corresponds to a given state |Ψ⟩ (e.g. the 4-MHCS).Without loss of generality, we can set λ min ( Ŵ ) = 0 as the identity operator is in the data list.Then, we need to find an operator Ŵ ∈ Span( D) such that it features |Ψ⟩ as its unique minimal eigenstate.Such problem can be cast as an SDP: where ker(|Ψ⟩⟨Ψ|) projects to the space orthogonal to |Ψ⟩⟨Ψ|.Intuitively, ∆ is the gap between the ground energy λ min ( Ŵ ) and the first excitation (with multiplicity), which is to be maximized.If the optimal ∆ is zero, it means that the ground state is degenerate, and a witness Ŵ cannot be found.
Conversely, ∆ > 0 implies that we succeed in finding a witness Ŵ based on D such that |Ψ⟩ is its unique extremal operator.In this case, we add the constraint ∆ ≤ cutoff to avoid unbounded solutions.With the above method, we find the existence of such witness for N ≤ 10 (even).Note that for N ≤ 4, up to fourth moments data are tomographically complete, so a witness could be simply Ŵ = I − |4MHCS⟩⟨4MHCS|.Here, we check that unexpectedly for N = 6, 8, 10 fourth moments are also sufficient to determine uniquely the 4-MHCS.

H Spin chain
In this last appendix, we go beyond permutationally-invariant states such as those generated in the OAT dynamics and turn to spatially-structured models.In particular, we consider the nearest-neighbour version of the OAT Hamiltonian studied in the main text, namely the Ising model: where the N spins form a one-dimensional chain with open boundary conditions.Because of the reduced interaction range, the evolution from the initial CSS along x (|π/2, 0⟩ = |→⟩ ⊗N ) will produce squeezing on a longer timescale than in the OAT model.The spatial structure of the model motivates the use of spatially-structured observables, beyond the collective observables considered in the rest of the paper.Hence, we use structure factor defined as: where Here, we will consider data of second moments {C (k) ab } k≤K only up to momentum K ≤ N − 1 as well as the mean spin ⟨ Ĵ⟩.There is no squeezing parameter based on such data as it is not algebraically closed.Indeed, two momenta k, k ′ couple to the sum k + k ′ mod N .Nevertheless, our algorithm is flexible in the data and we are able to verify that nonzero-momentum observables can improve the metrological bound.In Figure 10 we show the results for N = 4 taking as generator Ĵas .We observe that the inclusion of momentum K = 1 indeed improves over the QFI lower bound.This illustrates the potentiality of our approach to certify the metrological usefulness in generic multipartite systems, using spatially-structured data.Note however that in this general case, our algorithm requires computing memory resources scaling exponentially with N .

Figure 2 :
Figure2: Evolution of the metrological bounds from second moments of the collective spin during OAT dynamics for N = 10 particles.From bottom to top: in blue (dashed) conventional spin squeezing parameter and in orange (solid line), the SDP bound based on the same data.We also depict horizontal lines at QFI = 0 and QFI = N to represent no metrological usefulness and the coherent state limit respectively.

Figure 4 :
Figure4: Evolution of the metrological bounds from second moments of spin and parity measurements during OAT dynamics for N = 10 particles.From bottom to top: in blue (dashed) non-linear squeezing parameter derived from Ref.[18] using Ô = { Ĵx , Ĵy , Ĵz , Πx , Πy , Πz } as measurement operators; orange curve (solid line)the corresponding SDP bound using the same data.In green (dotted) parity squeezing parameter of Eq. (10) with a = y (generator direction) and b = x (measurement direction).In red (dashed-dotted) is plotted the maximal QFI of the OAT-generated states with respect to collective spin rotations.

Figure 5 :
Figure 5: Fidelity (left) and QFI (right) bounds as a function of δθ.In blue (solid line), the value computed from the fidelity F max (SDP output).The orange dashed line corresponds to the final inferred QFI, F Q [ρ * δθ=0.02 , Ĝ] = 0.2401.The green dotted curve is deduced from the QFI of the optimal state F Q [ρ * δθ , Ĝ] for each δθ using Eq.(11).

Figure 6 :
Figure 6: a) Computing time for both approaches for different values of N and δθ = 0.01 for the SDP used in the main text.b) Relative difference of the two bounds (F stands for fidelity and V for variance) for the same values of N .

Figure 7 :
Figure 7: Blue solid curve W µ (r) of Eq. (37), for the optimized r maximizing its global minimum (orange dashed line), r * = arg max r β(r).For the present example, we have r * = (1.31,−0.95).We check that the minimum, min µ W µ (r * ) = 0.2401, is compatible with our previous value introduced in Appendix A. In green (dotted), the same for a random (not optimal) setting r = (0.89, 0.08) whose global minimum constitutes generally a lower bound of the QFI.

Figure 8 :
Figure 8: a) Minimum QFI compatible with ⟨ Ĵ2 ⟩ = (N/2)(N/2 + 1) and ⟨ Ĵ2 x ⟩ = aN/4 for different values of a and number of particles N .In panel b) the same plot for nonzero values of a with the horizontal axis rescaled to ⟨ Ĵ2x ⟩.The dashed vertical line indicates ⟨ Ĵ2x ⟩ = 0.25.

Figure 8
Figure 8 suggests that the Heisenberg scaling is valid only when ⟨J 2x ⟩ ≲ 1/4, where the QFI bound shows a maximum.Beyond this point, the Heisenberg scaling is lost and the metrological usefulness is reduced.

Figure 9 :
Figure 9: Gap between our SDP-based bound and the analytical squeezing lower bound with the same data is considered: a) up to second moments where we depict also the QFI = 0.32N horizontal line, b) up to fourth moments, during the OAT dynamics and the line QFI = 0.The squeezing lower bound for second moments corresponds to the conventional spin squeezing parameter.If, in addition, up to fourth moments are taken into account, such bound corresponds to the generalized spin squeezing parameter introduced in Ref. [18].

⟩ − N δ ab / 4 .
a, b ∈ {x, y, z}.The structure factor captures the spin fluctuations at momentum k ∈ {0, 1, .., N − 1}.Indeed, it can be expressed as the second moment of the collective spin at momentum k, Such operators generalize to arbitrary momentum k the uniform collective spin used throughout this work, Ĵ = Ĵ(k=0) .