# Many-Body Entropies and Entanglement from Polynomially Many Local Measurements

Benoît Vermersch, Marko Ljubotina, J. Ignacio Cirac, Peter Zoller, Maksym Serbyn, Lorenzo Piroli

► To cite this version:

Benoît Vermersch, Marko Ljubotina, J. Ignacio Cirac, Peter Zoller, Maksym Serbyn, et al., Many-Body Entropies and Entanglement from Polynomially Many Local Measurements. Physical Review X, 2024, 14 (3), pp.031035. 10.1103/PhysRevX.14.031035. hal-04918172

HAL Id: hal-04918172

https://hal.science/hal-04918172v1

Submitted on 29 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Many-body entropies and entanglement from polynomially-many local measurements

Benoît Vermersch,¹,²,³ Marko Ljubotina,⁴ J. Ignacio Cirac,⁵,⁶ Peter Zoller,²,³ Maksym Serbyn,⁴ and Lorenzo Piroli⁷

¹Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France

²Institute for Theoretical Physics, University of Innsbruck, 6020 Innsbruck, Austria

³Institute for Quantum Optics and Quantum Information of the Austrian Academy of Sciences, 6020 Innsbruck, Austria

⁴Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria

⁵Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Straße 1, D-85748 Garching, Germany

⁶Munich Center for Quantum Science and Technology (MCQST), Schellingstraße 4, D-80799 München, Germany

⁷Dipartimento di Fisica e Astronomia, Università di Bologna and INFN, Sezione di Bologna, via Irnerio 46, I-40126 Bologna, Italy

(Dated: July 1, 2024)

Estimating global properties of many-body quantum systems such as entropy or bipartite entanglement is a notoriously difficult task, typically requiring a number of measurements or classical post-processing resources growing exponentially in the system size. In this work, we address the problem of estimating global entropies and mixed-state entanglement via partial-transposed (PT) moments, and show that efficient estimation strategies exist under the assumption that all the spatial correlation lengths are finite. Focusing on one-dimensional systems, we identify a set of approximate factorization conditions (AFCs) on the system density matrix which allow us to reconstruct entropies and PT moments from information on local subsystems. This yields a simple and efficient strategy for entropy and entanglement estimation. Our method could be implemented in different ways, depending on how information on local subsystems is extracted. Focusing on randomized measurements (RMs), providing a practical and common measurement scheme, we prove that our protocol only requires polynomially-many measurements and post-processing operations, assuming that the state to be measured satisfies the AFCs. We prove that the AFCs hold for finite-depth quantum-circuit states and translation-invariant matrix-product density operators, and provide numerical evidence that they are satisfied in more general, physically-interesting cases, including thermal states of local Hamiltonians. We argue that our method could be practically useful to detect bipartite mixed-state entanglement for large numbers of qubits available in today's quantum platforms.

## CONTENTS

I. Introduction 1

II. Entanglement and randomized measurements 3
A. Mixed-state entanglement and PPT conditions 3

B. Randomized measurements and classical shadows 3

III. A toy model for the estimation protocol 5

A. Purity estimation: Factorization formula 5

B. The normalized PT-moment estimation 7

IV. Finite-range correlated states 8

A. The approximate factorization conditions 8

B. Estimation protocols for states satisfying the AFCs 8

C. Proof of AFCs in Matrix Product Density Operators 9

D. Numerical study of AFCs in Gibbs states 10

E. Practical feasibility of purity estimation 12

V. Mixed-state entanglement detection 14

A. The $p_n$-PPT conditions and mixed-state entanglement probes 14

B. Mixed-state entanglement detection in FDQC and Gibbs states 15

VI. Outlook 16

Acknowledgments 17

A. Factorization condition for FDQC states 17

1. The purity 17

2. The PT moments 18

B. Statistical-error analysis in FDQC states 19

1. Statistical-error analysis for the purity 19

2. Statistical-error analysis for the PT moments 20

C. Statistical-error analysis for states satisfying the AFCs 21

D. AFCs and MPDOs 23

E. Details on the numerical computations 24

F. Technical details on the PPT conditions 25

References 26

## I. INTRODUCTION

In the context of today's digital quantum technologies [1–9], an outstanding challenge is to devise measure-
ment schemes of many-qubit states which are efficient and yet simple enough to be performed in current noisy intermediate-scale quantum (NISQ) devices [10, 11]. This problem has motivated new ideas and protocols to improve our ability to characterize complex quantum states. A notable example is that of the so-called randomized-measurement (RM) toolbox [12–15], which has provided us with novel opportunities to experimentally investigate entanglement, a cornerstone in both quantum-information [16, 17] and quantum many-body theory [18, 19].

For instance, paralleling earlier experiments based on quantum interference [20–22], cf. also Refs. [23–26], pure-state entanglement can be detected by exploiting the two-copy representation of subsystem purities [13, 15, 27–29], as is now routinely done in various experimental platforms [27, 30–33]. RMs have been also applied to study both mixed-state bipartite entanglement, based on the estimation of the so-called partial-transposed (PT) moments [34, 35], and multipartite entanglement, as characterized by the quantum Fisher information [36–39].

Despite these developments, estimating global properties of quantum systems with a very large number of qubits $N$ remains a difficult task. In particular, while RM approaches to estimate any local observable are more efficient than performing full state tomography [40–55], estimating global properties requires performing exponentially-many measurements or post-processing operations. Therefore, it is unfeasible to significantly scale up existing protocols to probe global purities and bipartite entanglement [13, 15, 27–29, 35], raising the question of whether these quantities will be experimentally accessible at all as larger NISQ devices become available.

In this work, we address precisely the problem of estimating global entropies and PT moments in many-qubit systems, and show that efficient strategies exist under the assumption that all spatial correlation lengths in the system are finite. This condition encompasses a large class of physically interesting cases, including ground and thermal states of local Hamiltonians, and is thus very natural when considering NISQ devices from the point of view of quantum simulation [56–58].

In more detail, we focus on one-dimensional (1D) systems and put forward a strategy for the estimation of entropies and PT moments, requiring only polynomially-many measurements and post-processing operations. Our protocol is provably accurate for states satisfying a set of approximate factorization conditions (AFCs), which express the absence of long-range correlations and which are shown to be a general feature of short-range correlated states. The basic idea, which is conveyed in Fig. 1, is that the AFCs allow one to reconstruct global purities and PT moments from local information. The efficiency of our method is independent of how local information is extracted (for instance, this could be done via full tomography of certain reduced density matrices or via protocols using multiple physical copies, as discussed later). However, we will focus on an implementation

<figure><img src="image_2.png" /><figcaption>Figure 1. A 1D quantum system $\mathcal{S}$ is partitioned into adjacent intervals $I_j$ of size $k$, $|I_j| = k$. In this work, we consider extracting the purities over $I_j$ and $I_j \cup I_{j+1}$ using the classical-shadow approach [15]. Performing $M$ local measurements with respect to randomly chosen bases and classical post-processing operations, they yield faithful estimators for $\text{Tr}[\rho_I^2]$, denoted by $\mathcal{P}_2^{(e)}[I]$. Such subsystem-purity estimators are then combined to obtain a prediction $\mathcal{P}_2^{(e)}$ for the global purity. The accuracy of the method is controlled by $k$ and the state correlation lengths.</figcaption></figure>

The proposed approach is different from full tomographic methods relying on prior assumptions on the system state. There, a common strategy is to assume that the latter is described by a sufficiently simple ansatz wavefunction, and estimate the values of its free parameters. This logic was first put forward in the context of matrix-product states (MPSs) [59–62], while recently extended to matrix-product density operators (MPDOs) [63–67], Gibbs states [68–73], permutation-invariant states [74–76], low-rank states [77], stabilizers [78], tensor- and neural-network wavefunctions [79–84]. While these methods may be practically very useful, they also face drawbacks. For instance, as we discuss in more detail in Sec. III A, an accurate estimation of, say, global purities might require reconstructing the state up to exponential precision, thus leading to unpractical overheads in $N$. On the contrary, our approach does not learn the state wavefunction, targeting purities and PT moments directly, cf. Fig. 1.

Finally, we mention that our ideas could be extended, in some cases straightforwardly, to probe other types of quantities generally requiring exponentially many measurements, including participation entropies [85–91] or stabilizer Rényi entropies [92], which were recently con-
sidered in the many-body setting [92–99]. In addition, while we will focus on 1D systems, where analytic and numerical analyses are simpler, we expect that our approach could be generalized to higher spatial dimensions. Therefore, our work also opens up a number of important directions for future research.

The rest of this manuscript is organized as follows. After reviewing a few preliminary notions and tools in Sec. II, we start by introducing the main ideas underlying our approach in Sec. III. To this end, we consider the class of so-called finite-depth quantum-circuit (FDQC) states, which provide an ideal toy model for short-range correlated many-body quantum states. Their minimal structure allows us to remove unnecessary technical complications from the discussion, and present the logic of our method in the simplest possible setting. We consider both purities (Sec. III A) and PT moments (Sec. III B), working out efficiency performance guarantees for their estimation.

The most general form of our protocol is presented in Sec. IV. After introducing the AFCs in Sec. IV A, we rigorously derive performance guarantees for the accurate estimation of the purity and PT moments (Sec. IV B) under the assumption that the state to be measured satisfies the AFCs. We then discuss the generality of the AFCs, proving that they are satisfied in translation-invariant MPDOS (Sec. IV C), and presenting numerical evidence for their validity in thermal states of local Hamiltonian, cf. Sec. IV D. We also present a full classical simulation of the measurement protocol (Sec. IV E), studying in a concrete example the typical number of measurements required to estimate the bipartite purity of area-law pure states.

Next, in Sec. V we discuss a few natural examples of highly mixed states where bipartite entanglement can be detected for large system sizes by estimating just the first few PT moments. Since the latter can be very simply extracted by our approach, we argue that our results could be practically useful to probe mixed-state entanglement in experimentally available noisy quantum platforms. Finally, we report our conclusions in Sec. VI, while the most technical parts of our work are consigned to several appendices.

## II. ENTANGLEMENT AND RANDOMIZED MEASUREMENTS

### A. Mixed-state entanglement and PPT conditions

We consider a system of $N$ qubits, denoted by $S$. Given a region $I \subset S$, the associated Hilbert space is $\mathcal{H}_I \simeq \mathbb{C}^{\otimes 2|I|}$, where we denoted by $|I|$ the number of qubits in $I$. We will be interested in the Rényi entropies of the region $I$

$$S_n[I] = \frac{1}{1-n} \ln \mathcal{P}_n[I], \quad (1)$$

where

$$\mathcal{P}_n[I] = \mathrm{Tr}[\rho_I^n], \quad (2)$$

and $\rho_I$ is the reduced density matrix on the region $I$. For $n=2$, $\mathcal{P}_2[I]$ coincides with the purity, which is a simple probe for the subsystem entropy, with $\mathcal{P}_2[I] = 1$ and $\mathcal{P}_2[I] = 2^{-|I|}$ for pure and maximally mixed states, respectively.

Consider now a partition of $S$ into two disjoint sets $S = A \cup B$, yielding $\mathcal{H}_S = \mathcal{H}_A \otimes \mathcal{H}_B$. If the system is in a pure state $|\psi\rangle_{AB} \in \mathcal{H}_S$, its bipartite entanglement is quantified by the Rényi entropies $S_n(\rho_A)$ [17]. Conversely, when the state of the system is mixed, its entanglement can be quantified by the logarithmic negativity [100, 101]

$$\mathcal{E}(\rho) = \ln \sum_j |\lambda_j|, \quad (3)$$

where the sum is over all eigenvalues $\{\lambda_j\}_j$ of the operator $\rho_{AB}^{T_A}$, and $(\cdot)^{T_A}$ denotes partial transpose with respect to subsystem $A$. The spectrum of the PT density matrix, and hence the logarithmic negativity, is completely fixed by the PT moments

$$p_n = \mathrm{Tr} \left[ \left( \rho_{AB}^{T_A} \right)^n \right], \quad (4)$$

for $n=1, 2, \dots, \dim(\mathcal{H}_A \otimes \mathcal{H}_B)$. Note that $p_1 = 1$, while the second PT moment coincides with the purity, $p_2 = \mathcal{P}_2$ [35].

The importance of the PT moments is two-fold. On the one hand, they can be accessed directly via RMs [34, 35] (or using quantum interference [102, 103]), see Sec. II B. On the other hand, the knowledge of the first few PT moments is enough to certify bipartite entanglement based on the non-positivity of the partial-transposed density matrix [34, 35, 104–106], or to detect different types of entanglement structures [107]. In this work, we will consider a particular set of conditions on the PT moments to certify bipartite entanglement, which were derived in Refs. [35, 106] and which we call the $p_n$-PPT conditions. Denoting by SEP the set of separable, i.e. not entangled, states in $S$, the $p_n$-PPT conditions take the form

$$\rho \in \mathrm{SEP} \Rightarrow p_n p_{n-2} \ge p_{n-1}^2. \quad (5)$$

Therefore, when the state of the system $\rho$ violates the $p_n$-PPT conditions, it is entangled, and the difference $p_{n-1}^2 - p_n p_{n-2}$ is a probe for mixed-state entanglement. Note that, for $n=3$, Eq. (5) coincides with the relation first derived in Ref. [35]. These conditions are in general not optimal and are a strict subset of those derived in Ref. [105, 106]. Still, their simplicity makes them particularly convenient for our purposes.

### B. Randomized measurements and classical shadows

The power of RM schemes lies in the fact that they need not be tailored to a specific property of the sys-
tem. Rather, one performs measurements which are randomly sampled from a fixed ensemble independent of the observable of interest. Subsequently, the outcomes are processed differently depending on the quantity to be estimated [28, 29, 108, 109]. Denoting by $\rho$ the system density matrix, this approach gives us access to all observable expectation values $\text{Tr}[\rho O]$ and, more generally, to multi-copy objects of the form $\text{Tr}[\rho^{\otimes n}O]$, where the integer $n \ge 1$ is called the copy (or replica) index.

In this section, we recall the basic aspects of RMs used in our work. While the logic explained in the next sections may be implemented in different ways, we will focus on a set of protocols making use the (local) classical shadows introduced in Ref. [15], a prominent element in the RM toolbox [12]. We briefly recall the main aspects of the formalism, while we refer to Refs. [12, 15] for a thorough introduction.

In what follows, we denote by $|0\rangle_j$ and $|1\rangle_j$ the basis elements of the local computational basis corresponding to qubit $j$, spanning $\mathcal{H}_j \simeq \mathbb{C}^2$. In the classical-shadow framework, one performs a set of $M$ measurements (one per experimental run, each labeled by an integer $r$), consisting of local unitary operations $\prod_j u_j^{(r)}$ followed by a projective measurement onto the computational basis

$$|k_1, \dots, k_N\rangle = \otimes_{j=1}^N |k_j\rangle_j, \quad (6)$$

with $k_j = 0, 1$. The unitaries are sampled from a Haar-random ensemble, identically and independently for each qubit $j$ and experimental run $r$, see Fig. 2. Denoting by $\{k_j^{(r)}\}$ the set of outcomes of this two-step process, the values $\{k_j^{(r)}\}$ and the unitaries $\{u_j^{(r)}\}$ are used to define the so-called classical shadows

$$\rho_S^{(r)} = \bigotimes_{i \in S} \left[ 3 (u_i^{(r)})^\dagger |k_i^{(r)}\rangle \langle k_i^{(r)}| u_i^{(r)} - 1 \right]_2, \quad (7)$$

which can be classically stored in an efficient way.

As mentioned, the measurement protocol does not depend on the observable of interest. Rather, one adapts the post-processing operations on the classical shadows based on the quantity to be estimated. For instance, given any observable $O$, an estimator for its expectation value is

$$\hat{o} = \frac{1}{M} \sum_{r=1}^{M} \text{Tr}[O \rho_S^{(r)}]. \quad (8)$$

It is easy to see that $\hat{o}$ is faithful, i.e. unbiased, while different bounds for the statistical variance of this estimator may be derived depending on the locality properties of $O$ [12, 15].

It is important to recall that classical shadows also give access to entropy and PT moments. Let us consider in particular the purity, which is given by Eq. (2) for $n = 2$. For this quantity, the classical shadows allow us to write the estimator [15]

$$\mathcal{P}_2^{(e)}[I] = \frac{1}{M(M-1)} \sum_{r \neq r'} \text{Tr} \left( \rho_I^{(r)} \rho_I^{(r')} \right), \quad (9)$$

<figure><img src="image_3.png" /><figcaption>Figure 2. (a): Within the classical-shadow approach, each measurement process consists of random on-site unitaries followed by local projective measurements. Their outcomes are stored and later post-processed to construct the classical shadows. (b): Pictorial representation of a local finite-depth quantum circuit. The gates are arranged in a brickwork pattern, forming layers of mutually commuting unitary operations acting on pairs of neighboring qubits. Lower (upper) dangling legs correspond to the input (output) qubits. The depth of the circuit is the number of applied layers. In this picture, the depth is $\ell = 3$.</figcaption></figure>

where $\rho_I$ is defined as in Eq. (7). $\mathcal{P}^{(e)}[I]$ is a faithful estimator. Note that, alternatively, one can use another estimator for the purity [13, 27, 29] that provides robustness against miscalibration errors using the same data. The statistical errors associated with $\mathcal{P}_2^{(e)}[I]$ are quantified by its variance, which can be bounded by [15, 35, 38]

$$\text{Var} \left[ \mathcal{P}_2^{(e)}[I] \right] \le 4 \left( \frac{2^{|I|} \mathcal{P}_2[I]}{M} \right) + 2 \left( \frac{2^{2|I|}}{M-1} \right)^2, \quad (10)$$

where $|I|$ denotes the number of qubits in $I$. This bound is known to be essentially optimal [35, 38], telling us that an exponentially large number of measurements is needed to estimate the purity.

The PT moments (4) can be treated similarly. In this case, one constructs the estimator [35]

$$p_n^{(e)}[AB] = \frac{1}{n!} \binom{M}{n}^{-1} \times \sum_{r_1 \neq r_2 \neq \dots \neq r_n} \text{Tr} \left[ [\rho_{AB}^{(r_1)}]^{T_A} \dots [\rho_{AB}^{(r_n)}]^{T_A} \right]. \quad (11)$$

Once again, $p_n^{(e)}[AB]$ is faithful and it is possible to derive explicit bounds on its variance, although it becomes increasingly involved for higher $n$. For instance, for $n = 3$ one has [35, 38]

$$\text{Var} \left[ p_3^{(e)}[AB] \right] \le 9 \frac{2^{|AB|}}{M} \text{Tr} \left( \rho_{AB}^4 \right) + 18 \frac{2^3 |AB|}{(M-1)^2} p_2[AB] + 6 \frac{2^6 |AB|}{(M-2)^3}. \quad (12)$$

Estimating the statistical errors by the variance, Eqs. (10) and (12) imply that, in order to guarantee an accurate reconstruction of the purity and PT moments of the system, exponentially-many measurements in its size are needed. As mentioned, this makes it unfeasible to significantly scale up previous experiments making use
of this strategy [13, 15, 27-29, 35]. The goal of this work is to show that these limitations may be overcome under assumptions which are very common in the context of many-body physics, and thus also natural from the point of view of quantum simulation. Namely, we will put forward a set of protocols for entropy and entanglement estimation which are provably efficient assuming that all spatial correlation lengths of the state to be measured are finite. We will focus on 1D systems, where analytic and numerical analyses are simpler, although we expect that our approach could be generalized to higher spatial dimensions, see Sec. VI.

### III. A TOY MODEL FOR THE ESTIMATION PROTOCOL

In this section we introduce the main ideas underlying our approach, focusing on a simplified setting where we can get rid of unnecessary technical complications. We analyze the case of FDQC states, where the state to be measured is prepared by a shallow (local) quantum circuit

$$ \rho = U^{(\ell)} \left( \bigotimes_{j=1}^{L} \sigma_j \right) [U^{(\ell)}]^{\dagger}, \quad (13) $$

where L is length of the system, i.e. the number of qubits (we use the letter L instead of N, as in the previous section, to emphasize that we are focusing on the 1D case). Here $\sigma_j$ are arbitrary single-qubit density matrices, while $U^{(\ell)}$ is a local circuit of depth $\ell$, namely $U^{(\ell)} = V_\ell \cdots V_2 V_1$, where $V_j$ contains quantum gates acting on disjoint pairs of nearest-neighbor qubits, cf. Fig. 2. We will assume that $\ell$ is fixed, i.e. not increasing with the system size. We do not ask for translation symmetry and, unless specified otherwise, assume open boundary conditions. The gates making up $U^{(\ell)}$ can be arbitrary, i.e. they need not be taken out of some finite gate set. The FDQC states (13) have a very simple structure from the point of view of many-body physics, but are known to approximate physically-interesting states such as MPSs [110, 111] and, more generally, ground states of gapped local Hamiltonians [112-115].

Throughout this section, we will assume that we know that the state of the system is exactly of the form (13) for some finite $\ell$. We develop a protocol to efficiently estimate the purity and the PT moments of such a state, based on this knowledge. The protocol only takes the depth of the circuit, $\ell$, as an input and does not make use of state tomography. In fact, it will be later generalized replacing the assumption of the FDQC structure with the AFCs, which make no explicit reference to the form of the state wavefunction.

<figure><img src="image_4.png" /><figcaption>Figure 3. (a): Partition considered in Eq. (14): A 1D interval $I$ is divided into three disjoint regions $A, B$, and $C$, where $B$ separates $A$ and $C$. (b): Partition considered in Eq. (26): A 1D interval $I$ is divided into two halves, $A$ and $B$. Each half is partitioned into two regions $A = A_1 \cup A_2$, $B = B_1 \cup B_2$</figcaption></figure>

#### A. Purity estimation: Factorization formula

The starting point of our method is a factorization property for powers of the system density matrix. Let $I$ be any interval of adjacent qubits ($I$ can coincide with the full system $S$). Consider a partition $I = A \cup B \cup C$, where $B$ separates $A$ and $C$ and denote by $|B|$ the number of qubits in $B$, cf. Fig. 3(a). Based on the fact that state is prepared by a finite depth circuit, Eq. (13), one can prove

$$ \mathrm{Tr}(\rho_I^2) = \frac{\mathrm{Tr}_{AB}(\rho_{AB}^2)\mathrm{Tr}_{BC}(\rho_{BC}^2)}{\mathrm{Tr}_B(\rho_B^2)}, \quad (14) $$

for any partition with $|B| \ge 2\ell - 1$, where $\rho_X$ denotes the density matrix reduced to the interval $X$ and $XY$ is a short-hand notation for $X \cup Y$. Now, let $\{I_j\}_{j=1}^R$ be a collection of adjacent intervals covering the system $S$, with $|I_j| = k \ge 2\ell - 1$, and assume without loss of generality that $R = L/k$ is an integer, cf. Fig. 1. By applying Eq. (14) iteratively, we arrive at

$$ \mathrm{Tr}(\rho^2) = \frac{\prod_{j=1}^{R-1} \mathrm{Tr}_{I_j \cup I_{j+1}}(\rho_{I_j \cup I_{j+1}}^2)}{\prod_{j=2}^{R-1} \mathrm{Tr}_{I_j}[\rho_{I_j}^2]}, \quad (15) $$

where $\rho$ is the system density matrix (13). A proof of Eqs. (14) and (15) is given in Appendix A.

Equation (15) is remarkable because its right hand side (RHS) only depends on the purities of subsystems containing up to $2k$ qubits, and thus gives us a natural basis for an efficient estimation strategy of the global purity. The idea is to first reconstruct the “local” purities $\mathrm{Tr}[\rho_I^2]$ (with either $I = I_j$ or $I = I_j \cup I_{j+1}$), and subsequently plug their estimated values in the RHS of Eq. (15), yielding the global purity estimate. This method is more efficient than directly targeting $\mathrm{Tr}(\rho^2)$, as the local purities can be reconstructed from a number of measurements growing exponentially in $k$, not in $L$, see Eq. (10). Below, we give a more precise description of the protocol, proving that an accurate estimation of the purity only requires a polynomial (in $L$) number of measurements and post-processing operations.

We consider performing $M_I$ measurements (i.e. experimental runs) to estimate the purity of each interval $I$, with $I = I_j$ or $I = I_j \cup I_{j+1}$, see Eq. (15). The set of
measurements $M_I$ is only used to determine the purity over $I$ and for simplicity we perform the same number of measurement $M_I = M$ for each interval. Although this increases the total number of experiments $M_T$ by a factor $L$, since

$$M_T = \sum_I M_I \propto LM, \qquad (16)$$

this procedure guarantees that statistical errors associated with distinct regions are independent and facilitate their rigorous analysis. For each interval, we then construct the estimators $\mathcal{P}^{(e)}[I]$ given in Eq. (9), and define

$$r_2^{(e)} = \frac{\prod_{j=1}^{R-1} \mathcal{P}_2^{(e)}[I_j \cup I_{j+1}]}{\prod_{j=2}^{R-1} \mathcal{P}_2^{(e)}[I_j]}, \qquad (17)$$

which is our estimator for the global purity. Since the latter is typically exponentially small in $L$, we quantify the accuracy of $r_2^{(e)}$ by the relative error

$$\varepsilon_r = \left| \frac{r_2^{(e)}}{\mathcal{P}_2} - 1 \right| . \qquad (18)$$

Next, we seek to bound the number of measurements required to guarantee that $\varepsilon_r$ is small with high probability. This problem is solved in Appendix B 1, where we prove the following result: for any arbitrarily small $\delta > 0$, choosing

$$M \geq \max \left\{ 2^{8k}, L^2 \frac{2^{4k+10}}{k^2 \delta^2} \right\}, \qquad (19)$$

the probability that $|\varepsilon_r| \geq \delta$ satisfies

$$\Pr[|\varepsilon_r| \geq \delta] \leq \frac{2^{4k+11} L^3}{\delta^2 k^3 M}, \qquad (20)$$

where we recall that $k = |I_j|$. The proof of Eq. (20) is based on a careful analysis of how the error on each factor in (17) affects the global purity, making use of the statistical independence of $\mathcal{P}_2^{(e)}[I]$ for different $I$ and the so-called Chebyshev's inequality. Eqs. (16), (19), and (20) then imply that polynomially-many (in $L$) measurements and post-processing operations are enough to accurately estimate the purity, as anticipated.<fnref n="1" /> We can rephrase this result in a more transparent way. For a given confidence level $\gamma = \Pr[|\varepsilon_r| < \delta]$, Eq. (20) implies that it is enough to take a number of measurements

$$M \geq \frac{2^{4k+11} L^3}{\delta^2 k^3 (1 - \gamma)}. \qquad (21)$$

Before concluding this section, a few remarks are in order. First, it is important to note that the individual factors in the formula (17) could also be extracted using different RM schemes or even standard tomography for the density matrices $\rho_I$. Indeed, the possibility of estimating the global purity from polynomially-many measurements does not depend on the fact that we are using classical shadows, but rather on the factorization property (17) (since its RHS only involves local subsystems). Still, classical shadows and RMs in general offer several advantages compared to state tomography. For instance, while RM schemes to estimate $\mathrm{Tr}[\rho_I^2]$ require a number of measurements scaling exponentially in $|I|$, the exponents are typically favorable compared to tomography [27, 31, 116]. In addition, classical shadows make it very easy to rigorously bound statistical errors, facilitating the analyses presented throughout this work.

Second, taking the logarithm of Eq. (15), we obtain

$$S^{(2)}(\rho) = \sum_{j=1}^{R-1} S^{(2)}(\rho_{I_j \cup I_{j+1}}) - \sum_{j=2}^{R-1} S^{(2)}(\rho_{I_j}), \qquad (22)$$

where $S^{(2)}(\rho) = -\ln \mathrm{Tr} \rho^2$ is the second Rényi entropy.<fnref n="2" /> Therefore, our results can equivalently be formulated in terms of Rényi entropies, rather than purities (note that a small relative error on the latter implies a small additive error on the Rényi entropy).

Finally, as the state (13) can be represented exactly as a matrix-product operator (MPO) [129], it is instructive to compare our strategy with those based on MPDO tomography [63–67]. In many cases, these methods can efficiently provide an estimate of the system density matrix $\rho^{(e)}$, satisfying

$$\|\rho^{(e)} - \rho\|_1 \leq \delta, \qquad (23)$$

where $\| \cdot \|_1$ denotes the trace norm, while $\delta = O(L^{-\alpha})$ is a small parameter vanishing polynomially in $L$. While this approximation allows us to accurately estimate the expectation value of any local observable, it might not be enough to extract the global Rényi-2 entropy. Indeed, given $\rho$ and $\sigma$ with $\delta = \|\rho - \sigma\|_1$, we have the following bound, which is known to be tight in general [130, 131]

$$|S^{(2)}(\rho) - S^{(2)}(\sigma)| \leq 2^L \left[1 - (1 - 2\delta)^2 - \frac{4\delta^2}{(2^L - 1)}\right] \\ \sim 2^L \delta. \qquad (24)$$

Therefore, a precise estimation of the Rényi-2 entropy requires an exponentially accurate reconstruction of the system density matrix, typically leading to unpractical

<footnote n="1">¹ The polynomial scaling of the post-processing operations follows straightforwardly from the definition (9), because the number of elements in each sum is $\sim M^2$, where $M = O(L^3)$.</footnote>

<footnote n="2">² We note that similar formulas previously appeared (albeit for the von Neumann, rather than for the Rényi entropy) in the context of approximate Markov-chain states [117–120], see also [121–128].</footnote>
overheads in $L^3$. Instead, our method gets around this technical issue, as it does not rely on state tomography.

## B. The normalized PT-moment estimation

The ideas presented in the previous section may be applied to the PT moments, although a few subtleties must be taken into account. Denoting by $p_n^{(e)}$ our estimate for $p_n$, one is tempted to ask for a protocol that makes the relative error $|p_n^{(e)}/p_n - 1|$ sufficiently small. This is, however, problematic: contrary to the moments $\mathcal{P}_n$, it is non-trivial to bound $|p_n|$ from below by a positive number. In order to get around this issue, we define the normalized PT moment

$$ \tilde{p}_n[AB] = \frac{p_n[AB]}{\mathcal{P}_n[A]\mathcal{P}_n[B]}, \quad (25) $$

and, instead of focusing on the relative error, we ask for an accurate estimate of $\tilde{p}_n$ up to a small additive error. The motivation for this choice is two-fold. First, as it will be clear later, for FDQC states one can show that $\tilde{p}_n[AB]$ is independent of the system size. Therefore, contrary to the purity, one does not have to deal with numbers which are exponentially small in $L$. Second, we will see in Sec. V that entanglement certification based on the $p_n$-PPT conditions requires approximating $\tilde{p}_n[AB]$ up to a small additive error.

In order to estimate the normalized PT moment (25), we once again rely on certain factorization properties of the density matrix. Consider the FDQC state (13) and take a partition of the system as in Fig. 3(b), with $|B_1| = |A_2| = k \ge 2\ell - 1$, where $\ell$ is the depth of the circuit. As we show in Appendix A 2, one can prove

$$ \mathrm{Tr} \left[ \left( \rho_{AB}^{T_A} \right)^n \right] = \mathrm{Tr} \left[ \left( \rho_{A_2 B_1}^{T_{A_2}} \right)^n \right] \frac{\mathrm{Tr}_A[\rho_A^n] \mathrm{Tr}_B[\rho_B^n]}{\mathrm{Tr}_{A_2}[\rho_{A_2}^n] \mathrm{Tr}_{B_1}[\rho_{B_1}^n]} \quad (26) $$

This equation implies

$$ \tilde{p}_n[AB] = s_n[A_2 B_1], \quad (27) $$

where we defined

$$ s_n[XY] = \frac{\mathrm{Tr} \left[ \left( \rho_{XY}^{T_X} \right)^n \right]}{\mathrm{Tr}_X[\rho_X^n] \mathrm{Tr}_Y[\rho_Y^n]} \quad (28) $$

Note that $s_n[A_2B_1]$ only depends on the density matrix of a local subsystem and does not scale with the system

size, as anticipated.<fnref n="4" /> Therefore, it is possible to obtain an accurate estimate from a number of measurements and post-processing operations independent of $L$, for any arbitrary small additive error.

To see this explicitly, we consider the following protocol. We introduce an estimator for $s_n[A_2B_1]$, namely

$$ s_n^{(e)} = \frac{p_n^{(e)}[A_2 B_1]}{\mathcal{P}_n^{(e)}[A_2] \mathcal{P}_n^{(e)}[B_1]} \quad (29) $$

The estimators for the moments $\mathcal{P}_n^{(e)}[A_2]$, $\mathcal{P}_n^{(e)}[B_1]$ and the PT-moment $p_n^{(e)}[A_2B_1]$ are defined in Eq. (11). We compute each of them out of $M$ different classical shadows, so that the total number of experimental runs is $M_T = 3M$. This is done so that the statistical errors are independent, facilitating their analysis. After these steps, we obtain an estimate for $s_n[A_2B_1]$, and thus, due to Eq. (26), for $\tilde{p}_n[AB]$.

Importantly, we can bound the additive statistical error which affects our estimate, namely

$$ \varepsilon_a = \left| s_n^{(e)} - \tilde{p}_n[AB] \right| \quad (30) $$

For instance, focusing for simplicity on the case $n = 3$, we prove in Appendix B 2 the following result: for any small $\delta > 0$ and choosing

$$ M \ge 27 \frac{2^{11k+9}}{\delta^2}, \quad (31) $$

we have

$$ \mathrm{Pr} \left[ \left| s_3^{(e)} - \tilde{p}_3[AB] \right| \ge \delta \right] \le 81 \frac{2^{11k+9}}{M \delta^2}, \quad (32) $$

where $|A_2| = |B_1| = k \ge 2\ell - 1$, with $\ell$ being the circuit depth. This result is proved by the same techniques used to derive Eq. (20). Again, it is useful to rephrase our result in terms of the confidence level $\gamma$, yielding

$$ M \ge 81 \frac{2^{11k+9}}{(1-\gamma)\delta^2} \quad (33) $$

Equation (32) states that a number of measurements not scaling with the system size is enough to guarantee that $\tilde{p}_3$ is approximated with arbitrary precision and high probability. Explicit performance guarantees such as (32) for higher integer values of $n$ are more cumbersome to derive, and will be omitted. Still, based on the analysis presented in Appendix B 2, one can easily see that inequalities of the form (32) hold for higher $n$ too, where the RHS is always independent of $L$.

<footnote n="3">³ We note that a precise estimation of the von Neumann entanglement entropy, instead, only requires to reconstruct the target state up to an error which decays polynomially in the system size [132]. However, computing the von Neumann entanglement entropy for MPDOs is expected to be computationally hard in general.</footnote>

<footnote n="4">⁴ Because of the FDQC structure, the reduced density matrix over the region $A_2B_1$ is independent of the system size for $L \ge |A_2B_1| + 2\ell - 1$.</footnote>
## IV. FINITE-RANGE CORRELATED STATES

In this section we present and discuss the most general form of our protocols. First, in Sec. IV A we give the definition of the AFCs, which state that the n-th powers of the system density matrix effectively factorize over regions of size smaller than some length scales $\xi_n$. Our definition might appear arbitrary at first, but we show later that such AFCs hold in large classes of states (see Secs. IVC and IVD). Next, in Sec. IV B we describe a protocol for purity and PT moment estimation in states satisfying the AFCs. Such protocols take as an input the value of the length scales $\xi_n$ (assumed to be known) and the desired accuracy, yielding as an output the numerical estimations for the purity and PT moments. The protocols only require polynomially-many (in $L$) measurements and post-processing operations. In addition, we can rigorously prove that the estimation errors are smaller than the desired threshold with high probability. Finally, we discuss the generality of the AFCs, proving analytically that they hold for MPDOS (Sec. IV C) and presenting numerical evidence for their validity in Gibbs states of local Hamiltonians (Sec. IVD).

### A. The approximate factorization conditions

The strategies developed so far rely on exact factorization conditions on the system density matrix, which are related to the sharp light-cone structure of correlation functions in FDQC states [110, 133, 134]. Away from these ideal cases, Eqs. (14) and (26) do not hold, but it is natural to conjecture that they can be modified taking into account exponential corrections over scales governed by the system correlation lengths. This is done in this section, where we introduce a set of AFCs generalizing Eqs. (14) and (26), respectively.

Given a state $\rho$, we say that it satisfies the *purity AFC* if there exist $\alpha_2, \beta_2, k_c > 0$ such that, for any interval $I = A \cup B \cup C$ as in Fig. 3(a), with $|B| = k \ge k_c$, we have

$$ \left| \mathrm{Tr}(\rho_I^2)^{-1} \frac{\mathrm{Tr}_{AB}(\rho_{AB}^2) \mathrm{Tr}_{BC}(\rho_{BC}^2)}{\mathrm{Tr}_B(\rho_B^2)} - 1 \right| \le \alpha_2 e^{-\beta_2 |B|}. \quad (34) $$

Similarly, we say that $\rho$ satisfies the *PT-moment AFCs* if there exist $\alpha_n, \beta_n > 0$ and an integer $k_c$ such that

$$ |\tilde{p}_n[AB] - s_n[A_2B_1]| \le \alpha_n e^{-\beta_n(|A_2|+|B_1|)}, \quad (35) $$

for all partitions as in Fig. 3(b) with $|B_1| = |A_2| = k \ge k_c$. Here $\tilde{p}_n[AB]$ and $s_n[A_2B_1]$ are defined in Eqs. (25) and (28), respectively. Equations (34) and (35) obviously generalize (14) and (26), introducing exponential corrections over length scales $\xi_n = \beta_n^{-1}$.

While the definitions (34) and (35) might appear arbitrary at first, we show later that such AFCs hold in large classes of states (see Secs. IVC and IVD). This statement is very intuitive: on the one hand, we have shown

that they hold exactly for FDQC states; on the other hand, in the absence of topological order ground states are known to be well-approximated by the latter [112–115]. It is then natural to expect that the AFCs continue to hold as the temperature is increased or as the system is perturbed by incoherent noise, as these ingredients do not introduce long-range quantum correlations. In the next section, we assume that the AFCs hold, and describe a protocol for purity and PT moment estimation. We stress that we do not need any additional assumption on the state to be measured. For instance, we do not require that it can be represented efficiently by an MPDO.

Our protocols take as an input the values of $\alpha_n, \beta_n$, and the desired accuracy. For instance, the purity estimation protocol takes as an input $\alpha_2$, the correlation length $\xi_2 = \beta_2^{-1}$, and the target threshold relative error $\delta$. Importantly, it is not necessary to know the values of $\alpha_2$ and $\beta_2$ exactly: it is sufficient to have two estimates $\tilde{\alpha}_2, \tilde{\xi}_2$, which upper-bound them. This is because if $\alpha_2 \le \tilde{\alpha}_2$ and $\xi_2 \le \tilde{\xi}_2$, then Eq. (34) also holds replacing $\alpha_2$ and $\beta_2$ with $\tilde{\alpha}_2$ and $\tilde{\beta}_2 = \tilde{\xi}_2^{-1}$, respectively. We can choose $\tilde{\alpha}_2$ and $\tilde{\xi}_2$ arbitrarily, but the number of measurements and post-processing operations increase parametrically with $\tilde{\alpha}_2$ and $\tilde{\xi}_2$, see Eq. (37). Therefore, the more accurately we can estimate the values of $\alpha_2$ and $\beta_2$, i.e. the more information we have on the state to be measured, the more efficiently we can estimate the purity and, similarly, the PT moments. It is important to note that, in this respect, the situation is completely analogous to MPDO tomography. In that case, one assumes that the density matrix of the system $\rho$ can be described accurately by some MPO with a bond dimension $D$. Such bond dimension does not need to be known exactly, but only an upper bound for it needs to be known. In practice, one can take $\alpha_2$ and $\beta_2^{-1}$ larger than any expected length scale in the system.

Finally, suppose that, for given state $\rho$, we have an ansatz for $\alpha_2$ and $\xi_2$. From the experimental point of view, an important question is whether it is possible to efficiently certify that Eq. (34) holds, with the ansatz values $\alpha_2$ and $\xi_2$, for the state to be measured. While at the moment this is an open question, a simple consistency check consists in repeating the purity estimation protocol (explained in IV B) for increasing values of the input ansatz values $\alpha_2$ and $\xi_2$, and checking that the estimated value for the purity does not change, up to the expected precision. If the estimated purity does change, this is a “red flag” signalling the failure of the AFCs. We refer to Sec. VI for further discussions on the possibility to efficiently certify the validity of the AFCs.

### B. Estimation protocols for states satisfying the AFCs

First, we consider estimating the global purity of a state satisfying the AFCs (34). We assume the following:
* we know that the system is prepared in a state $\rho$ satisfying the AFCs (34);

* we know two numerical values $\alpha_2$ and $\beta_2$ for which Eq. (34) is satisfied.

Below, we detail our protocol to efficiently estimate the purity. The proof that the protocol works is non-trivial and is presented in detail in Appendix C. Intuitively, we use the fact that, because of the AFC (34), the estimator $r_2^{(e)}$ defined in Eq. (17) yields a good approximation for the purity, for values of $|I_j|$ which grow mildly (logarithmically) in $L$.

The protocol takes as an input the values of $\alpha_2$ and $\beta_2$, and consists of the following steps:

1. Choose the desired relative error $\delta$ on the purity. This can be any arbitrarily small number $\delta > 0$;

2. Take a partition of the system as in Fig. 1, where $|I_j| = k$ and choose

$$k \ge \xi_2 \ln(7\alpha_2 L/\delta), \quad (36)$$

where $\xi_2 = \beta_2^{-1}$. The reason for the functional form appearing in the RHS is technical and explained in Appendix C;

3. As in Sec. III A, perform $M_I = M$ classical-shadow measurements to estimate the purity of each interval $I$, with $I = I_j$ or $I = I_j \cup I_{j+1}$ [we use the set of measurements $M_I$ to only determine the purity over $I$, so that the total number of experimental runs is given by (16)].

4. From the classical shadows obtained out of the measurements, compute $r_2^{(e)}$ defined Eq. (17). This is the output of our protocol, giving us an estimate for the global purity.

As mentioned, it is intuitive that the output of this protocol $r_2^{(e)}$ is, with high probability, an accurate approximation for the purity. This is proven rigorously in Appendix C, where we show<fnref n="5" />

$$\mathrm{Pr} \left[ \left| (r_2^{(e)} / \mathcal{P}_2) - 1 \right| \ge \delta \right] \le \frac{7^2 2^{11} L^3}{\delta^2 M} \left( \frac{7 \alpha_2 L}{\delta} \right)^{4 \xi_2 \ln 2} . \quad (37)$$

Therefore, by a polynomial number of measurements $M$ we can guarantee that the probability that the relative error $|r_2^{(e)} / \mathcal{P}_2 - 1|$ is larger than $\delta$ is arbitrarily small. Rephrasing our result in terms of the confidence level $\gamma$, we thus obtain

$$M \ge \frac{7^2 2^{11} L^3}{\delta^2 (1 - \gamma)} \left( \frac{7 \alpha_2 L}{\delta} \right)^{4 \xi_2 \ln 2} . \quad (38)$$

This formalizes the anticipated result.

Note that $r_2^{(e)}$ defined above is a faithful estimator for the quantity

$$r_2 = \frac{\prod_{j=1}^{R-1} \mathrm{Tr}_{I_j \cup I_{j+1}} (\rho_{I_j \cup I_{j+1}}^2)}{\prod_{j=2}^{R-1} \mathrm{Tr}_{I_j} [\rho_{I_j}^2]}, \quad (39)$$

but because (15) does not hold exactly, it is not a faithful estimator of the global purity $\mathcal{P}_2$. Still, it is a good approximation, as it is clear from the relation

$$\left| \frac{r_2}{\mathcal{P}_2} - 1 \right| \le \frac{4\alpha_2 L}{k} e^{-\beta_2 k}, \quad (40)$$

which holds for

$$k \ge \xi_2 \ln(2\alpha_2 L), \quad (41)$$

and follows directly from (34), cf. Appendix C. In practice, the RHS of (40) makes it necessary to consider intervals $I_j$ whose length $k$ grows logarithmically in $L$ to keep the error below the desired threshold.

A similar protocol works for the PT moments. In this case, a number of measurements which does not scale with the system size is enough to guarantee arbitrary accuracy with high probability, cf. Appendix C. We note that the results of this section imply that the number of post-processing operations to estimate the purity and PT moments is also at most polynomial in $L$. This is because the estimators (9) and (11) are constructed summing a number of terms which is polynomial in $M$, and hence at most polynomial in $L$.

In the rest of this section, we discuss the generality of the AFCs. We prove analytically that they hold for the important class of translation-invariant MPDOs and provide numerical evidence that they are also satisfied by thermal states of local Hamiltonians. These results showcase the generality and versatility of our approach.

### C. Proof of AFCs in Matrix Product Density Operators

We now focus on MPDOs, a very general class of states providing accurate approximations for several short-range correlated density matrices, including thermal states of local Hamiltonians [60, 135]. We study the translation-invariant case, where we can provide analytic results. We recall that a translation-invariant MPDO $\sigma$ is defined by a four-index tensor $A_{a,b}^{j,k}$, with $j, k = 0, 1$ and $a, b = 0, 1 \dots \chi - 1$, where $\chi$ is its bond dimension [129]. For a system size $L$, its matrix elements read

$$\langle i_1, \dots, i_L | \sigma_L | j_1, \dots, j_L \rangle = A_{a_1, a_2}^{i_1, j_1} A_{a_2, a_3}^{i_2, j_2} \cdots A_{a_L, a_1}^{i_L, j_L}, \quad (42)$$

<footnote n="5">⁵ More precisely, Eq. (37) holds for $M$ sufficiently large, namely $M \ge \max \left\{ (7\alpha_2 L/\delta)^{8\xi_2 \ln 2}, (L^2 7^2 2^{10}/\delta^2) (7\alpha_2 L/\delta)^{4\xi_2 \ln 2} \right\}$, see Appendix C.</footnote>
where repeated indices are summed over. This defines a family of normalized states<fnref n="6" />

$$ \rho_L = \frac{\sigma_L}{\mathrm{Tr}[\sigma_L]} . \qquad (43) $$

Using the standard notation from tensor-network theory [129], we can represent Eq. (42) as

$$ \sigma_L = \cdots \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \node[circle, fill, inner sep=1.2pt] (0) at (0,0) {}; \node[circle, fill, inner sep=1.2pt] (1) at (1,0) {}; \node[circle, fill, inner sep=1.2pt] (2) at (2,0) {}; \node[circle, fill, inner sep=1.2pt] (3) at (3,0) {}; \node at (0.5,0) {$\cdots$}; \node at (1.5,0) {$\cdots$}; \node at (2.5,0) {$\cdots$}; \node at (3.5,0) {$\cdots$}; \end{tikzpicture} \cdots \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \node[circle, fill, inner sep=1.2pt] (0) at (0,0) {}; \node[circle, fill, inner sep=1.2pt] (1) at (1,0) {}; \node[circle, fill, inner sep=1.2pt] (2) at (2,0) {}; \node[circle, fill, inner sep=1.2pt] (3) at (3,0) {}; \node at (0.5,0) {$\cdots$}; \node at (1.5,0) {$\cdots$}; \node at (2.5,0) {$\cdots$}; \node at (3.5,0) {$\cdots$}; \end{tikzpicture} \cdots \quad (44) $$

Here each circle correspond to a tensor $A_{a,b}^{j,k}$. Lower and upper legs are associated with input and output degrees of freedom, respectively, while the remaining ones correspond to the “virtual” indices $a_j = 0, \dots, \chi - 1$. Note that contracted legs indicate pairs of indices which are summed over.

We first focus on the purity. We are able to show that MPDOs satisfy the purity AFCs under a few technical assumptions, which encode the fact that all correlation lengths are finite but are otherwise very general. Technically, we impose a few conditions on the spectrum and eigenstates of the transfer matrices<fnref n="7" />

$$ \tau_n = n \left\{ \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \node[circle, fill, inner sep=1.2pt] (0) at (0,0) {}; \node[circle, fill, inner sep=1.2pt] (1) at (1,0) {}; \node[circle, fill, inner sep=1.2pt] (2) at (2,0) {}; \node[circle, fill, inner sep=1.2pt] (3) at (3,0) {}; \node at (0.5,0) {$\cdots$}; \node at (1.5,0) {$\cdots$}; \node at (2.5,0) {$\cdots$}; \node at (3.5,0) {$\cdots$}; \end{tikzpicture}, \right. \quad (45) $$

where the top and bottom legs are contracted. Informally, we ask that $\tau_1$ and $\tau_2$ have non-degenerate eigenvalues with largest absolute value, and that the corresponding eigenstates are not orthogonal. These conditions are general in the sense that one needs to choose fine-tuned examples to violate them. We refer the reader to Appendix D for a precise statement and a detailed discussion.

Under these assumptions, we show that $\rho_L$ satisfies the purity AFCs. Note that, different from the rest of this work, here we are considering periodic boundary conditions. Since Eq. (34) was introduced for open boundary conditions, we consider directly its global version (40), which can immediately be generalized to the periodic case. Clearly, the expression for $r_2$ must be modified with respect to (39): taking into account periodic boundary conditions, we replace it with

$$ r_2 = \frac{\prod_{j=1}^{R} \mathrm{Tr}_{I_j \cup I_{j+1}} (\sigma_{I_j \cup I_{j+1}}^2)}{\prod_{j=1}^{R} \mathrm{Tr}_{I_j} [\sigma_{I_j}^2]}, \quad (46) $$

where we assume without loss of generality that $R = L/k$ is an integer, while the normalization factors $\mathrm{Tr}[\sigma^2]$ cancel out.

The task of proving Eq. (40) [using the definition (46)] is carried out in Appendix D. There, we derive the following statement: given the family of MPDOs $\rho_L$ as above, there exist $\zeta > 0$, and $C > 0$ (independent of $L$) such that, for any integer $k \ge k_{\min}$, with

$$ k_{\min} = \max \{1, \zeta \ln(20C\chi^2L)\}, \quad (47) $$

we have

$$ \left| \frac{r_2}{P_2} - 1 \right| \le \chi^2(80C + 32) \frac{L}{k} e^{-k/\zeta}, \quad (48) $$

for all $L \ge \max\{\zeta \ln(2^5\chi^2), 4k + \zeta \ln(2\chi)\}$. The correlation length $\zeta$ and the constant $C$ depend on both the eigenvalues and eigenvectors of the transfer matrices $\tau_1$ and $\tau_2$, cf. Appendix D. Therefore, the purity approximate factorization formula (40) holds with $\alpha_2 = \chi^2(20C + 8)$, and $\beta_2 = 1/\zeta$. In turn, using the results of the previous sections, this implies that the purity of MPDOs can be estimated efficiently using our approach.<fnref n="8" />

Similarly, the AFCs for PT moments could be verified analytically for bulk-translation-invariant MPDOs with suitable open boundary conditions. This is slightly more involved, as a few additional technical assumptions on the boundaries are needed. For this reason, we do not discuss this explicitly. Instead, the AFCs for PT moments will be numerically analyzed in detail in the next section, together with the purity AFCs, for the physically interesting case of Gibbs states of local Hamiltonians in a few concrete examples.

### D. Numerical study of AFCs in Gibbs states

As mentioned, it is natural to conjecture that the AFCs are very general, as they encode the fact that all the spatial correlation lengths in the system are finite. In this section, we provide numerical evidence supporting this claim, studying thermal states of two prototypical 1D local Hamiltonians. We focus in particular on the Ising model with transverse and longitudinal fields

$$ H_I = -\frac{1}{4} \sum_{j=1}^{L-1} [\sigma_j^z \sigma_{j+1}^z + h_x \sigma_j^x + h_z \sigma_j^z], \quad (49) $$

and the so-called XXZ Heisenberg spin chain

$$ H_{XXZ} = -\frac{1}{4} \sum_{j=1}^{L-1} [\sigma_j^x \sigma_{j+1}^x + \sigma_j^y \sigma_{j+1}^y + \Delta \sigma_j^z \sigma_{j+1}^z]. \quad (50) $$

<footnote n="6">⁶ We assume that the local tensors $A$ generate a positive operator for all system sizes $L$, namely $\sigma_L \ge 0$ for all $L$.</footnote>

<footnote n="7">⁷ Note that $\tau_n$ is not Hermitian and that, for $n = 1$, we recover the standard transfer matrix [135].</footnote>

<footnote n="8">⁸ Note that, with these identifications, Eq. (36) implies $k \ge k_{\min}$, where $k_{\min}$ is defined in Eq. (47).</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4. Relative error $\epsilon_r(k, L, \beta)$, defined in Eq. (52), as a function of $k$, for $\beta = 2$ and different values of $L$. In the left and right plots we report numerical data corresponding to the Ising and XXZ spin chains, respectively. In the former case, we chose $h_x = 1.1$, $h_z = -0.04$, while the maximal bond dimension used in the computation is $\chi = 32$. In the latter case, we set $\Delta = 2$, while $\chi = 64$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Numerical results for the minimal region size $k^*(\delta, L)$ required to achieve a specified precision $\delta$ of purity estimation, see Eq. (53). In the plots we chose $\delta = 0.01$ and $\beta = 2$, while the left and right figure correspond to the Ising and XXZ spin chains, respectively. The Hamiltonian parameters and bond dimensions are the same as in Fig. 4.</figcaption></figure>

We will study the corresponding Gibbs states

$$ \rho_{\beta} = \frac{\exp(-\beta H)}{Z_{\beta}}, \qquad (51) $$

where $Z_\beta$ is a normalization factor.

Note that the Hamiltonians $H_I$ and $H_{XXZ}$ are integrable for $h_z = 0$ [136] and all values of $\Delta$ [137], respectively, although we do not expect that integrability plays any role in the following discussion. Note also that both Hamiltonians feature second-order quantum phase transitions at zero temperature: For $H_I$ the critical line is $h_z = 1$, while $H_{XXZ}$ displays a critical phase for $|\Delta| < 1$. However, we will be interested in finite-temperature states, and since in 1D the temperature always introduces a finite length scale, we similarly expect that quantum criticality does not play any role.<fnref n="9" />

In order to assess the validity of the AFC for the purity, we test its global version (40). To this end, we numerically compute $r_2(k)$, as defined in Eq. (39), for increasing values of the interval size $k = |I_j|$, and the global purity $P_2$. Then, we evaluate the corresponding relative error

$$ \epsilon_r(k, L, \beta) = \left| \frac{r_2(k)}{P_2} - 1 \right|, \qquad (52) $$

where we have made the dependence on $k$, $\beta$, and $L$ explicit. The calculations are carried out using the iTensor

<footnote n="9">⁹ An interesting question is whether our approach could be ex-</footnote>

<footnote>tended to pure critical states. Although in this case there is no spatial length scale, it is known that they can be accurately approximated by MPSs with bond-dimension scaling polynomially in the system size [138], suggesting that our approach could be generalized to include that case as well. We leave this problem for future work.</footnote>
library [139], by first approximating the thermal states by MPOs of finite bond dimension, denoted by $\chi$, and subsequently taking powers and traces of such MPOs. For each quantity, this procedure is repeated for increasing $\chi$, until convergence is reached. We refer to Appendix E for further details.

For fixed values of $L$, we first test the exponential dependence in $k$ of $\epsilon_r(k, L, \beta)$ predicted by (40). Examples of our numerical data are reported in Fig. 4, showing a very clear exponential decay. For the Ising model, we chose a non-zero value of the longitudinal field, in order to break integrability. Note that, for very large $k$, we see an apparent plateau. We attribute this behavior to finite numerical precision of our computations (note the very small values of $\epsilon_r$ corresponding to these plateaus). It is interesting to note that the decay rate for the XXZ chain is faster than the Ising model, although higher bond dimensions are needed in order to accurately approximate the corresponding Gibbs state by an MPO. In any case, we see that quite small values of $k$ are enough to obtain a very good approximation of the purity. We have repeated the calculations for different values of the Hamiltonian parameters and $\beta$, finding consistent results.

Next, we study the functional form of $\epsilon_r(k, L, \beta)$. From the numerical point of view, it is not straightforward to test that it is asymptotically bounded by the RHS of Eq. (40). Therefore, we proceed differently. Given an arbitrary fixed $\delta > 0$, we define the minimal region size required to achieve a specified precision $\delta$,

$$k^*(\delta, L) = \min \{k : \epsilon_r(k, L, \beta) < \delta\}. \quad (53)$$

Assuming that $\epsilon_r$ satisfies an asymptotic bound of the form (40), it is easy to show that $k^*(\delta, L)$ grows at most logarithmically in $L$, namely there exist $\alpha_2$ and $\beta_2$ such that

$$k^*(\delta, L) \le \frac{1}{\beta_2} \ln(\alpha_2 L / \delta). \quad (54)$$

When this condition holds, the estimation protocol explained in Sec. IV B is efficient. Importantly, compared to Eq. (35), Eq. (54) is straightforward to test by our numerical computations.

We have verified that Eq. (54) is satisfied for several values of the temperature and the Hamiltonian parameters. An example of our numerical results is given in Fig. 5. For the Ising model, we see a very clear logarithmic scaling. For the Heisenberg chain, higher bond dimensions are required to achieve the same accuracy. Accordingly, we were not able to simulate sufficiently large system sizes to obtain a clear logarithmic behavior. Still, our results are always consistent with a growth of $k^*(\delta, L)$ which is at most logarithmic in $L$, cf. Fig. 5. Thus, our numerical data fully support the validity of Eq. (54). In fact, note that the values of $k^*(\delta, L)$ are found to be very small in practice. This is true even when the required bond dimension for MPO calculations is large, as in the case of the XXZ Heisenberg chain. This might make our

approach useful already for relatively small system sizes, even if compared against MPO tomographic methods.

Finally, we study the AFC for the PT moments (34). To this end, we define the additive error

$$\epsilon_a(k, L, \beta) = |\tilde{p}_n[AB] - s_n[A_2B_1]|, \quad (55)$$

where we have made the dependence on $k$, $\beta$, and $L$ explicit. Because the numerical cost of evaluating the trace of the $n$-th power of an MPO increases with $n$, we have restricted our analyses to $n = 3$.

In Fig. 6 we report an example of our data for $\epsilon_a(k, L, \beta)$ as a function of $k$, for increasing system sizes. From the plots, two things are apparent. First, for fixed $L$, the additive error displays a clear exponential decay. Second, its values are independent of $L$ (for values of $k$ much smaller than the system size). This behavior is perfectly consistent with Eq. (35). Note that, as in the previous plots, we observe irregular behavior for sufficiently large $k$, where $\epsilon_a$ approaches a plateau and displays a $L$-dependence. We interpret these effects as arising due to numerical inaccuracies (note the very small values of $\epsilon_a$). Similarly, when $k$ approaches $L/2$ we observe $L$-dependent behavior, which can be attributed to finite size effects. We have repeated our calculations for different $\beta$ and Hamiltonian parameters (in each case, we have checked that the bond dimension was large enough), always finding consistent results.

## E. Practical feasibility of purity estimation

In Sec. IV B, we have shown that, for states satisfying the AFCs, the purity and the PT moments may be estimated by a number of measurements growing only polynomially in the system size. Still, the degrees of the polynomials in our bounds are high, cf. for example Eq. (33). Accordingly, one may wonder whether our approach could be applied in practice in experiments with current repetition rates. It is important to note, however, that our analysis was mainly concerned with establishing the polynomial dependence rigorously, so that our bounds may likely be improved by more refined mathematical derivations. More importantly, as we have stressed, our protocol may be improved by choosing different ways of estimating the local purities and PT moments, a simple variant being the method of “common randomized measurements” [48]. These improvements are expected to lead to polynomial bounds with lower exponents.

Going beyond general bounds, it is important to illustrate how our estimation protocol performs in practice. This is done in this section, where we exemplify its feasibility by presenting a classical simulation of the full protocol in an explicit case. For simplicity, we focus on random MPSs [60, 135] of $N$ qubits and consider the purity of the reduced density matrix over half system. The choice of pure MPSs is motivated by the fact that it makes it easier to classically simulate the probability
<figure><img src="image_7.png" /><figcaption>Figure 6. Additive error $\epsilon_a(k, L, \beta)$, defined in Eq. (55), as a function of $k$, for $\beta = 2$ and different values of $L$. In the left and right plots we report numerical data corresponding to the Ising and XXZ spin chains, respectively. The Hamiltonian parameters and bond dimensions are chosen as in Fig. 4.</figcaption></figure>

distribution of the measurement outcomes. At the same time, random states model typical physical states without allowing for fine-tuning. Our simulations show that the number of measurements needed to reconstruct the purities is order of magnitudes smaller compared to our theoretical bounds, demonstrating that our protocol is experimentally feasible in current quantum platforms.

The details of our simulation are as follows. We initialize the system in a random MPS with complex coefficients and bond dimension $\chi$ [139]. Once the random MPS is constructed, we save a copy of it and simulate the RM process (in particular, we imagine that we have experimental access to many identical copies of the state). In our simulation, we consider a slightly more general measurement scheme compared to the one detailed in Sec. IIB. In particular, we make use of the RMs described in Ref. [29] where each random unitary $u_j$ is re-used for $n_M$ local measurements. Therefore, denoting by $n_U$ the number of distinct unitaries, the total number of experiments is $n_U n_M$ (so that the classical-shadow approach described in Ref. IIB is a special case with $n_M = 1, n_U = M$). Each simulated experimental run produces a bit-string $\mathbf{s}^{(m,k)} = (s_1^{m,k}, \dots, s_L^{m,k})$, with $m = 1, \dots, n_U$, $k = 1, \dots, n_M$. Note that, numerically, we can obtain the correct probability distribution for the measurement outcomes using the perfect-sampling algorithm for MPSs [140]. We have chosen to simulate the RM scheme with re-use of local unitaries as this scheme is more efficient to simulate and it has been used in a number of recent experiments [27, 30–33]

Finally, we collect the simulated measurement outcomes to construct the purity estimator (17). Instead of Eq. (9), we use the following estimator for the local

purities

$$ \mathcal{P}_2^{(e)}[I] = \frac{2^N}{n_U n_M (n_M - 1)} \sum_{m=1}^{n_U} \sum_{\substack{k, k'=1 \\ k \neq k'}}^{n_M} (-2)^{-D[\mathbf{s}_I^{(m,k)}, \mathbf{s}_I^{(m,k')}]} \quad (56) $$

where $D[\mathbf{s}_I^{(m,k)}, \mathbf{s}_I^{(m,k')}]$ is the Hamming distance between the bit-stings $\mathbf{s}_I^{(m,k)}$ and $\mathbf{s}_I^{(m,k')}$ on subsystem I. Compared to (9), Eq. (56) achieves a more accurate estimate for many repetitions $n_M$ using a few set of distinct unitaries $u_M$ and is expected to be more robust against miscalibration-errors [12]. We also checked numerically that our results are qualitatively unchanged when using estimators based on classical shadows.

Using the above simulation procedure, we have studied the estimator $r_2^{(e)}$ defined in Eq. (29) as a function of the number of measurements $n_U n_M$. We computed in particular the statistical error $\epsilon_s = r_2^{(e)}/r_2$, where $r_2$ is the r.h.s. of Eq. (15) and the full relative error $\epsilon = \epsilon_r$ defined in Eq. (18). In Fig. 7, we report our data for $\chi = 2, 3, k = 4, 6$ and system sizes up to $L = 96$ (each plot corresponds to an average error over 20 instances of the random MPS).

From the plots, we see that roughly $10^4$ measurements are enough to obtain a relative error below 10% for system sizes up to $L \sim 100$. $^{10}$ In addition, roughly $10^5$ measurements are needed to reach an error around 5% for system sizes up to $L \sim 50$ (note that, for these sizes, a direct application of RM schemes to estimate global purities is unfeasible). As announced, these numbers are order of magnitudes smaller than what predicted by

<footnote>$^{10}$ Note that in this range of the parameters, statistical and full errors are seen to be numerically close, signalling that the deviations from exact AFCs are small for $k = 6$, both for $\chi = 2$ and $\chi = 3$.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 7. Statistical and full relative error on the bipartite purity in random MPSs, as a function of the number of experimental runs $n_{UNM}$, with fixed $n_M = 1000$. For each plot, we report the bond dimension $\chi$ of the random MPS and the length of the intervals $|I_j| = k$ in which the system is partitioned in our estimation protocol.</figcaption></figure>

our general bounds, suggesting that these bounds are far from being tight.

We stress that the number of measurements performed in our simulations are of the order of those performed in current trapped-ion and superconducting qubit experiments, see e.g. Refs. [27, 30–33, 39]. The results reported in this section thus substantiate the claim that our protocol is practically useful for global purity and entanglement estimation in today’s quantum platforms.

## V. MIXED-STATE ENTANGLEMENT DETECTION

In this section we study the $p_n$-PPT conditions for states satisfying the AFCs. In Sec. V A we first introduce a set of quantities, which we call $f_n(\rho)$, displaying the following properties: (i) $f_n(\rho) < 0$ if and only if $\rho$ violates the $p_n$-PPT conditions and (ii) $f_n(\rho)$ can be estimated efficiently if $\rho$ satisfies the AFCs, representing natural probes for mixed-state entanglement. Next, in Sec. V B, we compute $f_n$ in a few concrete cases including Gibbs states of local Hamiltonians, and identify the values of the system parameters for which $f_2, f_3 < 0$. The results of this section show that our approach can be successful in detecting entanglement even for large numbers of qubits and highly mixed states, thus being practically useful in many situations routinely encountered in available quantum platforms.

## A. The $p_n$-PPT conditions and mixed-state entanglement probes

Given the density matrix $\rho_{AB}$ on a bipartite system, its logarithmic negativity (3) can be extracted from the knowledge of all PT moments $p_n$, with $n = 1, \dots, \dim(\mathcal{H}_A \otimes \mathcal{H}_B)$. Therefore, while our method is efficient to estimate the individual PT moments, reconstructing the logarithmic negativity is still hard, as it requires an exponential number of them.

Luckily, the $p_n$-PPT conditions introduced in Sec. II A allow us to detect mixed-state entanglement from only a few PT moments. For example, the first two non-trivial elements in the family (5) only involve the PT moments up to $n = 5$. Explicitly, they read

$$p_3-\text{PPT}: \quad p_3p_1 \ge p_2^2, \qquad (57a)$$

$$p_5-\text{PPT}: \quad p_3p_5 \ge p_4^2. \qquad (57b)$$

The conditions (57) can be equivalently rewritten in terms of the quantities

$$f_3(\rho_{AB}) = \tilde{p}_3\tilde{p}_1 - \tilde{p}_2^2 \frac{\mathcal{P}_2^2[A]\mathcal{P}_2^2[B]}{\mathcal{P}_1[A]\mathcal{P}_1[B]\mathcal{P}_3[A]\mathcal{P}_3[B]}, \qquad (58a)$$

$$f_5(\rho_{AB}) = \tilde{p}_5\tilde{p}_3 - \tilde{p}_4^2 \frac{\mathcal{P}_4^2[A]\mathcal{P}_4^2[B]}{\mathcal{P}_3[A]\mathcal{P}_3[B]\mathcal{P}_5[A]\mathcal{P}_5[B]}, \qquad (58b)$$

where $\tilde{p}_n$ are the normalized PT moments defined in Eq. (25). Indeed $f_3$ and $f_5$ are negative if and only if the inequalities (57) are violated.

Crucially, $f_3(\rho_{AB})$ and $f_5(\rho_{AB})$ can be estimated efficiently if $\rho_{AB}$ satisfies the AFCs. Indeed, they only depend on $\mathcal{P}_n$ and the normalized PT moments $\tilde{p}_n$, both of
which can be estimated using the protocol introduced in Sec. IV.<fnref n="11" /> We note that such estimation protocols allow us to reconstruct $f_n(\rho_{AB})$ up to any arbitrary small additive error. Namely, for any arbitrary positive constant $\delta$, one can construct an estimator $f_n^{(e)}$, with $|f_n^{(e)} - f_n| \le \delta$ by polynomially many (in $L$) measurements and post-processing operations. This is easily seen recalling that $\tilde{p}_n(\mathcal{P}_n)$ can be estimated up to any additive (relative) error, and using that

$$|\tilde{p}_n[I]| \le 1, \quad \frac{(\mathcal{P}_n[I])^2}{\mathcal{P}_{n-1}[I]\mathcal{P}_{n+1}[I]} \le 1, \quad (59)$$

where the first inequality from $|\tilde{p}_n[I]| \le \tilde{p}_2[I] = \mathcal{P}_2[I]$, while the second follows from Hölder's inequality.

At this point, it is important to ask whether the quantities $f_n$ are practically useful in the many-body context considered in this work. For example, one could be worried that the conditions (57) are never violated for large $L$ (even if the logarithmic negativity is non-zero), or that $f_n$ is exponentially small in $L$ for highly mixed states, thus requiring an exponential accuracy for its estimation. To address this question, consider a family of states $\rho_{AB}^{(L)}$, defined for increasing system sizes $L$ and satisfying the AFCs. Suppose that we extract $f_n(\rho_{AB}^{(L)})$ using our estimation protocol for a fixed approximation error $\delta_n$. Then, the $p_n$-PPT conditions (57) allow us to successfully detect entanglement if

$$f_3(\rho_{AB}^{(L)}) \le -C_3, \quad (60a)$$

$$f_5(\rho_{AB}^{(L)}) \le -C_5, \quad (60b)$$

where $C_3, C_5 > 0$ are constants (independent of $L$) with $C_n \ge 2\delta_n$. Indeed, in these cases the statistical inaccuracy is smaller (with high probability) than the amount by which $f_n(\rho_{AB}^{(L)})$ is negative, allowing us to certify entanglement.

In the next subsection, we study a few natural examples of states satisfying the AFCs showing that (60) hold either for all values of $L$ or up to very large system sizes. Therefore, we exhibit concrete examples where $f_n$ are provably useful for entanglement detection.

## B. Mixed-state entanglement detection in FDQC and Gibbs states

We start by studying the relations (60) in the simplest case of FDQC states (13). We consider a family of states

where the single-qubit density matrices in Eq. (13) are parametrized as

$$\sigma_i(\gamma) = \gamma |a_i\rangle\langle a_i| + (1-\gamma) |b_i\rangle\langle b_i|. \quad (61)$$

Here $|a_i\rangle, |b_i\rangle$ are arbitrary orthonormal states, while $0 \le \gamma \le 1/2$ plays the role of a depolarizing parameter, controlling both purity and bipartite entanglement $[\gamma$ is not to be confused with the symbol previously used for confidence levels]. We ask for which values of $\gamma$ the conditions (57) are violated.

This problem is analyzed in detail in Appendix F. We prove that, for generic choices of the unitary gates, there exist $\gamma_3, \gamma_5 > 0$ such that Eqs. (60a) and (60b) are satisfied, respectively, for all $\gamma \le \gamma_3/L$ and $\gamma \le \gamma_5/L^{1/3}$. Now, choosing in particular $\gamma_L = \gamma_5/L^{1/3}$, we have

$$\begin{aligned} \mathrm{Tr}[\rho^{(L)}(\gamma_L)^2] &= [\gamma_L^2 + (1 - \gamma_L)^2]^L \\ &\le \exp[-2(\gamma_L - \gamma_L^2)L] \\ &\le \exp[-\gamma_L L] = \exp[-\gamma_5 L^{2/3}], \end{aligned} \quad (62)$$

namely

$$S_2[\rho^{(L)}(\gamma_L)] \ge \gamma_5 L^{2/3}. \quad (63)$$

Therefore, $\rho^{(L)}(\gamma_L)$ provides an example where the $p_5$-PPT condition is violated by a finite amount $C_5$ and, at the same time, the global Rényi-2 entropy grows polynomially, albeit sub-linearly, in the system size.<fnref n="12" />

As a final interesting example, let us consider once again finite-temperature states of local Hamiltonians. In this case, the purity decays exponentially in $L$ and we expect that the conditions (58) fail to detect the presence of entanglement in the thermodynamic limit. In the following we give evidence that they can nevertheless be useful when considering systems of large but finite size.

To support this claim, we study thermal states in the Ising model (49), and numerically compute the coefficients (60) for different temperatures and system sizes. The computations are performed using tensor-network methods, following the same strategy explained in Sec. IV D. An example of our results is reported in Fig. 8. In these plots, we have fixed small but otherwise arbitrary positive numbers $C_3$ and $C_5$, and identified the values of $\beta$ and $L$ for which $f_3(\beta, L) \le -C_3$ and $f_5(\beta, L) \le -C_5$.

The data corresponding to the $p_3$-PPT condition are reported in the left plot, which clearly shows that the values of the temperature $T = \beta^{-1}$ for which the entanglement can be detected decrease with $L$, consistent

<footnote n="11">¹¹ Note that, while in Sec. IV B we only considered the purity $\mathcal{P}_2$, similar AFCs can be defined for $\mathcal{P}_n$ with $n > 2$, and hence similar estimation protocols work to estimate higher moments. Note also that, for the case of FDQC, the proof of exact factorization conditions for $\mathcal{P}_2$ can be trivially extended to any $\mathcal{P}_n$ with $n > 2$, cf. Appendix A.</footnote>

<footnote n="12">¹² It is interesting to note that the $p_3$-PPT condition only allows us to detect bipartite entanglement in states with $O(1)$ entropy. Therefore, the $p_5$-PPT condition improves this result by allowing the entropy to scale as $O(L^{2/3})$. It is then natural to conjecture that higher-$n$ PPT conditions lead to a further improvement of the maximal scaling to $O(L^\alpha)$, with $\alpha$ approaching 1 in the large-$n$ limit.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 8. Density plots for the coefficients $f_3(\beta, L)$, $f_5(\beta, L)$, defined in Eq. (58) as a function of $\beta$ and $L$. The data correspond to thermal states in the Ising model (49), with $h_x = 1.1$ and $h_z = -0.04$. Red (blue) regions indicate the values of $\beta$ and $L$ for which $f_3(\beta, L)$, $f_5(\beta, L)$ are larger (smaller) than some small values $-C_3 = 0.01$ and $-C_5 = 0.001$, respectively.</figcaption></figure>

with our previous analysis on FDQC states. On the other hand, the right plot shows that $f_5(\beta, L) < -C_5$ up to values of the temperature ($T \approx 0.1$) that do not significantly depend on $L$, at least up to very large system sizes ($L \approx 1024$). We expect that this is a finite-size effect, namely that increasing $L$ further, the values of $T$ for which entanglement can be detected will decrease. Still, this example shows how the $p_n$-PPT conditions may be able to detect entanglement even for highly mixed states and up to very large system sizes, making the estimation protocol for PT moments practically very useful.

## VI. OUTLOOK

We have studied the problem of estimating global entropies and entanglement in many-qubit states, assuming the validity of certain AFCs which encode the fact that all the spatial correlations lengths in the system are finite. We have shown that the AFCs allow one to reconstruct entropies and PT moments from local information which can be extracted efficiently using standard RM schemes. Exploiting this fact, we have devised a very simple strategy for entropy and entanglement estimation which is provably accurate, requiring polynomially-many local measurements and post-processing operations. We have discussed the generality of the AFCs, providing both analytic and numerical evidence for their validity in different classes of states, including MPDOs and thermal states of local Hamiltonians.

As we have argued, the protocols proposed in this work could be practically useful to detect bipartite mixed-state entanglement for large numbers of qubits in today's quantum platforms. This is especially true when considering

current NISQ devices from the point of view of quantum simulation, since the AFCs hold under assumptions that are very common in the context of many-body physics. In addition, our work also raises important questions and opens up a number of interesting directions.

First, our protocols assume that the state to be measured satisfies the AFCs. Therefore, an important question from the experimental perspective is whether it is possible to devise a complementary (and still efficient) strategy to certify their validity, similar to what can be done in the case of MPS [61] or entanglement-Hamiltonian tomography [69].

As we have discussed in Sec. IV A, a simple but effective consistency check can be carried out by repeating our protocols multiple times, each time with an increased ansatz value for the correlation lengths. If the estimated purities or PT moments change beyond the expected precision over different repetitions, we obtain a “red flag” signalling the failure of the AFCs in the state to be measured. This process, however, does not allow one to rigorously rule out the presence of long-range correlations, as this would require running the estimation protocols using, as an input, correlation lengths which are proportional to the system size. We envision that a way to tackle the certification problem could be to identify a hierarchy of properties which are strictly stronger than the AFCs and for which rigorous certification protocols can be found more easily. At the moment, however, this remains an open problem.

Next, while we have focused on 1D systems, a very natural direction is to extend our approach to higher spatial dimensions. In this case, additional technical complications arise, but we expect that efficient estimation strategies for the purity and PT moments exists under similar
$$ \mathrm{Tr}_{B^{\otimes n}}[\vec{\Pi}_B (\rho_{ABC}^{\otimes n})] = T_1 T_2 T_3, \text{ where} $$

<figure><img src="image_10.png" /></figure>

Note that $T_1$ and $T_3$ are operators, while $T_2$ is a number. Let us compare this graphical expression with the one in the RHS of Eq. (A1). The first and second terms in the numerator read, respectively

<figure><img src="image_11.png" /></figure>

Similarly, the denominator in Eq. (A1) reads

<figure><img src="image_12.png" /></figure>

Putting everything together, and simplifying the common factors in the numerator and denominator, we arrive at Eq. (A1). Finally, using the fact that $\mathrm{Tr}(\rho_X^n) = \mathrm{Tr}(\vec{\Pi}_X \rho_X^{\otimes n})$, and the relation $\vec{\Pi}_{ABC} = \vec{\Pi}_A \otimes \vec{\Pi}_B \otimes \vec{\Pi}_C$, we arrive at

$$ \mathrm{Tr}(\rho_{ABC}^n) = \frac{\mathrm{Tr}_{AB}(\rho_{AB}^n) \mathrm{Tr}_{BC}(\rho_{BC}^n)}{\mathrm{Tr}_B(\rho_B^n)}. \quad (A4) $$

For $n=2$, this proves Eq. (14).

Next, we prove Eq. (15). Let us consider a partition of the whole system $S = A \cup B \cup C$, where $|B| = |C| = k \ge 2\ell - 1$. We first make use of Eq. (14). Setting now $A^{(0)} := A$, $B^{(0)} := B$, $C^{(0)} := C$ we can iterate this operation, each time splitting the interval $A^{(j)} \cup B^{(j)}$ into three adjacent regions $A^{(j+1)}$, $B^{(j+1)}$, and $C^{(j+1)}$ with $|B^{(j+1)}| = |C^{(j+1)}| = k$. We continue this procedure until $|A^R| = k$, which gives us Eq. (15).

## 2. The PT moments

The goal of this Appendix is to prove Eq. (26).

We start by recalling the expression for the PT moments [102]

$$ \mathrm{Tr} \left[ \left( \rho_{AB}^{T_A} \right)^n \right] = \mathrm{Tr} \left[ \overleftarrow{\Pi}_A \overrightarrow{\Pi}_B (\rho_{AB}^{\otimes n}) \right], \quad (A5) $$

where we introduced the $n$-copy forward (backward) cyclic permutation operator in the $n$-replica space, defined by

$$ \begin{aligned} \overrightarrow{\Pi}_A |i_1, \dots, i_n\rangle &= |i_n, i_1, \dots, i_{n-1}\rangle, \\ \overleftarrow{\Pi}_B |i_1, \dots, i_n\rangle &= |i_2, i_3, \dots, i_1\rangle. \end{aligned} \quad (A6) $$

Using this notation, it is easy to see that Eq. (A1) implies

$$ \mathrm{Tr}_{B^{\otimes n}} \left[ \overleftarrow{\Pi}_B (\rho_{ABC}^{\otimes n}) \right] = \frac{\mathrm{Tr}_{B^{\otimes n}} \left[ \overleftarrow{\Pi}_B (\rho_{AB}^{\otimes n}) \right] \otimes \mathrm{Tr}_{B^{\otimes n}} \left[ \overrightarrow{\Pi}_B (\rho_{BC}^{\otimes n}) \right]}{\mathrm{Tr}_B(\rho_B^n)}, \quad (A7) $$

Eq. (A7) can be proved by taking the partial-transpose with respect to $A^{\otimes n}$ and $C^{\otimes n}$ in both sides of (A1), followed by complex conjugation. Next, let us take a bipartition of the system into A and B, each divided into two intervals $A = A_1 \cup A_2$ and $B = B_1 \cup B_2$ with $|A_2| = |B_1| = k$, cf. Fig. 3(b). Making use of Eq. (A1), we have

$$ \begin{aligned} \mathrm{Tr} \left[ (\rho_{AB}^{T_A})^n \right] &= \mathrm{Tr}_{A^{\otimes n}} \mathrm{Tr}_{B_2^{\otimes n}} \overleftarrow{\Pi}_A \overrightarrow{\Pi}_B \mathrm{Tr}_{B_1^{\otimes n}} \overrightarrow{\Pi}_B \rho^{\otimes n} \\ &= \frac{\mathrm{Tr}_{A^{\otimes n} \otimes B_1^{\otimes n}} (\overleftarrow{\Pi}_A \overrightarrow{\Pi}_B \rho_{AB_1}^{\otimes n})}{\mathrm{Tr}_{B_1}(\rho_{B_1}^n)} \\ &\quad \times \frac{\mathrm{Tr}_{B_1^{\otimes n} \otimes B_2^{\otimes n}} (\overrightarrow{\Pi}_{B_1} \overleftarrow{\Pi}_{B_2} \rho_{B_1 B_2}^{\otimes n})}{\mathrm{Tr}_{B_1}(\rho_{B_1}^n)}. \end{aligned} \quad (A8) $$
Using

$$
\begin{align*}
& \mathrm{Tr}_{A^{\otimes n} \otimes B_1^{\otimes n}} (\overleftarrow{\Pi}_A \overleftarrow{\Pi}_{B_1} \rho_{AB_1}^{\otimes n}) \\
&= \mathrm{Tr}_{A_1^{\otimes n}} \mathrm{Tr}_{B_1^{\otimes n}} \overleftarrow{\Pi}_{A_1} \overleftarrow{\Pi}_{B_1} \mathrm{Tr}_{A_2^{\otimes n}} \overleftarrow{\Pi}_{A_2} \rho_{A_1 A_2}^{\otimes n} \\
&= \mathrm{Tr}_{A_1^{\otimes n} \otimes A_2^{\otimes n}} (\overleftarrow{\Pi}_{A_1} \overleftarrow{\Pi}_{A_2} \rho_{A_1 A_2}^{\otimes n}) \\
& \mathrm{Tr}_{A_2^{\otimes n} \otimes B_1^{\otimes n}} (\overleftarrow{\Pi}_{A_2} \overleftarrow{\Pi}_{B_1} \rho_{A_2 B_1}^{\otimes n}) [\mathrm{Tr}_{A_2}(\rho_{A_2}^n)]^{-1}, \tag{A9}
\end{align*}
$$

we finally arrive at Eq. (26).

**Appendix B: Statistical-error analysis in FDQC states**

1. Statistical-error analysis for the purity

In this Appendix we prove Eq. (20). We consider the protocol explained in Sec. III A and denote by $\mathcal{P}_2^{(e)}[I]$ the estimates for $\mathcal{P}_2[I] = \mathrm{Tr}(\rho_I^2)$ obtained from the classical-shadow approach. These local purities are estimated with a non-zero relative error, i.e. $\mathcal{P}_2^{(e)}[I]/\mathcal{P}_2[I] = (1 + \varepsilon^I)$. We define

$$
\varepsilon = \max \{ |\varepsilon^I| : I \in \{I_j\}_j \cup \{I_j \cup I_{j+1}\}_j \}. \quad (\text{B1})
$$

We will take $r_2^{(e)}$, defined in Eq. (17), as the experimental
estimate for the global purity $\mathcal{P}_2 = \mathrm{Tr}(\rho^2)$, choosing the
intervals $I_j$ with $|I_j| = k \ge 2\ell - 1$, where $\ell$ is the circuit
depth, so that Eq. (15) holds.

We start our proof from a few preliminary lemmas.

**Lemma 1.** Suppose $\varepsilon \le 1/2$. Then

$$
\left| \frac{r_2^{(e)}}{\mathcal{P}_2} - 1 \right| \le e^{(4L/k)\varepsilon} - 1. \qquad (\text{B2})
$$

Proof. Setting $\mathcal{P}_2^{(e)}[I] = \mathcal{P}_2[I](1+\varepsilon^I)$, recalling $R=L/k$, where $k=|I_j|$ and using (15), we have

$$
\frac{r_2^{(e)}}{\mathcal{P}_2} = \frac{\prod_{j=1}^{R-1} (1 + \varepsilon^{I_j \cup I_{j+1}})}{\prod_{j=2}^{R-1} (1 + \varepsilon^{I_j})}, \quad (\text{B3})
$$

as so, using $|\epsilon^I| \le \epsilon$,

$$
\begin{align*}
\frac{r_2^{(e)}}{\mathcal{P}_2} &\le \frac{(1+\varepsilon)^{R-1}}{(1-\varepsilon)^{R-2}} \le \left[ \frac{1+\varepsilon}{1-\varepsilon} \right]^R, \\
\frac{r_2^{(e)}}{\mathcal{P}_2} &\ge \frac{(1-\varepsilon)^{R-1}}{(1+\varepsilon)^{R-2}} \ge \left[ \frac{1-\varepsilon}{1+\varepsilon} \right]^R. \tag{B4}
\end{align*}
$$

Since $\varepsilon \le 1/2$, and using $(1+x) \le e^x$, we have

$$
\frac{1 + \varepsilon}{1 - \varepsilon} = 1 + \frac{2\varepsilon}{1 - \varepsilon} \le 1 + 4\varepsilon \le e^{4\varepsilon}. \quad (\text{B5})
$$

Therefore,

$$
e^{-(4L/k)\varepsilon} \le \frac{r_2^{(e)}}{\mathcal{P}_2} \le e^{(4L/k)\varepsilon}. \quad (\text{B6})
$$

Finally, using

$$
(1 - e^{-a}) \leq e^{a} - 1, \qquad (B7)
$$

we obtain (B2)

<figure><img src="image_13.png" /></figure>

**Lemma 2.** Let $M \ge 2^{4|I|}$. Then

$$
\mathrm{Pr} \left( |\varepsilon^I| \ge \alpha \right) \le \frac{2^{2|I|+3}}{\alpha^2 M}. \quad (\mathrm{B8})
$$

Proof. We start from the bound (10). Since $M \ge 2^{4|I|}$,
we have in particular $M \ge 4$ and Eq. (10) implies

$$
\mathrm{Var} \left[ \mathcal{P}_2^{(e)}[I] \right] \le 4 \left( \frac{2^{|I|} \mathcal{P}_2[I]}{M} \right) + 4 \left( \frac{2^{2|I|}}{M} \right)^2 . \quad (\text{B9})
$$

Since by hypothesis $M \ge 2^{4|I|}$, while $\mathcal{P}_2[I] \ge 2^{-|I|}$, we
have

$$
4 \left( \frac{2^{|I|}}{M} \right)^2 \le 4 \left( \frac{2^{|I|} \mathcal{P}_2[I]}{M} \right), \quad (\text{B10})
$$

and so

$$
\mathrm{Var} \left[ \mathcal{P}_2^{(e)}[I] \right] \le 8 \frac{2^{|I|} \mathcal{P}_2[I]}{M}. \quad (\text{B11})
$$

Eq. (B8) then follows using Chebyshev's inequality and
$\mathcal{P}_2[I] \ge 2^{-|I|}$. $\square$

We are ready to state and prove the main result of this
section.

**Theorem 1.** Let $0 \le \delta \le 1$ and take

$$
M \ge \max \left\{ 2^{8k}, L^2 \frac{2^{4k+10}}{k^2 \delta^2} \right\}, \quad (B12)
$$

where $k = |I_j| \ge 2$ (the case $k = 1$ is trivial). Then

$$
\mathrm{Pr} \left[ \left| (r_2^{(e)} / \mathcal{P}_2) - 1 \right| \ge \delta \right] \le \frac{2^{4k+11} L^3}{\delta^2 k^3 M}. \quad (\text{B13})
$$

Proof. Recalling the definition (B1), note that if $\varepsilon < (k/8L)\delta$, then $\varepsilon < 1/2$ and also $e^{(4L/k)\varepsilon} - 1 < \delta$. Using Lemma 1, this implies

$$
\begin{align*}
\mathrm{Pr} \left[ \left| (r_2^{(e)} / \mathcal{P}_2) - 1 \right| \ge \delta \right] &\le \mathrm{Pr}[\varepsilon \ge (k/8L)\delta] \\
&= 1 - \mathrm{Pr}[\varepsilon < (k/8L)\delta]. \tag{B14}
\end{align*}
$$

In more detail, the validity of the first inequality can be
seen as follows: Suppose $\varepsilon < (k/8L)\delta$. Then necessarily
$\varepsilon < 1/2$ and so also $|r_2^{(e)}/\mathcal{P}_2 - 1| < \delta$ (by Lemma 1).
Therefore, the set of cases in which $|r_2^{(e)}/\mathcal{P}_2 - 1| \ge \delta$ must
be contained in the set of cases in which $\varepsilon \ge (k/8L)\delta$. In
formula, this is the first line of Eq. (B14).
From the definition of $\varepsilon$ and Lemma 2, we have

$$
\begin{align*}
\Pr[\varepsilon < x] &= \prod_j (1 - \Pr[|\varepsilon^{I_j}| \ge x]) \\
&\quad \times \prod_j (1 - \Pr[|\varepsilon^{I_j \cup I_{j+1}}| \ge x]) \\
&\ge \left[ \left(1 - \frac{2^{2k+3}}{x^2 M}\right) \left(1 - \frac{2^{4k+3}}{x^2 M}\right) \right]^R, \quad (\text{B15})
\end{align*}
$$

where we used that the random variables $\varepsilon^{I_j}$ are statis-
tically independent. Therefore

$$
\Pr[\varepsilon < x] \ge \left(1 - \frac{2^{4k+3}}{x^2 M}\right)^{2R} . \quad (\text{B16})
$$

Setting $x = \delta k/(8L)$, we have by hypothesis
$2^{4k+3}/x^2 M \le 1/2$. Therefore, using $1-z \ge e^{-2z}$ for
$0 \le z \le 1/2$, we obtain

$$
\Pr[\varepsilon < (k/8L)\delta] \ge \exp\left[-\frac{2^{4k+11}L^3}{M\delta^2k^3}\right]. \quad (\text{B17})
$$

Plugging this into (B14), and using $1-e^{-z} \le z$, we finally obtain (B13). $\square$

2. Statistical-error analysis for the PT moments

In this Appendix we prove Eq. (32). Considering the
same protocol and using the same notation as Sec. III B,
we set

$$
p_3^{(e)}[A_2 B_1] = p_3[A_2 B_1] + \varepsilon^{A_2 B_1}, \quad (\text{B18a})
$$

$$
\mathcal{P}_3^{(e)}[A_2] = \mathcal{P}_3[A_2](1 + \varepsilon^{A_2}), \quad (\text{B18b})
$$

$$
\mathcal{P}_3^{(e)}[B_1] = \mathcal{P}_3[B_1](1 + \varepsilon^{B_1}). \quad (\text{B18c})
$$

Note that $\varepsilon^{A_2 B_1}$ is an additive error, while $\varepsilon^{A_2}$, $\varepsilon^{B_1}$ are
relative errors. We also define

$$
\epsilon = 2^{|A_2|+2|B_1|} \max \{|\epsilon^{A_2 B_1}|, |\epsilon^{A_2}|, |\epsilon^{B_1}|\} . \quad (\text{B19})
$$

In the following, we will set $|A_2| = |B_1| = k$, where
$k \ge 2\ell - 1$ with $\ell$ being the circuit depth. With this
choice, Eq. (27) holds.

For clarity, we organize the proof into lemmas.

**Lemma 3.** *Suppose* $\epsilon \le 1$. *Then*

$$
|s_3^{(e)}[A_2 B_1] - s_3[A_2 B_1]| \le 16\epsilon, \quad (\text{B20})
$$

where $s_n[A_2 B_1]$ is defined in Eq. (27).

*Proof.* We start from

$$
s_3^{(e)}[A_2 B_1] = \frac{p_3[A_2 B_1] + \varepsilon^{A_2 B_1}}{\mathcal{P}_3[A_2](1 + \varepsilon^{A_2})\mathcal{P}_3[B_1](1 + \varepsilon^{B_1})}, \quad (\text{B21})
$$

which gives us

$$
\begin{align}
|s_3^{(e)}[A_2 B_1] - s_3[A_2 B_1]| &\le \frac{|p_3[A_2 B_1]|}{\mathcal{P}_3[A_2]\mathcal{P}_3[B_1]} \frac{|\varepsilon^{A_2}| + |\varepsilon^{B_1}| + |\varepsilon^{A_2} \varepsilon^{B_1}|}{(1 + \varepsilon^{A_2})(1 + \varepsilon^{B_1})} \nonumber \\
&\quad + \frac{|\varepsilon^{A_2 B_1}|}{\mathcal{P}_3[A_2]\mathcal{P}_3[B_1](1 + \varepsilon^{A_2})(1 + \varepsilon^{B_1})}. \tag{B22}
\end{align}
$$

Denoting by $\lambda_j$ the eigenvalues of $\rho_{A_2 B_1}^{T_{A_2}}$, we have

$$
|p_3[A_2 B_1]| \le \sum_j |\lambda_j|^3 \le \sum_j \lambda_j^2 = \mathcal{P}[A_2 B_1] \le 1, \quad (\text{B23})
$$

where we used $\lambda_j \in [-1/2, 1]$ [144] and $\mathrm{Tr}[(\rho_{AB}^{T_A})^2] =$
$\mathrm{Tr}[\rho_{AB}^2]$ [35]. In addition, since $\varepsilon \le 1$, clearly $|\varepsilon^{A_2}| \le$
$1/2$, $|\varepsilon^{B_1}| \le 1/2$ and so also $|\varepsilon^{A_2} \varepsilon^{B_1}| < |\varepsilon^{A_2}|$. Putting
everything together, and using $\mathcal{P}_3[I] \le 2^{-2|I|}$, we finally
arrive at Eq. (B20). $\square$

**Lemma 4.** Let $M \ge 3 \cdot 2^{8k}$. Then

$$
\Pr (|\varepsilon^{A_2}| \ge \alpha) \le 27 \frac{2^{3k}}{\alpha^2 M}, \quad (\text{B24a})
$$

$$
\Pr \left( |\varepsilon^{B_1}| \ge \alpha \right) \le 27 \frac{2^{3k}}{\alpha^2 M}, \quad (\text{B24b})
$$

$$
\Pr \left( |\varepsilon^{A_2 B_1}| \ge \alpha \right) \le 27 \frac{2^{2k}}{\alpha^2 M}. \quad (\text{B24c})
$$

*Proof.* We start from the bound [38]

$$
\begin{equation}
\begin{split}
\operatorname{Var} [p_3^{(e)}[AB]] \le {}& 9 \frac{2^{|AB|}}{M} \operatorname{Tr} (\rho_{AB}^4) \\
& + 18 \frac{2^{3|AB|}}{(M-1)^2} p_2[AB] + 6 \frac{2^{6|AB|}}{(M-2)^3}.
\end{split}
\tag{B25}
\end{equation}
$$

Since $M \ge 3 \cdot 2^{8k}$ (and $k \ge 1$), then

$$
\frac{18}{(M-1)^2} \le \frac{27}{M^2}, \quad \frac{6}{(M-2)^3} \le \frac{9}{M^3}. \qquad (\text{B26})
$$

Therefore

$$
\begin{equation}
\begin{split}
\operatorname{Var} [p_3^{(e)}[AB]] \le {}& 9 \frac{2^{|AB|}}{M} \operatorname{Tr} (\rho_{AB}^4) \\
& + 27 \frac{2^{3|AB|}}{M^2} p_2[AB] + 9 \frac{2^{6|AB|}}{M^3}.
\end{split}
\tag{B27}
\end{equation}
$$

Since $M \ge 3 \cdot 2^{8k}$, we have

$$
\begin{align}
& 9 \frac{2^{|AB|}}{M} \operatorname{Tr} (\rho_{AB}^4) \ge 27 \frac{2^{3|AB|}}{M^2} p_2[AB], \nonumber \\
& 9 \frac{2^{|AB|}}{M} \operatorname{Tr} (\rho_{AB}^4) \ge 9 \frac{2^{6|AB|}}{M^3}. \tag{B28}
\end{align}
$$

In the first line we have used that $p_2[I] = \mathrm{Tr}[\rho_I^2]$ and
Hölder's inequality, which guarantees

$$
\mathrm{Tr}[\rho_I^2]/\mathrm{Tr}[\rho_I^4] \le 2^{|I|}/\mathrm{Tr}[\rho_I^2] \le 2^{2|I|}, \quad (\text{B29})
$$
while in the second line we have used Tr(ρ<sub>AB</sub><sup>4</sup>) ≥ 2<sup>-3|AB|</sup> = 2<sup>-6k</sup>. Putting all together, we get

$$
\mathrm{Var} \left[ p_3^{(e)}[AB] \right] \le 27 \frac{2^{2k}}{M} \mathrm{Tr} \left( \rho_{AB}^4 \right). \quad (\mathrm{B30})
$$

Similarly, we have [38]

$$
\begin{align*}
\mathrm{Var}[\mathcal{P}_3[A_2]] \le {}& 9 \frac{2^{|A_2|}}{M} \mathrm{Tr}(\rho_{A_2}^4) \\
& + 18 \frac{2^{3|A_2|}}{(M-1)^2} \mathrm{Tr}(\rho_{A_2}^2) + 6 \frac{2^{|B|}}{(M-2)^3} \\
& \le 27 \frac{2^k}{M} \mathrm{Tr}(\rho_{A_2}^4). \tag{B31}
\end{align*}
$$

Using Tr (ρ<sub>A<sub>2</sub></sub><sup>4</sup>) / (Tr(ρ<sub>A<sub>2</sub></sub><sup>3</sup>))<sup>2</sup> ≤ 1/Tr(ρ<sub>A<sub>2</sub></sub><sup>3</sup>) ≤ 2<sup>2k</sup>, we obtain

$$
\mathrm{Var} \left[ \frac{\mathcal{P}_3^{(e)}[A_2]}{\mathcal{P}_3[A_2]} \right] \le 27 \frac{2^k}{M} \frac{\mathrm{Tr}(\rho_{A_2}^4)}{(\mathrm{Tr}(\rho_{A_2}^3))^2} \le 27 \frac{2^{3k}}{M}, \quad (\text{B32})
$$

and analogously for $\mathcal{P}_3^{(2)}(B_1)$. Eqs. (B24) follow using
Chebyshev's inequality. $\square$

We are ready to state and prove the main result of this
section, yielding Eq. (32).

**Theorem 2.** Let $0 \le \delta \le 1$ and take

$$
M \geq 27 \frac{2^{11k+9}}{\delta^2}. \qquad (\text{B33})
$$

Then

$$
\mathrm{Pr} \left[ \left| s_3^{(e)}[A_2 B_1] - s_3[A_2 B_1] \right| \ge \delta \right] \le 81 \frac{2^{11k+9}}{M \delta^2}. \quad (\text{B34})
$$

*Proof.* Recall the definition (B19) and note that if $\varepsilon \le \delta/16$, then trivially $\varepsilon < 1$ and so, by Lemma 3, $|s_3^{(e)} - s_3| \le \delta$. This implies

$$
\begin{align}
\mathrm{Pr}[|s_3^{(e)} - s_3| \ge \delta] &\le \mathrm{Pr}[\varepsilon \ge \delta/16] \nonumber \\
&= 1 - \mathrm{Pr}[\varepsilon < \delta/16]. \tag{B35}
\end{align}
$$

In more detail, the validity of the first inequality can be
seen as follows: Suppose $\varepsilon < \delta/16$. Then trivially $\varepsilon < 1$
and so necessarily also $|s_3^{(e)} - s_3| < \delta$ (by Lemma 3).
Therefore, the set of cases in which $|s_3^{(e)} - s_3| \ge \delta$ is
contained in the set of cases for which $\varepsilon \ge \delta/16$. In
formula, this is the first line of Eq. (B35).

From the definition of $\varepsilon$ and Lemma 4, we have

$$
\begin{align*}
\mathrm{Pr}[\varepsilon < x] &= (1 - \mathrm{Pr}[|\varepsilon^{A_2}| \ge x 2^{-4k}]) \\
&\quad \times (1 - \mathrm{Pr}[|\varepsilon^{B_1}| \ge x 2^{-4k}]) (1 - \mathrm{Pr}[|\varepsilon^{A_2 B_1}| \ge x 2^{-4k}]) \\
&\ge \left(1 - 27 \frac{2^{11k}}{x^2 M}\right)^3, \tag{B36}
\end{align*}
$$

where we used that the random variables $\varepsilon^{I_j}$ are statisti-
cally independent. Setting $x = \delta/16$, we have by hypoth-
esis $2 \cdot 2^{11k}/x^2 M < 1/2$. Therefore, using $1-z \ge e^{-2z}$
for $0 \le z \le 1/2$, we obtain

$$
\mathrm{Pr}[\varepsilon < \delta/16] \ge \exp \left[ -81 \frac{2^{11k+9}}{M\delta^2} \right]. \quad (\mathrm{B37})
$$

Plugging into (B35), and using $1 - e^{-z} \le z$, we finally obtain (B34). $\square$

**Appendix C: Statistical-error analysis for states**
**satisfying the AFCs**

In this Appendix we prove the main claims reported in
Sec. IV B. Let $\rho$ be a state satisfying (34) for all partitions
as in Fig. 3(a), with $|B| = k \ge k_c$ (and open boundary
conditions). We first prove Eq. (37).

We start with the following:

**Lemma 5.** Let

$$
k \ge \frac{\ln(\alpha_2 L / \delta)}{\beta_2}, \quad (C1)
$$

for $0 \le \delta \le 1/2$. Then, defining $r_2$ as in Eq. (39), we
have

$$
\left| \frac{r_2}{\mathrm{Tr}(\rho^2)} - 1 \right| \le 4\delta. \qquad (\mathrm{C2})
$$

Proof. First, applying iteratively (34), we get (recalling
that $R = L/k$)

$$
\frac{r_2}{\mathrm{Tr}(\rho^2)} = \prod_{j=2}^{R-1} (1 + \varepsilon_j), \quad (\mathrm{C3})
$$

with $|\epsilon_j| \le \alpha_2 e^{-\beta_2 k}$, and so

$$
\left| \frac{r_2}{\mathrm{Tr}(\rho^2)} - 1 \right| = \left| \prod_{j=2}^{R-1} (1 + \epsilon_j) - 1 \right|. \quad (\text{C4})
$$

We have

$$
\begin{align*}
& \prod_{j=2}^{R-1} (1 + \epsilon_j) - 1 \le (1 + \alpha_2 e^{-\beta_2 k})^{R-2} - 1 \\
& \le \exp[\alpha_2(L/k)e^{-\beta_2 k}] - 1. \tag{C5}
\end{align*}
$$

Analogously,

$$
\prod_{j=2}^{R-1} (1 + \epsilon_j) - 1 \ge (1 - \alpha_2 e^{-\beta_2 k})^{R-2} - 1. \quad (\text{C6})
$$

Due to (C1), we have $\alpha_2 e^{-\beta_2 k} \le 1/2$. Therefore, using $1-z \ge e^{-2z}$ for $0 \le z \le 1/2$, we have

$$
\prod_{j=2}^{R-1} (1 + \epsilon_j) - 1 \ge \exp[-2\alpha_2(L/k)e^{-\beta_2 k}] - 1. \quad (\text{C7})
$$
Therefore

$$e^{-y} - 1 \le \prod_{j=2}^{R-1} (1 + \epsilon_j) - 1 \le e^{y/2} - 1 \le e^y - 1, \quad (C8)$$

where $y = 2\alpha_2(L/k)e^{-\beta_2 k}$. Since $1 - e^{-y} \le e^y - 1$, this implies

$$\left| \prod_{j=2}^{R-1} (1 + \epsilon_j) - 1 \right| \le \exp[2\alpha_2(L/k)e^{-\beta_2 k}] - 1. \quad (C9)$$

Finally, due to (C1), we have $2\alpha_2(L/k)e^{-\beta_2 k} \le 1$. Therefore, using $e^z - 1 \le 2z$ for $0 \le z \le 1$, we arrive at

$$\left| \frac{r_2}{\text{Tr}(\rho^2)} - 1 \right| \le \frac{4\alpha_2 L}{k} e^{-\beta_2 k} \le 4\delta. \quad (C10)$$

<figure><img src="image_14.png" /></figure>

Next, we prove Eq. (37), via the following:

**Theorem 3.** Let $0 \le \delta \le 1/2$ and set

$$k = \frac{\ln(\alpha_2 L / \delta)}{\beta_2}. \quad (C11)$$

Choosing

$$M \ge \max \left\{ \left( \frac{\alpha_2 L}{\delta} \right)^{\frac{8 \ln 2}{\beta_2}}, \frac{L^2 2^{10}}{\delta^2} \left( \frac{\alpha_2 L}{\delta} \right)^{\frac{4 \ln 2}{\beta_2}} \right\}, \quad (C12)$$

and recalling the definition (17), we have

$$\Pr \left[ \left| (r_2^{(e)} / \mathcal{P}_2) - 1 \right| \ge 7\delta \right] \le \frac{2^{11} L^3}{\delta^2 M} \left( \frac{\alpha_2 L}{\delta} \right)^{\frac{4 \ln 2}{\beta_2}}. \quad (C13)$$

*Proof.* Set $r_2^{(e)}/r_2 = (1+\epsilon_1)$ and $r_2/\mathcal{P}_2 = (1+\epsilon_2)$, where $r_2^{(e)}$ is defined in Eq. (17). Using Lemma 5 (and that $\delta \le 1/2$), we have $|\epsilon_2| \le 4\delta$ and so

$$\left| \frac{r_2^{(e)}}{\mathcal{P}_2} - 1 \right| = \left| \frac{r_2^{(e)}}{r_2} \frac{r_2}{\mathcal{P}_2} - 1 \right| \le 4\delta + 3|\epsilon_1|. \quad (C14)$$

If $|\epsilon_1| < \delta$, then $|r_2^{(e)}/\mathcal{P}_2 - 1| < 7\delta$. Therefore, the set of cases in which $|r_2^{(e)}/\mathcal{P}_2 - 1| \ge 7\delta$ is contained in the set of cases in which $|\epsilon_1| \ge \delta$. In formula,

$$\Pr[|(r_2^{(e)}/\mathcal{P}_2) - 1| \ge 7\delta] \le \Pr[|(r_2^{(e)}/r_2) - 1| \ge \delta]. \quad (C15)$$

Thanks to Eq. (C12), we can use Theorem 1, yielding

$$\Pr[|(r_2^{(e)}/r_2) - 1| \ge \delta] \le \frac{2^{4k+11} L^3}{\delta^2 k^3 M} \le \frac{2^{4k+11} L^3}{\delta^2 M} \\ \le \frac{2^{11} L^3}{\delta^2 M} \left( \frac{\alpha_2 L}{\delta} \right)^{\frac{4 \ln 2}{\beta_2}}, \quad (C16)$$

which completes the proof. □

Note that, in order to recover Eqs. (36) and (37), we simply rename $\delta' = 7\delta$.

Finally, we present a technical result showing that the PT moments can be estimated efficiently, assuming that the state to be measured satisfies the AFCs. We focus for simplicity on the case $n = 3$ and follow a protocol similar to that of Sec. III B. Namely, we take $s_3^{(e)}[A_2 B_1]$, defined in (29), as our estimator for $\tilde{p}_3[AB]$, and, for any $0 \le \delta \le 1$, we choose

$$k = \frac{\ln(2\alpha_3/\delta)}{2\beta_3}. \quad (C17)$$

Finally, we perform $M$ measurements to estimate $p_3[A_2 B_1]$, $\mathcal{P}_3[A_2]$, and $\mathcal{P}_3[B_1]$ each (so the total number is $M_T = 3M$). We can then prove the following:

**Theorem 4.** For any $0 \le \delta \le 1/2$, set

$$k = \frac{\ln(2\alpha_3/\delta)}{2\beta_3}. \quad (C18)$$

If

$$M \ge 27 \frac{2^9}{\delta^2} \left( \frac{2\alpha_3}{\delta} \right)^{\frac{11 \ln 2}{2\beta_3}}, \quad (C19)$$

then

$$\Pr \left[ \left| s_3^{(e)} - \tilde{p}_3 \right| \ge \delta \right] \le \frac{81 \cdot 2^9}{\delta^2 M} \left( \frac{2\alpha_3}{\delta} \right)^{\frac{11 \ln 2}{2\beta_3}}. \quad (C20)$$

*Proof.* Our estimator for the normalized PT moment $\tilde{p}_3[AB]$ is

$$s_3^{(e)} = \frac{p_3^{(e)}[A_2 B_1]}{\mathcal{P}_3^{(e)}[A_2] \mathcal{P}_3^{(e)}[B_1]}. \quad (C21)$$

First, suppose $|s_3 - s_3^{(e)}| < \delta/2$. Then, using Eqs. (35) and (C18), we have

$$|s_3^{(e)} - \tilde{p}_3| \le |s_3^{(e)} - s_3| + |s_3 - \tilde{p}_3| \\ \le (\delta/2) + (\delta/2) = \delta. \quad (C22)$$

Therefore, the set of cases in which $|s_3^{(e)} - \tilde{p}_3| \ge \delta$ is contained in the set of cases in which $|s_3^{(e)} - s_3| \ge \delta/2$, namely

$$\Pr[|s_3^{(e)} - \tilde{p}_3| \ge \delta] \le \Pr[|s_3^{(e)} - s_3| \ge \delta]. \quad (C23)$$

Finally, thanks to Eq. (C19), we can apply Theorem 2, yielding

$$\Pr \left[ \left| s_3^{(e)} - \tilde{p}_3 \right| \ge \delta \right] \le \frac{81 \cdot 2^9}{\delta^2 M} \left( \frac{2\alpha_3}{\delta} \right)^{\frac{11 \ln 2}{2\beta_3}}. \quad (C24)$$

<figure><img src="image_15.png" /></figure>
## Appendix D: AFCs and MPDOs

In this Appendix we prove that the purity AFCs hold for MPDOs. To this end, we assume the following conditions on the transfer matrices (45):

(A) The matrices $\tau_1$ and $\tau_2$ admit the spectral decomposition

$$ \tau_1 = \sum_{j=0}^{\chi-1} \lambda_j |R_j^{(1)}\rangle\langle L_j^{(1)}|, \quad (\text{D1a}) $$

$$ \tau_2 = \sum_{j=0}^{\chi^2-1} \mu_j |R_j^{(2)}\rangle\langle L_j^{(2)}|, \quad (\text{D1b}) $$

where we assume $|\lambda_0| > |\lambda_j|$, $|\mu_0| > |\mu_j|$ for all $j \ge 1$, i.e. $\tau_1$, $\tau_2$ have a trivial Jordan form and a finite gap.<fnref n="13" /> Note that we can also assume without the loss of generality that $\lambda_0 = 1$ (which implies $\mu_0 > 0$, since $\rho_L > 0$ for all $L$). Finally, $|R_j^{(n)}\rangle$ and $\langle L_j^{(n)}|$ are the left and right eigenstates, i.e, they are vectors on the left/right virtual indices of $\tau_1$, and $\tau_2$ which are normalized such that

$$ \langle L_j^{(n)} | R_k^{(n)} \rangle = \delta_{j,k} \quad n = 1, 2. \quad (\text{D2}) $$

(B) We further need to assume the technical conditions

$$ (\langle L_0^{(1)} | \otimes \langle L_0^{(1)} | | R_0^{(2)} \rangle) \neq 0, \quad (\text{D3a}) $$

$$ \langle L_0^{(2)} | (|R_0^{(1)}\rangle \otimes |R_0^{(1)}\rangle) \neq 0, \quad (\text{D3b}) $$

which again are quite general, as orthogonality requires fine-tuning.

We use the same notations and assumptions as in Sec. IV C, so that

$$ \mathcal{P}_2 = \frac{\mathrm{Tr}[\sigma^2]}{(\mathrm{Tr}[\sigma])^2}, \quad (\text{D4}) $$

and

$$ r_2 = \frac{\prod_{j=1}^{R} \mathrm{Tr}_{I_j \cup I_{j+1}} (\sigma_{I_j \cup I_{j+1}}^2)}{\prod_{j=1}^{R} \mathrm{Tr}_{I_j} [\sigma_{I_j}^2]}, \quad (\text{D5}) $$

where $|I_j| = k$ and $R = L/k$ is an integer. First, we introduce the correlation lengths

$$ e^{-\frac{1}{\zeta_1}} = \frac{\max_{j>0}\{|\lambda_j|\}}{\lambda_0}, \quad e^{-\frac{1}{\zeta_2}} = \frac{\max_{j>0}\{|\mu_j|\}}{\mu_0}, \quad (\text{D6}) $$

and also

$$ \zeta = \max(\zeta_1, \zeta_2). \quad (\text{D7}) $$

Next, given the spectral decomposition in Eqs. (D1), we define

$$ C = \max_{\substack{(j,k,l) \\ \neq (0,0,0)}} \left\{ \frac{|\langle L_j^{(1)} | \langle L_k^{(1)} | R_l^{(2)} \rangle \langle L_l^{(2)} | R_j^{(1)} \rangle | R_k^{(1)}|}{|\langle L_0^{(1)} | \langle L_0^{(1)} | R_0^{(2)} \rangle \langle L_0^{(2)} | R_0^{(1)} \rangle | R_0^{(1)}|} \right\}. \quad (\text{D8}) $$

Note that the denominator is non-vanishing because of Eqs. (D3). We can now state the main result of this section.

**Theorem 5.** Take $|I_j| = k$ with

$$ k \ge k_{\min} = \max\{1, \zeta \ln(20C\chi^2L)\}. \quad (\text{D9}) $$

Then, for all

$$ L \ge \max\{\zeta \ln(2^5\chi^2), 4k + \zeta \ln(2\chi)\}, \quad (\text{D10}) $$

we have

$$ \left| \frac{r_2}{\mathcal{P}_2} - 1 \right| \le \chi^2(80C + 32) \frac{L}{k} e^{-k/\zeta}. \quad (\text{D11}) $$

*Proof.* Using Eqs. (D1), we have

$$ \mathcal{P}_2 = \frac{\sum_{j=0}^{\chi^2-1} \mu_j^L}{\left(\sum_{j=0}^{\chi^2-1} \lambda_j^L\right)^2} = \frac{\mu_0^L}{\lambda_0^{2L}} (1 + \tilde{\varepsilon}), \quad (\text{D12}) $$

where

$$ \tilde{\varepsilon} = \frac{1 + \sum_{j=1}^{\chi^2-1} (\mu_j/\mu_0)^L}{\left(1 + \sum_{j=1}^{\chi^2-1} (\lambda_j/\lambda_0)^L\right)^2} - 1. \quad (\text{D13}) $$

Since by hypothesis $L \ge \zeta_1 \ln(2\chi)$, we have $\chi e^{-L/\zeta_1} \le 1/2$, and so

$$ \begin{aligned} |\tilde{\varepsilon}| &\le 4 \left( \chi^2 e^{-L/\zeta_2} + 2\chi e^{-L/\zeta_1} + \chi^2 e^{-2L/\zeta_1} \right) \\ &\le 2^4 \chi^2 e^{-L/\zeta}. \end{aligned} \quad (\text{D14}) $$

On the other hand

$$ \mathrm{Tr}[\sigma_I^2] = \sum_{j,k=0}^{\chi-1} \lambda_j^{L-|I|} \lambda_k^{L-|I|} (\langle L_j^{(1)} | \otimes \langle L_k^{(1)} |) \times \langle L_0^{(2)} | R_0^{(1)} \rangle | R_0^{(1)} \rangle (1 + \varepsilon^I), \quad (\text{D15}) $$

Therefore, recalling the definition (D8), we have

$$ \mathrm{Tr}[\sigma_I^2] = \lambda_0^{2(L-|I|)} \mu_0^{|I|} \langle L_0^{(1)} | \langle L_0^{(1)} | R_0^{(2)} \rangle \times \langle L_0^{(2)} | R_0^{(1)} \rangle | R_0^{(1)} \rangle (1 + \varepsilon^I), \quad (\text{D16}) $$

<footnote n="13">¹³ The assumption that $\tau_1$ and $\tau_2$ can be diagonalized is purely technical and not necessary. However, we keep it here as it makes the analysis simpler and it is in any case quite general.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 9. Scaling of the relative error $\epsilon_r(k, L, \beta)$, defined in Eq. (52), as a function of the bond dimension $\chi$ used to approximate the thermal state. In the left panel we show data for the Ising chain while the right panel shows data for the XXZ chain. In both cases we use the same parameters as in Fig. 4. Data shown for $k = 4$, other values of $k$ show similar behavior with $\chi = 16$ being sufficient for convergence in the Ising model, while the XXZ model appears to be well converged from $\chi = 64$.</figcaption></figure>

where

$$ |\epsilon^I| \le 2C\chi e^{-(L-|I|)/\zeta_1} + C\chi^2 e^{-|I|/\zeta_2} \\ +C\chi^2 e^{-2(L-|I|)/\zeta_1} + 2C\chi^3 e^{-(L-|I|)/\zeta_1} e^{-|I|/\zeta_2} \\ +C\chi^4 e^{-2(L-|I|)/\zeta_1} e^{-|I|/\zeta_2}. \qquad (D17) $$

Since by hypothesis $L \ge 4k + \zeta \ln(2\chi)$, it is easy to verify that all the five terms above are upper bounded by $C\chi^2 e^{-|I|/\zeta}$ (recall that $|I| \le 2k$ for all $I$), and so

$$ |\epsilon^I| \le 5C\chi^2 e^{-|I|/\zeta} \le 5C\chi^2 e^{-k/\zeta}, \qquad (D18) $$

where we used that either $|I| = k$ or $|I| = 2k$, and so $|I| \ge k$. Therefore

$$ \frac{\prod_{j=1}^{R} \mathrm{Tr}_{I_j \cup I_{j+1}} (\sigma_{I_j \cup I_{j+1}}^2)}{\prod_{j=1}^{R} \mathrm{Tr}_{I_j} [\sigma_{I_j}^2]} \\ = \frac{\mu_0^L}{\lambda_0^{2L}} \left( \frac{\prod_{j=1}^{R} (1 + \varepsilon^{I_j \cup I_{j+1}})}{\prod_{j=1}^{R} (1 + \varepsilon^{I_j})} \right), \qquad (D19) $$

with $|\epsilon^I| \le 5C\chi^2 e^{-k/\zeta}$. Combining this with Eq. (D12), we arrive at

$$ \frac{r_2}{\mathcal{P}_2} - 1 = \left( \frac{\prod_{j=1}^{R} (1 + \varepsilon^{I_j \cup I_{j+1}})}{\prod_{j=1}^{R} (1 + \varepsilon^{I_j})} \right) (1 + \tilde{\varepsilon})^{-1} - 1 \\ = (1 + \tilde{\varepsilon})^{-1} \left[ \frac{\prod_{j=1}^{R} (1 + \varepsilon^{I_j \cup I_{j+1}})}{\prod_{j=1}^{R} (1 + \varepsilon^{I_j})} - 1 - \tilde{\varepsilon} \right], \qquad (D20) $$

where $\tilde{\varepsilon}$ is given in Eq. (D13).

Next, we define

$$ \varepsilon = \max \{ |\varepsilon^I| : I \in \{I_j\}_j \cup \{I_j \cup I_{j+1}\}_j \}. \qquad (D21) $$

By hypothesis $k \ge k_{\min}$ [$k_{\min}$ is given in (D9)]. Therefore, using (D18) we have $|\varepsilon^I| \le 1/2$ for all $|I|$ and so

also $\varepsilon \le 1/2$. Therefore, we can apply the derivation in Lemma 1 to show

$$ \left| \frac{\prod_{j=1}^{R} (1 + \varepsilon^{I_j \cup I_{j+1}})}{\prod_{j=1}^{R} (1 + \varepsilon^{I_j})} - 1 \right| \le e^{(4L/k)\varepsilon} - 1. \qquad (D22) $$

Since by hypothesis $L \ge \xi \ln(2^5\chi^2)$, we also have $\tilde{\varepsilon} \le 1/2$ [cf. Eq. (D14)], and so (D20) yields

$$ \left| \frac{r_2}{\mathcal{P}_2} - 1 \right| \le 2(e^{(4L/k)\varepsilon} - 1) + 2|\tilde{\varepsilon}|. \qquad (D23) $$

Finally, we note that Eq. (D9) implies that $(4L/k)\varepsilon \le 1$, and using $e^z - 1 \le 2z$ for $0 \le z \le 1$, we arrive at

$$ \begin{aligned} \left| \frac{r_2}{\mathcal{P}_2} - 1 \right| & \le 80\chi^2 C(L/k) e^{-k/\zeta} + 2^5 \chi^2 e^{-L/\zeta} \\ & \le \chi^2 (80C + 32) \frac{L}{k} e^{-k/\zeta}. \end{aligned} \qquad (D24) \qquad \square $$

This theorem proves Eq. (48), under the condition (47). Therefore, the approximate factorization property (40) holds for MPDOs, with the identification

$$ \alpha_2 = (20C + 8)\chi^2, \quad \beta_2 = 1/\zeta. \qquad (D25) $$

## Appendix E: Details on the numerical computations

In this Appendix we provide further details on the numerical computations performed to obtain the data presented in Secs. IVD and V. As mentioned, the calculations are carried out using the iTensor library [139], by first approximating the thermal states by MPOs of bond dimension $\chi$ and subsequently taking powers and traces of the density matrices represented in this way. For each quantity, we have always verified that the results were
stable upon increasing the bond dimension $\chi$. An example of our data is reported in Fig. 9, where we study $\epsilon_r(k, L, \beta)$, defined in Eq. (52), as a function of the bond dimension $\chi$ used to approximate the thermal state. In general, we have found that relatively small bond dimensions are enough in the quantum Ising chain, while larger bond dimensions are required in order to observe convergence in the Heseinberg model.

## Appendix F: Technical details on the PPT conditions

The goal of this Appendix is to identify the conditions under which the class of states $\rho^{(L)}(\gamma)$ introduced in Sec. VB satisfy Eqs. (60) for all system sizes $L$ and some suitable constants $C_3, C_5$ independent of $L$.

For concreteness, let the circuit depth $\ell$ be even, $\ell = 2k$ with $k \ge 1$ (a similar discussion holds for $\ell = 2k+1$) and take a bipartition of the system as in Fig. 10. We group neighboring sites into blocks containing $q = \ell - 1$ qubits, forming a new super qudit associated with a Hilbert-space of dimension $d = 2^{\ell-1}$. It is easy to show that the circuit can be rewritten as a depth-2 quantum circuit acting on the super qudits, cf. Fig. 10. Therefore (assuming without loss of generality $R = L/q$ is an integer)

$$ \rho_L(\gamma) = V^{(2)} \left( \bigotimes_{j=1}^{L/q} \omega_j \right)^{\dagger}, \quad (F1) $$

where $\omega_j = \bigotimes_{i=1}^q \sigma_{jq+i}$, while

$$ V^{(2)} = \left( \prod_j V_{2j,2j+1} \right) \left( \prod_j V_{2j+1,2j+2} \right), \quad (F2) $$

is the depth-2 FDQC acting on the super lattice.

As it is manifest from Fig. 10, the region A contains all super qudits with labels from 1 to $L/2q$, while B contains all those with labels from $L/2q+1$ to $R = L/q$. Define the sets of super qudits $S_1 = \{R/2 - 1, R/2\}$, $S_2 = \{R/2 + 1, R/2 + 2\}$ and

$$ \tilde{\rho}_{S_1 S_2} = W \omega_{R/2-1} \otimes \omega_{R/2} \otimes \omega_{R/2+1} \otimes \omega_{R/2+2} W^{\dagger}, \quad (F3) $$

where $W = V_{R/2,R/2+1}V_{R/2-1,R/2}V_{R/2+1,R/2+2}$. Note that $\tilde{\rho}_{S_1 S_2}$ is different from the reduced density matrix over $S_1 \cup S_2$. Using the graphical representation for the blocked circuits and Eq. (27), it is simple to show

$$ \tilde{p}_n[AB] = \frac{\mathrm{Tr} \left[ \left( \tilde{\rho}_{S_1 S_2}^{T_{S_1}} \right)^n \right]}{\mathrm{Tr}_{S_1} \left[ \tilde{\rho}_{S_1}^n \right] \mathrm{Tr}_{S_2} \left[ \tilde{\rho}_{S_2}^n \right]} =: \tilde{s}_n, \quad (F4) $$

where $\tilde{\rho}_{S_1} = \mathrm{Tr}_{S_2}[\tilde{\rho}_{S_1 S_2}]$ and $\tilde{\rho}_{S_2} = \mathrm{Tr}_{S_1}[\tilde{\rho}_{S_1 S_2}]$. Therefore, $\tilde{p}_n[AB]$ coincides with the normalized PT moments of the state $\tilde{\rho}_{S_1 S_2}$, supported on four super qudits. Next, it is also easy to compute

$$ \begin{aligned} \mathcal{P}_n[A] &= \mathrm{Tr}_{S_1}[\tilde{\rho}_{S_1}^n][\gamma^n + (1-\gamma)^n]^{L/2-2q}, \\ \mathcal{P}_n[B] &= \mathrm{Tr}_{S_2}[\tilde{\rho}_{S_2}^n][\gamma^n + (1-\gamma)^n]^{L/2-2q}. \end{aligned} \quad (F5) $$

Finally, setting

$$ \tilde{t}_3 = \frac{\mathrm{Tr}_{S_1}[\tilde{\rho}_{S_1}^2]^2 \mathrm{Tr}_{S_2}[\tilde{\rho}_{S_2}^2]^2}{\mathrm{Tr}_{S_1}[\tilde{\rho}_{S_1}^3] \mathrm{Tr}_{S_2}[\tilde{\rho}_{S_2}^3]}, \quad (F6) $$

$$ \tilde{t}_5 = \frac{\mathrm{Tr}_{S_1}[\tilde{\rho}_{S_1}^4]^2 \mathrm{Tr}_{S_2}[\tilde{\rho}_{S_2}^4]^2}{\mathrm{Tr}_{S_1}[\tilde{\rho}_{S_1}^3] \mathrm{Tr}_{S_2}[\tilde{\rho}_{S_2}^3] \mathrm{Tr}_{S_1}[\tilde{\rho}_{S_1}^5] \mathrm{Tr}_{S_2}[\tilde{\rho}_{S_2}^5]}, \quad (F7) $$

we arrive at

$$ f_3 = \tilde{s}_3 - \tilde{s}_2^2 \tilde{t}_3 \frac{[\gamma^2 + (1-\gamma)^2]^{2L-8q}}{[\gamma^3 + (1-\gamma)^3]^{L-4q}}, \quad (F8a) $$

$$ f_5 = \tilde{s}_5 \tilde{s}_3 - \tilde{s}_4^2 \tilde{t}_5 \frac{[\gamma^4 + (1-\gamma)^4]^{2L-8q}}{[\gamma^3 + (1-\gamma)^3]^{L-4q} [\gamma^5 + (1-\gamma)^5]^{L-4q}}. \quad (F8b) $$

Now, define

$$ K_3 = \max\{\tilde{s}_3(\gamma) - \tilde{s}_2^2(\gamma)\tilde{t}_3(\gamma) : 0 \le \gamma \le 1/4\}, \quad (F9) $$

$$ K_5 = \max\{\tilde{s}_5(\gamma)\tilde{s}_3(\gamma) - \tilde{s}_4^2(\gamma)\tilde{t}_5(\gamma) : 0 \le \gamma \le 1/4\}, \quad (F10) $$

and

$$ H_3 = \max\{|\tilde{s}_2^2(\gamma)\tilde{t}_3(\gamma)| : 0 \le \gamma \le 1/4\}, \quad (F11) $$

$$ H_5 = \max\{|\tilde{s}_4^2(\gamma)\tilde{t}_5(\gamma)| : 0 \le \gamma \le 1/4\}. \quad (F12) $$

The constants $K_3, K_5, H_3$, and $H_5$ depend on the specific choices of the unitary gates forming $U^{(\ell)}$. For non-entangling gates, $K_3$ and $K_5$ are positive, but for generic choices of gates one has $K_3 < 0$ and $K_5 < 0$. We are ready to state our main result.

**Theorem 6.** Suppose $K_3, K_5 < 0$. Then, assuming without loss of generality $\gamma \le 1/4$ and defining

$$ C_l = \frac{|K_l|}{2}, \quad \gamma_3 = -\frac{K_3}{4H_3}, \quad \gamma_5 = \left(-\frac{K_5}{8H_5}\right)^{1/3}, \quad (F13) $$

we have

$$ 0 \le \gamma \le \gamma_3/L \Rightarrow f_3 \le -C_3, \quad (F14) $$

while

$$ 0 \le \gamma \le \gamma_5/L^{1/3} \Rightarrow f_5 \le -C_5, \quad (F15) $$

for all system sizes $L$.

*Proof.* We start from Eq. (F8a) and note

$$ \frac{[\gamma^2 + (1-\gamma)^2]^2}{\gamma^3 + (1-\gamma)^3} = 1 - \epsilon_3, \quad (F16) $$

with $0 \le \epsilon_3 \le \gamma$ for $0 \le \gamma \le 1/4$. Therefore

$$ \begin{aligned} f_3 &\le K_3 + \tilde{s}_2^2 \tilde{t}_3 [1 - (1 - \epsilon_3)^{L-4q}] \\ &\le K_3 + H_3 [1 - (1 - \gamma)^{L-4q}]. \end{aligned} \quad (F17) $$
<figure><img src="image_17.png" /><figcaption>Figure 10. Any FDQC of depth $ℓ$ can be rewritten as a depth-2 quantum circuit after grouping $ℓ-1$ qubits into a single qudit. In the middle picture, we highlight with different colors sets of gates defining the two-qudit gates in the grouped lattice (colored rectangles in the right figure). In all the figures, a dashed orange line separates the regions A and B, defining the bipartition of the system.</figcaption></figure>

Using $1 - \gamma \ge e^{-2\gamma}$ for $0 \le \gamma \le 1/4$, we get

$$ \begin{aligned} & f_3 \le K_3 + H_3[1 - e^{-2\gamma(L-4q)}] \\ & \le K_3 + 2H_3\gamma L. \end{aligned} \quad (F18) $$

Hence, if $\gamma \le \gamma_3/L$, we finally arrive at

$$ f_3 \le K_3/2 = -C_3. \quad (F19) $$

Analogously, we have

$$ \frac{[\gamma^4 + (1-\gamma)^4]^2}{[\gamma^3 + (1-\gamma)^3][\gamma^5 + (1-\gamma)^5]} = 1 - \epsilon_5, \quad (F20) $$

with $0 \le \epsilon_3 \le 2\gamma^3$ for $0 \le \gamma \le 1/4$. Therefore

$$ \begin{aligned} & f_5 \le K_5 + \tilde{s}_4^2 \tilde{t}_5 [1 - (1 - \epsilon_5)^{L-4q}] \\ & \le K_3 + H_5 [1 - (1 - 2\gamma^3)^{L-4q}]. \end{aligned} \quad (F21) $$

Using $1 - x \ge e^{-2x}$ for $0 \le x \le 1/4$, we get

$$ \begin{aligned} & f_5 \le K_5 + H_5[1 - e^{-4\gamma^3(L-4q)}] \\ & \le K_5 + 4H_5\gamma^3 L. \end{aligned} \quad (F22) $$

Hence, if $\gamma \le \gamma_5/L^{1/3}$, we arrive at

$$ f_5 \le K_5/2 = -C_5. \quad (F23) $$

<figure><img src="image_18.png" /></figure>

[1] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nature Phys. 8, 277 (2012).

[2] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).

[3] F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu, and Y. Takahashi, Tools for quantum simulation with ultracold atoms in optical lattices, Nature Rev. Phys. 2, 411 (2020).

[4] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Ann. Rev. Cond. Matt. Phys. 11, 369 (2020).

[5] M. Morgado and S. Whitlock, Quantum simulation and computing with rydberg-interacting qubits, AVS Quantum Science 3, 023501 (2021).

[6] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).

[7] Y. Alexeev, D. Bacon, K. R. Brown, R. Calderbank, L. D. Carr, F. T. Chong, B. DeMarco, D. Englund, E. Farhi, B. Fefferman, et al., Quantum computer systems for scientific discovery, PRX Quantum 2, 017001 (2021).

[8] E. Pelucchi, G. Fagas, I. Aharonovich, D. Englund, E. Figueroa, Q. Gong, H. Hannes, J. Liu, C.-Y. Lu, N. Matsuda, et al., The potential and global outlook of integrated photonics for quantum technologies, Nature Rev. Phys. 4, 194 (2022).

[9] G. Burkard, T. D. Ladd, A. Pan, J. M. Nichol, and J. R. Petta, Semiconductor spin qubits, Rev. Mod. Phys. 95, 025003 (2023).

[10] J. Preskill, Quantum computing in the nisq era and beyond, Quantum 2, 79 (2018).

[11] E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksen, K.-M. C. Fu, et al., Quantum simulators: Architectures and opportunities, PRX Quantum 2, 017003 (2021).

[12] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox, Nature Rev. Phys. 5, 9 (2023).

[13] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Statistical correlations between locally randomized mea-
surements: A toolbox for probing entanglement in many-body quantum states, Phys. Rev. A **99**, 052323 (2019).

[14] P. Cieśliński, S. Imai, J. Dziewior, O. Gühne, L. Knips, W. Laskowski, J. Meinecke, T. Paterek, and T. Vértesi, Analysing quantum systems with randomised measurements, arXiv:2307.01251 (2023).

[15] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Phys. **16**, 1050 (2020).

[16] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. **81**, 865 (2009).

[17] M. A. Nielsen and I. L. Chuang, *Quantum computation and quantum information* (Cambridge university press, Cambridge, UK, 2010).

[18] P. Calabrese and J. Cardy, Entanglement entropy and conformal field theory, J. Phys. A: Math. Theor. **42**, 504005 (2009).

[19] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Rev. Mod. Phys. **82**, 277 (2010).

[20] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature **528**, 77 (2015).

[21] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science **353**, 794 (2016).

[22] N. M. Linke, S. Johri, C. Figgatt, K. A. Landsman, A. Y. Matsuura, and C. Monroe, Measuring the rényi entropy of a two-site fermi-hubbard model on a trapped ion quantum computer, Phys. Rev. A **98**, 052334 (2018).

[23] C. Moura Alves and D. Jaksch, Multipartite entanglement detection in bosons, Phys. Rev. Lett. **93**, 110501 (2004).

[24] J. Cardy, Measuring entanglement using quantum quenches, Phys. Rev. Lett. **106**, 150404 (2011).

[25] A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, Measuring entanglement growth in quench dynamics of bosons in an optical lattice, Phys. Rev. Lett. **109**, 020505 (2012).

[26] D. A. Abanin and E. Demler, Measuring entanglement entropy of a generic many-body system with a quantum switch, Phys. Rev. Lett. **109**, 020504 (2012).

[27] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing rényi entanglement entropy via randomized measurements, Science **364**, 260 (2019).

[28] S. J. van Enk and C. W. J. Beenakker, Measuring $Trρ^n$ on single copies of $ρ$ using random measurements, Phys. Rev. Lett. **108**, 110503 (2012).

[29] A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac, and P. Zoller, Rényi entropies from random quenches in atomic hubbard and spin models, Phys. Rev. Lett. **120**, 050406 (2018).

[30] K. J. Satzinger, Y.-J. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, et al., Realizing topologically ordered states on a quantum processor, Science **374**, 1237-1241 (2021).

[31] R. Stricker, M. Meth, L. Postler, C. Edmunds, C. Ferrie, R. Blatt, P. Schindler, T. Monz, R. Kueng, and M. Ring-

bauer, Experimental single-setting quantum state tomography, PRX Quantum **3**, 040310 (2022).

[32] D. Zhu, Z. P. Cian, C. Noel, A. Risinger, D. Biswas, L. Egan, Y. Zhu, A. M. Green, C. H. Alderete, N. H. Nguyen, et al., Cross-platform comparison of arbitrary quantum states, Nature Comm. **13** (2022).

[33] J. C. Hoke, M. Ippoliti, E. Rosenberg, D. Abanin, R. Acharya, T. I. Andersen, M. Ansmann, F. Arute, K. Arya, A. Asfaw, et al., Measurement-induced entanglement and teleportation on a noisy quantum processor, Nature **622**, 481 (2023).

[34] Y. Zhou, P. Zeng, and Z. Liu, Single-copies estimation of entanglement negativity, Phys. Rev. Lett. **125**, 200502 (2020).

[35] A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, and B. Vermersch, Mixed-state entanglement from local randomized measurements, Phys. Rev. Lett. **125**, 200501 (2020).

[36] M. Cerezo, A. Sone, J. L. Beckey, and P. J. Coles, Sub-quantum fisher information, Quantum Sci. Tech. **6**, 035008 (2021).

[37] M. Yu, D. Li, J. Wang, Y. Chu, P. Yang, M. Gong, N. Goldman, and J. Cai, Experimental estimation of the quantum fisher information from randomized measurements, Phys. Rev. Res. **3**, 043122 (2021).

[38] A. Rath, C. Branciard, A. Minguzzi, and B. Vermersch, Quantum fisher information from randomized measurements, Phys. Rev. Lett. **127**, 260501 (2021).

[39] V. Vitale, A. Rath, P. Jurcevic, A. Elben, C. Branciard, and B. Vermersch, Estimation of the quantum fisher information on a quantum processor, arXiv:2307.16882 (2023).

[40] S. T. Flammia, D. Gross, Y.-K. Liu, and J. Eisert, Quantum tomography via compressed sensing: error bounds, sample complexity and efficient estimators, New J. Phys. **14**, 095022 (2012).

[41] J. Haah, A. W. Harrow, Z. Ji, X. Wu, and N. Yu, Sample-optimal tomography of quantum states, in *Proceedings of the forty-eighth annual ACM symposium on Theory of Computing* (2016) pp. 913–925.

[42] R. O'Donnell and J. Wright, Efficient quantum tomography, in *Proceedings of the forty-eighth annual ACM symposium on Theory of Computing* (2016) pp. 899–912.

[43] M. Paini and A. Kalev, An approximate description of quantum states, arXiv:1910.10543 (2019).

[44] C. Hadfield, Adaptive pauli shadows for energy estimation, arXiv:2105.12207 (2021).

[45] H.-Y. Huang, R. Kueng, and J. Preskill, Efficient estimation of pauli observables by derandomization, Phys. Rev. Lett. **127**, 030503 (2021).

[46] C. Hadfield, S. Bravyi, R. Raymond, and A. Mezzacapo, Measurements of quantum hamiltonians with locally-biased classical shadows, Comm. Math. Phys. **391**, 951 (2022).

[47] A. Rath, R. van Bijnen, A. Elben, P. Zoller, and B. Vermersch, Importance sampling of randomized measurements for probing entanglement, Phys. Rev. Lett. **127** (2021).

[48] B. Vermersch, A. Rath, B. Sundar, C. Branciard, J. Preskill, and A. Elben, Enhanced estimation of quantum properties with common randomized measurements, PRX Quantum **5**, 010352 (2024).
[49] T.-C. Yen, A. Ganeshram, and A. F. Izmaylov, Deterministic improvements of quantum measurements with grouping of compatible operators, non-local transformations, and covariance estimates, npj Quantum Inf. 9, 14 (2023).

[50] M. Arienzo, M. Heinrich, I. Roth, and M. Kliesch, Closed-form analytic expressions for shadow estimation with brickwork circuits, Quantum Inf. Comput 23, 961 (2023).

[51] A. A. Akhtar, H.-Y. Hu, and Y.-Z. You, Scalable and flexible classical shadow tomography with tensor networks, Quantum 7, 1026 (2023).

[52] C. Bertoni, J. Haferkamp, M. Hinsche, M. Ioannou, J. Eisert, and H. Pashayan, Shallow shadows: Expectation estimation using low-depth random clifford circuits, arXiv:2209.12924 (2022).

[53] M. Ippoliti, Y. Li, T. Rakovszky, and V. Khemani, Operator relaxation and the optimal depth of classical shadows, Phys. Rev. Lett. 130, 230403 (2023).

[54] M. Ippoliti, Classical shadows based on locally-entangled measurements, arXiv:2305.10723 (2023).

[55] G. García-Pérez, M. A. Rossi, B. Sokolov, F. Tacchino, P. K. Barkoutsos, G. Mazzola, I. Tavernelli, and S. Maniscalco, Learning to measure: Adaptive informationally complete generalized measurements for quantum algorithms, PRX Quantum 2, 040342 (2021).

[56] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).

[57] F. Tacchino, A. Chiesa, S. Carretta, and D. Gerace, Quantum computers as universal quantum simulators: state-of-the-art and perspectives, Ad. Quantum Tech. 3, 1900052 (2020).

[58] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature 607, 667 (2022).

[59] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Matrix product state representations, Quantum Inf. Comput. 7, 401 (2007).

[60] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).

[61] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Efficient quantum state tomography, Nature Comm. 1, 149 (2010).

[62] B. Lanyon, C. Maier, M. Holzäpfel, T. Baumgratz, C. Hempel, P. Jurcevic, I. Dhand, A. Buyskikh, A. Daley, M. Cramer, et al., Efficient tomography of a quantum many-body system, Nature Phys. 13, 1158 (2017).

[63] Z. Qin, C. Jameson, Z. Gong, M. B. Wakin, and Z. Zhu, Stable tomography for structured quantum states, arXiv:2306.09432 (2023).

[64] T. Baumgratz, D. Gross, M. Cramer, and M. B. Plenio, Scalable reconstruction of density matrices, Phys. Rev. Lett. 111, 020401 (2013).

[65] T. Baumgratz, A. Nüßeler, M. Cramer, and M. B. Plenio, A scalable maximum likelihood method for quantum state tomography, New J. Phys. 15, 125004 (2013).

[66] M. Holzäpfel, M. Cramer, N. Datta, and M. B. Plenio, Petz recovery versus matrix reconstruction, J. Math. Phys. 59, 042201 (2018).

[67] A. Lidiak, C. Jameson, Z. Qin, G. Tang, M. B. Wakin, Z. Zhu, and Z. Gong, Quantum state tomography with tensor train cross approximation, arXiv:2207.06397 (2022).

[68] C. Kokail, R. van Bijnen, A. Elben, B. Vermersch, and P. Zoller, Entanglement hamiltonian tomography in quantum simulation, Nature Phys. 17, 936 (2021).

[69] M. K. Joshi, C. Kokail, R. van Bijnen, F. Kranzl, T. V. Zache, R. Blatt, C. F. Roos, and P. Zoller, Exploring large-scale entanglement in quantum simulation, arXiv:2306.00057 (2023).

[70] A. Anshu, S. Arunachalam, T. Kuwahara, and M. Soleimanifar, Sample-efficient learning of interacting quantum systems, Nature Phys. 17, 931 (2021).

[71] A. Anshu and S. Arunachalam, A survey on the complexity of learning quantum states, arXiv:2305.20069 (2023).

[72] C. Rouzé and D. S. França, Learning quantum many-body systems from a few copies, arXiv:2107.03333 (2021).

[73] E. Onorati, C. Rouzé, D. S. França, and J. D. Watson, Efficient learning of ground & thermal states within phases of matter, arXiv:2301.12946 (2023).

[74] G. Tóth, W. Wieczorek, D. Gross, R. Krischek, C. Schwemmer, and H. Weinfurter, Permutationally invariant quantum tomography, Phys. Rev. Lett. 105, 250403 (2010).

[75] T. Moroder, P. Hyllus, G. Tóth, C. Schwemmer, A. Niggebaum, S. Gaile, O. Gühne, and H. Weinfurter, Permutationally invariant state reconstruction, New J. Phys. 14, 105001 (2012).

[76] C. Schwemmer, G. Tóth, A. Niggebaum, T. Moroder, D. Gross, O. Gühne, and H. Weinfurter, Experimental comparison of efficient tomography schemes for a six-qubit state, Phys. Rev. Lett. 113, 040503 (2014).

[77] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, Quantum state tomography via compressed sensing, Phys. Rev. Lett. 105, 150401 (2010).

[78] A. Montanaro, Learning stabilizer states by bell sampling, arXiv:1707.04012 (2017).

[79] S. S. Kuzmin, V. I. Mikhailova, I. V. Dyakonov, and S. S. Straupe, Learning a few tensor network model of a quantum state using a few single-qubit measurements, arXiv:2309.00397 (2023).

[80] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Neural-network quantum state tomography, Nature Phys. 14, 447 (2018).

[81] L. Zhao, N. Guo, M.-X. Luo, and P. Rebentrost, Probable learning of quantum states with graphical models circuits, arXiv:2309.09235 (2023).

[82] M. Rieger, M. Reh, and M. Gärttner, Sample-efficient estimation of entanglement entropy through supervised learning, Phys. Rev. A 109, 012403 (2024).

[83] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, Reconstructing quantum states with generative models, Nature Machine Intell. 1, 155 (2019).

[84] T. Schmale, M. Reh, and M. Gärttner, Efficient quantum state tomography with convolutional neural networks, npj Quantum Inf. 8, 115 (2022).

[85] D. J. Luitz, N. Laflorencie, and F. Alet, Participation spectroscopy and entanglement hamiltonian of quantum spin models, J. Stat. Mech. 2014, P08007 (2014).

[86] J.-M. Stéphan, S. Furukawa, G. Misguich, and V. Pasquier, Shannon and entanglement entropies of one- and two-dimensional critical wave functions, Phys. Rev. B 80, 184421 (2009).
[87] J.-M. Stéphan, G. Misguich, and V. Pasquier, Rényi entropy of a line in two-dimensional ising models, Phys. Rev. B 82, 125455 (2010).

[88] F. C. Alcaraz and M. A. Rajabpour, Universal behavior of the shannon mutual information of critical quantum chains, Phys. Rev. Lett. 111, 017201 (2013).

[89] J.-M. Stéphan, Shannon and rényi mutual information in quantum critical spin chains, Phys. Rev. B 90, 045424 (2014).

[90] P. Sierant and X. Turkeshi, Universal behavior beyond multifractality of wave functions at measurement-induced phase transitions, Phys. Rev. Lett. 128, 130605 (2022).

[91] X. Turkeshi and P. Sierant, Error-resilience phase transitions in encoding-decoding quantum circuits, arXiv:2308.06321 (2023).

[92] L. Leone, S. F. E. Oliviero, and A. Hamma, Stabilizer rényi entropy, Phys. Rev. Lett. 128, 050402 (2022).

[93] S. F. E. Oliviero, L. Leone, and A. Hamma, Magic-state resource theory for the ground state of the transverse-field ising model, Phys. Rev. A 106, 042426 (2022).

[94] T. Haug and L. Piroli, Quantifying nonstabilizerness of matrix product states, Phys. Rev. B 107, 035148 (2023).

[95] T. Haug and L. Piroli, Stabilizer entropies and nonstabilizerness monotones, Quantum 7, 1092 (2023).

[96] G. Lami and M. Collura, Nonstabilizerness via perfect pauli sampling of matrix product states, Phys. Rev. Lett. 131, 180401 (2023).

[97] P. S. Tarabunga, E. Tirrito, T. Chanda, and M. Dalmonte, Many-body magic via pauli-markov chains-from criticality to gauge theories, arXiv:2305.18541 (2023).

[98] P. Niroula, C. D. White, Q. Wang, S. Johri, D. Zhu, C. Monroe, C. Noel, and M. J. Gullans, Phase transition in magic with random quantum circuits, arXiv:2304.10481 (2023).

[99] X. Turkeshi, M. Schirò, and P. Sierant, Measuring nonstabilizerness via multifractal flatness, Phys. Rev. A 108, 042408 (2023).

[100] G. Vidal and R. F. Werner, Computable measure of entanglement, Phys. Rev. A 65, 032314 (2002).

[101] M. B. Plenio, Logarithmic negativity: A full entanglement monotone that is not convex, Phys. Rev. Lett. 95, 090503 (2005).

[102] H. A. Carteret, Noiseless quantum circuits for the peres separability criterion, Phys. Rev. Lett. 94, 040502 (2005).

[103] J. Gray, L. Banchi, A. Bayat, and S. Bose, Machine-learning-assisted many-body entanglement measurement, Phys. Rev. Lett. 121, 150503 (2018).

[104] H. A. Carteret, Estimating the entanglement negativity from low-order moments of the partially transposed density matrix, arXiv:1605.08751 (2016).

[105] A. Neven, J. Carrasco, V. Vitale, C. Kokail, A. Elben, M. Dalmonte, P. Calabrese, P. Zoller, B. Vermersch, R. Kueng, et al., Symmetry-resolved entanglement detection using partial transpose moments, npj Quantum Inf. 7, 152 (2021).

[106] X.-D. Yu, S. Imai, and O. Gühne, Optimal entanglement certification from moments of the partial transpose, Phys. Rev. Lett. 127, 060504 (2021).

[107] J. Carrasco, M. Votto, V. Vitale, C. Kokail, A. Neven, P. Zoller, B. Vermersch, and B. Kraus, Entanglement phase diagrams from partial transpose moments, arXiv:2212.10181 (2022).

[108] L. Knips, J. Dziewior, W. Kłobus, W. Laskowski, T. Paterek, P. J. Shadbolt, H. Weinfurter, and J. D. Meinecke, Multipartite entanglement analysis from random correlations, npj Quantum Inf. 6, 51 (2020).

[109] A. Ketterer, N. Wyderka, and O. Gühne, Characterizing multipartite entanglement with moments of random correlations, Phys. Rev. Lett. 122, 120505 (2019).

[110] L. Piroli, G. Styliaris, and J. I. Cirac, Quantum circuits assisted by local operations and classical communication: Transformations and phases of matter, Phys. Rev. Lett. 127, 220503 (2021).

[111] D. Malz, G. Styliaris, Z.-Y. Wei, and J. I. Cirac, Preparation of matrix product states with log-depth quantum circuits, arXiv:2307.01696 (2023).

[112] X. Chen, Z.-C. Gu, and X.-G. Wen, Local unitary transformation, long-range quantum entanglement, wave function renormalization, and topological order, Phys. Rev. B 82, 155138 (2010).

[113] M. B. Hastings, Classifying quantum phases with the kirby torus trick, Phys. Rev. B 88, 165114 (2013).

[114] B. Zeng, X. Chen, D.-L. Zhou, and X.-G. Wen, Quantum information meets quantum matter-from quantum entanglement to topological phase in many-body systems, arXiv:1508.02595 (2015).

[115] B. Zeng and X.-G. Wen, Gapped quantum liquids and topological order, stochastic local transformations and emergence of unitarity, Phys. Rev. B 91, 125121 (2015).

[116] A. Rath, V. Vitale, S. Murciano, M. Votto, J. Dubail, R. Kueng, C. Branciard, P. Calabrese, and B. Vermersch, Entanglement barrier and its symmetry resolution: Theory and experimental observation, PRX Quantum 4, 010318 (2023).

[117] D. Poulin and M. B. Hastings, Markov entropy decomposition: A variational dual for quantum belief propagation, Phys. Rev. Lett. 106, 080403 (2011).

[118] K. Kato and F. G. Brandao, Quantum approximate markov chains are thermal, Comm. Math. Phys. 370, 117 (2019).

[119] F. G. Brandão and M. J. Kastoryano, Finite correlation length implies efficient preparation of quantum thermal states, Comm. Math. Phys. 365, 1 (2019).

[120] I. H. Kim, Markovian matrix product density operators: Efficient computation of global entropy, arXiv:1709.07828 (2017).

[121] P. Svetlichnyy, S. Mittal, and T. Kennedy, Matrix product states and the decay of quantum conditional mutual information, arXiv:2211.06794 (2022).

[122] P. Svetlichnyy and T. Kennedy, Decay of quantum conditional mutual information for purely generated finitely correlated states, J. Math. Phys. 63, 072201 (2022).

[123] D. Haag, F. Baccari, and G. Styliaris, Typical correlation length of sequentially generated tensor network states, PRX Quantum 4, 030330 (2023).

[124] C.-F. Chen, K. Kato, and F. G. Brandão, Matrix product density operators: When do they have a local parent hamiltonian?, arXiv:2010.14682 (2020).

[125] I. H. Kim, Entropy scaling law and the quantum marginal problem, Phys. Rev. X 11, 021039 (2021).

[126] I. H. Kim, Entropy scaling law and the quantum marginal problem: simplification and generalization, arXiv:2109.11688 (2021).

[127] T. Kuwahara, K. Kato, and F. G. S. L. Brandão, Clustering of conditional mutual information for quantum gibbs states above a threshold temperature, Phys. Rev.
Lett. 124, 220601 (2020).

[128] T. Kuwahara, A. M. Alhambra, and A. Anshu, Improved thermal area law and quasilinear time algorithm for quantum gibbs states, Phys. Rev. X 11, 011047 (2021).

[129] P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, The tensor networks anthology: Simulation techniques for many-body quantum lattice systems, SciPost Phys. Lect. Notes, 8 (2019).

[130] K. M. Audenaert, A sharp continuity estimate for the von neumann entropy, J. Phys. A: Math. Theor. 40, 8127 (2007).

[131] Z. Chen, Z. Ma, I. Nikoufar, and S. Fei, Sharp continuity bounds for entropy and conditional entropy, Science China, Phys. Mech. Astr. 60, 020321 (2017).

[132] M. A. Nielsen, Continuity bounds for entanglement, Phys. Rev. A 61, 064301 (2000).

[133] T. Farrelly, A review of quantum cellular automata, Quantum 4, 368 (2020).

[134] L. Piroli and J. I. Cirac, Quantum cellular automata, tensor networks, and area laws, Phys. Rev. Lett. 125, 190402 (2020).

[135] J. I. Cirac, D. Perez-Garcia, N. Schuch, and F. Verstraete, Matrix product density operators: Renormalization fixed points and boundary theories, Ann. Phys. 378, 100 (2017).

[136] S. Sachdev, *Quantum phase transitions*, Vol. 12 (Cambridge University Press, Cambridge (UK), 1999).

[137] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, *Quantum inverse scattering method and correlation functions*, Vol. 3 (Cambridge University Press, Cambridge (UK), 1997).

[138] F. Verstraete and J. I. Cirac, Matrix product states represent ground states faithfully, Phys. Rev. B 73, 094423 (2006).

[139] M. Fishman, S. R. White, and E. M. Stoudenmire, The itensor software library for tensor network calculations, SciPost Phys. Codebases, 4 (2022).

[140] A. J. Ferris and G. Vidal, Perfect sampling with unitary tensor networks, Phys. Rev. B 85, 165146 (2012).

[141] A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, Cross-platform verification of intermediate scale quantum devices, Phys. Rev. Lett. 124, 010504 (2020).

[142] H. T. Quan, Z. Song, X. F. Liu, P. Zanardi, and C. P. Sun, Decay of loschmidt echo enhanced by quantum criticality, Phys. Rev. Lett. 96, 140604 (2006).

[143] P. Gondolf, S. O. Scalet, A. Ruiz-de Alarcon, A. M. Alhambra, and A. Capel, Conditional independence of 1d gibbs states with applications to efficient learning, arXiv:2402.18500 (2024).

[144] S. Rana, Negative eigenvalues of partial transposition of arbitrary bipartite states, Phys. Rev. A 87, 054301 (2013).