Symmetry-resolved entanglement detection using partial transpose moments

We propose an ordered set of experimentally accessible conditions for detecting entanglement in mixed states. The k-th condition involves comparing moments of the partially transposed density operator up to order k. Remarkably, the union of all moment inequalities reproduces the Peres-Horodecki criterion for detecting entanglement. Our empirical studies highlight that the first four conditions already detect mixed state entanglement reliably in a variety of quantum architectures. Exploiting symmetries can help to further improve their detection capabilities. We also show how to estimate moment inequalities based on local random measurements of single state copies (classical shadows) and derive statistically sound confidence intervals as a function of the number of performed measurements. Our analysis includes the experimentally relevant situation of drifting sources, i.e. non-identical, but independent, state copies.


INTRODUCTION
In the past years, a considerable effort led to the building of larger and larger Noisy Intermediate-Scale Quantum (NISQ) devices [1][2][3] . For the benchmarking of such devices comes the need for more scalable tools in order to characterize the underlying many-body quantum state (see e.g. 4 and references therein). For instance, characterizing the entanglement properties of these quantum states is, besides the intrinsic theoretical interest, essential to gauge the performance and verify the proper working of the NISQ devices.
As a first prominent example among the tools to characterize entanglement, there is the concept of entanglement witness 5,6 . An entanglement witness is a functional of the quantum density matrix that separates a specific entangled state from the set of all separable states (for examples of entanglement witnesses in various types of systems, see e.g. ref. 7 and references therein). When this function is linear, it can be identified with an observable whose expectation value can be used to decide whether the target state is entangled or not. By contrast, in this work, we shall focus on a superset of the set of separable states: the set of states with positive partial transpose. In other words, we will focus on sufficient conditions for entanglement (equivalently, necessary conditions for separability).
From the numerous theoretical sufficient conditions for entanglement that have been developed in the literature, many cannot be straightforwardly implemented experimentally, mainly because they require the (exponentially expensive) knowledge of the full density matrix [8][9][10] . This is for instance the case of the celebrated PPT condition 5,11 , which states that a separable state ρ always has a positive semi-definite (PSD) partial transpose (PT) ρ Γ for any bipartite splitting of its subsystems. Thus, if ρ Γ has (at least) a single negative eigenvalue, then ρ is entangled. The negativity, which resulted from this condition, is a highly used entanglement measure for mixed states [12][13][14] .
This powerful entanglement condition, which found many applications in theoretical works [15][16][17][18][19][20][21][22] , is difficult to apply in experimental conditions as the PT spectrum is difficult to access. To overcome this challenge, it was shown in ref. 23 that valuable information about the PT spectrum can be obtained from a few PT moments trðρ Γ Þ k only. Using the first three PT moments, an entanglement condition, called p 3 -PPT, was proposed and shown to be useful for detecting entangled states in several different contexts. Moments trðρ Γ Þ k have the advantage that they can be estimated using shadow tomography 23 in a more efficient way than if one had to reconstruct ρ via full quantum state tomography. Indeed, the k-th order PT moment of a state ρ can be expressed as the expectation value of some k-copy permutation operator O 23 , i.e. trðρ Γ Þ k ¼ tr Oρ k À Á . As shown in ref. 24 , the classical shadow formalism allows to estimate such an expectation value from independent classical snapshotsρ 1 ; ;ρ k (which can each be obtained from post-processing single-qubit measurement outcomes) through the U-statistics estimatorô k ¼ tr Oρ 1 Á Á Á ρ k ð Þ . Therefore, as in other randomized measurements protocols probing entanglement [23][24][25][26][27][28][29][30][31][32][33] , the classical shadows formalism only requires (randomized) single-qubit measurements in experiments realizing the single-copy state ρ. In this paper, we follow this idea of using PT moments to build experimentally computable entanglement conditions, and extend the p 3 -PPT condition in two directions.
On the one hand, we propose different entanglement detection strategies depending on how many PT moments can be estimated. Starting from the third-order moment, we show that the estimation of each higher-order moment gives access to an independent entanglement condition. Interestingly, if all the PT moments can be estimated, this set of conditions is then necessary and sufficient for the state to be PPT (i.e. to have a positive semi-definite partial transpose). Of course, the higher the moment, the larger the number of experimental runs needed. In case of higher-order moments cannot be accessed, we show how to obtain an optimal entanglement condition using PT moments of order up to three.
On the other hand, we investigate the effect of symmetries on this entanglement detection method. As shown in ref. 34 for the case of dynamical purification, taking symmetries into account to define symmetry-resolved (SR) versions of the tools usually used to characterize quantum states can enable a finer characterization of some quantum features and even reveal phenomena that cannot be observed without symmetry-resolution. For states preserving an extensive quantity, we define SR versions of the PT-moment inequalities mentioned previously and show that these are indeed better suited to detect the entanglement of such states. Furthermore, we also show that these SR inequalities provide a sufficient entanglement condition for states that do not possess any symmetry.
The conditions derived here are particularly interesting from an experimental (and numerical) point of view, as low moments of (partially transposed) density operators are accessible. We show how source drifts in an experiment can be taken into account and how the quantities which are of interest here can be accurately estimated via local measurements on single copies of the state.

Definitions and summary of results
In this section, we introduce the basic definitions needed for the entanglement detection criteria below, and summarize in a succinct manner our main results. Given a bipartite state ρ = ρ AB , we denote by ρ Γ its partial transpose with respect to subsystem B. We say that ρ is PPT if ρ Γ is positive semi-definite, and NPT otherwise. All NPT states are entangled, however there are entangled states, known as bound-entangled states, that are not NPT. We focus here on the detection of NPT entangled states.
We denote the k-th order moment of a matrix M by p k ðMÞ tr M k : We will mostly consider moments of ρ Γ , and sometimes use the short-hand notation p k ≡ p k (ρ Γ ). In the presence of symmetries, the partial transpose can be cast in block diagonal form: we denote as ρ Γ ðqÞ the resulting blocks, where q indicates a quantum number, and define the corresponding moments p k ðρ Γ ðqÞ Þ as from Eq. (1). We start by recalling the p 3 -PPT condition of ref. 23 , i.e. that any PPT state satisfies Any state violating this condition is NPT and therefore entangled.
The p 3 -PPT condition will serve as a reference point below in accessing the predictive power of the proposed relations. Throughout this paper, we will establish and analyse different sets of necessary (and sometimes also sufficient) PPT conditions based on PT moments. An exhaustive list summarizing these conditions is given below, together with two important results regarding the experimental estimation of the PPT conditions using classical shadows.
i) the first set of conditions, which we dub D n conditions, also contains polynomial inequalities in the moments p k of order up to k ⩽ n. The first non-trivial such a condition is D 3 , and reads: Knowing only the first three moments p 1 (ρ Γ ), p 2 (ρ Γ ) and p 3 (ρ Γ ), this condition is optimal for detecting entanglement if 1/2 ⩽ p 2 (ρ Γ ) ⩽ 1.
Knowing moments of order up to the dimension of ρ Γ , the set of D n conditions provides a necessary and sufficient condition for NPT entanglement; ii) the second set of conditions, dubbed Stieltjes n , involves inequalities among the moments p k of order up to n. The condition Stieltjes 3 is equivalent to p 3 -PPT, while Stieltjes 5 reads: and similar conditions are obtained including higher moments; iii) in case high-order moments are difficult or too expensive to access, we also show how to obtain an optimized, necessary condition for PPT using only PT moments of order up to three. We call this condition D opt 3 ; iv) all of the above conditions can be cast in terms of ρ Γ ðqÞ , in which case we add the prefix SR (for symmetry-resolved). For instance, the SR-p 3 -PPT condition for sector q reads Since these conditions are sensitive to the presence of negative eigenvalues in a specific symmetry sector, they are typically much more sensitive than their non-SR counterparts, as illustrated in Fig. 1.
In the SR case, it is worth mentioning that also the SR-D 2 condition, is non-trivial; v) we show how SR conditions can, in fact, be applied to arbitrary states, via the application of a proper transformation on the density matrix of interest. In practice, this transformation is Fig. 1 An illustration of the proposed method for entanglement detection. We assume the experimentally relevant situation of a source producing non-identical but independent copies {ρ 1 , …, ρ N } ("drift''). Randomly chosen unitaries U i are applied to the qubits of each copy and then measured in the standard basis. Using classical shadows 24 , these measurement outcomes are post-processed to obtain the moments p j ¼ trðρ Γ avg Þ j . As explained in the main text, we combine those moments to derive inequalities whose violation implies that the state ρ avg is NPT, thus showing that at least one of the states ρ k produced by the source is entangled. We also show how symmetry-resolution techniques can be used to enhance entanglement detection capabilities.
A. Neven et al. effectively carried out in the post-processing step of the classical shadows; vi) as illustrated in Fig. 1, the uncertainty in estimating the moments can be bounded, in principle, using the classical shadows formalism. Here, we show how to combine those bounds to provide rigorous confidence intervals for SR-D 2 , which considerably strengthens the impact of our results in real experiments.
Entanglement detection from partial transpose moments In this section, we present entanglement conditions based on PT moments. We extend the idea behind the p 3 -PPT condition (c.f. Eq. (2)) of ref. 23 in two directions. On the one hand, we present a set of inequalities involving all the PT moments which provides a necessary and sufficient condition for the underlying state to be PPT. In addition, each condition of this set is itself a necessary PPT condition. On the other hand, we show how to optimize such entanglement conditions when only a few low-order PT moments are accessible.
The idea behind this set of conditions is to use Descartes' rule of signs on the characteristic polynomial of a Hermitian matrix to obtain a set of moment inequalities that has to be satisfied by any PSD matrix. Applied to the partially transposed matrix ρ Γ , such conditions can then be used to detect the entanglement of NPT quantum states. More precisely, using the definition of the elementary symmetric polynomials on d variables, for i = 1, …, d, and e 0 (x 1 , …, x d ) = 1, we derive in the Methods section the following lemma.  7)).
Using Newton's identities, which relate the elementary symmetric polynomials, e k , in the eigenvalues of A to the moments of A through the recursive formula ðÀ1Þ iÀ1 e kÀi p i ðAÞ; (8) each inequality e i ⩾ 0 can be transformed into an inequality involving moments of A of order up to i. We denote by D i these moments inequalities. As an illustration, the conditions D 1 to D 4 read respectively. One has p 1 (ρ Γ ) = 1 for any quantum state ρ, implying that D 1 is trivially satisfied. Similarly, since p 2 (ρ Γ ) is equal to p 2 (ρ) (i.e., to the purity of ρ) for any quantum state ρ, the inequality D 2 is also trivially satisfied. Therefore, when ρ is a quantum state, the first non-trivial inequality for ρ Γ is D 3 . As will be shown in the next section, it is sometimes more efficient (in order to detect entanglement) to apply these inequalities to projections of ρ Γ onto specific subspaces, rather than to ρ Γ itself. We would like to stress here that, in that case, the argument above does not hold, so that the inequality D 2 is not trivially satisfied and can already reveal the presence of entanglement. When applied to ρ Γ , Lemma 1 and Newton's identities (8) can thus be used to detect NPT entangled states from PT moments only. From an experimental point of view, this is an important aspect of this entanglement detection scheme, as PT moments are experimentally more affordable to estimate than, for instance, the whole spectrum of ρ Γ . As PT moments are more expensive to be estimated the higher the order, these inequalities should be considered starting from those involving the lowest moment orders. Even though showing that a state is NPT with this method can in principle require the knowledge of all the PT moments, we will provide many experimentally relevant instances where entanglement can be effectively detected from low-order moments even in the presence of errors. To this end, we provide confidence intervals for the quantities of interest granting that a certain inequality is violated with high probability (see Theorem 1 and the Supplementary Information).
Similarly, let us mention here that necessary and sufficient conditions for a matrix to be PSD can be expressed as different sets of polynomial inequalities in its moments. One of such sets can be deduced from the well-known (truncated) Stieltjes moment problem (see Methods). In the part dedicated to applications, we illustrate the usefulness of these inequalities by applying them to the entanglement detection of the ground state of the XXZ model (c.f. Fig. 5). Let us finally also mention that, from a few moments of a Hermitian matrix, one can also obtain bounds on the distance between this matrix and the PSD cone 35 .
Optimized condition for low-order moments. Due to (experimental) constraints, it might not be possible to determine all, but only a few PT moments. This is why, we show here how to optimize necessary PPT conditions using only PT moments of order up to three. From the previous sections, we already have two examples of such conditions, namely the p 3 -PPT and D 3 conditions. As illustrated in Fig. 2, the p 3 -PPT (D 3 ) condition is tighter than D 3 (p 3 - Fig. 2 Entanglement conditions of order three. The plot of the value of the third moment p 3 saturating the p 3 -PPT (dashed orange curve), the D 3 (dashed green line), and the optimal D opt 3 (thick black curve) conditions as a function of the second moment p 2 for a normalized Hermitian matrix. According to the p 3 -PPT condition, any state ρ with a value of p 3 (ρ Γ ) below the dashed orange curve is entangled. Similarly, the condition D 3 shows that any state ρ with a value of p 3 (ρ Γ ) below the dashed green line is entangled. From this plot, it is clear that, for p 2 (ρ Γ ) > 1/2, all entangled states detected by the p 3 -PPT condition are also detected by D 3 , which coincides with D opt 3 in this case. When p 2 (ρ Γ ) < 1/2, the p 3 -PPT condition is then stronger than D 3 , and D opt 3 represents a slight improvement over the p 3 -PPT condition. As illustrated in the applications, this slight improvement can nevertheless be important for the detection of physically relevant states.
PPT) for states with purity larger (smaller) than 1/2, respectively. As the low-order moments are easier to access experimentally, we now address the question about the optimal inequality involving PT moments of order up to three.
To answer this question, we use the following approach. For fixed values p 1 and p 2 of the first two moments, we determine the minimal value p min 3 that the third moment can reach for any PSD matrix. From this bound, we know that any Hermitian matrix with a smaller third moment is necessarily not PSD. Note that we want here to minimize p 3 because it is an odd moment (for which negative eigenvalues would have the tendency to decrease the value of the moment). For an even moment, we would instead maximize the value of this moment over PSD matrices. This is also reflected in the D n conditions (9)- (12), where the inequality sign alternates between even and odd values of n. Naturally, we restrict ourselves to values of p 1 and p 2 which are compatible with a PSD matrix, and therefore satisfy Eqs. (9) and (10). Recall that for the partial transpose of a density operator, this is always fulfilled.
Given a d × d PSD matrix A, with non-zero eigenvalues λ 1 , …, λ r , for some r ∈ [1, d], this optimization can be performed with the help of Lagrange multipliers. As shown in the Methods section, this leads to solutions with only two distinct eigenvalues λ a , λ b with multiplicity r a , r − r a , respectively, for r a ∈ [1, r]. Assuming, without loss of generality, that λ a ⩾ λ b , the optimization of p 3 leads then to r a = r − 1. For each value of r, the optimal value of p 3 can be easily determined in the interval [1/r, 1/(r − 1)]. For r = 2 this leads to D 3 whereas for r > 2 one obtains an optimal value of p 3 which is slightly better than p 3 -PPT. Observe that p min 3 ðp 2 Þ is a piece-wise function and the derivative ∂p min 3 =∂p 2 is discontinuous at points p 2 = 1/r (see Fig. 2).
In case higher-order PT moments can be accessed experimentally, one may consider optimizing higher-order PPT conditions. In principle, the Lagrange multipliers method we used to optimize the entanglement condition using PT moments of order at most three can be straightforwardly extended to optimize higher-order PPT conditions. In practice, however, the number of variables would grow with the moment order, making the analytical resolution of the optimization problem more involved. In the simultaneous but independent work 36 , the authors also consider this optimization problem and explicitly carry out the optimization of the PPT condition using PT moments of order up to four (though with a slightly different method). In addition, they provide a numerical code to perform the optimization with PT moments of order up to five.

Symmetry-resolved entanglement detection
Symmetries, as they often occur in physical situations, can be exploited to observe relevant phenomena (see e.g. refs. 34,[37][38][39][40][41][42][43][44][45][46][47] ). Here, we use symmetries to ease the detection of entanglement. More precisely, we apply the previously developed tools to symmetric states, which will lead to conditions of entanglement involving much lower moments of the partial transpose projected onto certain subspaces. Despite the fact that these quantities differ significantly from the moments of ρ Γ , we will show later on that they can nevertheless be estimated using the framework of classical shadows.
We consider a bipartite state ρ = ρ AB , with subsystems A and B containing n and m qubits, respectively. We assume that this state commutes with P nþm i¼1 Z i , or similarly with the total number operator N ¼ N A þ N B . Here and in the following, we denote by X, Y, Z the Pauli operators. Obviously, such a state has a block diagonal form, i.e., where each block (or sector) is labeled by an eigenvalue q ∈ {0, 1, …, n + m} of the operator N and has support in the corresponding eigenspace. Here, with and similarly for B. It has been shown 48 that, for this type of symmetry, the partial transpose ρ Γ is also block diagonal, but in a different basis. In fact, This can be easily seen as follows. Consider a matrix element This shows that matrix elements within a block of ρ are mapped, via partial transposition, to matrix elements within a block of ρ Γ .
The size of the sector corresponding to the eigenvalue q in the block-decomposition of ρ Γ is given by When the partial transpose of a density matrix has a block structure, it is naturally PSD iff each block is itself a PSD matrix. Therefore, one can apply the conditions of the previous section directly to the blocks ρ Γ ðqÞ of the partial transpose. For the p 3 -PPT condition, the corresponding symmetry-resolved (SR) inequalities are simply for all q = − m, − m + 1, …, n. Any violation of a PSD condition in one of the blocks is then sufficient to show that ρ Γ has at least one negative eigenvalue and that ρ is therefore entangled. When using the D i conditions, symmetry-resolution can be a significant advantage (see e.g. the applications section). First, the necessary and sufficient PSD conditions involve moments of order at most equal to the dimension of the largest block, that is n þ m bðn þ mÞ=2c , which is necessarily smaller than the dimension of the density matrix itself. Second, since a block ρ Γ ðqÞ of ρ Γ is (in general) not the partial transpose of any density matrix (i.e. there could be no σ > 0 such that ρ Γ ðqÞ ¼ σ Γ ), the inequality: is not necessarily satisfied. This implies that moments of order two can already be sufficient to detect entanglement. As stressed in the introduction, using PT-moment inequalities to detect entanglement is particularly interesting from an experimental point of view, because such PT moments can be estimated, for instance using shadow tomography 23 . As we show in the following lemma, the shadow tomography protocol used in ref. 23 can also be used to estimate moments of blocks of the partial transpose (which differ significantly from the PT moments) by slightly modifying the non-linear observable that has to be measured.
A. Neven et al.
where the operators L ðkÞ i are given by Here, the sum over each a i (b i ) runs from 1 to 2 n (2 m ) respectively.
Proof. First, it is straightforward to see that Then, we have that Here, the first equality follows from the definition of L ðkÞ i and Eq. (22); the second, from tr ðRS Γ Þ ¼ tr ðR Γ SÞ for any two matrices R, S; and the third, from P Γ i ¼ P i . Finally, using that P i ¼ P 2 i are orthogonal projectors, the cyclic property of the trace, and the block structure of ρ Γ , we have which completes the proof.
Classical shadows. Classical shadows are a convenient formalism to reason about predicting properties of a quantum system based on randomized single-qubit measurements performed in singlecopy experiments 24 . We refer the reader to 24 or to the Supplementary Information for an introduction to classical shadows. The original classical shadow formalism is contingent on noiseless measurements and sources that produce i.i.d. states. Subsequently, it was shown that classical shadows can also handle noisy measurements 49,50 . As we show in detail in the Supplementary Information, the formalism can also be extended to take nonidentical, but independent, state preparations ("drifts") into account, i.e. the source produces the states {ρ 1 , ρ 2 , …, ρ N }. In this case, each snapshotρ i will have a different expectation value and i.e., the average of the snapshots converges to the average state. Since different snapshots are statistically independent, it turns out that one can estimate linear functions, say tr ðOρ avg Þ. These ideas regarding the prediction of linear observables do extend to higher-order polynomials. Here, we restrict ourselves to the quadratic case 24 , involving up to second-order moments of a block of the partial transpose. An extension to higher-order polynomials is conceptually straightforward, but can become somewhat tedious to analyze 23 . Let us fix a block label i and consider the second-order moment inequality restricted to this block: where, using Lemma 2, Recall that D ðiÞ 2 ðρÞ < 0 implies that ρ is entangled. Here, we will provide confidence intervals in the estimation of D ðiÞ 2 ðρ avg Þ given a fixed number of measurements. Observe that D ðiÞ 2 ðρ avg Þ < 0 implies that there is at least a value k ∈ {1, 2, …, N} such that D ðiÞ 2 ðρ k Þ < 0. Thus, D ðiÞ 2 ðρ avg Þ < 0 implies that the source is able to produce entangled states.
Let us finally introduce (see the Supplementary Information for details) the empirical average, over N(N − 1) pairs of independent snapshots, b D We can fix the desired approximation accuracy ϵ and a probabilityof-error threshold δ to obtain a lower bound on the measurement budget N. For simplicity, let us consider a non-trivial sector (i.e. q ≠ − m, n) and assume n + m ⩾ 4. Then one has the following theorem (see the Supplementary Information for a slightly better bound).
2 ðρ avg Þ ⩽ ϵ with prob: ðat leastÞ1 À δ: It should be noted that this bound corresponds to randomized measurements performed using local Clifford unitaries (see Supplementary Information). Using a different measurement strategy, for instance with global Clifford unitaries, we could expect a much better scaling of the measurement budget 24 , but such global unitary gates are more difficult to implement, and are still out of reach for most of the current experimental platforms.
Finally, let us stress that this error bound addresses the estimation of D ðiÞ 2 ðρ avg Þ in terms of a single U-statistics estimator. The poor scaling in 1/δ can be exponentially improved by dividing the classical shadow into equally-sized batches and performing a median-of-U-statistics estimation instead 24 : 1=δ ! const log ð1=δÞ. However, numerical experiments conducted in ref. 23 suggest that this trade-off is only worthwhile if one attempts to predict many properties with the same data set.
SR inequalities applied to states without symmetries. In the last part of this section, we show that the SR inequalities can also be used to detect the entanglement of arbitrary states, including those that do not have any symmetry.
The reason for that is that there exists a local channel C that transforms any state ρ into a state σ CðρÞ that has the desired block structure. The channel can be realized with local operations assisted by classical communication, and can thus not generate entanglement. Therefore, the initial state ρ must be at least as entangled as the final block diagonal state σ. This statement holds for any entanglement measure. As a consequence, if entanglement is detected in σ (which can be investigated using the symmetry-resolved tools), then ρ is necessarily also entangled. In other words, looking at the entanglement of σ, the "blockdiagonalized" version of ρ, provides a sufficient condition of entanglement for ρ. This condition is not necessary as it could be A. Neven et al. that the channel C destroys all the entanglement of ρ.
The local channel that can be used for this approach is the following: (30) where k ¼ blog ðn þ mÞc þ 1 and U i ¼ Z i=2 k . The fact that this channel maps ρ to a state σ that is block-diagonal in the numberof-excitations basis can easily be seen as follows. First, observe that for any j ∈ {0, …, 2 (n+m) − 1}, the computational basis state j j i is an eigenvector of U i , associated to an eigenvalue, ðÀ1Þ jjji=2 k , that essentially depends on |j|, the number of excitations of j j i. Therefore, we have iðjjjÀjj 0 jÞ 2 k ρ j;j 0 j j i j 0 h j: For any j and j 0 having different number of excitations, i.e. such that jjj ≠ jj 0 j, the sum over i in Eq. (31) vanishes, explaining why σ is diagonal in the number-of-excitations basis. As can be seen from the argument above, the non-zero elements of σ are all equal to the corresponding elements of ρ. This implies that the channel can effectively be replaced by a sum of projectors onto all the number-of-excitations sectors. From an experimental point of view, the practical implementation of this channel can thus be circumvented by using the observables of Lemma 2 in the post-processing of the classical shadows.

Applications
In this section, we apply and compare the entanglement conditions presented in the previous sections on various physical systems. For the systems possessing a symmetry as discussed above, we highlight some of the advantages that can result from considering symmetry-resolved entanglement detection tools.
Entanglement detection in quench dynamics. We begin by considering the situation of quench dynamics, where entanglement emerges from the dynamics of a many-body Hamiltonian. We consider the model presented in ref. 34 , where the interplay between coherent dynamics with U(1) symmetry and dissipation leads to a dynamics of "purification". Here, we will use the same formalism to show how entanglement is generated at short times, and can be detected via the symmetry-resolved versions of the D 2 and p 3 -PPT conditions. In the next subsection, we will consider an analogous experimental situation obtained with trapped ions 28 .
Our model is described by a master equation with the lowering and raising operators σ À j ¼ ðX j À iY j Þ=2, σ þ j ¼ ðX j þ iY j Þ=2, and the Hamiltonian and where γ is the spontaneous emission rate. Here, we consider open boundary conditions. The hopping between spins i and j is described by the coefficient J ij and we consider in this subsection nearest-neighbor hopping J ij = Jδ i+1,j with strength J. The initial state is the Néel state ρð0Þ ¼ ψð0Þ j i ψð0Þ h j, with ψð0Þ j i¼ #" j i N=2 . As shown in ref. 34 , the time evolved state ρ(t) of the N spin system has the block diagonal form of Eq. (13). Moreover, the partially transposed matrix w.r.t a partition A, ρ Γ is also block diagonal with blocks ρ Γ ðqÞ . Here, the index q represents the difference between the number of spin excitations in A and the one in the complement partition B.
As we are interested in short-time dynamics, we can solve Eq. (32) in first-order in perturbation theory, which is valid for t ≪ 1/ J, 1/γ. Considering for concreteness a half-partition A, made of the first N A sites, we obtain a block with a negative eigenvalue 34 . Assuming for simplicity N A = N/2, N A even, we obtain ρ Γ ðÀ1Þ ðtÞ ¼ γt The presence of a negative eigenvalue in this sector can be detected from the value of the moments in leading order in J ≫ γN A . In particular, the p 3 -PPT ratio and the D 2 condition can be used to reveal the presence of entanglement at short times. We show in Fig. 3 a numerical confirmation of these results for various values of γ/J and N = 8, which was obtained by simulating Eq. (32). We note that, in the present context, utilizing symmetryresolution is fundamental to detect entanglement: this is due to the fact that the negative eigenvalues in ρ Γ appear in sectors that are not macroscopically populated 34 , so that moments without symmetry resolution would not be able to detect them.

Experimental demonstration in a trapped-ion quantum simulator.
In the previous section, we showed in an idealized theoretical setting that entanglement is generated -and can be revealed via SR-entanglement conditions-at early times after a quantum quench. Here, we demonstrate this effect experimentally via the measurement of the SR-D 2 and SR-p 3 -ppT condition using randomized measurement data taken at early times after quantum Fig. 3 Symmetry resolved entanglement detection in quench dynamics with spin excitation loss. We study SR-entanglement in quench dynamics in a system consisting of N = 8 spins initialized in a Néel state #" j i N=2 and evolved with H XX subject to spin excitation loss with various rates γ (γ/J increases with the darkness of the color, see insets). We take A = [1, 2, 3, 4] and B = [5,6,7,8]. In panels (a) and (b), the D 2 ratio and p 3 -PPT ratio of sector q = − 1 are shown, respectively. Entanglement is detected for values below unity in the shaded gray areas. The insets in (a) and (b) show the early time value at t = 0 + of the D 2 -ratio a) and p 3 -PPT ratio (b), respectively, as a function of the decoherence rate γ/J. Black lines are the perturbation theory results displayed in Eqs. (40) and (39).
quench in a trapped-ion quantum simulator (c.f. ref. 28 ). In particular, we show that the SR-D 2 condition and SR-p 3 -PPT condition allow for a fine-grained detection of bipartite entanglement, in regimes where the corresponding global conditions 23 and conditions relying on the purities of different subsystems 28 are not conclusive.
In the experiment reported in ref. 28 , a one-dimensional spin-1/ 2-chain, consisting of N = 10 spins, was initialized in the Néel state "# j i 5 and time-evolved with the Hamiltonian H XY [Eq. (33)] where the coupling parameter J ij follows the approximate power-law decay J ij ≈ J 0 /|i − j| α , with α ≈ 1.24, J 0 = 420s −1 . The Hamiltonian evolution exhibits a global U(1)-symmetry conserving the total magnetization of the system (i.e., [H, ∑ i Z i ] = 0). Symmetry-breaking terms (such as σ þ i σ þ j þ h.c. ) are strongly suppressed due to a large effective magnetic field 28 . As detailed in refs. 28,34 weak decoherence effects are present in the experiment, including imperfect initial state preparation, local spin-flips and spontaneous emission during the dynamics, and measurement errors model as local depolarization. Note that coherent spin-flips do not preserve the global magnetization and block-diagonal form of the (reduced) density matrix. On the timescales accessed in the experiment, these effects are however very weak (causing in numerical simulations including the above decoherence model a purity mismatch of the order of 10 −5 of the full 10-spin density matrix ρ vs. the projected one ρ Q = ∑ q Q q ρQ q at t = 5 ms).
In ref. 28 randomized measurements were performed at various times (t = 0ms, …, 5ms) after the quantum quench. As described in detail in ref. 34 (see also the Supplementary Information), we can use this data to estimate SR-PT moments and the SR entanglement conditions via classical shadow formalism 24 . In Fig. 4, we present the SR D 2 and p 3 -PPT conditions in the different sectors, for a subsystem consisting of the neighbouring spins A, B = [4,5], [6,7] and where the partial transpose is taken in the subsystem A = [4,5]. Similar to the results of the previous subsection, both conditions detect entanglement at short times in the sector q = −1. The corresponding global conditions, in particular the global p 3 -PPT condition, do not reveal the presence of entanglement in this regime [see Fig. 4 b)].
The fact that the SR-D 2 condition can reveal the presence of entanglement is particularly interesting from an experimental point of view as it implies that entanglement can be detected from the estimation of only two moments of the partial transpose (in a sector). For the shadow estimation of D 2 ( − 1), our rigorous bound from Theorem 1 ensures that~1.3 × 10 6 measurements would be sufficient to guarantee entanglement detection with a probability of 95%. While this represents an upper bound, valid irrespective of the quantum state in question, for the specific states in the experiment only 8 × 10 5 have been performed. The errorbars of the experimental are then drawn at 1.96σ where the standard error of the mean σ has been estimated for each data point using Jackknife resampling. For normally distributed data with empirical mean μ, μ ± 1.96σ defines a 95% confidence interval. While normal distribution is here not guaranteed a priori, we checked through additional numerical simulations of many experiments (with fixed number of runs per experiment) that errorbars of length 1.96σ indeed approximate a confidence interval with confidence level 95%.
While the SR-D 2 condition requires only the estimation of first and second PT-moment, the third order SR-p 3 -PPT condition [panel b)], allows to detect entanglement in an even wider time window. In comparison to the global p 3 -PPT condition [red curve in panel b)], this clearly demonstrates the benefit of taking symmetry-resolution into account.
Entanglement detection in the ground state of the XXZ model. The XXZ spin chain is a generalization of the Heisenberg chain including an anisotropy in the interaction along the z direction, whose Hamiltonian reads: We will fix J = 1 as energy unit: J z sets the strength of the anisotropy along the z-axis. The phase diagram at zero temperature is known 51 : the system hosts an antiferromagnetic phase when J z < − 1, a Luttinger liquid for J z ∈ [ − 1, 1], and a ferromagnetic one for J z > 1. We might expect that the entanglement conditions we described in the previous sections will detect that the state is not PPT in the range J z 2 À 1; À1. Since the XXZ spin chain exhibits a U(1) symmetry related to magnetization conservation, we can exploit the symmetry-resolved counterpart of the D k conditions, the p 3 -PPT and their optimized version D opt 3 . The simulation results are shown in Fig. 5. We consider the ground state of an open chain of L = 14 sites. In a) and b), we select ℓ = 10 sites in the middle as subsystem A and divide it in two parts A = A 1 ∪ A 2 . We use the negativity as a reference to benchmark the efficiency of some entanglement conditions to detect entanglement between A 1 and A 2 .
In Fig. 5a), we calculate the p 3 -PPT, the D 3 , the optimal D opt 3 condition, and the Stieltjes condition using moments up to order five (see Methods). The convention we choose in the plot is that entanglement is detected whenever the value is positive. All the conditions work in most of the interval J z ∈ [ − 4, 1], where we expect entanglement to be sizeable, except for D 3 failing in the vicinity of J z = 1. The presence of entanglement is confirmed by the calculation of the negativity (red line). In Fig. 5b) we focus on the q = 1 sector. In this case, we observe that all conditions indicate the presence of at least a negative eigenvalue in the sector ρ Γ A ðq ¼ 1Þ -that is, they are informative about which sector for the reduced density matrix contributes to violating PPT. In this specific instance, SR is however not fundamental in detecting entanglement beyond what non-SR conditions can.
In Fig. 5c) and d), we carry out the same analysis for disconnected partitions. We consider A = A 1 ∪ A 2 , where A 1 consists of the first l/2 sites and A 2 of the last l/2, and L = 14, l = 10. In Fig. 5c) for J z~− 1.9 all the quantities except Stieltjes 5 are below zero, thus not revealing entanglement even though the negativity is positive. In this plot, one can also see that, for J z < − 2, the optimized condition D opt 3 detects entanglement whereas both p 3 -PPT and D 3 fail. This illustrates that the slight improvement obtained from the optimization (see Fig. 2) can be decisive to detect the entanglement of physically relevant states from the first three moments only. Fig. 4 Experimental SR-entanglement detection in a trapped ion quantum simulator using data obtained in ref. 28 . For a total system of N = 10 spins and subsystem A, B = [4,5], [6,7], we present in a) the SR-D 2 ratio and b) SR-p 3 -PPT ratio as a function of time for various symmetry sectors . In both panels, entanglement is detected in regimes where the corresponding global conditions do not reveal entanglement, as indicated in the shaded grey areas (values below unity). The points with error bars correspond to the values and uncertainties extracted from the experimental data from ref. 28 , whereas the dashed (solid) lines are theoretical simulations of unitary dynamics (taking decoherence effects into account), as detailed in refs. 28,34 .
Entanglement detection under constrained dynamics. As a third example, we study the detection of mixed-state entanglement in subsystems of constrained spin models after a global quantum quench. Such models have been realized experimentally with neutral atoms in optical tweezer arrays coupled to Rydberg states 52,53 . Below we simulate an experiment, in which moments of the partially transposed density matrix are obtained from a classical shadow involving global random unitaries available in current experimental setups. In particular, we demonstrate that periodic revivals of mixed state entanglement can be detected from the conditions D 3 and D 4 (Eqs. (11) and (12)) requiring only a small number of experimental runs.
We consider the Fibonacci chain with open boundary conditions described by the Hamiltonian where P i ¼ 0 j i i 0 h j are local projectors. As can be seen from Eq. (42), each spin undergoes independent Rabi-oscillations as long as the neighbouring spins are in their ground state 0 j i. This constraint breaks the tensor product structure of the Hilbert space (as it is the case in a lattice gauge theory 54 ). The model effectively resembles the experimental situation in 52 if the Rydberg atoms are driven close to resonance and neighbouring atoms cannot be simultaneously in the state 1 j i due to the Rydberg blockade mechanism. The Hamiltonian (42) has recently attracted great interest in the context of quantum many-body scarring 55,56 . In particular, performing a quantum quench on special unentangled product states results in long-lived periodic revivals which have been attributed to the existence of quantum scarred eigenstates in the many-body spectrum 55 .
In the following we study the conditions given in Eqs. (11) and (12), when a quench is performed from a product state that leads to kinetically constrained dynamics. To this end, the initial state Ψ 0 j i ¼ 10 j i N=2 is time evolved with the Hamiltonian (42) up to t = 50/Ω. Figure 6a) shows the local Z iexpectation values exhibiting long-lived persistent oscillations. This striking departure from a thermalizing behaviour is also reflected in the slow growth of entanglement entropy (Panel b). We now analyse the time-resolved behaviour of mixed state entanglement for a subsystem depicted in the inset of Fig. 6c). The revivals in the negativity indicate that spins in the subsystem get periodically entangled and disentangled with each other. Interestingly, the p 3 -PPT condition is unable to detect the revivals, while D 3 yields positive values at the first 3 peaks in the negativity. At later times, the D 3 fails to detect the entanglement present in the system, but this entanglement is still captured by D 4 .
Finally, we investigate the required number of experimental runs in order to measure the conditions up to the given error bar. The classical shadow is constructed by sampling bit strings from the quantum state after applying a global random unitary on subsystem A. At each point in Fig. 6c), we collect 5000-bit strings in different random basis. Such global random unitaries in Rydberg systems can be implemented via random quenches with local disorder potentials 57 . The entire estimation of the conditions is repeated 20 times in order to obtain statistical uncertainties. Note that statistical covariances among the measured moments trðρ Γ Þ n can give rise to nonuniform sizes for the error bars. In Fig. 6c)

DISCUSSION
The study of entanglement has a long and prominent history in a variety of disciplines. And with the advent of serious quantum technologies, reliable entanglement generation is more important than ever. This work provides a principled approach to reliably detect bipartite entanglement between subsystem A and subsystem B. We have presented a set of inequality conditions D k (1 ⩽ k ⩽ 2 |AB| ). Each D k is an inequality that involves the first k moments of the partially transposed density operator. Violation of a single inequality implies that the underlying density operator cannot have a PSD partial transpose. This in turn implies that the state must be entangled. Conversely, if the underlying state is not PPT, then, there must exist at least one D k that is violated. This motivates a sequence of one-sided entanglement tests. Start with D 3the lowest non-trivial conditionand check whether it is violated. If this is the case, we are done. If not, we move on to the next higher condition (D 4 ) and repeat until we find a violation. For states having an extensive conserved quantity (such as total magnetization, in the case of spin systems), both the density matrix and its partial transpose have a block-diagonal structure 48 . In this case, it is advisable to apply these conditions directly to individual symmetry sectors of the partial transpose. The resulting sequence of symmetry-resolved conditions is stronger in the sense that lower-order moments (of blocks of the partial transpose) suffice to detect entanglement. Importantly, this approach is not only conceptually sound, but also tractable from an experimental perspective. The classical shadows formalism 24 allows for reliably estimating moments of the partial transpose from randomized single-qubit measurements. We demonstrated how to include the experimentally relevant situation of nonidentical (however independent) copies in the analysis and derived error bounds and confidence intervals for D 2 , with a natural extension to quantities involving higher-order moments. Empirical evaluations complement our theoretical findings. Applications to several theoretical models, as well as experimental data, demonstrate both tractability and viability of our approach.
We are confident that this work opens up several interesting future research directions. Firstly, the sequence of D k 's is designed to detect bipartite entanglement in a reliable and experimentally accessible fashion. A natural next step is to try to extend similar ideas to multipartite entanglement detection, e.g. using non-linear entanglement witnesses 58,59 . Secondly, the complete sequence of D k 's is used to answer a binary question: is the partial transpose positive or not? Entanglement measures, like the negativity, address entanglement in a quantitative fashion, but are also harder to estimate. Is it possible to use moments (or other density matrix functionals) to define entanglement measures that are experimentally tractable? The statistical analysis of the estimation procedure is also far from complete. We have shown that independence between the states that are produced in each iteration of an experiment is enough to derive statistically sound confidence intervals for estimating matrix moments with classical shadows. This addresses the practically relevant case of drifting sources, but further extensions to correlated states would also be appealing. In this context, the quantum de Finetti theorem seems highly relevant. In future work, we will also investigate how importance sampling 60,61 and/or derandomization 24,62 can further improve moment estimation based on classical shadows. Finally, another promising direction of research would be to try to detect and characterize phase transitions in quantum mechanical Hamiltonians at finite (non-zero) temperatures. Quantum phase transitions at zero temperature originate from quantum fluctuations, whereas quantum phase transitions at finite temperature are due to thermal fluctuations. Following ref. 63 , quantum phase transitions at finite temperature can be studied using negativity. It would be interesting to investigate whether low-order PT moments, intimately related to the negativity, can also be used to this end.

Note
While completing the writing of the present work, we became aware of a work by Yu et al. 36 , in which similar questions have been addressed.

Descartes' rule of signs
Let A be a Hermitian matrix of dimension d. Its eigenvalues λ 1 , …, λ d are the roots of the characteristic polynomial For convenience, let us now consider the polynomial P( − t), which effectively replaces the positive eigenvalues of A by negative ones and vice versa. The coefficients of this polynomial can be expressed using the elementary symmetric polynomials in its roots, e i (λ 1 , …, λ d ), defined as for i = 1, …, d and with e 0 (λ 1 , …, λ d ) = 1. This yields For a polynomial with real roots (as it is the case here), Descartes' rule of sign states that the number of positive roots is given by the number of sign changes between consecutive elements in the ordered list of its nonzero coefficients (see ref. 64 and references therein). The matrix A is PSD iff P( − t) has only negative roots, which by Descartes' rule is the case iff there is no sign change in the ordered list of its non-zero coefficients, i.e., iff e i (λ 1 , …, λ d ) ⩾ 0 for all i = 1, …, d, since e 0 (λ 1 , …, λ d ) = 1.

Stieltjes moment problem
Given a sequence of moments, ðm n Þ d n¼0 , the (truncated) Stieltjes moment problem consists in finding necessary and sufficient conditions for the existence of a measure μ on the half-line [0, ∞) such that x n dμðxÞ; 8n 2 f0; ; dg: If such a measure exists, one may wonder whether it is unique or not. For our purposes, it will be enough to discuss only its existence. Defining the matrices and a solution to this problem can be stated as follows 65 . If d is oddsuch that d = 2k + 1 for some integer kthere exists such a measure μ if and only if AðkÞ ⩾ 0; BðkÞ ⩾ 0 and ðm k ; ; m 2kþ1 Þ T 2 R½AðkÞ; where, given a matrix M, the notation M ⩾ 0 indicates that M is PSD and RðMÞ denotes the range of M. If d is evensuch that d = 2k for some integer kthere exists such a measure μ if and only if AðkÞ ⩾ 0; Bðk À 1Þ ⩾ 0 and ðm kþ1 ; ; m 2k Þ T 2 R½Bðk À 1Þ: These solutions to the Stiltjes moment problem can be used to obtain entanglement conditions. Given λ 1 , …, λ r the eigenvalues of ρ Γ for some density matrix ρ, let us define the (atomic) eigenvalue distribution function dμðxÞ ¼ where δ is the Dirac delta distribution. If ρ is PPT, this density function has support on [0, ∞) and reproduces the moments of ρ Γ , as x n dμðxÞ ¼ p n ðρ Γ Þ: Therefore, according to the solution of the Stieltjes moment problem mentioned above, the moments of any PPT state necessary satisfy either condition (49) or (50), depending on the value of r. The violation of any of these conditions for a set of PT moments thus reveals that the corresponding state must be entangled. The range condition may require the knowledge of all the moments to be checked, but the PSD conditions can be broken into sets of simpler conditions. And some of them only involve low order moments. Indeed, it is well known (see e.g. 66 ) that a matrix is PSD if and only if all its principal minors are non-negative. For instance, looking at the principal minor at the intersection of the first two rows and columns of B(k), one obtains the condition m 1 m 3 À ðm 2 Þ 2 ⩾ 0. This condition is nothing but the p 3 -PPT condition (which we know is useful to detect entanglement 23 ). Extending this principal minor to the third row and column, one gets another PPT condition: We call this condition Stieltjes 5 . We illustrate in the Applications section (c.f. Fig. 5) that this condition is also useful for entanglement detection. Numerical computations suggest that this condition is a powerful tool to detect the entanglement of random mixed states, in the sense that it detects more random entangled states than either p opt 3 or D 5 . Not all Stieltjes moment conditions are this powerful, though. For instance, the principal minor condition for the first two rows and columns of A(k) is trivial.
Note that, because we consider here an atomic density function, we have m 0 = r. We could naturally renormalize the density function so that m 0 = 1, but it would imply a re-scaling of the first moment, i.e., the trace of ρ Γ would be 1/r. Since the partial transpose and the density function cannot be normalized at the same time, we chose to keep normalized partial transposes.
Optimizing conditions involving moments up to degree three Given a PSD matrix A, with non-zero eigenvalues λ 1 , …, λ r , for some r 2 ½1; dim A, consider the Lagrangian function Lðλ 1 ; ; λ r ; C 1 ; (54) where C 1 and C 2 are Lagrange multipliers.
Here we show that, for all 1⩽r⩽ dim A, the stationary points (λ 1 , …, λ r ) of the Lagrangian function (54) are such that the variables λ i can take at most two distinct values. These stationary points, for which the derivatives of the Lagrangian (54) with respect to each variable vanish, satisfy the set of equations We first sum up Eq. (55) for all values of i and then insert Eqs. (56) and (57) into it. This yields Inserting this relation into Eq. (55) and considering this equation for two distinct values of i, say 1 and k ≠ 1, one can eliminate the variable C 1 to get a relation between λ 1 and λ k . After some algebra, one finds λ k ¼ λ 1 or λ k ¼ λ 1 p 1 À p 2 λ 1 r À p 1 : Since this argument holds for any k ≠ 1, it must hold that the eigenvalues λ i are either all equal or can only take two distinct values. In the first case, in which all the eigenvalues are equal, one obtains the isolated points (p 2 , p 3 ) = (1/r, 1/r 2 ) in Fig. 2 in the main text. In the second case, the rank r PSD matrices corresponding to the stationary points of the Lagrangian (54) have a spectrum with r a degenerate eigenvalues λ a and r − r a eigenvalues λ b . Assuming, without loss of generality, λ a > λ b , one can show that the minimal value of the third moment is obtained when r a = r − 1.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.