{'fraction_non_alphanumeric': 0.03428996714407595, 'fraction_numerical': 0.013814010056454102, 'mean_word_length': 4.391466040321002, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Entanglement estimation in tensor network states via sampling

We introduce a method for extracting meaningful entanglement measures of tensor network states in general dimensions. Current methods require the explicit reconstruction of the density matrix, which is highly demanding, or the contraction of replicas, which requires an effort exponential in the number of replicas and which is costly in terms of memory. In contrast, our method requires the stochastic sampling of matrix elements of the classically represented reduced states with respect to random states drawn from simple product probability measures constituting frames. Even though not corresponding to physical operations, such matrix elements are straightforward to calculate for tensor network states, and their moments provide the R\'enyi entropies and negativities as well as their symmetry-resolved components. We test our method on the one-dimensional critical XX chain and the two-dimensional toric code in a checkerboard geometry. Although the cost is exponential in the subsystem size, it is sufficiently moderate so that - in contrast with other approaches - accurate results can be obtained on a personal computer for relatively large subsystem sizes.


I. INTRODUCTION
Entanglement is the key feature of quantum mechanics that renders it different from classical theories. It takes centre stage in quantum information processing where it plays the role of a resource. The significance of notions of entanglement for capturing properties of condensed matter systems has also long been noted and appreciated 1,2 . The observation that ground states of gapped phases of matter are expected to feature little entanglement -in fact, they feature what are called area laws for entanglement entropies 3 -is at the basis of tensor networks (TN) methods 4,5 accurately describing interacting quantum many-body systems. It has been noted that certain scalings of entanglement measures can indicate the presence of quantum phase transitions 6,7 . Indeed, the very fact that locally interacting quantum many body systems tend to be much less entangled than they could possibly be renders TN methods a powerful technique to capture their properties [3][4][5]8 . Maybe most prominently, topologically ordered systems can be regarded as long ranged entangled systems 9 . In addition, detailed information about the scaling of entanglement properties can provide substantial diagnostic information about properties of condensed matter systems.
Accepting that tensor network states often provide an accurate and efficient classical description of interacting quantum systems, the question arises how one can meaningfully read off entanglement properties from tensor network states. This, however, constitutes a challenge. Current entanglement calculation methods in tensor network states in two and higher dimensions are highly impractical even for moderate-size systems, since they require a full reconstruction of the quantum states at heavy computational costs. For Rényi entropies one may instead employ the replica trick, which uses several copies of the RDM (as explained in Section II B below); this, however, comes with an exponential scaling of the computational effort in the number of copies, often making the calculation unfeasible.
Here, we bring these ideas to a new level by applying them to quantum states that are classically represented in the first place by tensor networks. The core idea of these methods is basically the following: While the entanglement moments naively require several copies of the system, we can refrain from this requirement by resorting to random sampling. The general protocol is to evolve the system under a random unitary drawn from the Haar measure followed by a measurement of a suitable projector. The process is repeated and moments of the results are averaged over different unitaries, giving as a result entanglement moments or other density-matrixbased properties. The effectiveness of this mindset has been demonstrated experimentally in a number of platforms, including that of cold atoms for Rényi entropies 16,22 and Rényi negativities 24,26 .
While these ideas have been further developed into estimation techniques 23 giving rise to classical representations in their own right, we turn these ideas upside down by applying them to quantum systems that are already classically represented by tensor networks. There are some crucial differences that arise in classical representations compared to quantum experiments: First, they are much more suitable for a direct calculation of expectation values, rather than estimating them from sampling from measurements. Second, and importantly, when performing a classical simulation, we are not limited to physically-allowed actions, and specifically, we are not constrained to the application of unitary operators and measurements. This feature is to an extent reminiscent of shadow estimation in that also there, unphysical maps are made use of. It is the estimation procedure itself that is not physical here, however. The method we develop allows for having only a single copy of the simulated state, and at the same time for estimating the entanglement moments based on matrix elements that are naturally calculated. Instead of sampling operators from the Haar measure or some unitary n-design 27 , we only need to sample from a simple, finite set of tensor products of independent d-dimensional vectors -specifically from what are called frames or spherical complex 1-designs 28 , where d is the Hilbert space of a single site in the system. Furthermore, this simple structure allows our method to be applied to arbitrary system geometries.
The remainder of this work is organized as follows. Section II includes preliminary theoretical background. The Rényi moments we aim to estimate are defined and their relation to standard entanglement measures is discussed in Section II A. Section II B covers the basic ideas of the TN ansatzes we use in our work: For one-dimensional systems, the matrix product state (MPS) ansatz, and in higher dimensions, projectedentangled-paired-states (PEPS) and its infinite system size version known as iPEPS. We discuss the algorithms we used for extracting the reduced density matrix and the naive method for estimating entanglement moments of states represented by these ansatzes. The solvable models used as benchmarks for testing our method are presented in Section II D. In Section III we explain our proposed algorithm for using random variables for estimating the entanglement moments of TN in general dimension, and study the variance of the estimate in Section IV B, from which arises the complexity of an estimation up to a chosen error. We benchmark the algorithm with the ground states of the exactly solvable toric code model, Eq. (13), using iPEPS, and the XX chain, Eq. (16), using MPS, in Section IV. Finally, we discuss the results and future steps in the conclusions, Section V. In the appendix, we present variance estimations of the Rényi moments in the general case (Appendix A), as well as specifically in the toric code model, which is used as a benchmark (Appendix B).

II. PRELIMINARIES
A. Entanglement measures For a quantum system in a pure state ρ = |ψ ψ|, we define for a subsystem the reduced quantum state or reduced density matrix (RDM) as The entanglement of the subsystem A with its environment (consitutung its complement) A is encoded in the RDM. We introduce the moments of the RDM. For a positive integer n, the n-th RDM moment is defined to be also referred to as the Rényi moments. On top of being entanglement monotones 29 and hence measures of entanglement in their own right, these moments are used for defining various entanglement measures with useful mathematical properties [30][31][32][33][34][35][36][37] . The RDM moments can -under mild mathematical conditions -be analytically continued to the entanglement measure featuring the strongest interpretation for pure bi-partite quantum states, the von Neumann entanglement entropy 38 defined as for RDMs ρ A , as the 1-Rényi entropy. The von Neumann entropy is obtained in the limit S(ρ A ) = lim n→1 (1 − n) −1 log (p n (ρ A )). The Rényi moments are especially popular as entanglement indicators since they do not require a full reconstruction of the RDM spectrum. Therefore, they are often easier to either calculate theoretically or measure experimentally than other entanglement measures [12][13][14]16,20,22,[39][40][41][42] . The measures above are appropriate when quantifying the entanglement between a subsystem A and its environment A when A ∪ A is in a pure state. When characterizing the entanglement between two subsystems A 1 and A 2 whose union A is not necessarily pure, the quantities above will no longer be suitable to quantify entanglement, as they cannot distinguish between the quantum entanglement between A 1 and A 2 from their entanglement with the environment. One of the best known measures for the entanglement between two subsystems labeled as A 1 and A 2 is the entanglement negativity 35,43,44 , based on the positive partial transpose (PPT) criterion 2,45,46 where · 1 denotes the trace norm, and ρ T2 A the partial transposition of the degrees of freedom corresponding to A 2 in ρ A , for all vectors (|i , |k ), (|j , |l ) in an orthonormal basis of the Hilbert spaces of A 1 , A 2 , respectively. The usefulness of the negativity as an entanglement measure for two subsystems in a mixed state 35,44 leads us to define the moments of the partially-traced RDM, further referred to as PT moments. The n-th PT moment is defined to be for a positive integer n. The negativity can be obtained by an analytic continuation of the PT even integer moments by ρ T2 A 1 = lim n→1/2 R 2n (ρ A ). The PT moments are not entanglement monotones, but they can be used to detect entanglement between A 1 and A 2 24,26 , as well as for estimating the negativity 47 . The popularity of the PT moments as entanglement indicators stems from the fact that they too do not require a full reconstruction of the partially transposed RDM, and are therefore easier to calculate and measure experimentally 24,42 .
For systems with a conserved charge Q, the quantum state of the full system typically commutes with the charge operator, A partial trace can be applied to the permutation relation above to give [ρ A ,Q A ] = 0, whereQ A is the charge operator on subsystem A. The RDM is thus composed of blocks, each corresponding to a charge value q in subsystem A, as illustrated in the inset of Fig. 6. We denote the RDM block corresponding to charge q by ρ A (q). The entanglement measures, and specifically the RDM moments, can then be decomposed into suitable contributions from the different blocks called charge-resolved or symmetry-resolved moments [48][49][50][51] , again for positive integers n. This definition could be extended to the negativity as well 52 . The study of symmetryresolved entanglement has drawn much interest lately, both analytically and numerically 26,[53][54][55][56][57][58][59][60][61][62][63][64][65][66] as well as in the development of experimental measurement protocols 42,49,52,67 . It reveals the relation between entanglement and charge and can point to effects such as topological phase transitions 53,55,61 or to instances of dissipation in open systems dynamics 65 .
The estimation of symmetry-resolved entanglement can be done based on the analysis in Ref. 49: We introduce the socalled flux-resolved RDM moments as where ϕ ∈ [0, 2π) can be thought of as an Aharonov-Bohm flux inserted in the replica trick. The symmetry-resolved moments can be extracted from the flux-resolved moments by a Fourier transform according to where ϕ = 2πk/N A and k = 0, . . . , N A − 1.

B. Tensor networks
We will now briefly review the TN tools that are being used to compute the quantities presented above in the remainder of this work by means of sampling techniques.

Matrix product states
We here consider a one-dimensional spin system, so of a finite local dimension, featuring N lattice sites. The state vector of the system can be written as where Ψ is the rank N tensor of the coefficients of |ψ . Ψ has d N complex amplitudes, where d is the Hilbert space size of a single spin. In an MPS representation we decompose Ψ into N different tensors, each corresponding to a single site, as illustrated in Fig. 1a. Each such tensor will have a single index corresponding to the indices of the original tensor, often called the 'physical leg' or the 'physical index', and two additional indices connecting with the tensors corresponding to the site's neighbours, often called 'entanglement legs' or 'bond indices' (with only one bond index for the sites at the edges). Contracting all the bond indices will result in the original tensor Ψ. In many cases one can limit the dimension of the bond index to be some chosen constant D, also known as the bond dimension, and discard the least significant variables. In this way, the number of real parameters will be scaling as O(N D 2 d), at the cost of getting an approximate representation for the state. States that are expected to be well approximated by such limited tensors obey an entanglement area-law [3][4][5]8 (in fact, this is provably the case for area laws of suitable Rényi entropies 69 ). This MPS decomposition is a widely used method for the simulation of ground states 70,71 , thermal states [72][73][74] and states undergoing a time evolution [75][76][77][78] generated by local Hamiltonians of one-dimensional systems.
For a system partitioned into two contiguous subsystems, the extraction of the spectrum of the RDM, also called the entanglement spectrum, is very natural 79 and can often be useful in classifying phases of matter in one spatial dimension [80][81][82][83] . We note that a decomposition of the system into two tensors, one corresponding to subsystem A and one to A, is built into the decomposition of the systems into site tensors, and that this decomposition can be transformed into the Schmidt decomposition of the state vector The RDM eigenvalues are thus the squared absolute values of the Schmidt values, and by obtaining them, we can extract the RDM moment in all ranks n, as well as the von Neumann entropy. Specific techniques have been developed for the extraction of entanglement measures in some additional cases, such as the entanglement of a contiguous subsystem of an infinite system 84 or the negativity of two contiguous subsystems 85 .

Projected entangled pair states
For two or higher dimensional lattice systems, the MPS formalism is extended to an ansatz called projected entangled pair states (PEPS) 86,87 . The tensor capturing the state vector of the entire lattice is then decomposed into N site tensors, each with a single physical index and a bond index for each neighbour of the site in the system. An example for a square lattice is depicted in Fig. 1b, and the generalization to other lattice configurations is straightforward.
The infinite version of PEPS, known as iPEPS 87 , can be used to represent states in the thermodynamic limit in 2D. They are defined by a finite set of tensors repeated all over the lattice with some periodicity. iPEPS have found numerous applications in studying ground states [88][89][90][91][92][93] , thermal states [94][95][96] and non-equilibrium problems [97][98][99][100][101][102] in two spatial dimensions, and have become state of the art numerical technique for studying strongly correlated two-dimensional problems. The technique does not suffer from the infamous sign 2 Figure 1. (a-b) A general pure quantum state can be thought of as a rank N tensor represented by a single tensor with N legs (indices), each with degree d. The full tensor capturing the entire quantum state can be decomposed into N tensors, each representing a single site. The legs of the original tensor (in pink) are now divided between all tensors, and will be referred to as the physical index. Each tensor has additional legs connecting it to its neighbours, which will be referred to as the bond index. Contracting all the bond indices will result in the original tensor. (a) and (b) represent states in a one and two spatial dimensions, respectively. (c) Representing the density matrix of a pure state by taking two copies and applying complex conjugation to one of them, which is indicated by the physical index facing downwards. (d) The RDM of the sites in blue is obtained by tracing out the degrees of freedom of the sites in grey, which in turn is obtained by contracting the physical legs of both copies. (e) The RDM can be represented by the site tensors of the sites in A from both copies, which is here represented by a single tensor with two physical legs, and boundary tensors, which have only bond indices. In order to obtain the boundary tensors, we use the boundary MPS method 68 . The dimension of the bond indices in black is D 2 , while that of the legs in green, connecting the boundary tensors to each other, is χ, the bond dimension of the environment. problem 103,104 and can go to very large system sizes, thus allowing access to regimes where techniques like Quantum Monte Carlo and exact diagonalization fail.
The pure quantum state ρ = |ψ ψ| of the quantum system can be obtained by taking a PEPS state vector and its Hermitian conjugate and placing them back to back as a tensor product, as depicted for PEPS in Fig. 1c. We now examine a rectangular subsystem A with N A = w 1 × w 2 sites, where w 1 ≤ w 2 (as will be the notation throughout this work). In order to get the RDM of A as defined in Eq. (1), the degrees of freedom of A need to be traced out. This can be obtained by contracting the physical legs of all tensors corresponding to sites in A with the physical legs of the same tensor in the complex conjugate. We get a RDM composed of site tensors for the tensors in A, and boundary tensors resulting from the tensors in A, as depicted in Fig. 1e. Such boundary tensors can be approximately computed for an infinite system. We remark here that exactly contracting PEPS tensors is a computationally hard problem (in worst case complexity and for meaningful probability measures even in average case) 105,106 and therefore, we have to rely on approximation algorithms such as the corner transfer matrix renormalization group algorithm 107,108 boundary MPS techniques 68 , higher order tensor renormalization group methods 109 or others. It is also known that those PEPS that are ground states of uniformly gapped parent Hamiltonians -which are interesting in the condensed matter context -can actually be contracted in quasi-polynomial time 110 . In this work, we make use of the boundary MPS technique: We create a one-dimensional TN representing the boundary of the (supposedly infinite) system, and multiply it by the 'traced out' tensors indicated in Fig. 1d. The boundary bond dimension is limited to a constant dimension χ. This process is then repeated until the one-dimensional boundary tensors are converged, resulting in a one-dimensional boundary as depicted in Fig. 1e.

C. Entanglement measures computed from reduced states
The entanglement measures presented in Section II A can be extracted for two contiguous systems in an MPS as presented in Section II B 1, as well as in additional specific cases in one 84,85 and two 111 spatial dimensions. However, for a general dimension and partition, there is no efficient way known to quantify the entanglement. A straightforward method can be contracting the tensors such that the RDM is obtained explicitly to then obtain its spectral decomposition. However, the explicit RDM is of dimension d N A × d N A , which comes along with substantial computational effort and which imposes a strong restriction on the accessible system sizes.
That being said, the n-th RDM or PT moments defined in Eqs. (2,5) can be calculated in polynomial time in N A using n copies of the system tensors as depicted in Fig. 2 for a two-dimensional PEPS. The space complexity required for performing this multiplication for an MPS scales as O d 2n D 2n + D 4n . The space complexity for a twodimensional PEPS is given by O χ 2n D 2(w1+1)n + D 8n d , where χ is the bond dimension of the environment as depicted in Fig. 1e, and w 1 is the short edge of the rectangular system, as define in Section II B 2 above. For n > 1, the exponential dependence of the cost on n quickly makes it prohibitively large, despite the fact that for a narrow system (constant w 1 ), The exponential dependence of the dimension of ρA on NA turns into a linear dependence. However, a prohibitive exponential dependence on n emerges instead. the time complexity is linear in N A and the space complexity does not depend on N A (except for the possible dependence of D on N A , as can sometimes happen in finite systems).

D. Benchmark models
Before turning to presenting the actual sampling method for computing entanglement measures in quantum systems captured by tensor networks, we here first present the models we used for benchmarking our method: The two-dimensional gapped toric code model on a square lattice and the onedimensional gapless XX model.

The toric code model
The first benchmark model we elaborate on is the analytically solvable toric code model on a square lattice. The toric code, introduced by Kitaev 112 , transferring insights from topological quantum field theory to the realm of quantum spin systems, is a model of spins on a square lattice with local dimension d = 2. The spins live on the edges of the lattice rather than its nodes. The Hamiltonian of the model is given by s in the equation above represents the set of edges around a single node in the lattice (a star) and p represents the set of edges forming a plaquette in the lattice, as shown in Fig. 3. The ground state of toric code model displays several important properties, among which are topological order, which leads to robustness to local errors, making it an important candidate for fault-tolerant error correction code. For J s , J p > 0, the ground state vector of the model with open boundary conditions is known and can be written as where N is the number of all sites in the system. In the limit N → ∞, the iPEPS representation of the infinite toric code ground state is given in Refs. 5 and 113. A set of two site tensors, T A and T B , are repeated infinitely such that all of the nearest-neighbours of a site represented by T A are of the form T B and vice versa. The bond dimension of all bond indices of For a subsystem of the infinite system in the state defined in Eq. (14), the density-matrix-based measures can be analytically calculated 114 . This relies on the symmetry of the ground state under the application of ⊗ i∈s σ x i for all stars s and of ⊗ i∈p σ z i for all plaquettes p. Due to this symmetry, the RDM is block diagonal, where the size of each block equals the order of the group generated by each operator, which is 2 for the operators above. Considering the fact that all non-zero blocks are identical, as can be seen from Eq. (14), the eigenvalues of the RDM can be extracted analytically. Note that the symmetry mentioned above is not utilized in the numerical method, so as to make our performance results applicable to general analytically-unsolvable models, which do not posses such local symmetries. Here, we estimate the 2nd, 3rd and 4th RDM moments, as well as the 3rd PT moment, for a checkerboardlike partition of a square subsystem (w 1 = w 2 = w), as shown for w = 6 in Fig. 3. We study the cases of w = 4, 6, 8. Note that while the toric code displays an area law type entanglement structure, here the entire system is in the area and therefore a volume law is reached. Such extensive partitions were shown to be interesting for the study of topological phases in Refs. 115,116 and following works (a more traditional geometry is studied in Appendix B). For such systems, the n-th RDM moment is shown to be 114 The log of the moment deviates from an area law by an additive constant term, reflecting the the topological order of the model 117,118 . Note that for the toric code ground state p 3 = R 3 due to the structure of the RDM discussed above. However, we compute an estimate for R 3 based on the generallyapplicable estimator defined in Eq. (25) for completeness.

The XX model
While suited for high dimensions, we note that our method is blind to the dimensionality of the system, and will apply to one-dimensional systems in precisely the same way as it would for higher dimensions. Therefore, we can use one-dimensional models as benchmark models for testing the system. We test our model on the ground state of the onedimensional XX model captured by the local Hamiltonian where i stands for a site in the system. The Hamiltonian can be seen as a Hamiltonian of non-interacting fermions by virtue of the Jordan-Wigner transformation 119 and is thus analytically solvable 120 . We compute the ground state of a system of length 2l and extract the 2nd, 3rd and 4th RDM moments of a contiguous half of the system. In contrast to the toric code, the XX model is gapless in the absence of a large magnetic field and can be well approximated as a conformal system. For such systems, the RDM moments of a subsystem when the total system is in the ground state is to a good approximation predicted to be 121,122 , where c = 1 is the conformal charge. As opposed to the toric code ground state RDM, which is composed of 2 × 2 blocks, the XX ground state is not as structured, and the performance of the method is harder to predict. As such, the XX model ground state is a good complement to the toric code ground state in the study of the method's performance.

III. METHOD
We now turn to describing the method that is at the heart of this work. The core idea is that with suitable stochastic sampling techniques, one can more resource-efficiently estimate entanglement properties of systems captured by tensor networks. Inspired by the growing body of methods based on random unitaries, 10-23 and described in Section I, we now turn to describe our random-variables-based method for estimating the entanglement contained in a TN state. As mentioned in the introduction, our method differs from the protocols that are routinely implemented in experiments by two key aspects: First, we do not base the protocol on sampling local measurement results, which are cumbersome to extract from TNs, but on sampling of expectation values, which can be naturally calculated in TNs. Note that while actual sampling from MPS can be performed exactly 123 , sampling from PEPS is shown to be computationally hard, in the worst as well as average case 124,125 . The second difference between our method and the experimental protocols is that we do not have to limit ourselves to physically allowed processes, and specifically, our random operations are neither unitary nor quantum channels, which allows for a significant simplification of the protocol.

A. Sampling random vectors
In what follows, random vectors |v ∈ C d drawn from appropriate probability measures will feature with the property such that where E refers to the average over the chosen probability measure. This is up to the normalization that is only different by a factor of d 1/2 than what is commonly called a frame or a spherical complex 1-design 27,28 . This convention is helpful in what follows. The set of vectors can be a discrete or a continuous set. We use the random variable to get simple estimators for the entanglement quantifiers based on Rényi moments of RDM of subsystems A, each site α ∈ A of which corresponds to a system of local dimension d. In this setting, we consider random vectors where |v [α] ∈ C d are vectors drawn in an i.i.d. fashion as in Eq. (18), one for each site α. Naturally still holds true in this multi-partite setting.

B. Estimators of entanglement measures
By applying these random vectors to the reduced density matrix, we obtain an estimator of the second entanglement moment, also referred to as the purity, from expressions of the formp Indeed, averaging over the (independent) random vectors |V , |V drawn from a product probability measure as defined in Eq. (19), we consistently obtain the second moment as defined in Eq. (2) as the expectation . Note again that these quantities can be readily computed at hand of the classical description of the quantum state, but cannot be natively measured in a quantum system. In this sense, the random sampling technique proposed here resorts to 'unphysical operations'.
The n-th RDM moment can be obtained by a generalization of Eq. (21) as the expectation E(p n (ρ A )) = p n (ρ A ) of where |V (1) , . . . , |V (n) are drawn in an i.i.d. fashion from the same probability measure. For the PT moments, perform the partial transposition with respect to subsystem A 2 . Specifically, we define, for i, j = 1, . . . , n, pairs of product vectors as so that the correct ordering of random product vectors can be reflected. The estimator of the negativity moment is obtained byR so that Computing such an estimator on a system represented by TN is pursued by separately computing each element The calculation of a single element is illustrated in Fig. 4, and is equivalent in terms of complexity to an expectation value calculation: For example, for a two-dimensional PEPS, the space complexity is O χ 2 D 2(w1+1) + D 4 d , and the time complexity is O nw 2 χD 2(w1+1) χ + D 2 + D 4 d . The calculation is repeated M times over realizations of the respective random vectors and the outcomes are averaged in order to get an estimate for the desired quantity. Thus the cost of the calculation of a single density matrix element given above, times the number of repetitions M , which is discussed in Secs. III D and IV B, as well as in Appendix A.

C. Candidate probability measures
The required property of the random vectors, captured in Eq. (18), can be naturally obtained in a wealth of ways: After all, all that is required is to have up to normalization a spherical 1-design property. Still, since we do not require the vectors to necessarily constitute a spherical complex 2-design, the second moments will depend on the specific choice of the probability measure. For example, this can be done by choosing the vectors randomly out of some orthogonal basis, or several orthogonal bases. For prime dimension d, the clock and shift operators, Weyl operators, or simply d−dimensional Pauli matrices, are defined to be the operators (27) where ω := e 2πi/d . The j-th normalized eigenvector of the i-th d-dimensional Pauli matrix is denoted by |p (i,j) , and we note that for prime d, the number of non-commuting Pauli matrices is d + 1. We then compare two possible distributions that are particularly practical in the context given: First is the 'full-basis' distribution, The vectors are normalized such that Eq. (18) For non-prime d, the Pauli matrices can be defined to be tensor products of the matrices in Eq. (27) in the dimensions of the factors of d. In this case, the vector distributions defined as in Eqs. (28,29), but with the eigenbases of the independent products of the clock and shift operators. We compare the two probability measures (and discuss why it is sufficient to only consider these distributions) in Secs. III D and IV B and in Appendix A below. For states represented efficiently by TN, the partial-basis distribution (with an optimized basis choice, as detailed in Section IV C) turns out to be more efficient, and therefore most of the presented results were obtained using this method.

D. Required number of repetitions
The method suggested makes use of random quantum states for the estimation of entanglement measures. When drawing random vectors from the probability measures indicated above, one finds for the probability of deviating from the expectation to be bounded by for real k and σ(p n ) 2 := Var(p n ). This is true by virtue of Chebychev's inequality, a large deviation bound. We here and in the following suppress the dependence on ρ A . For M repetitions, the variance Var M of the estimator of the mean E(p n ) is given by Since then for a given ε > 0, the number of required repetitions scales as M = Var(p n )/(p 2 n ε 2 ). The characteristics of the method lead us to expect an exponential dependence of the required number of repetitions M on system size (whilst as can be seen below, a weak one). We thus define the scaling factor ξ n by We use the notation ξ n for the scaling factor ofR n as well, since the scaling factors for both properties are expected to behave similarly.
In Appendix A, we show that the variance of the estimators defined in Eqs. (23) and (25) is given by and where Given the product structure of the probability measure, this expression is found to be In a coordinate representation, this is found to be for the full-basis distribution, and i,j;k,l = δ i,k δ j,l + δ i,l δ j,k − δ i,k δ j,l δ i,l for the partial-basis distribution. To give an even more specific example in a coordinate dependent form, for d = 2, we have for the full-basis and partial-basis distributions, respectively. The variance of the symmetry-resolved moments estimator is shown to be bounded from above by Eq. (34) in Appendix A. We first focus on discussing the full-basis distribution. As shown in Appendix A, the squared coefficient of variation in a product state ρ = ⊗ α∈A |ψ α ψ| α and in a maximally mixed Here, we use the partial-basis distribution, Eq. (29). The normalization of M by the scaling factor ξn in the horizontal axis is extracted from Fig. 8 and presented in Table I. The presented plots are averaged over 20 different permutations of the estimation results. case (i.e., subsystem A being maximally entangled with the rest of the system), and respectively. The same results apply for the PT moments. For example, for d = 2, for a product state and for a maximally mixed state. Note that the scaling factors obtained below for the benchmark models, as displayed in Table  I, are in-between these two extreme cases.
As for the partial-basis distribution, as explained in detail in Appendix A, the highest variance, hence the largest additive sampling error, for both the RDM moments and PT moments arises for where |ψ E stands for a state vector of the form with an equal magnitude of the amplitude for each state in the computational basis. In this case, the contribution of each site to the first term of the variance is and the overall variance is given by For example, for d = 2, the scaling factor in such a case is ξ n = (3/2) n . The best case is ρ A = |0 0| ⊗N A (or any other basis state in the computational basis), in which the variance is 0. Both cases are completely disentangled, and the moments equal 1. Therefore, the former of these is not maximal in terms of the squared coefficient of variation, which is determined by the ratio between the standard deviation and the expected value. However, the analysis of the variance itself already serves to demonstrate that that the choice of basis for the vector |v can have a significant impact on the variance and due to that on the performance of the algorithm, as discussed in Section IV C below. The coefficient of variation in the maximally nixed case in this method is calculated in Appendix A and equals The two cases represented in Eqs. (42), (43), and (48) above, in which the variance can be calculated exactly, are not promised to be the best or worst case for the two distribution methods. For a two-qubit system we have performed a gradient descent search for the extreme cases in both distributions, where a basis optimization (as described in Section IV C) has been included in the partial-basis method. Fig. 5 presents the results, which strongly support the hypothesis that the maximally-and minimally-mixed cases are indeed the extreme cases for the method's performance.
The best distribution choice is therefore case-dependent: For moderate or large n and highly mixed cases, the fullbasis distribution is advantageous. For small n or weakly entangled cases, the partial-basis method with a smart basis choice (as discuss in Section IV C) is more beneficial. Based on Eq. (50), the partial-basis seems to have a poor performance on highly entangled (mixed) states. However, states represented efficiently by TN feature an entanglement that exhibits an area law, which restrains the entanglement of relevant states to begin with. Note that even in this worst case our method is still favourable compared to the the time required for an exact diagonalization of the RDM, which scales as O(d 3N A ), for intermediate n. The variances calculated above compare favorably with the variances in the experimental sampling-based protocols [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] , as calculated in Ref. 24: For d = 2, the relative variances for n = 2, 3 are shown to be Var(p 2 )/p 2 2 ≥ (8/p 2 2 ) max 2 N A p 2 , 2 1.5N A and Var(p 3 )/p 2 3 ≥ 39/p 2 3 max 2 N A p 2 2 , 2 1.5N A p 2 , 2 2N A . It would be interesting to use the analysis above as a basis for a study regarding the number of local samples required for the estimation of Rényi moments in general. Naively, the moments are defined as a function of the full RDM, which has d 2N A elements, hence should require a comparable number of samples. However, an extraction of a single degree of freedom is expected to require a smaller number of samples, as is the case in small ns in our method. Such analysis may point to the amount of information contained in Rényi moments of different ranks.

E. Symmetry-resolved RDM moments
Using the locality of the phase operator from Eq. (8), e iϕQ A = ⊗ α∈A e iϕQα , we can extract the flux-resolved moments from the TN, and substitute them into Eq. (9) to get the symmetry-resolved moments. The estimator of the fluxresolved moment is similar to the estimator of the full moments in Eq. (23) and is obtained from The estimator of the charge-resolved moment is obtained from A similar analysis of symmetry-resolved PT moments has been done in Ref. 52, and the extension to their estimation is natural. Below we estimate the symmetry-resolved RDM moments for the XX model and its conserved total S z . For this model, the symmetry-resolved moments can be obtained exactly following Ref. 49. The expected and extracted results for q → p n (ρ A , q) for n = 2, 3 are displayed in Fig. 6.

A. Specific tests
We have tested the model against the exactly solvable twodimensional toric code model and one-dimensional XX model as detailed in Section II D, employing the TensorNetwork library 126 . The precision of the estimation for both models as a function of M , the number of samples of the expressions in Eqs. (21), (23), (25), and (51) is shown in Fig. 7. The results are obtained using the partial-basis distribution, Eq. (29), and are optimized based on the analysis in Section IV C below. In order to reduce the numerical noise in the dependence of the precision in M , we averaged this dependence over several permutations of the M repetitions. While the required number of repetitions M (for a given allowed error ε > 0) is exponential in system size (as discussed above in Section III D), it has a relatively small base ξ n . When considering the significant decrease in required memory space, our method can become advantageous for systems around N A = 20, for which the method described in Section II B can become too heavy in memory demands for a standard computer workstation.

B. Variance estimation
Here we follow the analysis of the scaling factor ξ n in Section III D and estimate the scaling factors in the benchmark models, in order to get some idea regarding the variance in the general case. The scaling factors ξ n obtained for the toric code and XX ground states, for n = 2, 3, 4 are estimated numerically. In Appendix B, we demonstrate an exact calculation of the expressions in Eqs. (34) and (35) for a narrow strip-like system in the toric code model, and show that the resulting expressions agree with the numerically estimated scaling factors.
We emphasize that the models used are not specifically suitable for the method, and are not expected to have a low scaling factors based on Eqs. (34) and (35). The estimated scaling factor ξ n can therefore be considered typical. Fig. 8 presents the estimated variances of the benchmark models in the full-basis distribution and partial-basis distri-  Table I. Note that for the toric code ground state shown in (a), the behaviours of the 3rd RDM moment and the 3rd PT moment are expected to be identical, and are slightly different only due to the random nature of our protocol.
bution after the basis-choice optimization detailed in Section IV C. The scaling factor can be extracted from the dependence of the variance on N A . We see that the scaling factors are smaller than the worst case presented above. In the XX model, for which the entanglement is log-dependent in system size, we get a better performance with the basis-optimized partialbasis distribution. In the more strongly entangled (and therefore less basis-sensitive) toric code model, the difference in performance between the two methods is clearly less significant.

C. Dependence on basis choice
In the partial-basis distribution, as shown above, the largest and smallest additive variance values both correspond to a completely disentangled case, and the difference between the two stems from the single-particle basis choice alone. This can be understood by the decomposition as can be seen from Eqs. (40) and (41), that demonstrates the orientation dependence of E in this case (in contrast with in the full-basis case). For a translationally invariant system, one can expect that the optimal basis choice for each site will be the same. We can now attempt to decrease the variance by finding the basis for which Var(p 1 ) is minimal, and use this basis for the estimation of higher moments. We have tested this idea against the first moment of the two benchmark models, by rotating the random vectors and finding the best basis, as demonstrated in Fig. 9. We plot the scaling factor for n = 1, ξ 1 , for the two benchmark models as a function of the basis choice. We then compare the best and worst choices of φ, θ ∈ [0, 2π) and extract the variance of the higher moments in the corresponding bases, as summarized in Table I. We can see that the p 1 case acts as a good indicator for the basis choice of the moments in higher ns, and allows for a smart basis choice which decreases the variance.     Table I. The extracted scaling factor defined in Eq. (33), ξn, for n = 2, 3, 4 for the best and worst case in the toric code ground state extracted from Fig. 9 and the XX model ground state. The scaling factors obtained in the full-basis method are displayed for comparison. The worst and best case in ξ1 are shown to predict well the dependence of ξn on the basis choice.

V. CONCLUSIONS AND OUTLOOK
In high-dimensional TN states, the naive computation of entanglement is highly sensitive to the size of the system (when explicitly extracting the RDM) or the bond dimension of the site tensors and boundaries (when performing the replica trick). We developed a method for estimating RDM moments in Eq. (2) and PT moments in Eq. (5) of such systems, as well as their symmetry-resolved components in Eqs. (7) and (8), without fully reconstructing the density matrix or contracting several copies of the state. The method uses randomization in order to correlate separate copies of the TN state, allowing for the estimation of properties that are defined using more than one copy of the RDM. Though we are inspired by recent experimental protocols 10-23 , we developed a completely new algorithm which is suitable to classical simulations, takes advantage of their strengths such as ability to estimate the expectation value of non-Hermitian operators, and avoids their weakness in sampling the outcomes of random measurements.
We have demonstrated our method with the iPEPS representation of the toric code ground state and the MPS representation of the XX ground state, and compared the results with analytical calculations. The method can be readily used for any tensor network ansatz representing a spin or bosonic system, and provide information on the entanglement of systems that were formerly unreachable by today's computers due to a strong exponential dependence of the memory space in the moment degree n. Additionally, our method is advantageous for nontrivial partitions in one or higher spatial dimensions, such as the checkerboard partition 115 or random partition 116 , for which the moments are hard to calculate even for onedimensional MPS.
We compare two options for the random distribution, where each of the methods turns out to be suitable for different cases. For small ns, the scaling of required samples number M with system size N A turns out to be lower than the scaling of RDM size, and can have implications regarding the information contained in these moments. It would be interesting to try and develop a sampling-based method which incorporates nonphysical operations and compare its performance with ours. The analysis of such a protocol may shed more light on the power of Rényi moments.
The method should be suitable to fermionic PEPS 103,104 , and can be generalized to additional Rényi measures, such as participation entropies, used for the detection of many-body localization 127 . Exploring the possibility of derandomizing the algorithm, similarly to the recent results of Huang et al. 128 would also be interesting. In contrast to setting of shadow estimation, the very quantum state is already classically efficiently represented, and computing overlaps with suitable random vectors gives rise to an effective estimation of entanglement properties. Now that the option to use non-physical sampling has been opened, it can be expanded to various platforms, including experimental setting with a vectorized density matrix. It is the hope that this work contributes to the program of exploiting the power of random measurements in quantum physics, even in situations where the sampling scheme itself is not reflected by physical operations.
Below we derive Eqs. (34), (40) and (41) of the main text, and use them to find density matrices which extremize the variance. First, we write the expression for the RDM moments estimator explicitly. In the above coordinate independent fashion, this derivation is straightforward. The n-th RDM moment is obtained as the expectation of where again, the product state vectors |V (1) , . . . , |V (n) are drawn in an i.i.d. fashion from the same probability measure. In expectation, we find Similarly, for the PT moments, one can make use of random vectors of the form for i, j = 1, . . . , n, so that the estimator of the negativity moment iŝ since one simply finds by performing partial transposes in all terms so that indeed the correct moment of the partially transposed operator is recovered. The variance can then be calculated from the expectation ofp from which the square of p n (ρ A ) is subtracted. The subtle point is now that projections appear twice rather than once. This can be reflected by making use of two tensor factors. Upon reordering the tensor entries, one immediately finds the expression In expectation, this is E(p n (ρ A ) 2 ) = Tr((ρ ⊗2 and In this way, one finds the expression for the variance Here it is relevant that the frames made use of do not necessarily have to constitute complex spherical 2-designs for estimating the above entanglement measures in an unbiased fashion, so that the average does not necessarily resemble that of the Haar average, and may depend on the ensemble. We now see why it is meaningful to consider the two probability measures specified above: Drawing vectors from the eigenbasis of a single Pauli matrix, i.e., randomly sampling RDM elements in this basis, would give, for example if we drew from the eigenbasis of Z, i,k;j,l = δ i,j δ k,l + δ i,k δ j,l + δ i,l δ j,k − 2δ i,k δ j,l δ i,l . (A11) This matrix equals the matrix obtained above, with added positive terms, hence, can only increase the variance. The same argument can be made for any number of Pauli matrices between 2 and d − 1. The Pauli matrices constitute a unitary 1-design 27 , which is the requirement for them to be a universal measure for the estimators in Eqs. (21), (23), (25), and (52). Therefore, it is unnecessary to consider additional distributions. In the maximally mixed case, ρ A = (I d /d) ⊗N A , the normalized first term is for the full-basis and for the partial-basis, respectively. In a product state, ρ A = ⊗ α∈A |ψ α ψ| α in the full-basis distribution, one finds (A14) The performance of the partial-basis distribution for a product state can be analyzed as follows: E [α] is a block diagonal matrix, with d blocks of the form B 1 = (1) and d 2 blocks of the form The largest eigenvalue of α∈A E [α] is thus 2 N A , and corresponds to an eigenvector of the form ⊗ α∈A |i, j , where i and j correspond to the spin of the site α in the two copies and i, j = 1, . . . , d, for i = j. However, since ρ ⊗2 A is constructed from two identical copies of ρ A , vectors of the form above cannot be the only contributors to the RDM. The RDM with the largest possible variance has an equal weight to all vectors of the form above, which means it will have the form Then, the contribution of each site to the first term of the variance is The RDM with smallest possible variance will be, for example, or any other product state in the computational basis. In this case, the first term of the variance sums up to 1 and Var(p n ) = 0. In both cases A is disentangled from its environment, which demonstrates that the variance of the estimated value depends on the basis choice for the vectors |v . The variance of the flux-resolved moment estimators of the complex valued random variable defined in Eq. (51) can similarly be computed from  (34) and (35), for n = 2 and a subsystem composed of l = 3 stars. Depicted in pink and blue are stabilizers under which the subsystem RDM is invariant, and therefore each of them can be applied to any of the copies and contribute to the trace (Note that σ y = iσ z σ x ). In orange and purple is a pair of operators which commute with the star operator ⊗i∈sσ x i due to having an even number of σ y matrices. Such a pair of operators which commute with the stabilizer operator will also contribute to the trace, provided they are applied on both sides of the same copy of ρ ⊗2 A factors. (c) An example element of the transfer matrix T (2) . Applying the operators in orange to the (l + 1)-th star, given that the operator configuration on the l-th operator is the one in pink, is allowed, and the expression will be multiplied by 2 −4 due to the four 1 2 σ y ⊗ σ y applied to the (l + 1)-th star. (d) An additional example element in the transfer matrix T (2) . Applying the operators in orange to the (l + 1)-th star will not contribute to the trace, since this operator configuration do not commute with the plaquette or star terms.
which is bounded from above by the variance for the nonresolved case, as the sum of the first term is composed of terms with the same absolute values, but with added phases. The estimator for the variance of the symmetry-resolved moments is therefore also bounded by Var(p n (q)) = 1 N A ϕ Var(p n (ϕ)) ≤ 1 N A ϕ Var(p n ) = Var(p n ). (A20) N A 9DU(p)/p 2 p 1 p 2 p 3 Figure 11. The estimated variances of a strip-like system in the toric code model.

Appendix B: Explicit variance calculation for the toric code
Here, we demonstrate how the variance can be calculated exactly in the toric code model for a subsystem shaped as a narrow strip (Fig. 10a for the partial-basis distribution, and compare it to the extracted variance. We start from Eqs. (34), (40) and (41), and the decomposition which applies for the partial-basis method in d = 2. The full E matrix can be written as where S is an N A -long configuration of the operators I, 1 2 σ x ⊗ σ x , 1 2 σ y ⊗ σ y . We use the local symmetry of the toric code ground state ⊗ i∈s σ x i ρ = ⊗ i∈p σ z i ρ = ρ ⊗ i∈s σ x i = ρ ⊗ i∈p σ z i = ρ, (B2) in order to distinguish configurations S which will contribute to Eq. (34). An allowed configuration S will contain any number of star (⊗ i∈s σ x i ) or plaquette (⊗ i∈p σ z i ) operators, also called the stabilizers. Pairs of some operator which commutes with the stabilizers and act on two sides of the same copy of ρ ⊗2 are also allowed, as well as combinations of operators which can be transformed into such pairs by a multiplication of the operators by stabilizers. This is illustrated in Fig. 10b.
We calculate the variance for a narrow system of dimension 2×2l, as depicted in Fig. 10a. Such a system is composed of a chain of l contiguous stars. For a moment of rank n, we think of n copies of the subsystem, and write a 3 4n -dimensional vector of combinations of operators [I, σ x ⊗ σ x , σ y ⊗ σ y ] on a single star in all n copies. We now write a transfer matrix T (n) which takes the operator combinations on the l-th star to the contributing combinations on the (l + 1)-th star: T (n) i,j Exact variance Estimated variance Table II. Exact dependence of the variance on system size for the toric code model as calculated in by the transfer matrix method compared to the variance estimated from sampling of the expressions in Eqs. (34) and (35). equals the contribution of an operator combination with operators j on the (l + 1)-th star to the variance in Eq. (34), given that the combination on the l-th star is i, where the symmetries in Eq. (B2) are considered, as well as the 1 2 factors in Eq. (B1). For clarity, we give a specific example in Fig. 10c-d. With these definitions, the first term in the left side of Eq. (34) for a subsystem of 2 × 2l sites is c 0 | (T (n) ) l−1 |c 0 , where |c 0 is the vector of allowed contributions for the edges of the subsystem as can be deduced from the ground states of the toric code model in Eq. (14). The dependence of this term in the system size is thus O λ 4l max , where is the largest eigenvalue of the Hermitian T (n) . One may, in fact, work with equivalent but much smaller transfer matrices, by considering only the two left sites of a star rather than the whole star. This allows decreasing the transfer matrix to dimension 3 2n × 3 2n . For n ≤ 3 the matrix T (n) can be extracted and diagonalized exactly. We have performed an estimation of the variance for such strip-like systems, similarly to the one done in Fig. 8 in the main text. The estimated variances are displayed in Fig. 11. We have compared the results obtained exactly using the transfer matrix to the numerical variances and got a good agreement, demonstrating the accuracy of our PEPS calculations, as can be seen in the Table II.